Open64: osprey/libF77/d_cis.c Source File

00001 /*
00002 
00003   Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
00004 
00005   This program is free software; you can redistribute it and/or modify it
00006   under the terms of version 2.1 of the GNU Lesser General Public License 
00007   as published by the Free Software Foundation.
00008 
00009   This program is distributed in the hope that it would be useful, but
00010   WITHOUT ANY WARRANTY; without even the implied warranty of
00011   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
00012 
00013   Further, this software is distributed without any warranty that it is
00014   free of the rightful claim of any third person regarding infringement 
00015   or the like.  Any license provided herein, whether implied or 
00016   otherwise, applies only to this software file.  Patent licenses, if
00017   any, provided herein do not apply to combinations of this program with 
00018   other software, or any other product whatsoever.  
00019 
00020   You should have received a copy of the GNU Lesser General Public 
00021   License along with this program; if not, write the Free Software 
00022   Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307, 
00023   USA.
00024 
00025   Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
00026   Mountain View, CA 94043, or:
00027 
00028   http://www.sgi.com
00029 
00030   For further information regarding this notice, see:
00031 
00032   http://oss.sgi.com/projects/GenInfo/NoticeExplan
00033 
00034 */
00035 
00036 
00037 
00038 #include "cmplx.h"
00039 #include <errno.h>
00040 #include "moremath.h"
00041 
00042 extern  double  __libm_qnan_d;
00043 extern  int  *__errnoaddr;
00044 
00045 double  fabs(double);
00046 #pragma intrinsic (fabs)
00047 
00048 int  round(double);
00049 #pragma intrinsic (round)
00050 
00051 #define  ROUND(d)  round(d)
00052 
00053 typedef union
00054 {
00055   struct
00056   {
00057     unsigned int hi;
00058     unsigned int lo;
00059   } word;
00060 
00061   double  d;
00062 } du;
00063 
00064 /* tables used to do argument reduction for args between +/- 16 radians;
00065    the sum of the high and low values of the kth entry is (k - 10)*pi/2
00066 */
00067 
00068 static const du  tblh[] =
00069 {
00070 {D(0xc02f6a7a,  0x2955385e)},
00071 {D(0xc02c463a,  0xbeccb2bb)},
00072 {D(0xc02921fb,  0x54442d18)},
00073 {D(0xc025fdbb,  0xe9bba775)},
00074 {D(0xc022d97c,  0x7f3321d2)},
00075 {D(0xc01f6a7a,  0x2955385e)},
00076 {D(0xc01921fb,  0x54442d18)},
00077 {D(0xc012d97c,  0x7f3321d2)},
00078 {D(0xc00921fb,  0x54442d18)},
00079 {D(0xbff921fb,  0x54442d18)},
00080 {D(0x00000000,  0x00000000)},
00081 {D(0x3ff921fb,  0x54442d18)},
00082 {D(0x400921fb,  0x54442d18)},
00083 {D(0x4012d97c,  0x7f3321d2)},
00084 {D(0x401921fb,  0x54442d18)},
00085 {D(0x401f6a7a,  0x2955385e)},
00086 {D(0x4022d97c,  0x7f3321d2)},
00087 {D(0x4025fdbb,  0xe9bba775)},
00088 {D(0x402921fb,  0x54442d18)},
00089 {D(0x402c463a,  0xbeccb2bb)},
00090 {D(0x402f6a7a,  0x2955385e)},
00091 };
00092 
00093 static const du  tbll[] =
00094 {
00095 {D(0xbcc60faf,  0xbfd97309)},
00096 {D(0xbcc3daea,  0xf976e788)},
00097 {D(0xbcc1a626,  0x33145c07)},
00098 {D(0xbcbee2c2,  0xd963a10c)},
00099 {D(0xbcba7939,  0x4c9e8a0a)},
00100 {D(0xbcb60faf,  0xbfd97309)},
00101 {D(0xbcb1a626,  0x33145c07)},
00102 {D(0xbcaa7939,  0x4c9e8a0a)},
00103 {D(0xbca1a626,  0x33145c07)},
00104 {D(0xbc91a626,  0x33145c07)},
00105 {D(0x00000000,  0x00000000)},
00106 {D(0x3c91a626,  0x33145c07)},
00107 {D(0x3ca1a626,  0x33145c07)},
00108 {D(0x3caa7939,  0x4c9e8a0a)},
00109 {D(0x3cb1a626,  0x33145c07)},
00110 {D(0x3cb60faf,  0xbfd97309)},
00111 {D(0x3cba7939,  0x4c9e8a0a)},
00112 {D(0x3cbee2c2,  0xd963a10c)},
00113 {D(0x3cc1a626,  0x33145c07)},
00114 {D(0x3cc3daea,  0xf976e788)},
00115 {D(0x3cc60faf,  0xbfd97309)},
00116 };
00117 
00118 static const du  rpiby2 =
00119 {D(0x3fe45f30,  0x6dc9c883)};
00120 
00121 static const du  twopm23 =
00122 {D(0x3e800000,  0x00000000)};
00123 
00124 static const du  zero =
00125 {D(0x00000000,  0x00000000)};
00126 
00127 static const du  half =
00128 {D(0x3fe00000,  0x00000000)};
00129 
00130 static const du  one =
00131 {D(0x3ff00000,  0x00000000)};
00132 
00133 static const du  ph =
00134 {D(0x3ff921fb,  0x50000000)};
00135 
00136 static const du  pl =
00137 {D(0x3e5110b4,  0x60000000)};
00138 
00139 static const du  pt =
00140 {D(0x3c91a626,  0x30000000)};
00141 
00142 static const du  pe =
00143 {D(0x3ae8a2e0,  0x30000000)};
00144 
00145 static const du  pe2 =
00146 {D(0x394c1cd1,  0x29024e09)};
00147 
00148 static const du  Pt =
00149 {D(0x3c91a626,  0x33145c07)};
00150 
00151 static const du  piby2low =
00152 {D(0x3e5110b4,  0x611a6263)};
00153 
00154 /* coefficients for polynomial approximation of sin on +/- pi/2 */
00155 
00156 static const du  S[] =
00157 {
00158 {D(0x3ff00000,  0x00000000)},
00159 {D(0xbfc55555,  0x55555548)},
00160 {D(0x3f811111,  0x1110f7d0)},
00161 {D(0xbf2a01a0,  0x19bfdf03)},
00162 {D(0x3ec71de3,  0x567d4896)},
00163 {D(0xbe5ae5e5,  0xa9291691)},
00164 {D(0x3de5d8fd,  0x1fcf0ec1)},
00165 };
00166 
00167 /* coefficients for polynomial approximation of cos on +/- pi/2 */
00168 
00169 static const du  C[] =
00170 {
00171 {D(0x3ff00000,  0x00000000)},
00172 {D(0xbfdfffff,  0xffffff96)},
00173 {D(0x3fa55555,  0x5554f0ab)},
00174 {D(0xbf56c16c,  0x1640aaca)},
00175 {D(0x3efa019f,  0x81cb6a1d)},
00176 {D(0xbe927df4,  0x609cb202)},
00177 {D(0x3e21b8b9,  0x947ab5c8)},
00178 };
00179 
00180 /* ====================================================================
00181  *
00182  * FunctionName    __dcis
00183  *
00184  * Description    computes cos(arg) + i*sin(arg)
00185  *
00186  * ====================================================================
00187  */
00188 
00189 dcomplex __dcis(double x)
00190 {
00191   double  xsq;
00192   double  cospoly, sinpoly;
00193   int  m, n;
00194   dcomplex  result;
00195   double  absx;
00196   double  y, z, dn;
00197   double  dn1, dn2;
00198   double  zz;
00199   int  ix, xpt;
00200 
00201   /* extract exponent of x and 1 bit of mantissa for some quick screening */
00202 
00203   ix = *(int *)(&x);
00204   xpt = (ix >> 19);
00205   xpt &= 0xfff;
00206 
00207   if (xpt < 0x7fd) {
00208     /*   |x| < 0.75   */
00209     n = 0;
00210     if (xpt >= 0x7c2)
00211     {
00212       /*   |x| >= 2^(-30)   */
00213       xsq = x*x;
00214       sinpoly = (((((S[6].d*xsq + S[5].d)*xsq + S[4].d)*xsq +
00215         S[3].d)*xsq + S[2].d)*xsq + S[1].d)*(xsq*x) +
00216         x;
00217       cospoly = (((((C[6].d*xsq + C[5].d)*xsq + C[4].d)*xsq +
00218         C[3].d)*xsq + C[2].d)*xsq + C[1].d)*xsq +
00219         one.d;
00220     } else {
00221       sinpoly = x;
00222       cospoly = one.d;
00223     }
00224 L:
00225     if (n&1) {
00226       if (n&2) {
00227         /*
00228          *  n%4 = 3
00229          *  result is sin(x) - i*cos(x)
00230          */
00231         result.dreal = sinpoly;
00232         result.dimag = -cospoly;
00233       } else {
00234         /*
00235          *  n%4 = 1
00236          *  result is -sin(x) + i*cos(x)
00237          */
00238         result.dreal = -sinpoly;
00239         result.dimag = cospoly;
00240       }
00241       return (result);
00242     }
00243     if (n&2) {
00244       /*
00245        *  n%4 = 2
00246        *  result is -cos(x) - i*sin(x)
00247        */
00248       result.dreal = -cospoly;
00249       result.dimag = -sinpoly;
00250     } else {
00251       /*
00252        *  n%4 = 0
00253        *  result is cos(x) + i*sin(x)
00254        */
00255       result.dreal = cospoly;
00256       result.dimag = sinpoly;
00257     }
00258     return(result);
00259   }
00260 
00261   if (xpt < 0x806) {
00262     /*   |x| < 16.0   */
00263     /*  do a table based argument reduction to +/- pi/4  */
00264     dn = x*rpiby2.d;
00265     n = ROUND(dn);
00266     /*  compute z = x - n*pi/2  */
00267     x = x - tblh[n+10].d;
00268     y = tbll[n+10].d;
00269     z = x - y;
00270     zz = (x - z) - y;
00271     /* z is reduced argument; zz is correction term */
00272 cont:
00273     xsq = z*z;
00274     cospoly = (((((C[6].d*xsq + C[5].d)*xsq + C[4].d)*xsq +
00275         C[3].d)*xsq + C[2].d)*xsq + C[1].d)*xsq -
00276         zz*z + one.d;
00277     sinpoly = (((((S[6].d*xsq + S[5].d)*xsq + S[4].d)*xsq +
00278       S[3].d)*xsq + S[2].d)*xsq + S[1].d)*(xsq*z) + zz + z;
00279     goto L;
00280   }
00281 
00282   if (xpt < 0x836) {
00283     /*  |x| < 2^28  */
00284     dn = x*rpiby2.d;
00285     n = ROUND(dn);
00286     dn = n;
00287     x = x - dn*ph.d;
00288     x = x - dn*pl.d;
00289     if (fabs(x) < twopm23.d)
00290       goto fix;
00291     y = dn*Pt.d;
00292     z = x - y;
00293     zz = (x - z) - y;
00294     goto cont;
00295 fix:
00296     y = dn*pt.d;
00297     z = x - y;
00298     zz = (x - z) - y;
00299     x = z;
00300     y = dn*pe.d;
00301     z = x - y;
00302     zz += (x - z) - y;
00303     x = z;
00304     y = dn*pe2.d;
00305     z = x - y;
00306     zz += (x - z) - y;
00307     goto cont;
00308   }
00309 
00310   if (xpt < 0x862) {
00311     /*  |x| < 2^50  */
00312     absx = fabs(x);
00313     dn = z = absx*rpiby2.d;
00314     /* round dn to the nearest integer */
00315     m = *(int *)&dn;
00316     m >>= 20;
00317     m &= 0x7ff;
00318     /* shift off fractional bits of dn */
00319     n = *((int *)&dn + 1);
00320     n >>= (0x433 - m);
00321     *((int *)&dn + 1) = (n << (0x433 - m));
00322     /* adjust dn and n if the fractional part of dn
00323        was >= 0.5
00324     */
00325     n &= 3;
00326     if ((z - dn) >= half.d) {
00327       dn += one.d;
00328       n += 1;
00329     }
00330     /* split dn into 2 parts */
00331     dn1 = dn;
00332     m = *((int *)&dn1 + 1);
00333     m >>= 27;
00334     m <<= 27;
00335     *((int *)&dn1 + 1) = m;
00336     dn2 = dn - dn1;
00337     z = absx - dn1*ph.d - dn2*ph.d - dn*piby2low.d;
00338     zz = zero.d;
00339     if (x < 0.0) {
00340       z = -z;
00341       n = -n;
00342     }
00343     goto cont;
00344   }
00345 
00346   if (x != x) {
00347     /* x is a NaN; return a pair of quiet NaNs */
00348 #ifdef _IP_NAN_SETS_ERRNO
00349     *__errnoaddr = EDOM;
00350 #endif
00351     result.dreal = __libm_qnan_d;
00352     result.dimag = __libm_qnan_d;
00353     return (result);
00354   }
00355 
00356   /* just give up and return 0.0 */
00357   result.dreal = 0.0;
00358   result.dimag = 0.0;
00359   return (result);
00360 }