GMP 给出了一系列扩展的错误结果
GMP gives wrong result of a series expansion
我正在使用 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 一致,后两个显然不一致。
如有任何帮助,我们将不胜感激,谢谢!
这是灾难性取消的教科书示例。
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 并且前四个参数是固定的,那么这个渐近展开会有一个标准公式,但是随着参数的增加,我不太清楚确定如何去做。由于上参数和下参数几乎 "cancel out",将它们视为常量并使用关于参数的标准渐近级数可能会起作用,但我没有检查这一点。在渐近分析方面具有更多专业知识的人将不得不发表评论。
您也可以使用连续关系将 2F2 函数简化为参数较小的函数,但我不确定这是否会在实践中有所改进。
我正在使用 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 一致,后两个显然不一致。
如有任何帮助,我们将不胜感激,谢谢!
这是灾难性取消的教科书示例。
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 并且前四个参数是固定的,那么这个渐近展开会有一个标准公式,但是随着参数的增加,我不太清楚确定如何去做。由于上参数和下参数几乎 "cancel out",将它们视为常量并使用关于参数的标准渐近级数可能会起作用,但我没有检查这一点。在渐近分析方面具有更多专业知识的人将不得不发表评论。
您也可以使用连续关系将 2F2 函数简化为参数较小的函数,但我不确定这是否会在实践中有所改进。