c++ - GMP 给出了一系列扩展的错误结果

Question

我正在使用 GMP 包来实现两个函数的乘积，我表示为两个收敛级数的柯西乘积。

更详细地说:我正在寻找一种方法来计算 f(x)=g(x)*h(x)其中 g(x)是指数函数，h(x)是一个特殊的超几何函数(见下文)，均表示为级数。

我的问题是这行得通并且与我自己的近似值和 wolframalpha 对 x<29 的结果一致, 但因 x>29 而失败.在实践中，我需要大约 x=10^6 的值.

---

我使用的 3 个公式如下图所示:

---

代码

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);

}

函数f(x)应该大约为 f(x)=1.05x+1还有，f(x)>0对于 x>0 .

但是实现给出了这个:

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

前两个值与 Wolframalpha 一致，后两个显然不一致。

如有任何帮助，我们将不胜感激，谢谢!

Answer 1

这是灾难性取消的教科书示例。

2F2 系列中的项增长到最大大小约为 exp(x)，但 2F2 函数的大小约为 exp(-x)。这意味着您至少需要使用 log_2(exp(2x)) ~= 2.886*x 额外的精度位来准确计算 2F2 级数，并且可能略多一些，具体取决于项的计算方式。

例如，如果 x = 29，则需要大约 83 位的精度。您的代码使用 MPF 类型的默认精度，我认为它类似于 64 位。您需要更改代码以将所有 MPF 变量的精度设置为 64 + 2.886*x 位以获得 64 位精确位(有关如何执行此操作，请参阅 GMP 手册)。

在实践中，您实现的系列评估效率不高，并且对于 x = 1e6 可能太慢。

一种可能是使用 Arb库(我开发的)。这支持开箱即用地计算广义超几何函数，并使用更高效的级数评估策略(在本例中，使用二进制拆分)。它还使用区间算法，因此您可以免费获得误差范围，并且可以自动设置精度，而不是提前预测所需的精度(但在这种情况下，预测精度很容易，而且速度更快)。

下面是演示如何使用它的代码:

#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);
}

这是输出:

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]

在 x = 1e6 时，由于使用了 290 万位，计算大约需要 20 秒(您的代码会花费很长时间)。如果这仍然太慢，您需要找到一个更好的公式来计算 f(x)，理想情况下是一个在 x -> 无穷大时有效的渐近展开，或者可能是一个积分表示(如果没有取消问题)。

现在，如果你的 2F2 函数只依赖于最后一个参数中的 x 并且前四个参数是固定的，那么这个渐近展开会有一个标准公式，但是随着涉及的参数越来越多，我不确定如何去做吧。由于上参数和下参数几乎“抵消”，因此将它们视为常量并使用关于参数的标准渐近级数可能会起作用，但我没有检查这一点。在渐近分析方面具有更多专业知识的人必须发表评论。

您也可以使用连续关系将 2F2 函数简化为参数较小的函数，但我不确定这是否会在实践中有所改进。