计算totient函数的总和高达10^16答案

【问题标题】：Computing the sum of the totient function up to 10^16计算totient函数的总和高达10^16
【发布时间】：2019-09-18 12:25:44
【问题描述】：

以清晰易懂的方式重新发布，没有任何无法正确显示的复杂 MathJax：

为了好玩，我探索了一些计算机科学/数论挑战网站，他们提出了以下问题，具体如下：

让P(n) = sum{1<=k<=n} φ(k)

查找P(10^16)

我为此搜索了很长时间并尝试了不同的方法：

使用φ(n)= n * product{1<=i<=k} (Pi-1)/Pi 的公式，我尝试计算范围内的φ(n)，但是对于大n，这变得非常低效。通过这种方法，我可以达到10^7。除此之外，它变得太慢了。
我尝试了另一种更直接的方法。 Wikipedia 和 Wolfram Alpha 建议直接计算 P(n) 的类似公式：

P(n) = sum {1<=k<=n} φ(k)= 0.5⋅(1+∑{1<=k<=n} μ(k)⋅⌊n/k⌋^2)

这个公式似乎更有希望。我试过了，并设法比10^7 走得更远，但离目标还很远。通过预先计算 Moebius 函数的筛子，我可以得到比10^9 少一点的结果。我的内存不足，无法再计算筛子中的值。而且即使可以，也需要很长时间，离10^16还很远。

下面是我用 Java 编写的第二种方法的部分代码：

public static BigInteger PhiSummatoryFunction (long limit)
{
    BigInteger sum = BigInteger.ZERO;
    int [] m = MoebiusSieve(limit);
    for (int i=1;i<m.length;i++)
        sum=sum.add(BigInteger.valueOf((long) (m[i]*Math.floor(limit/i)*Math.floor(limit/i))));
    return sum.add(BigInteger.ONE).divide(BigInteger.ONE.add(BigInteger.ONE));
}

其中 MoebiusSieve 是一个函数，它使用类似 Eratosthenes 的方法计算 Moebius 函数值直到筛子中的某个限制。

在了解并实现了我在 Internet 上找到的用于该计算的新递归方法后：

P(n)=n(n+1)/2−∑{2

我可以计算高达 P(10^11) 的值，并且在最大内存分配的情况下，预先计算尽可能多的 φ(n) 以及因此我可以用于记忆的所有 P(n)，我可以只计算 P(10^12)超过20分钟。一项重大改进，但距离P(10^16) 仍有一点距离。如果计算时间稍长一点也没关系，但我担心P(10^16) 会花费指数级更长的时间，从P(10^11) 和P(10^12) 之间计算时间的“跳跃”来看。我的记忆允许我“保存”最多350,000,000 φ(n) values或最多700,000,000 μ(k) values。也许有一种方法可以使用 μ(k) 值而不是 φ(n) 来执行求和？

我所有的计算都表明并表明我的递归是主要的时间消耗者。这很明显，但我相信它需要的时间比它应该的要长。我在递归代码下方发布了一些文档。在我看来，这是进行此计算的正确方法，但我的实现并不是最佳的。

public static BigInteger phiR (long limit, long [] s) // limit is 10^t, s is the sieve of precomputed values of `P(n)`. Can store maximum 350,000,000 values
    {                                                                                                                                                       
        if (limit<s.length)                                 
            return BigInteger.valueOf(s[(int) limit]);
        BigInteger sum = BigInteger.valueOf(limit).multiply(BigInteger.valueOf(limit).add(BigInteger.ONE)).divide(BigInteger.valueOf(2)); // this corresponds to the n'th triangular number
        BigInteger midsum1=BigInteger.ZERO; // the first sum
        BigInteger midsum2=BigInteger.ZERO; // the second sum
        long m = 2;
        while (limit/m != limit/(m+1) && m*m<=limit) // computing the first sum, first for changing floor(limit/m) values
        {
            midsum1=midsum1.add(phiR((long) Math.floor(limit/m),s));
            m++;
        }
        for (long k = m;k*k<=limit;k++) // once the floors become constant for some values,-->
        {                               //  can check how many times the value appears, and multiply accordingly,--> 
            BigInteger midPhi = phiR((long) Math.floor(limit/k),s);  // rather than compute the Phi every time
            long q = 1;
            while (limit/k==limit/(k+1)&&k*k<=limit)
            {
                q++;
                k++;
            }
            k--;
            midPhi=midPhi.multiply(BigInteger.valueOf(q));
            midsum1=midsum1.add(midPhi);
        }
        for (long d=1;d*d<=limit;d++) // computing the second sum
            if ((double)d!=Math.floor(limit/d))
                midsum2=midsum2.add(BigInteger.valueOf((long) (Math.floor(limit/d)-Math.floor(limit/(d+1)))).multiply(phiR(d,s)));
        sum=sum.subtract(midsum1).subtract(midsum2);
        return sum;
    }

除了数组之外，有人建议我使用 dictinaries 来获取 n 的大值，但我对此一无所知。是否可以进行另一项改进以使时间框架为一天左右？

【问题讨论】：

您能否分享方法 3 的链接和/或解释其背后的理论？
@גלעדברקן 当然，请在此处查看qwr 的回答：math.stackexchange.com/questions/316376/…
谢谢。有趣的是，那里的算法是由 Daniel Fischer 提供的，他是写 SO answer 我在我的（现已删除）答案中分享的同一个人 :)
是的，我也注意到了。我认为我的递归由于某种原因比它应该的慢。我不知道如何使它更快。有人建议使用字典来存储P(n) 值。
我很困惑：不使用 P(floor(n/i))（在方法 3 中）是否意味着我们从一开始就需要 P(5*10^15) 的值？

标签： java arrays algorithm function optimization

【解决方案1】：

如果你想知道单个数n的总和，最好的方法是对n进行因式分解，并取每个因数减去1的乘积;例如，30 = 2 * 3 * 5，从每个因子中减去 1，然后相乘，得到总和 1 * 2 * 4 = 8。但是如果你想找到小于给定 n，比将它们中的每一个分解更好的方法是筛分。思路很简单：设置一个从0到n的数组X，在每个X_i中存储i，然后从 0 开始遍历数组，并且每当 X_i = i 循环遍历 i 的倍数，将每个乘以 (i em> - 1)/i。您可以在最后计算总和，也可以随时累积。由于您的筛子很大，因此您需要对其进行分段。

以下是我博客中的一些有用页面：Sieving For Totients 和 Segmented Sieve of Eratosthenes。如果您在那里四处寻找，您可能还会发现其他一些有趣的东西。

【讨论】：

正如我所提到的，我尝试了这种方法，但与其他方法相比它相当慢。但即使我这样做了，我该如何应用细分？接下来的每个片段都依赖于前一个片段，所以在我看来，一旦我“抛出”前一个片段，我就无法真正产生后续片段。

【解决方案2】：

对于较低的订单：

public class Solution {

    static Map<Long,Long> X2 = new HashMap<>();

    static long F2(long N){
        return N*(N-1)/2;
    }

    static long R2(long N){
        if(N==1) return 0;
        if(X2.containsKey(N)) return X2.get(N);
        long sum = F2(N);
        long m=2;
        while(true){
            long x = N/m;
            long nxt = N/x;
            if(nxt >= N) {
                X2.put(N, sum - (N-m+1)*R2(N/m));
                return X2.get(N);
            }
            sum -= (nxt-m+1) * R2(N/m);
            m = nxt+1;
        }
    }

    public long odradi(int N){
        return R2(N)+1;
    }

    public static void main(String[] args) {
        System.out.println(R2(100000)+1);
    }
}

【讨论】：