00001 /* 00002 * Copyright (C) 1999,2000,2001 FSF 00003 * 00004 * This program is free software; you can redistribute it and/or modify 00005 * it under the terms of the GNU General Public License as published by 00006 * the Free Software Foundation; either version 2 of the License, or 00007 * (at your option) any later version. 00008 * 00009 * This program is distributed in the hope that it will be useful, 00010 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00012 * GNU General Public License for more details. 00013 * 00014 * You should have received a copy of the GNU General Public License 00015 * along with this program; if not, write to the Free Software 00016 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. 00017 */ 00018 00019 #include "mathvalues.h" 00020 00021 /****************************************************************************/ 00022 /* SIN() - sine */ 00023 /* */ 00024 /* Based on the algorithm from "Software Manual for the Elementary */ 00025 /* Functions", Cody and Waite, Prentice Hall 1980, chapter 8. */ 00026 /* */ 00027 /* N = round(x / PI) */ 00028 /* f = x - N * PI */ 00029 /* g = f * f */ 00030 /* R = polynomial expansion */ 00031 /* */ 00032 /* result = f + f * R */ 00033 /* */ 00034 /* if x < 0, result = - result */ 00035 /* if N is even, result = - result */ 00036 /* */ 00037 /* This will return the wrong result for x >= MAXINT * PI */ 00038 /****************************************************************************/ 00039 00040 double sin(double x) 00041 { 00042 double xn, f, g, rg; 00043 double sgn = (x < 0) ? -1.0 : 1.0; 00044 int n; 00045 00046 x = fabs(x); 00047 n = (int) ((x * INVSPI) + 0.5); 00048 xn = (double) n; 00049 00050 /*************************************************************************/ 00051 /* if n is odd, negate the sign */ 00052 /*************************************************************************/ 00053 if (n % 2) sgn = -sgn; 00054 00055 /*************************************************************************/ 00056 /* f = x - xn * PI (but mathematically more stable) */ 00057 /*************************************************************************/ 00058 f = (x - xn * C1) - xn * C2; 00059 00060 /*************************************************************************/ 00061 /* determine polynomial expression */ 00062 /*************************************************************************/ 00063 g = f * f; 00064 00065 #if BITS<=24 00066 rg = (((R4 * g + R3) * g + R2) * g + R1) * g; 00067 #elif BITS>=25 && BITS<=32 00068 rg = ((((R5 * g + R4) * g + R3) * g + R2) * g + R1) * g; 00069 #elif BITS>=33 && BITS<=50 00070 rg = ((((((R7*g+R6)*g+R5)*g+R4)*g+R3)*g+R2)*g+R1)*g; 00071 #else 00072 rg = (((((((R8*g+R7)*g+R6)*g+R5)*g+R4)*g+R3)*g+R2)*g+R1)*g; 00073 #endif 00074 00075 return (sgn * (f + f * rg)); 00076 } 00077 00078 00079 /****************************************************************************/ 00080 /* COS() - Cosine */ 00081 /* */ 00082 /* Based on the algorithm from "Software Manual for the Elementary */ 00083 /* Functions", Cody and Waite, Prentice Hall 1980, chapter 8. */ 00084 /* */ 00085 /* N = round(x / PI + 1/2) - 0.5 */ 00086 /* f = x - N * PI */ 00087 /* g = f * f */ 00088 /* R = polynomial expression */ 00089 /* */ 00090 /* result = f + f * R */ 00091 /* if N is even, result = - result */ 00092 /* */ 00093 /* This will return the wrong result for x >= MAXINT * PI */ 00094 /****************************************************************************/ 00095 double cos(double x) 00096 { 00097 double sgn; /* the sign of the result */ 00098 double xn, f, g, rg; 00099 int n; 00100 00101 /**************************************************************************/ 00102 /* cos(x) = cos(-x) */ 00103 /**************************************************************************/ 00104 x = fabs(x); 00105 00106 /**************************************************************************/ 00107 /* n = round(x/PI + 1/2) (can be rounded this way, since positive number) */ 00108 /**************************************************************************/ 00109 n = (int) (((x + HALFPI) * INVSPI) + 0.5); 00110 xn = (double) n - 0.5; 00111 00112 /**************************************************************************/ 00113 /* if n is odd, negate the sign */ 00114 /**************************************************************************/ 00115 sgn = (n % 2) ? -1.0 : 1.0; 00116 00117 /**************************************************************************/ 00118 /* f = x - xn * PI (but more mathematically stable) */ 00119 /**************************************************************************/ 00120 f = (x - xn * C1) - xn * C2; 00121 00122 /**************************************************************************/ 00123 /* determine polynomial expression */ 00124 /**************************************************************************/ 00125 g = f * f; 00126 00127 #if BITS<=24 00128 rg = (((R4 * g + R3) * g + R2) * g + R1) * g; 00129 #elif BITS>=25 && BITS<=32 00130 rg = ((((R5 * g + R4) * g + R3) * g + R2) * g + R1) * g; 00131 #elif BITS>=33 && BITS<=50 00132 rg = ((((((R7*g+R6)*g+R5)*g+R4)*g+R3)*g+R2)*g+R1)*g; 00133 #else 00134 rg = (((((((R8*g+R7)*g+R6)*g+R5)*g+R4)*g+R3)*g+R2)*g+R1)*g; 00135 #endif 00136 00137 return (sgn * (f + f * rg)); 00138 } 00139 00140 00141 00142 /****************************************************************************/ 00143 /* ATAN() - Arctangent */ 00144 /* */ 00145 /* Based on the algorithm from "Software Manual for the Elementary */ 00146 /* Functions", Cody and Waite, Prentice Hall 1980, chapter 11. */ 00147 /* */ 00148 /* if x > 1, x = 1 / x */ 00149 /* if x > 2 - sqrt(3), x = (x * sqrt(3) - 1) / (sqrt(3) + x) */ 00150 /* g = x * x */ 00151 /* R = polynomial expression */ 00152 /* */ 00153 /* result = (t * (x + x * R) + an) * s */ 00154 /****************************************************************************/ 00155 double atan(double x) 00156 { 00157 double g, p, q; 00158 double s = (x < 0.0) ? -1.0F : 1.0F; /* sign */ 00159 double t = 1.0; 00160 int n = 0; 00161 00162 static double a[4] = {0.0, 0.52359877559829887308, 1.57079632679489661923, 00163 1.04719755119659774615}; 00164 00165 if ((x = fabs(x)) > 1.0) 00166 { 00167 x = 1.0 / x; 00168 n = 2; 00169 00170 /******************************************************************/ 00171 /* the partial result needs to be negated */ 00172 /******************************************************************/ 00173 t = -1.0; 00174 } 00175 00176 /**********************************************************************/ 00177 /* for x > (2 - sqrt(3) ) */ 00178 /**********************************************************************/ 00179 if (x > TWO_SQRT3) 00180 { 00181 /******************************************************************/ 00182 /* x = (x * sqrt(3) -1) / (sqroot(3) + x) */ 00183 /******************************************************************/ 00184 x = (x * SQRTTHREE - 1.0) / (SQRTTHREE + x); 00185 ++n; 00186 } 00187 00188 /*********************************************************************/ 00189 /* determine polynomial expression */ 00190 /*********************************************************************/ 00191 g = x * x; 00192 00193 #if BITS<=24 00194 p = (ATP1 * g + ATP0) * g; 00195 q = g + ATQ0; 00196 #elif BITS>=25 && BITS<=32 00197 p = (ATP1 * g + ATP0) * g; 00198 q = (g + ATQ1) * g + ATQ0; 00199 #elif BITS>=33 && BITS<=50 00200 p = ((ATP2 * g + ATP1) * g + ATP0) * g; 00201 q = ((g + ATQ2) * g + ATQ1) * g + ATQ0; 00202 #else 00203 p = (((ATP3 * g + ATP2) * g + ATP1) * g + ATP0) * g; 00204 q = (((g + ATQ3) * g + ATQ2) * g + ATQ1) * g + ATQ0; 00205 #endif 00206 00207 /**********************************************************************/ 00208 /* calculate the result multiplied by the correct sign */ 00209 /**********************************************************************/ 00210 return ((((p / q) * x + x) * t + a[n]) * s); 00211 } 00212 00213 00214 /****************************************************************************/ 00215 /* FMOD() - Doubleing point remainder */ 00216 /* */ 00217 /* Returns the remainder after dividing x by y an integral number of times.*/ 00218 /* */ 00219 /****************************************************************************/ 00220 double fmod(double x, double y) 00221 { 00222 double d = fabs(x); 00223 double e = fabs(y); 00224 00225 /*************************************************************************/ 00226 /* if y is too small or y == x, any remainder is negligible. */ 00227 /*************************************************************************/ 00228 if (d - e == d || d == e) return (0.0); 00229 00230 /*************************************************************************/ 00231 /* if x and y are integers, just do a %. */ 00232 /*************************************************************************/ 00233 if (((double)(int)d == d) && ((double)(int)e == e)) 00234 return ((double)((int)x % (int)y)); 00235 00236 /*************************************************************************/ 00237 /* otherwise, divide; result = dividend - (quotient * divisor) */ 00238 /*************************************************************************/ 00239 /* modf(x/y, &d); */ 00240 return (x - (double)((int)(x/y)) * y); 00241 } 00242 00243 00244 00245 /****************************************************************************/ 00246 /* EXP() - e ^ x */ 00247 /* */ 00248 /* Based on the algorithm from "Software Manual for the Elementary */ 00249 /* Functions", Cody and Waite, Prentice Hall 1980, chapter 6. */ 00250 /* */ 00251 /* N = round(x / ln(2)) */ 00252 /* g = x - N * ln(2) */ 00253 /* z = g * g */ 00254 /* */ 00255 /* R = polynomial expansion */ 00256 /* */ 00257 /* result = R * 2 ^ (N + 1) */ 00258 /****************************************************************************/ 00259 double exp(double x) 00260 { 00261 double g, z, q, p; 00262 int n; 00263 00264 /*************************************************************************/ 00265 /* check if input would produce output out of the range of this function */ 00266 /*************************************************************************/ 00267 if (x > MAXX) { return (HUGE_VAL); } 00268 00269 if (x < 0) n = (int) (x * INVLOGe2 - 0.5); /* since (int) -1.5 = -1.0 */ 00270 else n = (int) (x * INVLOGe2 + 0.5); 00271 00272 /*************************************************************************/ 00273 /* g = x - n * ln(2) (but more mathematically stable) */ 00274 /*************************************************************************/ 00275 g = (x - n * LNC3) - n * LNC4; 00276 00277 /*************************************************************************/ 00278 /* determine polynomial expression */ 00279 /*************************************************************************/ 00280 z = g * g; 00281 00282 #if BITS <=29 00283 p = (EXP1 * z + EXP0) * g; 00284 q = EXQ1 * z + EXQ0; 00285 #elif BITS>=30 && BITS<=42 00286 p = (EXP1 * z + EXP0) * g; 00287 q = (EXQ2 * z + EXQ1) * z + EXQ0; 00288 #elif BITS>=43 && BITS<=56 00289 p = ((EXP2 * z + EXP1) * z + EXP0) * g; 00290 q = (EXQ2 * z + EXQ1) * z + EXQ0; 00291 #else 00292 p = ((EXP2 * z + EXP1) * z + EXP0) * g; 00293 q = ((EXQ3 * z + EXQ2) * z + EXQ1) * z + EXQ0; 00294 #endif 00295 00296 /*************************************************************************/ 00297 /* exp(x) = exp(g) * 2 ^ (n + 1) */ 00298 /*************************************************************************/ 00299 /* return ldexp(0.5 + p / (q - p), n + 1); */ 00300 return (0.5 + p / (q - p)) * (1 << ( n + 1)); 00301 } 00302 00303 double sqrt(double x) { 00304 #define DBITS2 (BITS/2+1) 00305 double a=x; 00306 double y=x; 00307 double diff; 00308 int i; 00309 // b = 64; 00310 // a = 2^ceil(b/2); //b is the wordlength you are using 00311 // y = a; //a is the first estimation value for sqrt(x) 00312 // for (i=0; i<ceil(b/2); i++) //each iteration achieves one bit 00313 for (i=0; i<DBITS2; i++) //each iteration achieves one bit 00314 { //of accuracy 00315 a *= 0.5; 00316 diff = x - y*y; 00317 if (diff > 0) 00318 { 00319 y = y + a; 00320 } 00321 else if (diff < 0) 00322 { 00323 y = y - a; 00324 } 00325 } 00326 return y; 00327 } 00328 00329 #if 0 // did not work for PPC 00330 00331 #define itable ((double *)xtable) 00332 00333 static int xtable[16] = { 00334 0x540bcb0d,0x3fe56936, 0x415a86d3,0x3fe35800, 0xd9ac3519,0x3fe1c80d, 00335 0x34f91569,0x3fe08be9, 0x8f3386d8,0x3fee4794, 0x9ea02719,0x3feb5b28, 00336 0xe4ff9edc,0x3fe92589, 0x1c52539d,0x3fe76672 }; 00337 00338 static int norm2(double *t) { 00339 unsigned int e,f,g; 00340 f = ((((unsigned int *)t)[1])>>1); 00341 e = ((unsigned int *)t)[1]; 00342 f += 0x1FF80000; 00343 g = e&0x000FFFFF; 00344 f &= 0xFFF00000; 00345 ((int *)t)[1] = g+0x40000000-(e&0x00100000); 00346 00347 return f; 00348 } 00349 00350 double sqrt(double y) { 00351 double a; 00352 int e,c; 00353 00354 e = norm2(&y); 00355 c = (((int *)&y)[1])>>(18)&(7); 00356 a = itable[c]; 00357 00358 for(c=0;c<6;c++) a = 0.5*a*(3.0-y*a*a); 00359 a*=y; 00360 00361 ((int *)&a)[1] &= 0x000FFFFF; 00362 ((int *)&a)[1] |= e; 00363 00364 return a; 00365 } 00366 #endif