Partial sums of Pillai's function Pillai’s arithmetical function (gcd-sum function) is defined by
$$
P(n) = \sum_{k=1}^n\gcd(k,n)
$$
Let $\sum_{n\leq x}P(n)$ be summation of all values of P for all $n$ up to given $x$.
I dervied that
$$
\sum_{n\leq x}P(n) = \sum_{d\leq x}\mu(d)\sum_{k\leq\frac{x}{d}}k\tau(k)
$$
where $\mu$ is a Moebius function and $\tau$ is a divisor function.
But it has at least an O(n) complexity.
I wonder is it possible to get the summation in sublinear time?

There is an  approximation that may serve your  purposes. First derive
an alternate expression for $P(n):$
$$P(n) = \sum_{k=1}^n \gcd(k, n)
= \sum_{d|n} \sum_{k=1 \atop \gcd(k, n)=d}^n d
= \sum_{d|n} d \sum_{\gcd(k/d, n/d)=1} 1
= \sum_{d|n} d \times \varphi\left(\frac{n}{d}\right).$$
Now recall that
$$\sum_{d|n} \varphi(d) = n$$
so that
$$\sum_{n\ge 1} \frac{\varphi(n)}{n^s} \times \zeta(s) =
\sum_{n\ge 1} \frac{n}{n^s}$$
which implies
$$\sum_{n\ge 1} \frac{\varphi(n)}{n^s} = \frac{\zeta(s-1)}{\zeta(s)}.$$
Introduce
$$L(s) = \sum_{n\ge 1} \frac{P(n)}{n^s}.$$
Using the convolution formula for $P(n)$ we thus have
$$L(s) = \frac{\zeta(s-1)^2}{\zeta(s)}.$$
We can  now use  the Mellin-Perron formula  to predict, but  not quite
prove, the asymptotics of the partial sums of $P(n).$ We have that
$$\sum_{k=1}^n P(n)
= \frac{1}{2} P(n) + 
\frac{1}{2\pi i} \int_{c-i\infty}^{c+i\infty} L(s) \frac{n^s}{s} ds$$
where $c>2.$
We shift this integral to the left.
But the dominant pole at $s=2$ has residue
$$\mathrm{Res}\left(L(s)\frac{n^s}{s}; s=2\right) =
\frac{3}{\pi^2} n^2 \log n
+ n^2
\left(\frac{6}{\pi^2}\gamma
- \frac{3}{2\pi^2}
- \frac{18}{\pi^4} \zeta'(2) \right).$$
This yields the following approximation:
$$\sum_{k=1}^n P(n) \sim
\frac{1}{2} P(n) + 
\frac{3}{\pi^2} n^2 \log n
+ n^2
\left(\frac{6}{\pi^2}\gamma
- \frac{3}{2\pi^2}
- \frac{18}{\pi^4} \zeta'(2) \right).$$
Shifting the integral further to the left would pick up a contribution
from the non-trivial zeros of the zeta function and hence is best left
to professionals.
For your information the exact value for the sum up to $P(1000)$ is
$$ 2475190 $$
while the approximation yields
$$ 2476126.$$
The exact value for the sum up to $P(1200)$ is
$$ 3649500$$
while the approximation yields
$$ 3647105.$$
Remark. As the  question text does not prove  the equation for the
sum that is quoted there I will  include it here, so that it counts as
verified.
Start with
$$P(n) = \sum_{d|n} d \times \varphi\left(\frac{n}{d}\right)
= n \sum_{d|n} \frac{\varphi(d)}{d}.$$
Now observe that
$$\frac{\varphi(d)}{d} = \prod_p \left(1-\frac{1}{p}\right)$$
where the product ranges over the prime divisors $p$ of $n.$
This immediately implies that
$$ \frac{\varphi(d)}{d} = \sum_{q|d} \frac{\mu(q)}{q}$$
by definition of the Mobius function.
Returning to the sum we obtain
$$\sum_{n\le x} P(n)
= \sum_{n\le x} n \sum_{d|n} \sum_{q|d} \frac{\mu(q)}{q}.$$
Switching the two inner sums and putting $d=qf$ yields
$$\sum_{n\le x} n
\sum_{q|n} \frac{\mu(q)}{q} \sum_{f|n/q} 1
= \sum_{n\le x} n
\sum_{q|n} \frac{\mu(q)}{q} \tau(n/q)$$
which yields
$$\sum_{n\le x}
\sum_{q|n} \mu(q) \times n/q \times \tau(n/q).$$
Summing over all $q$ and putting $n=kq$ we finally obtain
$$\sum_{q\le x} \mu(q) \sum_{k\le\lfloor x/q\rfloor} k \; \tau(k),$$
confirming the statement in the question.
