$$$\overset{n}{\underset{x = 1}{\LARGE \Sigma}} \tau(x)^2 = \overset{n}{\underset{x = 1}{\LARGE \Sigma}} \tau(x^2) \times \left \lfloor \frac{n}{x} \right \rfloor$$$

You can notice that the series $$$S = $$$ { $$$\left \lfloor \frac{n}{1} \right \rfloor, \left \lfloor \frac{n}{2} \right \rfloor, \dots, \left \lfloor \frac{n}{n} \right \rfloor$$$ } have $$$O(S) = O(2\sqrt{n})$$$

Therefore you can split into $$$O(\sqrt{n})$$$ interval of $$$\tau(x^2)$$$ to calculate

Now you need to calculate this fast enough $$$\overset{r}{\underset{x = l}{\LARGE \Sigma}} \tau(x^2)$$$

We also have $$$\overset{n}{\underset{x = 1}{\LARGE \Sigma}} \tau(x^2) = \overset{n}{\underset{x = 1}{\LARGE \Sigma}} S(\left \lfloor \frac{n}{x} \right \rfloor)$$$ where $$$S(n) = \overset{n}{\underset{x = 1}{\LARGE \Sigma}} \mu(x)^2 \times \left \lfloor \frac{n}{x} \right \rfloor$$$

But I am not sure if we can calculate this fast enough

Update: There is a also a formula that $$$T(n) = \overset{n}{\underset{x = 1}{\LARGE \Sigma}} 2^{\omega(x)} \left \lfloor \frac{n}{x} \right \rfloor$$$

Still, I am not sure if we can calculate $$$\overset{r}{\underset{x = l}{\LARGE \Sigma}} 2^{\omega(x)}$$$ fast enough

Editted: Fixed the barrack that was invisble on latex though it should not be like that

Let's iterate through every number $$$d \in [1, n]$$$, and let's see its contribution to the answer. Of course, $$$d$$$ contributes to the Sum of squared divisors of a single number in $$$d^2$$$. How many times does $$$d$$$ do that? Well, there are $$$\left\lfloor\frac{n}{d}\right\rfloor$$$ multiples of $$$d$$$, so the contribution of $$$d$$$ to the answer is $$$\left\lfloor\frac{n}{d}\right\rfloor\cdot d^2$$$.

Hence, the answer is $$$\displaystyle\sum_{d = 1}^n \left(\left\lfloor\frac{n}{d}\right\rfloor\cdot d^2\right)$$$.

Obviously, computing it directly is not good, because $$$n \le 10^{12}$$$. But we can optimize it and make it work in $$$\mathcal{O}(\sqrt{n})$$$.

Fact: The sequence $$$a_1, a_2, \dots, a_n$$$, defined by $$$a_i = \left\lfloor\frac{n}{i}\right\rfloor$$$, has $$$\mathcal{O}(\sqrt{n})$$$ distinct values.

• Proof: For every $$$i > \sqrt{n}$$$, it holds that $$$\left\lfloor\frac{n}{i}\right\rfloor < \sqrt{n}$$$. From that, it follows that there are at most $$$\sqrt{n}$$$ distinct values of $$$a_i$$$ for $$$i >\sqrt{n}$$$. Now, for $$$i \le \sqrt{n}$$$, obviously there are $$$\sqrt{n}$$$ values, and therefore there are at most $$$\sqrt{n}$$$ distinct values of $$$a_i$$$. Thus, there are at most $$$2\sqrt{n}$$$ distinct values in the sequence $$$a_1, a_2, \dots, a_n$$$. This concludes the proof.

Another fact: $$$a_i \ge a_{i+1}$$$.

• Proof: $$$i < i+1\implies \frac{n}{i} > \frac{n}{i+1}\implies \left\lfloor\frac{n}{i}\right\rfloor \ge \left\lfloor\frac{n}{i+1}\right\rfloor$$$.

Now, we could split the sequence $$$a_1, a_2, \dots, a_n$$$ in maximal buckets where all $$$a_i$$$'s in the same bucket are equal. From the previous facts, it follows that there are $$$\mathcal{O}(\sqrt{n})$$$ buckets, and now for each bucket $$$[l, r]$$$ we need to add to the answer $$$\left\lfloor\frac{n}{l}\right\rfloor\cdot (l^2 + (l+1)^2 + \dots + r^2)$$$.

For computing the sum of the squares, we can use a well-known formula: $$$1^2 + 2^2 + 3^2 + \dots + n^2 = \frac{n\cdot(n + 1)\cdot(2n + 1)}{6}$$$.

Now, some sample code:

def sum_of_squares(n):
    return n * (n + 1) * (2*n + 1)//6

n = int(input())
l = 1
ans = 0
while l <= n:
    r = n//(n//l)
    ans += (n//l) * (sum_of_squares(r) - sum_of_squares(l-1))
    l = r + 1
print(ans)