WARNING: The text below provides a guidance for a Project Euler problem.

  • If you pride yourself in tackling problems entirely on your own, you are encouraged to exit this post right now.

  • If you are simply looking for a hint for the problem, please visit the »Project Euler’s official forum«. You may find some hints here, but if you don’t want your problem spoiled, scroll cautiously.

  • If you are looking for a readily available copy/paste answer, you will not find it here. All of the code snippets in the post have been stripped of some crucial lines (these are clearly marked). It is my aim that anyone able to replicate the answer based on such snippets will have to understand the solution conceptually. At that point, by virtue of learning something new, I believe he deserves it anyway.

  • Guidance for future problems will not be published before 100 people have solved the problem and there are at least two more recent problems.

Gcd sum

The problem requires us to compute the double sum . We need to find for a large value , so we better find a sublinear algorithm. We will start by trying to rewrite the sum into a more enlightening form.

Notice that we may rewrite the internal sum by equating the condition with for any :

Where is the Euler’s totient function (Given , counts the number of coprime integers to ). The sum is by the way also known as the »Pillai’s arithmetical function«. From here, we obtain being:

As for any positive numbers such that , the term will appear exactly once in the summation (corresponding ), we may conveniently rewrite this sum as:

where we define to be the summatory function of the euler totients up to :

As the next step in finding a sublinear algorithm we need to realize that there is only possible values . Moreover, for any given value , holds precisely with ’s for which . Let be the sum of the first natural numbers. We continue the derivation of :

What are now the possible values of ? Well, defining for clarity, any will appear in the summation exactly once, at for some . Any will appear in the summation times. We thus divide our sum into two:

Now, provided that we can compute , all is good - we only have summands. Cracking the problem now turns into the following goal:

Find a sublinear way to compute .

For this we will actually use a similar trick as we did in the first time (»Dirichlet hyperbola method«). Recall from elementary number theory (or Wikipedia) that we can write . We may trickily obtain a sum similar to (1):

Analogously, for any positive numbers such that , the term will appear exactly once in the summation (corresponding ) and we conveniently rewrite:

Exposing again the fact that there are only possible values of , in the exactly same way as we did to obtain (2), we find:

And finally:

Using the recurrence (3) to compute will take time; with precomputing some number of the first values of this can be made faster. The following Julia code gives answer in about 80 seconds.

using Nemo
R = ResidueRing(ZZ, 998244353)
n = 10^11

sqrtn = floor(Int, n^0.5)
precompute_bound = floor(Int, n^0.7)

# Sieves Euler's totient up to n
function phis(n)
    phis = collect(1:n)
    for p=2:2:n
        phis[p] >>= 1
    for p=3:2:n
        if phis[p] == p
            phis[p] -= 1
            for j=2p:p:n
                phis[j] -= div(phis[j], p)
    return phis

# Sum of the first n numbers
S(n) = n % 2 == 0 ? R(div(n, 2)) * R(n + 1) : R(n) * R(div(n+1, 2))

# firstΦs[i] = Φ(i)
println("Precomputing first Phis.")
@time firstΦs = map(x -> R(x), accumulate( #<REMOVED>

# secondΦs[i] = Φ(n//i)            Iteratively applies the recurrence (3)
println("Precomputing second Phis.")
secondΦs = [R(0) for i=1:div(n, precompute_bound)]
@time for j = div(n, precompute_bound):-1:1
    k = div(n, j)
    ksqrt = floor(Int, k^0.5)
    for i = ksqrt:-1:2
        if div(k, i) > precompute_bound
            secondΦs[j] -= secondΦs[div(n, div(k, i))]
            secondΦs[j] -= firstΦs[div(k, i)]
    secondΦs[j] -= R(k - div(k, 2)) * firstΦs[1]
    if div(k,ksqrt) == ksqrt
        secondΦs[j] += firstΦs[ksqrt]
    secondΦs[j] += #<REMOVED>

G(n) = #<REMOVED>

# Computing G(n)
println("Computing G(n)")
@time G(n)