Divisor summatory function for squares The Divisor summatory function is a function that is a sum over the divisor function.
$$D(x)=\sum_{n\le x} d(n) = 2 \sum_{k=1}^u \lfloor\frac{x}{k}\rfloor - u^2, \;\;\text{with}\; u = \lfloor \sqrt{x}\rfloor$$
http://en.wikipedia.org/wiki/Divisor_summatory_function#Dirichlet.27s_divisor_problem
I am looking for a formula or an efficient algorithm (complexity less than $O(x)$) to calculate the sum of the dividers of the squares. 
$$E(x)=\sum_{n\le x} d(n^2)$$
e.g. $$E(3)=d(1)+d(4)+d(9)=1+3+3=7$$
 A: The number of divisors of a square is the divisor function convolved with the square of the Möbius function (see $g(n)$ here)
$$e(n)=d(n^2)=(d\star\mu^2)(n)$$
by combining these three identities from the linked table
$$e=1\star1\star Q_2,\space d=1\star1,\space Q_2=\mu^2$$
and since
$$\mu^2(n)=\sum_{d^2|n}\mu(d)$$
therefore
$$e(n)=\sum_{a \left| n \right.} d \left( \frac{n}{a} \right) \sum_{b^2 \left| a
\right.} \mu \left( b \right)$$
which can be simplified and rewritten as
$$e(n)=\sum_{b^2 \left| n \right.} \mu \left( b \right) d_3 \left( \frac{n}{b^2}
\right)$$
where $d_3(n)$ is the number of ways that a given number can be written as a product of three integers.  This identity can be verified by noting that $e(n)$ is multiplicative and checking at prime powers which yields $e(p^a)={2a+1}$ and can be compared with $d(p^a)={a+1}$.  In particular note that $d_3(p^a)={\binom{a+2}{2}}$ (see $d_k$ here).
Then the summation of the number of divisors of the square numbers can be computed as:
$$E(x)=\sum_{n\le x} d(n^2)
=\sum_{n \leq x} \sum_{b^2 \left| n \right.} \mu \left( b \right) d_3 \left( \frac{n}{b^2}
\right)$$
which can be reorganized as:
$$E(x)=\sum_{b \leq \sqrt{x}} \mu \left( b \right)  \sum_{n \leq x / b^2} d_3
\left( n \right)$$
$$E(x)=\sum_{a \leq \sqrt{x}} \mu \left( a \right) D_3 \left( \frac{x}{a^2} \right)$$
where $D_3$ is the summatory function for $d_3$.  Since $D_3(x)$ can be computed in $O(x^{2/3})$ time complexity using the three dimensional analogue of the hyperbola method, $E(x)$ can also therefore be computed in
$$O(\int_{a=1}^\sqrt{x}{(\frac{x}{a^2})}^{2/3}da)=O(x^{2/3})$$
which is better than $O(x)$ as desired.
By taking an $O(x^{1/3})$ algorithm to compute $D(x)$ and using it for $D_3(x)$, this bound can be reduced to $O(x^{5/9})$.  You can find such an algorithm and one formula for $D_3(x)$ in my article here which uses the notation $T(n)$ and $T_3(n)$.
