mathemu.c

Go to the documentation of this file.
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

Generated on Sat Apr 1 09:04:16 2006 for GXSM by  doxygen 1.4.6