Changeset 1026


Ignore:
Timestamp:
Oct 15, 2011, 2:03:35 PM (9 years ago)
Author:
sam
Message:

core: implement the gamma function for reals using Spouge's formula.

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r1024 r1026  
    606606        return real::R_0;
    607607    }
     608}
     609
     610real gamma(real const &x)
     611{
     612    /* We use Spouge's formula. FIXME: precision is far from acceptable,
     613     * especially with large values. We need to compute this with higher
     614     * precision values in order to attain the desired accuracy. It might
     615     * also be useful to sort the ck values by decreasing absolute value
     616     * and do the addition in this order. */
     617    int a = ceilf(logf(2) / logf(2 * M_PI) * real::BIGITS * real::BIGIT_BITS);
     618
     619    real ret = sqrt(real::R_PI << 1);
     620    real fact_k_1 = real::R_1;
     621
     622    for (int k = 1; k < a; k++)
     623    {
     624        real a_k = (real)(a - k);
     625        real ck = pow(a_k, (real)((float)k - 0.5)) * exp(a_k)
     626                / (fact_k_1 * (x + (real)(k - 1)));
     627        ret += ck;
     628        fact_k_1 *= (real)-k;
     629    }
     630
     631    ret *= pow(x + (real)(a - 1), x - (real::R_1 >> 1));
     632    ret *= exp(-x - (real)(a - 1));
     633
     634    return ret;
    608635}
    609636
  • trunk/src/real.h

    r1024 r1026  
    9292    friend real cbrt(real const &x);
    9393    friend real pow(real const &x, real const &y);
     94    friend real gamma(real const &x);
    9495
    9596    /* Rounding, absolute value, remainder etc. */
Note: See TracChangeset for help on using the changeset viewer.