直到 n 与 n 的所有数字的最大公约数之和答案

【问题标题】：Sum of Greatest Common Divisor of all numbers till n with n直到 n 与 n 的所有数字的最大公约数之和
【发布时间】：2016-02-07 16:26:26
【问题描述】：

从 1 到 n 有 n 个数。我需要找到 ∑gcd(i,n) 其中 i=1 到 i=n 对于 10^7 范围内的 n。我使用 euclid 的 gcd 算法，但它给了 TLE。有没有什么有效的方法可以求出上面的总和？

#include<bits/stdc++.h>
using namespace std;
typedef long long int ll;
int gcd(int a, int b)
{
    return b == 0 ? a : gcd(b, a % b);
}
int main()
{
   ll n,sum=0;
   scanf("%lld",&n);
   for(int i=1;i<=n;i++)
   {
       sum+=gcd(i,n);
   }
   printf("%lld\n",sum);
   return 0;
}

【问题讨论】：

你能贴出你试过的代码吗？
@AndyM 我已经添加了代码
我没有看到任何名为 A 的数组...您的 i 是 from 0 to n 而不是 from 1 to n...n 是从 cin 读取的，而您写的是 @ 987654327@...我猜t是测试用例的数量？也许你应该让你的问题与你的代码相匹配......
@SimonKraemer 对错误感到抱歉。 i 从 1 到 n 并且 n 的范围是 10^7 不完全是 10^7。
bits/stdc++ 请不要。

标签： c++ arrays algorithm greatest-common-divisor

【解决方案1】：

您可以通过批量 GCD 计算来完成。你应该找到所有简单的除数和这些除数的幂。这可以在 Sqtr(N) 复杂度中完成。之后需要编写 GCD 表。

可以在C#上写sn-p，转成C++也不难

int[] gcd = new int[x + 1];
for (int i = 1; i <= x; i++) gcd[i] = 1;
for (int i = 0; i < p.Length; i++)
    for (int j = 0, h = p[i]; j < c[i]; j++, h *= p[i])
        for (long k = h; k <= x; k += h)
            gcd[k] *= p[i];

long sum = 0;
for (int i = 1; i <= x; i++) sum += gcd[i];

p 它是简单除数的数组和这个除数的 c 次幂。

例如，如果 n = 125

p = [5]
c = [3]
125 = 5^3

如果 n = 12

p = [2,3]
c = [2,1]
12 = 2^2 * 3^1

【讨论】：

在这样的“在线裁判”中使用 ~80MB 是不可能的，而且乍一看这实际上看起来更慢
更新到 40Mb。 40 和 80 Mb 对于当前的在线评委来说是可以的。复杂性非常接近 NLogN。
啊，由于编辑前的缩进，我误读了循环的嵌套方式...
更新：复杂度 O(N) 不是 O(NlogN)
如果你有 x 的因式分解，那么就不需要所有这些循环和用数字填充巨大的数组 - 只需应用 deviantfan 提到的MathOverflow article 中的公式 (S(m*n) = S(m) * S(n) 表示 m,n 互质，S(p^a) = (a + 1) * p^a + a * p^(a-1) 表示 p 素数。这介于 O(log N) 和 O(sqrt(N)) 之间，具体取决于您如何进行分解。

【解决方案2】：

我刚刚在两个数字之间实现了 GCD 算法，这很容易，但我无法得到你想要在那里做的事情。我在那里读到的是你试图总结一系列 GCD；但是 GCD 是两个或多个数字之间一系列数学运算的结果，这些运算产生一个值。我不是数学家，但我认为你写的“sigma”意味着你试图总结 1 到 10.000.000 之间数字的 GCD；这对我来说根本没有意义。

您试图找到 GCD 的值是什么？ 1 到 10.000.000 之间的所有数字？我怀疑是这样的。

不管怎样，这里是欧几里得 GCD 算法的一个非常基本（而且很匆忙）的实现：

int num1=0, num2=0;
cout << "Insert the first number: ";
cin >> num1;

cout << "\n\nInsert the second number: ";
cin >> num2;
cout << "\n\n";
fflush(stdin);

while ((num1 > 0) && (num2 > 0))
{
    if ((num1 - num2) > 0)
    {
        //cout << "..case1\n";
        num1 -= num2;
    }
    else if ((num2 - num1) > 0)
    {
        //cout << "..case2\n";
        num2 -= num1;
    }
    else if (num1 = num2)
    {
        cout << ">>GCD = " << num1 << "\n\n";
        break;
    }
}

【讨论】：

【解决方案3】：

开始研究这个问题的一个好地方是整数序列在线百科全书上的here，因为您要做的是计算 1 和 N 之间的序列 A018804 的总和。正如您发现的方法尝试使用简单的 Euclid GCD 函数太慢，所以您需要一种更有效的方法来计算结果。

根据从 OEIS 链接的一个paper，可以根据欧拉函数重写总和。这将问题变成了质因数分解之一 - 仍然不容易，但可能比蛮力要快得多。

【讨论】：

【解决方案4】：

您可以使用 Seive 存储小于等于 10^7 的所有数字的最低素因数和给定数字的素数分解直接计算你的答案..

【讨论】：

【解决方案5】：

我有机会研究 GCD 和的计算，因为这个问题出现在一个名为 GCD Sum 的 HackerEarth 教程中。谷歌搜索出现了 some academic papers 有用的公式，我在这里报告，因为它们在 deviantfan 链接的 MathOverflow article 中没有提到。

对于互质 m 和 n（即 gcd(m, n) == 1），函数是乘法的：

gcd_sum[m * n] = gcd_sum[m] * gcd_sum[n]

素数 p 的 e 次幂：

gcd_sum[p^e] = (e + 1) * p^e - e * p^(e - 1)

如果只计算一个总和，那么这些公式可以应用于分解相关数字的结果，这仍然比重复 gcd() 调用或通过 rigmarole proposed by Толя 更快。

但是，这些公式可以很容易地用于有效地计算整个函数表。基本上，您所要做的就是将它们插入linear time Euler totient calculation 的算法中，您就完成了 - 这计算 all GCD 总和比计算单个 GCD 总和要快得多数字 10^6 通过调用 gcd() 函数。基本上，该算法有效地枚举了直到 n 的数字的最小因子分解，从而可以轻松计算任何乘法函数 - Euler totient（又名 phi）、sigmas，或者实际上是 GCD 和。

这里有一段哈希代码，用于计算 GCD 和的表以用于较小的限制 - “小”是指 sqrt(N) * N 不会溢出 32 位有符号整数。 IOW，它适用于 10^6 的限制（对于限制为 5 * 10^5 的 HackerEarth 任务来说已经足够了），但 10^7 的限制需要在几个战略位置坚持 (long) 演员表。然而，这种在更高范围内操作的功能强化留给读者作为众所周知的练习...... ;-)

static int[] precompute_Pillai (int limit)
{
    var small_primes = new List<ushort>();
    var result = new int[1 + limit];

    result[1] = 1;

    int n = 2, small_prime_limit = (int)Math.Sqrt(limit);

    for (int half = limit / 2; n <= half; ++n)
    {
        int f_n = result[n];

        if (f_n == 0)
        {
            f_n = result[n] = 2 * n - 1;

            if (n <= small_prime_limit)
            {
                small_primes.Add((ushort)n);
            }
        }

        foreach (int prime in small_primes)
        {
            int nth_multiple = n * prime, e = 1, p = 1;  // 1e6 * 1e3 < INT_MAX

            if (nth_multiple > limit)
                break;

            if (n % prime == 0)
            {
                if (n == prime)
                {
                    f_n = 1;
                    e = 2;
                    p = prime;
                }
                else break;
            }

            for (int q; ; ++e, p = q)
            {
                result[nth_multiple] = f_n * ((e + 1) * (q = p * prime) - e * p);

                if ((nth_multiple *= prime) > limit)
                    break;
            }
        }
    }

    for ( ; n <= limit; ++n)
        if (result[n] == 0)
            result[n] = 2 * n - 1;

    return result;
}

正如所承诺的，这将在 12.4 毫秒内计算所有 GCD 总和，达到 500,000，而在同一台机器上通过 gcd() 调用计算 500,000 的单个总和需要 48.1 毫秒。该代码已针对 OEIS list of the Pillai function (A018804) 进行了高达 2000 次验证，并针对基于 gcd 的函数进行了高达 500,000 次验证——这项工作耗时整整 4 小时。

可以应用一系列优化来显着加快代码速度，例如将模除法替换为乘法（使用逆运算）和比较，或者通过步进“素数”来减少更多毫秒Cleaner-upper' 循环模 6。但是，我想以基本的、未优化的形式展示该算法，因为 (a) 它的速度非常快，并且 (b) 它可能对其他乘法函数有用，而不仅仅是 GCD总和。

PS：Granlund/Montgomery 论文Division by Invariant Integers using Multiplication 的第 9 节描述了通过与逆相乘进行的模测试，但很难找到有关高效计算 2 的逆模幂的信息。大多数来源使用扩展欧几里得算法或类似的矫枉过正。所以这里有一个计算乘法逆模 2^32 的函数：

static uint ModularInverse (uint n)
{
    uint x = 2 - n;

    x *= 2 - x * n;
    x *= 2 - x * n;
    x *= 2 - x * n;
    x *= 2 - x * n;

    return x;
}

这实际上是Newton-Raphson 的五次迭代，以防万一。 ;-)

【讨论】：