Open64 (mfef90, whirl2f, and IR tools)
TAG: version-openad; SVN changeset: 916
|
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 }