51nod1220 约数之和

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1220

$G(n)=\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{k=1}^{n^2}[k|ij]\cdot k$
$ \small{[k|ij]\to[\frac{k}{gcd(i,k)}|j]}$
$G(n)=\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{k=1}^{n^2}[\frac{k}{gcd(i,k)}|j]\cdot k$
$ \small{gcd(i,k)\to g,i/g\to i,\sum\limits_{j=1}^X[Y|j]\to \lfloor\frac XY\rfloor}$
$G(n)=\sum\limits_{k=1}^{n^2}\sum\limits_{g|k}\sum\limits_{i=1}^{n/g}[gcd(i,k/g)=1]\lfloor\frac{n}{k/g}\rfloor\cdot k$
$ \small{k/g\to k}$
$G(n)=\sum\limits_{g=1}^{n}\sum\limits_{k=1}^{n}\sum\limits_{i=1}^{n/g}[gcd(i,k)=1]\lfloor\frac{n}{k}\rfloor\cdot k\cdot g$
$\small{[X=1]\to\sum\limits_{d|X}\mu(d)}$
$G(n)=\sum\limits_{g=1}^{n}\sum\limits_{k=1}^{n}\sum\limits_{i=1}^{n/g}\sum\limits_{d|i\wedge d|k}\mu(d)\lfloor\frac{n}{k}\rfloor\cdot k\cdot g$
$\small{i/d\to i,k/d\to k}$
$G(n)=\sum\limits_{d=1}^n\sum\limits_{g=1}^{n/d}\sum\limits_{k=1}^{n/d}\sum\limits_{i=1}^{n/gd}\mu(d)\lfloor\frac{n}{kd}\rfloor\cdot k\cdot g\cdot d$
$\small{\sum\limits_{i=1}^X1\to X}$
$G(n)=\sum\limits_{d=1}^n\sum\limits_{g=1}^{n/d}\sum\limits_{k=1}^{n/d}\mu(d)\lfloor\frac{n}{gd}\rfloor\lfloor\frac{n}{kd}\rfloor\cdot k\cdot g\cdot d$
$\small{\big(\sum\limits_{i=1}^{n}i\lfloor\frac{n}{d}\rfloor\big)^2\to F(n)}$
$G(n)=\sum\limits_{d=1}^n\mu(d)d\cdot F(\lfloor\frac{n}{d}\rfloor)$
$\small{G\to \sum g,F\to \sum f}$
$g(n)=\sum \limits_{d|n}\mu(d)d\cdot f(n/d)$
$f(n)=\sum \limits_{d|n}d\cdot g(n/d)$
$\small{\sum f\to F,\sum g\to G}$
$F(n)=\sum\limits_{d=1}^nd\cdot G(\lfloor\frac{n}{d}\rfloor)$
$G(n)=F(n)-\sum\limits_{d=2}^nd\cdot G(\lfloor\frac{n}{d}\rfloor)$

记录一下推公式过程。

最后用记忆化搜索就可以$O(n^{3/4})$过了，如果加上线性预处理可以做到$O(n^{2/3})$

#include<bits/stdc++.h>
typedef unsigned long long i64;
const int P=1e9+7;
int n,vG[32007],vg[32007],B;
int G(int n){
    int&ans=n<=B?vg[n]:vG[::n/n];
    if(ans)return ans;
    i64 s=0,sx=n*2;
    for(int l=1,r,c;l<n;l=r){
        r=n/(c=n/(l+1));
        i64 t=i64(r+l+1)*(r-l);
        sx+=t*c;
        s+=t%P*G(c);
        if(s>i64(1.5e19))s%=P;
    }
    sx=sx/2%P;
    s=(sx*sx%P-s%P*((P+1)/2)%P+P)%P;
    return ans=s;
}
int main(){
    scanf("%d",&n);
    B=sqrt(n);
    printf("%d\n",G(n));
    return 0;
}