Search code examples
c++seriesgmp

GMP gives wrong result of a series expansion

I am using the GMP package to implement the product of two functions which I express as a cauchy product of two convergent series.

In more detail: I am looking for a way to calculate f(x)=g(x)*h(x) where g(x) is the exponential function and h(x) is a special hypergeometric function (see below), both expressed as series.

My problem is that this works and agrees with my own approximation and wolframalpha's results for x<29, but fails for x>29. In practice I need values of about x=10^6.

The 3 formulas I used are depicted in the image below:

formulas

The Code

void invfak(mpf_t invn, unsigned int n) //Calculates inverse factorial, !/n!
{
        unsigned int i;
        mpf_t k;

        mpf_init_set_d(k,1.0);
        mpf_init_set_d(invn,0.0);
        i=2;

        for (i = 2; i <= n; ++i) {
            mpf_mul_ui(k, k, i);
        }

        mpf_ui_div(invn,1.0,k);
        mpf_clear(k);
}

 void E2F2E(mpf_t result, long double x, unsigned int N)
{
        mpf_t Q1,Q2;          ///gives Nth term in series expansion of exp(x)

        mpf_init_set_d(Q1,x);
        mpf_init_set_d(Q2,0.0);
        mpf_init_set_d(result,0.0);

        mpf_pow_ui(Q1,Q1,N); /// Q1=Q1^N=x^N
        invfak(Q2,N);       /// Q2=1/N!

        mpf_mul(result,Q1,Q2); ///result= Q1*Q2 = x^N/N!

        mpf_clear(Q1);
        mpf_clear(Q2);

}

void E2F2F(mpf_t result, long double x, unsigned int N)
{
        mpf_t Q1,Q2,Q3;         ///gives Nth term in series expansion of 2F2

        mpf_init_set_d(Q1,(N+x)*(N+x));
        mpf_init_set_d(Q2,-x);
        mpf_init_set_d(Q3,0.0);
        mpf_init_set_d(result,0.0);

        mpf_pow_ui(Q2,Q2,N+2);    /// Q2=Q2^(N+2)=(-x)^(N+2)
        invfak(Q3,N);             /// Q3=1/N!
        mpf_mul(Q2,Q2,Q3);        /// Q2=Q2*Q3

        mpf_div(result,Q2,Q1); ///result= Q2/Q1 = .../(N+x)^2

        mpf_clear(Q1);    mpf_clear(Q3);    mpf_clear(Q2)
    }

void Exp2F2gmp(mpf_t result, long double x, unsigned int N)
{
       mpf_t Q1,Qexp,Q2F2,Qsum;
       mpf_init_set_d(Q1,0.0);
       mpf_init_set_d(Qexp,0.0);
       mpf_init_set_d(Q2F2,0.0);
       mpf_init_set_d(Qsum,0.0);

       mpf_init_set_d(result,0.0);

       for(unsigned i = 0; i <= N; ++i){
               mpf_set_d(Qsum,0.0);
               mpf_set_d(Qexp,0.0);
               mpf_set_d(Q2F2,0.0);

               for(unsigned l = 0; l <= i; ++l){ /// a_l und b_i-l
                   E2F2E(Qexp,x,l);
                   E2F2F(Q2F2,x,i-l);
                   mpf_mul(Q1,Qexp,Q2F2);
                   mpf_add(Qsum,Qsum,Q1);
               }
               mpf_add(result,result,Qsum);
               mpf_abs(Qsum,Qsum);
              //if(mpf_cmp_d(Qsum,0.00000001)==-1){ cout << "reached precision at i="<<i; break;}
       }
       cout <<  "\n\n Result = " << result << endl; 
       mpf_clear(Q1);
       mpf_clear(Qexp);
       mpf_clear(Q2F2);
       mpf_clear(Qsum);

}

The function f(x) should approximately go as f(x)=1.05x+1 and also, f(x)>0 for x>0.

But the implementation gives this:

Exp2F2gmp(Q,10,1000) = 12.3707

Exp2F2gmp(Q,20,1000) = 23.1739

Exp2F2gmp(Q,30,1000) = -35195.1

Exp2F2gmp(Q,40,1000) = -2.92079e+13

The first two values agree with Wolframalpha, the second two obviously dont.

Any kind of help would be appreciated, thank you!

Solution

  • This is a textbook example of catastrophic cancellation.

    The terms in the 2F2 series grow to a maximum size of about exp(x), but the magnitude of the 2F2 function is about exp(-x). This means that you need to use at least log_2(exp(2x)) ~= 2.886*x extra bits of precision to compute the 2F2 series accurately, and possibly slightly more depending on how the terms are computed.

    So for example, if x = 29, you need about 83 bits of precision. Your code uses the default precision of the MPF type, which I believe is something like 64 bits. You need to change the code to set the precision of all MPF variables to 64 + 2.886*x bits to get 64 accurate bits (see the GMP manual for how to do this).

    In practice, the series evaluation as you have implemented is not very efficient, and will probably be too slow for x = 1e6.

    One possibility is to use the Arb library (which I develop). This supports computing generalized hypergeometric functions out of the box and uses a much more efficient series evaluation strategy (in this case, using binary splitting). It also uses interval arithmetic, so you get error bounds for free and can set the precision automatically instead of predicting the needed precision in advance (but in this case predicting the precision is easy, and faster).

    Here is code demonstrating how to use it:

    #include "acb_hypgeom.h"

int main()
{
    acb_t z, t;
    acb_struct a[2];
    acb_struct b[2];
    double x;

    acb_init(z); acb_init(t);
    acb_init(a + 0); acb_init(a + 1);
    acb_init(b + 0); acb_init(b + 1);

    for (x = 10.0; x <= 1000000; x *= 10)
    {
        acb_set_d(a + 0, x);
        acb_set_d(a + 1, x);
        acb_set_d(b + 0, x + 1);
        acb_set_d(b + 1, x + 1);
        acb_set_d(z, -x);
        acb_hypgeom_pfq(t, a, 2, b, 2, z, 0, 64 + 2.886 * x);
        acb_neg(z, z);
        acb_exp(z, z, 64);
        acb_mul(t, t, z, 64);
        printf("f(%f) = ", x); acb_printn(t, 20, 0); printf("\n");
    }

    acb_clear(z); acb_clear(t);
    acb_clear(a + 0); acb_clear(a + 1);
    acb_clear(b + 0); acb_clear(b + 1);
}

    Here is the output:

    f(10.000000) = [12.37067931727649929 +/- 5.38e-18]
f(100.000000) = [106.6161729468899444 +/- 4.93e-17]
f(1000.000000) = [1020.154983574938368 +/- 3.54e-16]
f(10000.000000) = [10063.00061277849954 +/- 2.57e-15]
f(100000.000000) = [100198.5001942224819 +/- 6.28e-14]
f(1000000.000000) = [1000626.990558714621 +/- 4.59e-13]

    At x = 1e6, the evaluation takes about 20 seconds (your code would take far, far longer) as a result of 2.9 million bits being used. If this is still too slow, you need to find a better formula to compute f(x), ideally an asymptotic expansion valid when x -> infinity, or perhaps an integral representation (if there is one without cancellation issues).

    Now, if your 2F2 function only depended on x in the final argument and had the first four parameters fixed, there would be a standard formula for this asymptotic expansion, but with the growing parameters involved, I'm not exactly sure how to do it. Since the upper and lower parameters almost "cancel out", it might work to treat them as constants and use the standard asymptotic series with respect to the argument, but I did not check this. Someone with more expertise in asymptotic analysis would have to comment.

    It's also possible that you could use contiguous relations to reduce the 2F2 function to something with small parameters, but I'm not sure if this would be an improvement in practice.