How to calculate the Pillai function

Question

I want to calculate

$$\sum_{i=1}^n \gcd(i,n)$$

(Pillai function) for large n ($\ 1 \leq n \leq 10^{15}\ $). There is one formula involving Euler's totient formula but that requires calculating values of Euler's totient function for all divisors of $\ n\ $. What is the best method to do so?

Answer 1

Even in python.

This script

from sympy.ntheory import totient, divisors
from random import sample

def pillai(n):
    return sum([n*totient(d)//d for d in divisors(n)])

pop = list(range(10**15-10**6,10**15))
test = sample(pop, 10)
for t in test:
    print(t, pillai(t))

prints

999999999867116 15528648791180712
999999999329163 12380949940095165
999999999556314 9978540768097875
999999999764414 22495561180844031
999999999634550 30747783543314175
999999999374335 14096008679619933
999999999569680 84064457099079216
999999999952264 38933431675509300
999999999387833 3996070724468241
999999999458806 5999960040520095

essentially instantaneously.