Open64: osprey/libm/mips/dcis.c Source File

00001 /*
00002  * Copyright 2003, 2004, 2005, 2006 PathScale, Inc.  All Rights Reserved.
00003  */
00004 
00005 
00006 /*
00007 
00008   Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
00009 
00010   This program is free software; you can redistribute it and/or modify it
00011   under the terms of version 2.1 of the GNU Lesser General Public License 
00012   as published by the Free Software Foundation.
00013 
00014   This program is distributed in the hope that it would be useful, but
00015   WITHOUT ANY WARRANTY; without even the implied warranty of
00016   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
00017 
00018   Further, this software is distributed without any warranty that it is
00019   free of the rightful claim of any third person regarding infringement 
00020   or the like.  Any license provided herein, whether implied or 
00021   otherwise, applies only to this software file.  Patent licenses, if
00022   any, provided herein do not apply to combinations of this program with 
00023   other software, or any other product whatsoever.  
00024 
00025   You should have received a copy of the GNU Lesser General Public 
00026   License along with this program; if not, write the Free Software 
00027   Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307, 
00028   USA.
00029 
00030   Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
00031   Mountain View, CA 94043, or:
00032 
00033   http://www.sgi.com
00034 
00035   For further information regarding this notice, see:
00036 
00037   http://oss.sgi.com/projects/GenInfo/NoticeExplan
00038 
00039 */
00040 
00041 /* ====================================================================
00042  * ====================================================================
00043  *
00044  * Module: dcis.c
00045  * $Revision: 1.5 $
00046  * $Date: 04/12/21 14:58:21-08:00 $
00047  * $Author: bos@eng-25.internal.keyresearch.com $
00048  * $Source: /home/bos/bk/kpro64-pending/libm/mips/SCCS/s.dcis.c $
00049  *
00050  * Revision history:
00051  *  05-Sep-96 - Original Version
00052  *
00053  * Description: source code for dcis function
00054  *
00055  * ====================================================================
00056  * ====================================================================
00057  */
00058 
00059 #ifdef _CALL_MATHERR
00060 #include <stdio.h>
00061 #include <math.h>
00062 #include <errno.h>
00063 #endif
00064 
00065 #include "libm.h"
00066 #include "complex.h"
00067 
00068 #if defined(mips) && !defined(__GNUC__)
00069 extern  dcomplex  __dcis(double);
00070 
00071 #pragma weak __dcis = __libm_dcis
00072 #endif
00073 
00074 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00075 extern dcomplex __libm_dcis(double);
00076 #pragma weak __dcis
00077 dcomplex __dcis(double x) {
00078   return __libm_dcis(x);
00079 }
00080 #elif defined(__GNUC__)
00081 extern  dcomplex  __libm_dcis(double);
00082 dcomplex    __dcis() __attribute__ ((weak, alias ("__libm_dcis")));
00083 #endif
00084 
00085 static const  du  Qnan =
00086 {D(QNANHI, QNANLO)};
00087 
00088 static const du half =
00089 {D(0x3fe00000, 0x00000000)};
00090 
00091 static const du one =
00092 {D(0x3ff00000, 0x00000000)};
00093 
00094 static const du Twop19xpi =
00095 {D(0x413921fb, 0x54442d18)};
00096 
00097 static const du rpiby2 =
00098 {D(0x3fe45f30, 0x6dc9c883)};
00099 
00100 static const du piby2hi =
00101 {D(0x3ff921fb, 0x54400000)};
00102 
00103 static const du piby2lo =
00104 {D(0x3dd0b461, 0x1a600000)};
00105 
00106 static const du piby2tiny =
00107 {D(0x3ba3198a, 0x2e037073)};
00108 
00109 static const du ph =
00110 {D(0x3ff921fb, 0x50000000)};
00111 
00112 static const du pl =
00113 {D(0x3e5110b4, 0x60000000)};
00114 
00115 static const du pt =
00116 {D(0x3c91a626, 0x30000000)};
00117 
00118 static const du pe =
00119 {D(0x3ae8a2e0, 0x30000000)};
00120 
00121 static const du pe2 =
00122 {D(0x394c1cd1, 0x29024e09)};
00123 
00124 static const du Ph =
00125 {D(0x3ff921fb, 0x54000000)};
00126 
00127 static const du Pl =
00128 {D(0x3e110b46, 0x10000000)};
00129 
00130 static const du Pt =
00131 {D(0x3c5a6263, 0x3145c06e)};
00132 
00133 /* coefficients for polynomial approximation of sin on +/- pi/4 */
00134 
00135 static const du S[] =
00136 {
00137 {D(0x3ff00000, 0x00000000)},
00138 {D(0xbfc55555, 0x55555548)},
00139 {D(0x3f811111, 0x1110f7d0)},
00140 {D(0xbf2a01a0, 0x19bfdf03)},
00141 {D(0x3ec71de3, 0x567d4896)},
00142 {D(0xbe5ae5e5, 0xa9291691)},
00143 {D(0x3de5d8fd, 0x1fcf0ec1)},
00144 };
00145 
00146 /* coefficients for polynomial approximation of cos on +/- pi/4 */
00147 
00148 static const du C[] =
00149 {
00150 {D(0x3ff00000, 0x00000000)},
00151 {D(0xbfdfffff, 0xffffff96)},
00152 {D(0x3fa55555, 0x5554f0ab)},
00153 {D(0xbf56c16c, 0x1640aaca)},
00154 {D(0x3efa019f, 0x81cb6a1d)},
00155 {D(0xbe927df4, 0x609cb202)},
00156 {D(0x3e21b8b9, 0x947ab5c8)},
00157 };
00158 
00159 #ifdef _32BIT_MACHINE
00160 
00161 /* tables used to do argument reduction for args between +/- 16 radians;
00162    the sum of the high and low values of the kth entry is (k - 10)*pi/2
00163 */
00164 
00165 static const du tblh[] =
00166 {
00167 {D(0xc02f6a7a, 0x2955385e)},
00168 {D(0xc02c463a, 0xbeccb2bb)},
00169 {D(0xc02921fb, 0x54442d18)},
00170 {D(0xc025fdbb, 0xe9bba775)},
00171 {D(0xc022d97c, 0x7f3321d2)},
00172 {D(0xc01f6a7a, 0x2955385e)},
00173 {D(0xc01921fb, 0x54442d18)},
00174 {D(0xc012d97c, 0x7f3321d2)},
00175 {D(0xc00921fb, 0x54442d18)},
00176 {D(0xbff921fb, 0x54442d18)},
00177 {D(0x00000000, 0x00000000)},
00178 {D(0x3ff921fb, 0x54442d18)},
00179 {D(0x400921fb, 0x54442d18)},
00180 {D(0x4012d97c, 0x7f3321d2)},
00181 {D(0x401921fb, 0x54442d18)},
00182 {D(0x401f6a7a, 0x2955385e)},
00183 {D(0x4022d97c, 0x7f3321d2)},
00184 {D(0x4025fdbb, 0xe9bba775)},
00185 {D(0x402921fb, 0x54442d18)},
00186 {D(0x402c463a, 0xbeccb2bb)},
00187 {D(0x402f6a7a, 0x2955385e)},
00188 };
00189 
00190 static const du tbll[] =
00191 {
00192 {D(0xbcc60faf, 0xbfd97309)},
00193 {D(0xbcc3daea, 0xf976e788)},
00194 {D(0xbcc1a626, 0x33145c07)},
00195 {D(0xbcbee2c2, 0xd963a10c)},
00196 {D(0xbcba7939, 0x4c9e8a0a)},
00197 {D(0xbcb60faf, 0xbfd97309)},
00198 {D(0xbcb1a626, 0x33145c07)},
00199 {D(0xbcaa7939, 0x4c9e8a0a)},
00200 {D(0xbca1a626, 0x33145c07)},
00201 {D(0xbc91a626, 0x33145c07)},
00202 {D(0x00000000, 0x00000000)},
00203 {D(0x3c91a626, 0x33145c07)},
00204 {D(0x3ca1a626, 0x33145c07)},
00205 {D(0x3caa7939, 0x4c9e8a0a)},
00206 {D(0x3cb1a626, 0x33145c07)},
00207 {D(0x3cb60faf, 0xbfd97309)},
00208 {D(0x3cba7939, 0x4c9e8a0a)},
00209 {D(0x3cbee2c2, 0xd963a10c)},
00210 {D(0x3cc1a626, 0x33145c07)},
00211 {D(0x3cc3daea, 0xf976e788)},
00212 {D(0x3cc60faf, 0xbfd97309)},
00213 };
00214 #endif
00215 
00216 
00217 /* ====================================================================
00218  *
00219  * FunctionName    __libm_dcis
00220  *
00221  * Description    computes cos(arg) + i*sin(arg)
00222  *
00223  * Note:  Routine __libm_dcis() computes cos(arg) + i*sin(arg) and should
00224  *    be kept in sync with sin() and cos() so that it returns 
00225  *    the same result as one gets by calling cos() and sin() 
00226  *    separately.
00227  *
00228  * ====================================================================
00229  */
00230 
00231 dcomplex
00232 __libm_dcis(double x)
00233 {
00234 #ifdef _32BIT_MACHINE
00235 int ix, xpt, m, l;
00236 #else
00237 long long ix, xpt, m, l;
00238 #endif
00239 
00240 double  xsq;
00241 double  cospoly, sinpoly;
00242 int  n;
00243 dcomplex  result;
00244 double  absx;
00245 double  z, dn;
00246 double  dn1, dn2;
00247 double  s, ss;
00248 double  t, w, ww;
00249 #ifdef _CALL_MATHERR
00250 struct exception  exstruct;
00251 #endif
00252 
00253   /* extract exponent of x and 1 bit of mantissa for some quick screening */
00254 
00255 #ifdef _32BIT_MACHINE
00256 
00257   DBLHI2INT(x, ix); /* copy MSW of x to ix  */
00258 #else
00259   DBL2LL(x, ix);    /* copy x to ix */
00260 #endif
00261   xpt = (ix >> (DMANTWIDTH-1));
00262   xpt &= 0xfff;
00263 
00264   if ( xpt < 0x7fd )
00265   {
00266     /*   |x| < 0.75   */
00267 
00268     if ( xpt >= 0x7c2 )
00269     {
00270             /*   |x| >= 2^(-30)   */
00271 
00272             xsq = x*x;
00273 
00274       result.dreal = (((((C[6].d*xsq + C[5].d)*xsq +
00275           C[4].d)*xsq + C[3].d)*xsq +
00276           C[2].d)*xsq + C[1].d)*xsq + C[0].d;
00277 
00278       result.dimag = (((((S[6].d*xsq + S[5].d)*xsq +
00279           S[4].d)*xsq + S[3].d)*xsq +
00280           S[2].d)*xsq + S[1].d)*(xsq*x) + x;
00281 
00282       return ( result );
00283     }
00284     else
00285     {
00286       result.dreal = 1.0;
00287       result.dimag = x;
00288 
00289       return ( result );
00290     }
00291   }
00292 
00293 #ifdef _32BIT_MACHINE
00294 
00295   if (xpt < 0x806)
00296   {
00297     /*   |x| < 16.0   */
00298 
00299     /*  do a table based argument reduction to +/- pi/4  */
00300 
00301     dn = x*rpiby2.d;
00302 
00303     n = ROUND(dn);
00304 
00305     /*  compute x - n*pi/2  */
00306 
00307     x = x - tblh[n+10].d;
00308     x = x - tbll[n+10].d;
00309 
00310     goto L;
00311 
00312   } else
00313 #endif
00314   if ( xpt < 0x827 )
00315   {
00316     /*  |x| < 1.5*2^20  */
00317 
00318     dn = x*rpiby2.d;
00319     n = ROUND(dn);
00320     dn = n;
00321 
00322     x = x - dn*piby2hi.d;
00323     x = x - dn*piby2lo.d;
00324     x = x - dn*piby2tiny.d; /* x = x - n*pi/2 */
00325 
00326 L:
00327     xsq = x*x;
00328 
00329     cospoly = (((((C[6].d*xsq + C[5].d)*xsq +
00330         C[4].d)*xsq + C[3].d)*xsq +
00331         C[2].d)*xsq + C[1].d)*xsq + C[0].d;
00332 
00333     sinpoly = (((((S[6].d*xsq + S[5].d)*xsq +
00334         S[4].d)*xsq + S[3].d)*xsq +
00335         S[2].d)*xsq + S[1].d)*(xsq*x) + x;
00336 
00337     if ( n&1 )
00338     {
00339       if ( n&2 )
00340       {
00341         /*
00342          *  n%4 = 3
00343          *  result is sin(x) - i*cos(x)
00344          */
00345 
00346         result.dreal = sinpoly;
00347         result.dimag = -cospoly;
00348       }
00349       else
00350       {
00351         /*
00352          *  n%4 = 1
00353          *  result is -sin(x) + i*cos(x)
00354          */
00355 
00356         result.dreal = -sinpoly;
00357         result.dimag = cospoly;
00358       }
00359 
00360       return (result);
00361     }
00362 
00363     if ( n&2 )
00364     {
00365       /*
00366        *  n%4 = 2
00367        *  result is -cos(x) - i*sin(x)
00368        */
00369 
00370       result.dreal = -cospoly;
00371       result.dimag = -sinpoly;
00372     }
00373     else
00374     {
00375       /*
00376        *  n%4 = 0
00377        *  result is cos(x) + i*sin(x)
00378        */
00379 
00380       result.dreal = cospoly;
00381       result.dimag = sinpoly;
00382     }
00383 
00384     return(result);
00385 
00386   }
00387   else if ( xpt < 0x836 )
00388   {
00389     /*  |x| < 2^28  */
00390 
00391     dn = x*rpiby2.d;
00392     n = ROUND(dn);
00393     dn = n;
00394 
00395     x = x - dn*ph.d;
00396     x = x - dn*pl.d;
00397     x = x - dn*pt.d;
00398     x = x - dn*pe.d;
00399     x = x - dn*pe2.d;
00400 
00401     goto L;
00402   }
00403   else if ( xpt < 0x862 )
00404   {
00405     /*  |x| < 2^50  */
00406 
00407     absx = fabs(x);
00408 
00409     dn = z = absx*rpiby2.d;
00410 
00411     /* round dn to the nearest integer */
00412 
00413 #ifdef _32BIT_MACHINE
00414 
00415     DBLHI2INT(dn, l);
00416     m = (l >> DMANTWIDTH);
00417     m &= 0x7ff;
00418 
00419     /* shift off fractional bits of dn */
00420 
00421     DBLLO2INT(dn, l);
00422 
00423     l >>= (0x433 - m);
00424     n = l;
00425     l <<= (0x433 - m);
00426     INT2DBLLO(l, dn);
00427 #else
00428     DBL2LL(dn, l);
00429     m = (l >> DMANTWIDTH);
00430     m &= 0x7ff;
00431 
00432     /* shift off fractional bits of dn */
00433 
00434     l >>= (0x433 - m);
00435     n = l;
00436     l <<= (0x433 - m);
00437     LL2DBL(l, dn);
00438 #endif
00439     /* adjust dn and n if the fractional part of dn 
00440        was >= 0.5
00441     */
00442 
00443     n &= 3;
00444 
00445     if ( (z - dn) >= half.d )
00446     {
00447       dn += one.d;
00448       n += 1;
00449     }
00450 
00451   /* compute x - dn*Ph - dn*Pl - dn*Pt by dividing dn into
00452      two parts and using the Kahan summation formula
00453   */
00454 
00455     /* split dn into 2 parts */
00456 
00457     dn1 = dn;
00458 
00459 #ifdef _32BIT_MACHINE
00460 
00461     DBLLO2INT(dn1, m);
00462     m >>= 28;
00463     m <<= 28;
00464     INT2DBLLO(m, dn1);
00465 #else
00466     DBL2LL(dn1, m);
00467     m >>= 28;
00468     m <<= 28;
00469     LL2DBL(m, dn1);
00470 #endif
00471     dn2 = dn - dn1;
00472 
00473     z = absx - dn1*Ph.d;  /* this operation is exact */
00474 
00475     t = dn2*Ph.d;
00476     s = z - t;
00477     ss = z - s - t;   /* correction term */
00478 
00479     t = ss - dn1*Pl.d;
00480     w = s + t;
00481     ww = s - w + t;   /* correction term */
00482 
00483     t = ww - dn2*Pl.d;
00484     s = w + t;
00485     ss = w - s + t;   /* correction term */
00486 
00487     t = ss - dn1*Pt.d;
00488     w = s + t;
00489     ww = s - w + t;   /* correction term */
00490 
00491     t = ww - dn2*Pt.d;
00492     z = w + t;    /* z = reduced arg */
00493 
00494     if ( x < 0.0 )
00495     { /* adjust for sign of arg */
00496 
00497       z = -z;
00498       n = -n;
00499     }
00500 
00501     x = z;
00502 
00503     goto L;
00504   }
00505 
00506 
00507   if (x != x)
00508   {
00509     /* x is a NaN; return a pair of quiet NaNs */
00510 
00511     NAN_SETERRNO(EDOM);
00512 
00513     result.dreal = Qnan.d;
00514     result.dimag = Qnan.d;
00515     return (result);
00516   }
00517 
00518   /* just give up and return 0.0 */
00519 
00520   result.dreal = 0.0;
00521   result.dimag = 0.0;
00522   return (result);
00523 }
00524