Sum of Greatest Common Divisor of all numbers till n with n

Question

There are n numbers from 1 to n. I need to find the ∑gcd(i,n) where i=1 to i=n for n of the range 10^7. I used euclid's algorithm for gcd but it gave TLE. Is there any efficient method for finding the above sum?

#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;
}

Answer 1

You can do it via bulk GCD calculation. You should found all simple divisors and powers of these divisors. This is possible done in Sqtr(N) complexity. After required compose GCD table.

May code snippet on C#, it is not difficult to convert into 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 it is array of simple divisors and c power of this divisor.

For example if n = 125

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

if n = 12

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

Answer 2

I've just implemented the GCD algorithm between two numbers, which is quite easy, but I cant get what you are trying to do there. What I read there is that you are trying to sum up a series of GCD; but a GCD is the result of a series of mathematical operations, between two or more numbers, which result in a single value. I'm no mathematician, but I think that "sigma" as you wrote it means that you are trying to sum up the GCD of the numbers between 1 and 10.000.000; which doesnt make sense at all, for me.

What are the values you are trying to find the GCD of? All the numbers between 1 and 10.000.000? I doubt that's it.

Anyway, here's a very basic (and hurried) implementation of Euclid's GCD algorithm:

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;
    }
}

Answer 3

A good place to start looking at this problem is here at the Online Encyclopedia of Integer Sequences as what you are trying to do is compute the sum of the sequence A018804 between 1 and N. As you've discovered approaches that try to use simple Euclid GCD function are too slow so what you need is a more efficient way to calculate the result.

According to one paper linked from the OEIS it's possible to rewrite the sum in terms of Euler's function. This changes the problem into one of prime factorisation - still not easy but likely to be much faster than brute force.

Answer 4

I had occasion to study the computation of GCD sums because the problem cropped up in a HackerEarth tutorial named GCD Sum. Googling turned up some academic papers with useful formulas, which I'm reporting here since they aren't mentioned in the MathOverflow article linked by deviantfan.

For coprime m and n (i.e. gcd(m, n) == 1) the function is multiplicative:

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

Powers e of primes p:

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

If only a single sum is to be computed then these formulas could be applied to the result of factoring the number in question, which would still be way faster than repeated gcd() calls or going through the rigmarole proposed by Толя.

However, the formulas could just as easily be used to compute whole tables of the function efficiently. Basically, all you have to do is plug them into the algorithm for linear time Euler totient calculation and you're done - this computes all GCD sums up to a million much faster than you can compute the single GCD sum for the number 10^6 by way of calls to a gcd() function. Basically, the algorithm efficiently enumerates the least factor decompositions of the numbers up to n in a way that makes it easy to compute any multiplicative function - Euler totient (a.k.a. phi), the sigmas or, in fact, GCD sums.

Here's a bit of hashish code that computes a table of GCD sums for smallish limits - ‘small’ in the sense that sqrt(N) * N does not overflow a 32-bit signed integer. IOW, it works for a limit of 10^6 (plenty enough for the HackerEarth task with its limit of 5 * 10^5) but a limit of 10^7 would require sticking (long) casts in a couple of strategic places. However, such hardening of the function for operation at higher ranges is left as the proverbial exercise for the reader... ;-)

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;
}

As promised, this computes all GCD sums up to 500,000 in 12.4 ms, whereas computing the single sum for 500,000 via gcd() calls takes 48.1 ms on the same machine. The code has been verified against an OEIS list of the Pillai function (A018804) up to 2000, and up to 500,000 against a gcd-based function - an undertaking that took a full 4 hours.

There's a whole range of optimisations that could be applied to make the code significantly faster, like replacing the modulo division with a multiplication (with the inverse) and a comparison, or to shave some more milliseconds by way of stepping the ‘prime cleaner-upper’ loop modulo 6. However, I wanted to show the algorithm in its basic, unoptimised form because (a) it is plenty fast as it is, and (b) it could be useful for other multiplicative functions, not just GCD sums.

P.S.: modulo testing via multiplication with the inverse is described in section 9 of the Granlund/Montgomery paper Division by Invariant Integers using Multiplication but it is hard to find info on efficient computation of inverses modulo powers of 2. Most sources use the Extended Euclid's algorithm or similar overkill. So here comes a function that computes multiplicative inverses modulo 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;
}

That's effectively five iterations of Newton-Raphson, in case anyone cares. ;-)

Answer 5

you can use Seive to store lowest prime Factor of all number less than equal to 10^7 and the by by prime factorization of given number calculate your answer directly..