【题意】给定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/dsum_{i=1}^{n}i*[(n,i)=d]$$
令i=i/d
$$g(n)=n*sum_{d|n}1/dsum_{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^2sum_{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^2sum_{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; }
1.F(n)的预处理：在i%prime[j]时，f[i*prime[j]]=f[i]+f[i/e[i]]*φ(e[i]*prime[j])*(e[i]*prime[j])，其中e[i]是 i 的所有最小素因子乘积（p1^k1)。
2. 1~3次幂前缀和
$$sum_{i=1}^{n}i=frac{n(n+1)}{2}$$
$$sum_{i=1}^{n}i^2=frac{n(n+1)(2n+1)}{6}$$
$$sum_{i=1}^{n}i^3=(frac{n(n+1)}{2})^2$$
3.x的逆元是(MOD+1)/x，如果这是个整数。
4.全程long long的话，define比较好。而且输入的n是long long的话，n*n会爆long long。
【另一种写法】
$$ans=sum_{i=1}^{n}sum_{i=1}^{n}frac{i*j}{(i,j)}$$
直接枚举gcd值
$$ans=sum_{d=1}^{n}1/dsum_{i=1}^{n/d}sum_{i=1}^{n/d}d^2*i*j*[(i,j)=1]$$
$$ans=sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{n/d}i*j*[(i,j)=1]$$
这里看到[(i,j)=1]很容易想到莫比乌斯反演，但是这个问题如果反演会变得相当复杂，应该留到后面用n*φ(n)/2化简。
从上式已经可以看出分块取值优化的形式了。
令
$$s(n)=sum_{i=1}^{n}sum_{j=1}^{n}i*j*[(i,j)=1]$$
那么
$$ans=sum_{d=1}^{n}d*s(n/d)$$
为了将s(n)表示成前缀和的形式，将矩形转化为上三角。
$$s(n)=2*(sum_{i=1}^{n}sum_{j=1}^{i}i*j*[(i,j)=1])-1$$
$$s(n)=2*P(n)-1$$
$$p(n)=n*sum_{i=1}^{n}i*[(n,i)=1]$$
然后就可以转化了
$$p(n)=n*frac{n*varphi(n)+[n==1]}{2}=frac{1}{2}(n^2varphi(n)+n)$$
令
$$f(n)=n^2varphi(n)$$
那么s(n)可以表示为
$$s(n)=F(n)+frac{n(n+1)}{2}-1$$
然后我们就可以用杜教筛求解f(n)的前缀和。
$$f=id^2 cdot varphi$$
$$g=id^2$$
$$f*g=id^3$$
因为分块和杜教筛都是求n/i，所以复杂度并列，最终O(n^2/3)。