我怎样才能用 PARI/GP 更快地计算这个主要产品？
How can I calculate this prime product faster with PARI/GP?
我想计算 1-1/p 的乘积，其中 p 在素数上运行高达 10^10
我知道近似值 exp(-gamma)/ln(10^10) ，其中 gamma 是 Euler-Mascheroni 常数，ln 是自然对数，但我想计算确切的乘积以查看与近似值是。
问题是 PARI/GP 需要很长时间来计算从大约 4.2 * 10^9 到 10^10 的素数。 prodeuler-command也需要很长时间。
有什么方法可以加快PARI/GP的计算速度吗？
我倾向于认为性能问题主要与有理数有关，而不是生成 10^10 以内的素数。
作为快速测试，我 运行
a(n)=my(t=0);forprime(p=1,n,t+=p);t
与 a(10^10) 计算并在几分钟内完成计算，这似乎是合理的。
您的要求对应的程序是：
a(n)=my(t=1);forprime(p=1,n,t*=(1-1/p));t
这比第一个程序运行得慢得多，所以我的问题是问是否有一种方法可以重新制定计算以避免有理数直到最后？我上面的表述是否如您所愿？ - 即使是 10^6 的数字也非常大，所以计算需要很长时间也就不足为奇了，也许这个问题与数字的合理性无关，而只是它们的大小。
我用来计算大乘积的一个技巧是拆分问题，以便在每个阶段乘法左右两侧的数字大小大致相同。例如计算一个大的阶乘，比如 8！计算 ((1*8)*(2*7))*((3*6)*(4*5)) 比明显的从左到右的方法更有效。
以下是使用精确算术快速尝试执行您想要的操作。达到 10^8 大约需要 8 分钟，但分子的大小已经是 190 万位数字，因此在 运行 内存不足之前这不太可能达到 10^10。 [即使对于这个计算，我也需要增加堆栈大小]。
xvecprod(v)={if(#v<=1, if(#v,v[1],1), xvecprod(v[1..#v]) * xvecprod(v[#v+1..#v]))}
faster(n)={my(b=10^6);xvecprod(apply(i->xvecprod(
apply(p->1-1/p, select(isprime, [i*b+1..min((i+1)*b,n)]))), [0..n\b]))}
使用小数肯定会加快速度。以下程序运行相当快，最多 10^8，精度为 1000 位。
xvecprod(v)={if(#v<=1, if(#v,v[1],1), xvecprod(v[1..#v]) * xvecprod(v[#v+1..#v]))}
fasterdec(n)={my(b=10^6);xvecprod(apply(i->xvecprod(
apply(p->1-1.0/p,select(isprime,[i*b+1..min((i+1)*b,n)]))),[0..n\b]))}
使用小数的最快方法也是最简单的：
a(n)=my(t=1);forprime(p=1,n,t*=(1-1.0/p));t
将精度设置为 100 位十进制数字，这将在 2 分钟内生成 (10^9)，在 22 分钟内生成 (10^10)。
10^9: 0.02709315486987096878842689330617424348105764850
10^10: 0.02438386113804076644782979967638833694491163817
处理小数时，拆分乘法的技巧不会提高性能，因为数字始终具有相同的位数。但是，我留下了代码，因为有可能提高准确性。 （至少在理论上是这样。）
我不确定我能否就所需的精度位数提供任何好的建议。 （我更像是一个程序员类型，倾向于使用整数）。但是，我的理解是每次乘法都有可能丢失 1 个二进制数字的精度，尽管由于舍入可以平均采用任何一种方式，所以不会那么糟糕。鉴于这是超过 4.5 亿项的产品，这意味着所有精度都将丢失。
但是，使用拆分计算的算法，每个值只经过大约 30 次乘法运算，因此最多只会损失 30 位二进制数字（10 位十进制数字）的精度，因此使用 100 位数字精度应该足够了。令人惊讶的是，无论哪种方式我都得到了相同的答案，所以简单的天真方法似乎有效。
在我的测试中，我注意到使用 forprime 比使用 isprime 快得多。 （例如，fasterdec 版本花费了将近 2 小时，而简单版本花费了 22 分钟才获得相同的结果。）。类似地，sum(p=1,10^9,isprime(p)) 大约需要 8 分钟，而 my(t=1);forprime(p=1,10^9,t++);t 只需要 11 秒。
我想计算 1-1/p 的乘积，其中 p 在素数上运行高达 10^10
我知道近似值 exp(-gamma)/ln(10^10) ，其中 gamma 是 Euler-Mascheroni 常数，ln 是自然对数，但我想计算确切的乘积以查看与近似值是。
问题是 PARI/GP 需要很长时间来计算从大约 4.2 * 10^9 到 10^10 的素数。 prodeuler-command也需要很长时间。
有什么方法可以加快PARI/GP的计算速度吗？
我倾向于认为性能问题主要与有理数有关，而不是生成 10^10 以内的素数。
作为快速测试，我 运行
a(n)=my(t=0);forprime(p=1,n,t+=p);t
与 a(10^10) 计算并在几分钟内完成计算，这似乎是合理的。
您的要求对应的程序是：
a(n)=my(t=1);forprime(p=1,n,t*=(1-1/p));t
这比第一个程序运行得慢得多，所以我的问题是问是否有一种方法可以重新制定计算以避免有理数直到最后？我上面的表述是否如您所愿？ - 即使是 10^6 的数字也非常大，所以计算需要很长时间也就不足为奇了，也许这个问题与数字的合理性无关，而只是它们的大小。
我用来计算大乘积的一个技巧是拆分问题，以便在每个阶段乘法左右两侧的数字大小大致相同。例如计算一个大的阶乘，比如 8！计算 ((1*8)*(2*7))*((3*6)*(4*5)) 比明显的从左到右的方法更有效。
以下是使用精确算术快速尝试执行您想要的操作。达到 10^8 大约需要 8 分钟，但分子的大小已经是 190 万位数字，因此在 运行 内存不足之前这不太可能达到 10^10。 [即使对于这个计算，我也需要增加堆栈大小]。
xvecprod(v)={if(#v<=1, if(#v,v[1],1), xvecprod(v[1..#v]) * xvecprod(v[#v+1..#v]))}
faster(n)={my(b=10^6);xvecprod(apply(i->xvecprod(
apply(p->1-1/p, select(isprime, [i*b+1..min((i+1)*b,n)]))), [0..n\b]))}
使用小数肯定会加快速度。以下程序运行相当快，最多 10^8，精度为 1000 位。
xvecprod(v)={if(#v<=1, if(#v,v[1],1), xvecprod(v[1..#v]) * xvecprod(v[#v+1..#v]))}
fasterdec(n)={my(b=10^6);xvecprod(apply(i->xvecprod(
apply(p->1-1.0/p,select(isprime,[i*b+1..min((i+1)*b,n)]))),[0..n\b]))}
使用小数的最快方法也是最简单的：
a(n)=my(t=1);forprime(p=1,n,t*=(1-1.0/p));t
将精度设置为 100 位十进制数字，这将在 2 分钟内生成 (10^9)，在 22 分钟内生成 (10^10)。
10^9: 0.02709315486987096878842689330617424348105764850
10^10: 0.02438386113804076644782979967638833694491163817
处理小数时，拆分乘法的技巧不会提高性能，因为数字始终具有相同的位数。但是，我留下了代码，因为有可能提高准确性。 （至少在理论上是这样。）
我不确定我能否就所需的精度位数提供任何好的建议。 （我更像是一个程序员类型，倾向于使用整数）。但是，我的理解是每次乘法都有可能丢失 1 个二进制数字的精度，尽管由于舍入可以平均采用任何一种方式，所以不会那么糟糕。鉴于这是超过 4.5 亿项的产品，这意味着所有精度都将丢失。
但是，使用拆分计算的算法，每个值只经过大约 30 次乘法运算，因此最多只会损失 30 位二进制数字（10 位十进制数字）的精度，因此使用 100 位数字精度应该足够了。令人惊讶的是，无论哪种方式我都得到了相同的答案，所以简单的天真方法似乎有效。
在我的测试中，我注意到使用 forprime 比使用 isprime 快得多。 （例如，fasterdec 版本花费了将近 2 小时，而简单版本花费了 22 分钟才获得相同的结果。）。类似地，sum(p=1,10^9,isprime(p)) 大约需要 8 分钟，而 my(t=1);forprime(p=1,10^9,t++);t 只需要 11 秒。