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 float __libm_qnan_f; 00043 extern int *__errnoaddr; 00044 00045 double fabs(double); 00046 #pragma intrinsic (fabs) 00047 00048 #if _COMPILER_VERSION >= 400 00049 int round(double); 00050 #pragma intrinsic (round) 00051 #define ROUND(d) round(d) 00052 #else 00053 #define ROUND(d) (int)(((d) >= 0.0) ? ((d) + 0.5) : ((d) - 0.5)) 00054 #endif 00055 00056 typedef union 00057 { 00058 struct 00059 { 00060 unsigned int hi; 00061 unsigned int lo; 00062 } word; 00063 00064 double d; 00065 } du; 00066 00067 /* coefficients for polynomial approximation of sin on +/- pi/2 */ 00068 00069 static const du S[] = 00070 { 00071 {D(0x3ff00000, 0x00000000)}, 00072 {D(0xbfc55554, 0x5268a030)}, 00073 {D(0x3f811073, 0xafd14db9)}, 00074 {D(0xbf29943e, 0x0fc79aa9)}, 00075 }; 00076 00077 static const du C[] = 00078 { 00079 {D(0x3ff00000, 0x00000000)}, 00080 {D(0xbfdffffb, 0x2a77e083)}, 00081 {D(0x3fa553e7, 0xf02ac8aa)}, 00082 {D(0xbf5644d6, 0x2993c4ad)}, 00083 }; 00084 00085 static const du rpiby2 = 00086 {D(0x3fe45f30, 0x6dc9c883)}; 00087 00088 static const du piby2hi = 00089 {D(0x3ff921fb, 0x50000000)}; 00090 00091 static const du piby2lo = 00092 {D(0x3e5110b4, 0x611a6263)}; 00093 00094 static const du zero = 00095 {D(0x00000000, 0x00000000)}; 00096 00097 static const du half = 00098 {D(0x3fe00000, 0x00000000)}; 00099 00100 static const du one = 00101 {D(0x3ff00000, 0x00000000)}; 00102 00103 static const du twopm12 = 00104 {D(0x3f300000, 0x00000000)}; 00105 00106 /* ==================================================================== 00107 * 00108 * FunctionName __rcis 00109 * 00110 * Description computes cos(arg) + i*sin(arg) 00111 * 00112 * ==================================================================== 00113 */ 00114 00115 complex __rcis(float x) 00116 { 00117 int m, n; 00118 int ix, xpt; 00119 double dx, xsq; 00120 double dn, dn1, dn2; 00121 double z, absdx; 00122 double sinpoly, cospoly; 00123 complex result; 00124 00125 ix = *(int *)(&x); 00126 xpt = (ix >> 22); 00127 xpt &= 0x1ff; 00128 00129 /* xpt is exponent(x) + 1 bit of mantissa */ 00130 00131 if (xpt < 0xfd) { 00132 /* |x| < .75 */ 00133 dx = x; 00134 n = 0; 00135 L: 00136 if (fabs(dx) >= twopm12.d) { 00137 /* |dx| >= 2^(-12) */ 00138 xsq = dx*dx; 00139 sinpoly = ((S[3].d*xsq + S[2].d)*xsq + S[1].d)*(xsq*dx) + dx; 00140 cospoly = ((C[3].d*xsq + C[2].d)*xsq + C[1].d)*xsq + one.d; 00141 } else { 00142 sinpoly = dx; 00143 cospoly = one.d; 00144 } 00145 00146 if (n&1) { 00147 if (n&2) { 00148 /* 00149 * n%4 = 3 00150 * result is sin(x) - i*cos(x) 00151 */ 00152 result.real = sinpoly; 00153 result.imag = -cospoly; 00154 } else { 00155 /* 00156 * n%4 = 1 00157 * result is -sin(x) + i*cos(x) 00158 */ 00159 result.real = -sinpoly; 00160 result.imag = cospoly; 00161 } 00162 return (result); 00163 } 00164 00165 if (n&2) { 00166 /* 00167 * n%4 = 2 00168 * result is -cos(x) - i*sin(x) 00169 */ 00170 result.real = -cospoly; 00171 result.imag = -sinpoly; 00172 } else { 00173 /* 00174 * n%4 = 0 00175 * result is cos(x) + i*sin(x) 00176 */ 00177 result.real = cospoly; 00178 result.imag = sinpoly; 00179 } 00180 00181 return(result); 00182 } 00183 00184 if (xpt < 0x136) { 00185 /* |x| < 2^28 */ 00186 dx = x; 00187 dn = dx*rpiby2.d; 00188 n = ROUND(dn); 00189 dn = n; 00190 dx = dx - dn*piby2hi.d; 00191 dx = dx - dn*piby2lo.d; /* dx = x - n*piby2 */ 00192 goto L; 00193 } 00194 00195 if (xpt < 0x162) { 00196 /* |x| < 2^50 */ 00197 dx = x; 00198 absdx = fabs(dx); 00199 dn = z = absdx*rpiby2.d; 00200 /* round dn to the nearest integer */ 00201 m = *(int *)&dn; 00202 m >>= 20; 00203 m &= 0x7ff; 00204 /* shift off fractional bits of dn */ 00205 n = *((int *)&dn + 1); 00206 n >>= (0x433 - m); 00207 *((int *)&dn + 1) = (n << (0x433 - m)); 00208 n &= 3; 00209 /* adjust dn and n if the fractional part of dn 00210 was >= 0.5 00211 */ 00212 if ((z - dn) >= half.d) { 00213 dn += one.d; 00214 n += 1; 00215 } 00216 /* split dn into 2 parts */ 00217 dn1 = dn; 00218 m = *((int *)&dn1 + 1); 00219 m >>= 27; 00220 m <<= 27; 00221 *((int *)&dn1 + 1) = m; 00222 dn2 = dn - dn1; 00223 z = absdx - dn1*piby2hi.d - dn2*piby2hi.d - dn*piby2lo.d; 00224 if (dx < zero.d) { 00225 z = -z; 00226 n = -n; 00227 } 00228 dx = z; 00229 goto L; 00230 } 00231 00232 if (x != x) { 00233 /* x is a NaN; return a pair of quiet NaNs */ 00234 #ifdef _IP_NAN_SETS_ERRNO 00235 *__errnoaddr = EDOM; 00236 #endif 00237 result.real = __libm_qnan_f; 00238 result.imag = __libm_qnan_f; 00239 return (result); 00240 } 00241 /* just give up and return 0.0 */ 00242 result.real = (float)0.0; 00243 result.imag = (float)0.0; 00244 return (result); 00245 }