Changeset 992


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

core: implement log() for real numbers, and start documenting our next
improved implementation of exp(), which relies on log().

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r990 r992  
    117117        return *this - (-x);
    118118
    119     /* Ensure *this is the larger exponent (no need to be strictly larger,
    120      * as in subtraction). Otherwise, switch. */
     119    /* Ensure *this has the larger exponent (no need for the mantissa to
     120     * be larger, as in subtraction). Otherwise, switch. */
    121121    if ((m_signexp << 1) < (x.m_signexp << 1))
    122122        return x + *this;
     
    515515}
    516516
     517static real fastlog(real const &x)
     518{
     519    /* This fast log method is tuned to work on the [1..2] range and
     520     * no effort whatsoever was made to improve convergence outside this
     521     * domain of validity. It can converge pretty fast, provided we use
     522     * the following variable substitutions:
     523     *    y = sqrt(x)
     524     *    z = (y - 1) / (y + 1)
     525     *
     526     * And the following identities:
     527     *    ln(x) = 2 ln(y)
     528     *          = 2 ln((1 + z) / (1 - z))
     529     *          = 4 z (1 + z^2 / 3 + z^4 / 5 + z^6 / 7...)
     530     *
     531     * Any additional sqrt() call would halve the convergence time, but
     532     * would also impact the final precision. For now we stick with one
     533     * sqrt() call. */
     534    real y = sqrt(x);
     535    real z = (y - (real)1) / (y + (real)1), z2 = z * z, zn = z2;
     536    real sum = 1.0;
     537
     538    for (int i = 3; i < 200; i += 2)
     539    {
     540        sum += zn / (real)i;
     541        zn *= z2;
     542    }
     543
     544    /* FIXME: sum.m_signexp += 2; (but needs private data access) */
     545    sum *= (real)4;
     546    return z * sum;
     547}
     548
     549static real LOG_2 = fastlog((real)2);
     550
     551real log(real const &x)
     552{
     553    /* Strategy for log(x): if x = 2^E*M then log(x) = E log(2) + log(M),
     554     * with the property that M is in [1..2[, so fastlog() applies here. */
     555    real tmp = x;
     556    if (x.m_signexp >> 31 || x.m_signexp == 0)
     557    {
     558        tmp.m_signexp = 0xffffffffu;
     559        tmp.m_mantissa[0] = 0xffffu;
     560        return tmp;
     561    }
     562    tmp.m_signexp = (1 << 30) - 1;
     563    return (real)(x.m_signexp - (1 << 30) + 1) * LOG_2 + fastlog(tmp);
     564}
     565
    517566real exp(real const &x)
    518567{
     568    /* Strategy for exp(x): the Taylor series does not converge very fast
     569     * with large positive or negative values.
     570     *
     571     * However, we know that the result is going to be in the form M*2^E,
     572     * where M is the mantissa and E the exponent. We first try to predict
     573     * a value for E, which is approximately log2(exp(x)) = x / log(2).
     574     *
     575     * 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).
     577     *
     578     * Thus the final algorithm:
     579     *  int E0 = x / log(2)
     580     *  real x0 = x - log(E0)
     581     *  real x1 = exp(x0)
     582     *  return x1 * 2^E0
     583     */
    519584    int square = 0;
    520585
  • trunk/src/real.h

    r990 r992  
    5858    friend real fres(real const &x);
    5959    friend real sqrt(real const &x);
     60    friend real log(real const &x);
    6061    friend real exp(real const &x);
    6162
Note: See TracChangeset for help on using the changeset viewer.