【题意】给定n，求Σi=1~nΣj=1~n lcm(i,j)，n<=10^10。
【算法】杜教筛
【题解】就因为写了这个非常规写法，我折腾了3天……
$$ans=\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i,j)$$
令
$$g(n)=n*\sum_{i=1}^{n}\frac{i}{(n,i)}$$
那么
$$ans(n)=2*g(n)-\sum_{i=1}^{n}i$$
枚举gcd，化简g(n)。
$$g(n)=n*\sum_{d|n}1/d\sum_{i=1}^{n}i*[(n,i)=d]$$
令i=i/d
$$g(n)=n*\sum_{d|n}1/d\sum_{i=1}^{n/d}id*[(n/d,i)=1]$$
$$g(n)=n*\sum_{d|n}\sum_{i=1}^{n/d}i*[(n/d,i)=1]$$
由于
$$\sum_{i=1}^{n}[(n,i)=1]*i=\frac{n*\varphi(n)+[n==1]}{2}$$
所以代入，得
$$g(n)=n*\sum_{d|n}\frac{d*\varphi(d)+[d==1]}{2}$$
这里需要注意取整的问题，当d>1时d*φ(d)一定是偶数，当d=1时d*φ(d)=1就必须结合[d==1]，于是可以化简成下面的形式。
$$g(n)=\frac{1}{2}n(1+\sum_{d|n}d*\varphi(d))$$
$$g(n)=\frac{1}{2}(n+n*\sum_{d|n}d*\varphi(d))$$
将上式代入ans，得
$$ans=2*\sum_{i=1}^{n}\frac{1}{2}(i+i*\sum_{d|i}d*\varphi(d)))-\sum_{i=1}^{n}i$$
$$ans=\sum_{i=1}^{n}i*\sum_{d|i}d*\varphi(d)$$
令
$$f(n)=n*\sum_{d|n}d*\varphi(d)$$
那么
$$ans=F(n)=\sum_{i=1}^{n}f(n)$$
★呼，经过上面一系列的化简，我们终于来到了杜教筛——求f(n)的前缀和。
$$f(n)=n*\sum_{d|n}d*\varphi(d)$$
将f表示成狄利克雷卷积的形式，根据点积的卷积分配律（乱起的名字）。
$$f=id \cdot (1*(id \cdot \varphi))=id*(id^2 \cdot \varphi)$$
发现其中有id^2的形式，我们知道同阶幂函数卷积有奇效www。（实际上应该是相同完全积性函数卷积有奇效）
$$g=id^2$$
$$f*g=(id^2 \cdot \varphi)*id*id^2=[(id^2 \cdot \varphi)*id^2]*id$$
双重卷积非常麻烦，考虑先化简方括号内的卷积
$$([(id^2 \cdot \varphi)*id^2])(i)=\sum_{d|n}d^2*\varphi(d)*\frac{n^2}{d^2}=n^2\sum_{d|n}\varphi(d)=n^3$$
成功化简！我们可以把f*g表示出来了！
$$(f*g)(i)=\sum_{d|n}d^3*\frac{n}{d}=n*\sum_{d|n}d^2$$
jiry_2在她的博客中表示这个柿子的前缀和是非常好求的。
我：？？？
假设h=f*g，那么
$$H(i)=\sum_{i=1}^{n}i*\sum_{d|i}d^2=\sum_{i=1}^{n}i^2\sum_{d=1}^{n/i}d*i=\sum_{i=1}^{n}i^3*\frac{\frac{n}{i}*(\frac{n}{i}+1)}{2}$$
然后就可以进行分块取值优化了，注意G和H的求解复杂度为O(1)~O(√n)对杜教筛的总复杂度都没有影响。
最终杜教筛的形式是
$$F(n)=H(n)-\sum_{i=2}^{n}i^2*F(\frac{n}{i})$$
#include<cstdio> #include<cstring> #include<algorithm> #define int long long using namespace std; const int MOD=1e9+7,preN=10000000,v=(MOD+1)/2; int N,f[preN+5],phi[preN+5],e[preN+5],a[200010],prime[1000010],tot; int A(int n){return n%MOD*(n%MOD+1)%MOD*v%MOD;} int B(int n){return n%MOD*(n%MOD+1)%MOD*(2*n%MOD+1)%MOD*((MOD+1)/6)%MOD;} int K(int x){return x*x%MOD;} int C(int n){return K(A(n));} int query(int n){ int pos=1,ans=0; for(int i=1;i<=n;i=pos+1){ pos=n/(n/i); ans=(ans+(C(pos)-C(i-1)+MOD)*A(n/i)%MOD)%MOD; } return ans; } int solve(int n){ if(n<=preN)return f[n]; if(~a[N/n])return a[N/n]; int ans=query(n); int pos=2; for(int i=2;i<=n;i=pos+1){ pos=n/(n/i); ans=(ans-(B(pos)-B(i-1)+MOD)*solve(n/i)%MOD+MOD)%MOD; } return a[N/n]=ans; } #undef int int main(){ #define int long long scanf("%lld",&N); phi[1]=1;f[1]=1;e[1]=1; for(int i=2;i<=preN;i++){ if(!e[i]){ phi[prime[++tot]=i]=i-1; e[i]=i; f[i]=(i*phi[i]+1)%MOD; } for(int j=1;j<=tot&&i*prime[j]<=preN;j++){ int k=i*prime[j]; if(i%prime[j]==0){ e[k]=e[i]*prime[j]; phi[k]=phi[i]*prime[j]; f[k]=(f[i]+f[i/e[i]]*phi[e[k]]%MOD*e[k]%MOD)%MOD; break; } e[k]=prime[j]; phi[k]=phi[i]*(prime[j]-1); f[k]=f[i]*f[prime[j]]%MOD; } } for(int i=2;i<=preN;i++)f[i]=(f[i-1]+f[i]*i%MOD)%MOD; memset(a,-1,sizeof(a)); printf("%lld",solve(N)); return 0; }