Changeset 1023


Ignore:
Timestamp:
Oct 14, 2011, 12:16:22 AM (8 years ago)
Author:
sam
Message:

core: implement log10, sinh and cosh for real numbers.

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r1022 r1023  
    679679}
    680680
    681 static real fast_exp(real const &x)
    682 {
    683     /* This fast exp method is tuned to work on the [0..1] range and
     681real log10(real const &x)
     682{
     683    return log(x) * real::R_LOG10E;
     684}
     685
     686static real fast_exp_sub(real const &x, real const &y)
     687{
     688    /* This fast exp method is tuned to work on the [-1..1] range and
    684689     * no effort whatsoever was made to improve convergence outside this
    685      * domain of validity. */
    686     real ret = real::R_1, fact = real::R_1, xn = x;
     690     * domain of validity. The argument y is used for cases where we
     691     * don't want the leading 1 in the Taylor series. */
     692    real ret = real::R_1 - y, fact = real::R_1, xn = x;
    687693
    688694    for (int i = 1; ; i++)
     
    722728    int e0 = x / real::R_LN2;
    723729    real x0 = x - (real)e0 * real::R_LN2;
    724     real x1 = fast_exp(x0);
     730    real x1 = fast_exp_sub(x0, real::R_0);
    725731    x1.m_signexp += e0;
    726732    return x1;
     
    732738    int e0 = x;
    733739    real x0 = x - (real)e0;
    734     real x1 = fast_exp(x0 * real::R_LN2);
     740    real x1 = fast_exp_sub(x0 * real::R_LN2, real::R_0);
    735741    x1.m_signexp += e0;
    736742    return x1;
     743}
     744
     745real sinh(real const &x)
     746{
     747    /* We cannot always use (exp(x)-exp(-x))/2 because we'll lose
     748     * accuracy near zero. We only use this identity for |x|>0.5. If
     749     * |x|<=0.5, we compute exp(x)-1 and exp(-x)-1 instead. */
     750    bool near_zero = (fabs(x) < real::R_1 >> 1);
     751    real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x);
     752    real x2 = near_zero ? fast_exp_sub(-x, real::R_1) : exp(-x);
     753    return (x1 - x2) >> 1;
     754}
     755
     756real tanh(real const &x)
     757{
     758    /* See sinh() for the strategy here */
     759    bool near_zero = (fabs(x) < real::R_1 >> 1);
     760    real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x);
     761    real x2 = near_zero ? fast_exp_sub(-x, real::R_1) : exp(-x);
     762    real x3 = near_zero ? x1 + x2 + real::R_2 : x1 + x2;
     763    return (x1 - x2) / x3;
     764}
     765
     766real cosh(real const &x)
     767{
     768    /* No need to worry about accuracy here; maybe the last bit is slightly
     769     * off, but that's about it. */
     770    return (exp(x) + exp(-x)) >> 1;
    737771}
    738772
  • trunk/src/real.h

    r1022 r1023  
    7373
    7474    /* Hyperbolic functions */
    75     /* FIXME: sinh cosh tanh */
     75    friend real sinh(real const &x);
     76    friend real cosh(real const &x);
     77    friend real tanh(real const &x);
    7678
    7779    /* Exponential and logarithmic functions */
     
    8082    friend real log(real const &x);
    8183    friend real log2(real const &x);
    82     /* FIXME: frexp ldexp log10 modf */
     84    friend real log10(real const &x);
     85    /* FIXME: frexp ldexp modf */
    8386
    8487    /* Power functions */
Note: See TracChangeset for help on using the changeset viewer.