Ignore:
Timestamp:
Sep 28, 2011, 7:04:37 PM (11 years ago)
Author:
sam
Message:

core: improve exp() on reals: faster (constant time) and a lot more
accurate.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r992 r993  
    574574     *
    575575     * Let E0 be an integer close to x / log(2). We need to find a value x0
    576      * such that exp(x) = 2^E0 * exp(x0). We get x0 = x - log(E0).
     576     * such that exp(x) = 2^E0 * exp(x0). We get x0 = x - E0 log(2).
    577577     *
    578578     * Thus the final algorithm:
    579579     *  int E0 = x / log(2)
    580      *  real x0 = x - log(E0)
     580     *  real x0 = x - E0 log(2)
    581581     *  real x1 = exp(x0)
    582582     *  return x1 * 2^E0
    583583     */
    584     int square = 0;
    585 
    586     /* FIXME: this is slow. Find a better way to approximate exp(x) for
    587      * large values. */
    588     real tmp = x, one = 1.0;
    589     while (tmp > one)
    590     {
    591         tmp.m_signexp--;
    592         square++;
    593     }
    594 
    595     real ret = 1.0, fact = 1.0, xn = tmp;
     584    int e0 = x / LOG_2;
     585    real x0 = x - (real)e0 * LOG_2;
     586    real x1 = 1.0, fact = 1.0, xn = x0;
    596587
    597588    for (int i = 1; i < 100; i++)
    598589    {
    599590        fact *= (real)i;
    600         ret += xn / fact;
    601         xn *= tmp;
    602     }
    603 
    604     for (int i = 0; i < square; i++)
    605         ret = ret * ret;
    606 
    607     return ret;
     591        x1 += xn / fact;
     592        xn *= x0;
     593    }
     594
     595    x1.m_signexp += e0;
     596    return x1;
    608597}
    609598
Note: See TracChangeset for help on using the changeset viewer.