simgear/math/fastmath.cxx

   1 /*
   2  * \file fastmath.cxx
   3  * fast mathematics routines.
   4  *
   5  * Refferences:
   6  *
   7  * A Fast, Compact Approximation of the Exponential Function
   8  * Nicol N. Schraudolph
   9  * IDSIA, Lugano, Switzerland
  10  * http://www.inf.ethz.ch/~schraudo/pubs/exp.pdf
  11  *
  12  * Base-2 exp, Laurent de Soras
  13  * http://www.musicdsp.org/archive.php?classid=5#106
  14  *
  15  * Fast log() Function, by Laurent de Soras:
  16  * http://www.flipcode.com/cgi-bin/msg.cgi?showThread=Tip-Fastlogfunction&forum=totd&id=-1
  17  *
  18  * Sin, Cos, Tan approximation
  19  * http://www.musicdsp.org/showArchiveComment.php?ArchiveID=115
  20  *
  21  * fast floating point power computation:
  22  * http://playstation2-linux.com/download/adam/power.c
  23  */
  24
  25 /*
  26  * $Id$
  27  */
  28
  29
  30 #include <simgear/constants.h>
  31
  32 #include "fastmath.hxx"
  33
  34 /**
  35  * This function is on avarage 9 times faster than the system exp() function
  36  * and has an error of about 1.5%
  37  */
  38 static union {
  39     double d;
  40     struct {
  41 #if BYTE_ORDER == BIG_ENDIAN
  42         int i, j;
  43 #else
  44         int j, i;
  45 #endif
  46     } n;
  47 } _eco;
  48
  49 double fast_exp(double val) {
  50     const double a = 1048576/M_LN2;
  51     const double b_c = 1072632447; /* 1072693248 - 60801 */
  52
  53     _eco.n.i = (int)(a*val + b_c);
  54
  55     return _eco.d;
  56 }
  57
  58 /*
  59  * Linear approx. between 2 integer values of val. Uses 32-bit integers.
  60  * Not very efficient but faster than exp()
  61  */
  62 double fast_exp2( const double val )
  63 {
  64     int e;
  65     double ret;
  66
  67     if (val >= 0) {
  68         e = int (val);
  69         ret = val - (e - 1);
  70
  71 #if BYTE_ORDER == BIG_ENDIAN
  72         ((*((int *) &ret)) &= ~(2047 << 20)) += (e + 1023) << 20;
  73 #else
  74         ((*(1 + (int *) &ret)) &= ~(2047 << 20)) += (e + 1023) << 20;
  75 #endif
  76     } else {
  77         e = int (val + 1023);
  78         ret = val - (e - 1024);
  79
  80 #if BYTE_ORDER == BIG_ENDIAN
  81         ((*((int *) &ret)) &= ~(2047 << 20)) += e << 20;
  82 #else
  83         ((*(1 + (int *) &ret)) &= ~(2047 << 20)) += e << 20;
  84 #endif
  85     }
  86
  87     return ret;
  88 }
  89
  90
  91 /*
  92  *
  93  */
  94 float _fast_log2(const float val)
  95 {
  96     float result, tmp;
  97     float mp = 0.346607f;
  98
  99     result = *(int*)&val;
 100     result *= 1.0/(1<<23);
 101     result = result - 127;
 102
 103     tmp = result - floor(result);
 104     tmp = (tmp - tmp*tmp) * mp;
 105     return tmp + result;
 106 }
 107
 108 float _fast_pow2(const float val)
 109 {
 110     float result;
 111
 112     float mp = 0.33971f;
 113     float tmp = val - floor(val);
 114     tmp = (tmp - tmp*tmp) * mp;
 115
 116     result = val + 127 - tmp;
 117     result *= (1<<23);
 118     *(int*)&result = (int)result;
 119     return result;
 120 }
 121
 122
 123
 124 /**
 125  * While we're on the subject, someone might have use for these as well?
 126  * Float Shift Left and Float Shift Right. Do what you want with this.
 127  */
 128 void fast_BSL(float &x, register unsigned long shiftAmount) {
 129
 130         *(unsigned long*)&x+=shiftAmount<<23;
 131
 132 }
 133
 134 void fast_BSR(float &x, register unsigned long shiftAmount) {
 135
 136         *(unsigned long*)&x-=shiftAmount<<23;
 137
 138 }
 139
 140
 141 /*
 142  * fastpow(f,n) gives a rather *rough* estimate of a float number f to the
 143  * power of an integer number n (y=f^n). It is fast but result can be quite a
 144  * bit off, since we directly mess with the floating point exponent.
 145  *
 146  * Use it only for getting rough estimates of the values and where precision
 147  * is not that important.
 148  */
 149 float fast_pow(const float f, const int n)
 150 {
 151     long *lp,l;
 152     lp=(long*)(&f);
 153     l=*lp;l-=0x3F800000l;l<<=(n-1);l+=0x3F800000l;
 154     *lp=l;
 155     return f;
 156 }
 157
 158 float fast_root(const float f, const int n)
 159 {
 160     long *lp,l;
 161     lp=(long*)(&f);
 162     l=*lp;l-=0x3F800000l;l>>=(n-1);l+=0x3F800000l;
 163     *lp=l;
 164     return f;
 165 }
 166
 167
 168 /*
 169  * Code for approximation of cos, sin, tan and inv sin, etc.
 170  * Surprisingly accurate and very usable.
 171  *
 172  * Domain:
 173  * Sin/Cos    [0, pi/2]
 174  * Tan        [0,pi/4]
 175  * InvSin/Cos [0, 1]
 176  * InvTan     [-1, 1]
 177  */
 178
 179 float fast_sin(const float val)
 180 {
 181     float fASqr = val*val;
 182     float fResult = -2.39e-08f;
 183     fResult *= fASqr;
 184     fResult += 2.7526e-06f;
 185     fResult *= fASqr;
 186     fResult -= 1.98409e-04f;
 187     fResult *= fASqr;
 188     fResult += 8.3333315e-03f;
 189     fResult *= fASqr;
 190     fResult -= 1.666666664e-01f;
 191     fResult *= fASqr;
 192     fResult += 1.0f;
 193     fResult *= val;
 194
 195     return fResult;
 196 }
 197
 198 float fast_cos(const float val)
 199 {
 200     float fASqr = val*val;
 201     float fResult = -2.605e-07f;
 202     fResult *= fASqr;
 203     fResult += 2.47609e-05f;
 204     fResult *= fASqr;
 205     fResult -= 1.3888397e-03f;
 206     fResult *= fASqr;
 207     fResult += 4.16666418e-02f;
 208     fResult *= fASqr;
 209     fResult -= 4.999999963e-01f;
 210     fResult *= fASqr;
 211     fResult += 1.0f;
 212
 213     return fResult;
 214 }
 215
 216 float fast_tan(const float val)
 217 {
 218     float fASqr = val*val;
 219     float fResult = 9.5168091e-03f;
 220     fResult *= fASqr;
 221     fResult += 2.900525e-03f;
 222     fResult *= fASqr;
 223     fResult += 2.45650893e-02f;
 224     fResult *= fASqr;
 225     fResult += 5.33740603e-02f;
 226     fResult *= fASqr;
 227     fResult += 1.333923995e-01f;
 228     fResult *= fASqr;
 229     fResult += 3.333314036e-01f;
 230     fResult *= fASqr;
 231     fResult += 1.0f;
 232     fResult *= val;
 233
 234     return fResult;
 235
 236 }
 237
 238 float fast_asin(float val)
 239 {
 240     float fRoot = sqrt(1.0f-val);
 241     float fResult = -0.0187293f;
 242     fResult *= val;
 243     fResult += 0.0742610f;
 244     fResult *= val;
 245     fResult -= 0.2121144f;
 246     fResult *= val;
 247     fResult += 1.5707288f;
 248     fResult = SGD_PI_2 - fRoot*fResult;
 249
 250     return fResult;
 251 }
 252
 253 float fast_acos(float val)
 254 {
 255     float fRoot = sqrt(1.0f-val);
 256     float fResult = -0.0187293f;
 257     fResult *= val;
 258     fResult += 0.0742610f;
 259     fResult *= val;
 260     fResult -= 0.2121144f;
 261     fResult *= val;
 262     fResult += 1.5707288f;
 263     fResult *= fRoot;
 264
 265     return fResult;
 266 }
 267
 268 float fast_atan(float val)
 269 {
 270     float fVSqr = val*val;
 271     float fResult = 0.0028662257f;
 272     fResult *= fVSqr;
 273     fResult -= 0.0161657367f;
 274     fResult *= fVSqr;
 275     fResult += 0.0429096138f;
 276     fResult *= fVSqr;
 277     fResult -= 0.0752896400f;
 278     fResult *= fVSqr;
 279     fResult += 0.1065626393f;
 280     fResult *= fVSqr;
 281     fResult -= 0.1420889944f;
 282     fResult *= fVSqr;
 283     fResult += 0.1999355085f;
 284     fResult *= fVSqr;
 285     fResult -= 0.3333314528f;
 286     fResult *= fVSqr;
 287     fResult += 1.0f;
 288     fResult *= val;
 289
 290     return fResult;
 291 }