自问自答。
网上流传着的做法都是洲阁筛做法，从ranking上看跑的都不是非常快，但是Min_25只跑了7.1s，想必正解是其它做法，我昨晚努力思考了一下，想出了一个$O(n^{2/3})$的别解，虽然不清楚是不是正解大抵是正解，但是应该常数上比洲阁筛快上不少。
先求$\sigma_0(n^3)$的值，令$n=\prod p_i^{e_i}$，考虑哪些因子只有$n^3$有，而$n$没有的。不妨假设$d=\prod p_i^{d_i}$，显然那些$d_i > e_i$的是$n^3$独有的。那么对于每个$n$的约数$d$，贡献就是$3^{\omega(d)}$（每个幂次可以加$e_i$，加$2e_i$或者保持不变），也就是说：$$\sigma_0(n^3)=\sum_{d|n}3^{\omega(d)}$$
然后，发挥了oeis的作用，得知$3^{\omega(n)} = \sum\limits_{d|n} \sigma_0(d) \mu^2(d)$。于是经过一系列化简，我们就得到了下面的式子：$$\sum\limits_{i=1}^{n} \sigma_0(i^3)=\sum_{i=1}^n\sigma_0(i)\mu^2(i)\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor} \sigma_0(j)$$
令$f(n)=\sum\limits_{i=1}^{n} \sigma_0(i)\mu^2(i)$，$g(n)=\sum\limits_{i=1}^{n} \sigma_0(i)$，只要我们能快速求出每个$f(\lfloor \frac{n}{i} \rfloor)$和$g(\lfloor \frac{n}{i} \rfloor)$，这题就算解决了。
$g(n)$可以参考DIVCNT1，我们有$O(n^{0.6})$的做法，或者暴力一些，有一个$O(n^{2/3})$的做法。
考虑求$f(n)$，回忆下$\sum\limits_{i=1}^{n}\mu^2(i)$的$O(\sqrt{n})$做法，把这个（其实是个容斥）套用过来，可以得到$$f(n)=\sum_{i=1}^{\lfloor \sqrt{n} \rfloor} \mu(i) \sum_{d=1}^{\lfloor \frac{n}{i^2} \rfloor}\sigma_0(i^2d)$$
考虑$\sigma_0(nm)$的值，可以简单推导（$d \vert nm$等价于$d=\frac{pm}{q}, p|n, q|m, (p,q)=1$）得到$$\sigma_0(nm)=\sum_{d|(n,m)}\mu(d)\sigma_0(\frac{n}{d})\sigma_0(\frac{m}{d})$$
带入上面的$f(n)$，可以得到$$f(n)=\sum_{i=1}^{\lfloor \sqrt{n} \rfloor} \mu(i) \sum_{d|i^2}\sigma_0(\frac{i^2}{d})\mu(d)g(\lfloor \frac{n}{i^2d} \rfloor)$$
可以观察到，只有$i$和$d$是squarefree number的时候才需要计算一些东西，本地跑了跑，假设需要的$g(n)$都算好了，这个式子计算量差不多是$O(\sqrt{n})$。
于是线性筛预处理前$n^{2/3}$的$f(n)$和$g(n)$，比$n^{2/3}$大的那些用$O(\sqrt{n})$的做法算，就可以$O(n^{2/3})$算出所有的$f(\lfloor \frac{n}{i} \rfloor)$和$g(\lfloor \frac{n}{i} \rfloor)$了。
附代码：还没写。写好了
update：和Min_25确认了下，他用的也是这个做法，经过一系列的优化，现在我的代码跑的比Min_25快了一点。