0

我想计算超过 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 的计算速度？

4

1 回答 1

1

我倾向于认为性能问题主要与有理数有关，而不是与高达 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\2]) * xvecprod(v[#v\2+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\2]) * xvecprod(v[#v\2+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 分钟内产生 a(10^9)，在 22 分钟内产生 a(10^10)。

10^9: 0.02709315486987096878842689330617424348105764850

10^10: 0.02438386113804076644782979967638833694491163817

使用小数时，拆分乘法的技巧不会提高性能，因为数字始终具有相同的位数。但是，我已经留下了代码，因为有可能获得更好的准确性。（至少在理论上。）

我不确定我能否就所需的精度位数给出任何好的建议。（我更像是程序员类型，并且倾向于使用整数）。但是，我的理解是，每次乘法都有可能丢失 1 个二进制数字的精度，尽管由于舍入平均可以采用任何一种方式，它不会那么糟糕。鉴于这是超过 4.5 亿个术语的乘积，这意味着所有精度都会丢失。

但是，使用拆分计算的算法，每个值只经过大约 30 次乘法，所以这只会导致最多 30 个二进制数字（10 个十进制数字）的精度损失，因此使用 100 位精度就足够了. 令人惊讶的是，无论哪种方式，我都得到了相同的答案，所以简单的天真方法似乎有效。

在我的测试中，我注意到 usingforprime比 using 快得多isprime。（例如，该fasterdec版本花费了将近 2 个小时，而简单版本则需要 22 分钟才能达到相同的结果。）。同样，sum(p=1,10^9,isprime(p))大约需要 8 分钟，相比之下，my(t=1);forprime(p=1,10^9,t++);t只需 11 秒。

于 2018-04-05T16:06:59.970 回答