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  * Fast log() Function, by Laurent de Soras:
  13  * http://www.flipcode.com/cgi-bin/msg.cgi?showThread=Tip-Fastlogfunction&forum=totd&id=-1
  14  *
  15  */
  16
  17 /*
  18  * $Id$
  19  */
  20
  21
  22
  23 #include <fastmath.hxx>
  24
  25 /**
  26  * This function is on avarage 9 times faster than the system exp() function
  27  * and has an error of about 1.5%
  28  */
  29 static union {
  30     double d;
  31     struct {
  32 #if BYTE_ORDER == BIG_ENDIAN
  33         int i, j;
  34 #else
  35         int j, i;
  36 #endif
  37     } n;
  38 } _eco;
  39
  40 double fast_exp(double val) {
  41     const double a = 1048576/M_LN2;
  42     const double b_c = 1072632447; /* 1072693248 - 60801 */
  43
  44     _eco.n.i = a*val + b_c;
  45
  46     return _eco.d;
  47 }
  48
  49
  50 double fast_log2 (double val)
  51 {
  52    int * const    exp_ptr = reinterpret_cast <int *> (&val);
  53    int            x = *exp_ptr;
  54    const int      log_2 = ((x >> 23) & 255) - 128;
  55    x &= ~(255 << 23);
  56    x += 127 << 23;
  57    *exp_ptr = x;
  58
  59    val = ((-1.0f/3) * val + 2) * val - 2.0f/3;   // (1)
  60
  61    return (val + log_2);
  62 }
  63
  64 /**
  65  * This function is about 3 times faster than the system log() function
  66  * and has an error of about 0.01%
  67  */
  68 double fast_log (double val)
  69 {
  70    return (fast_log2 (val) * 0.69314718f);
  71 }
  72
  73 double fast_log10 (double val)
  74 {
  75    return (fast_log (val) / fast_log (10));
  76 }
  77
  78
  79 /**
  80  * While we're on the subject, someone might have use for these as well?
  81  * Float Shift Left and Float Shift Right. Do what you want with this.
  82  */
  83 void fast_BSL(double &x, register unsigned long shiftAmount) {
  84
  85         *(unsigned long*)&x+=shiftAmount<<23;
  86
  87 }
  88
  89 void fast_BSR(double &x, register unsigned long shiftAmount) {
  90
  91         *(unsigned long*)&x-=shiftAmount<<23;
  92
  93 }
  94