java - 计算最大 10^16 的 totient 函数之和

Question

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

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

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

查找P(10^16)

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

1. 使用 φ(n)= n * product{1<=i<=k} (Pi-1)/Pi 的公式, 我试图计算 φ(n)在范围内，但这对于大 n 变得非常低效.我可以达到 10^7用这种方法。除此之外，它变得太慢了。

2. 我尝试了另一种方法，更直接。维基百科和 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 函数值。

1. 在理解并实现了我在互联网上找到的用于此计算的新递归方法之后:

P(n)=n(n+1)/2−∑{2<=i<=√n} P(⌊n/i⌋)−∑{1<=j<=√n} P( j)⋅(⌊n/j⌋−⌊n/(j+1)⌋)

我可以计算最大为 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 的大值。 ，但我对此一无所知。是否可以进行其他改进，将时间范围缩短为一天左右？

Answer 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 − 1)/i。您可以在最后计算总和，或边走边累加。由于您的筛子会很大，因此您需要将其分段。

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