acm-(数论、莫比乌斯反演、杜教筛)51nod 1227 平均最小公倍数

传送门 设函数$F(n)={\sum }_{k=1}^n\frac{1}{k}{\sum }_{i=1}^{k}lcm(i,k)$，题目就是让求$F(b)-F(a-1)$，于是对$F(n)$先套路反演一波，得到$F(n)=\frac{n+{\sum }_{i=1}^{n}f(i)}{2}$($f(n)={\sum }_{d\mid n}d\phi (d)$)，设$S(n)={\sum }_{i=1}^{n}f(i)$，看上去就像杜教筛的样子，由于$f=(id\cdot \phi )\ast I$，故我们可以设$g=f\ast id=((\phi \ast id)\ast id)\ast I={\sigma }_2$，又由于${\sum }_{i=1}^{n}g(i)={\sum }_{i=1}^n{i}^2\lfloor \frac{n}{i}\rfloor$，显然可以分块$O(\sqrt{n})$求解，不过我们可以预处理出$g$函数的$n^{\frac{2}{3}}$的前缀和，这样能加快速度，然后在此基础上套一个杜教筛式子，我们有$S(n)={\sum }_{i=1}^{n}g(i)-{\sum }_{i=2}^{n}iS(\lfloor \frac{n}{i}\rfloor )$，由于计算${\sum }_{i=1}^{n}g(i)$比较快，总复杂度应该还是$O(n^{\frac{2}{3}})$，注意$f$函数也要预处理一下前缀和，否则杜教筛复杂度会退化。

int nn,N,limt,prim[maxn],tot,s[maxn],e[maxn],f[maxn],g[maxm],dp[maxm],ss[maxn], dps[maxm],es[maxn],fs[maxn]; bool flag[maxn]; int id(int x){ return x<=limt?x:nn/x+limt; } void init(int n){ s[1]=1,ss[1]=1; FOR(i,2,n+1){ if(!flag[i])prim[tot++]=i,e[i]=1ll*i*(i-1)%mod,s[i]=e[i]+1,f[i]=1, es[i]=1ll*i*i%mod,ss[i]=es[i]+1,fs[i]=1; for(register int j=0;j<tot && prim[j]*i<=n;++j){ flag[i*prim[j]]=1; if(i%prim[j]==0){ e[i*prim[j]]=1ll*e[i]*sqr(prim[j])%mod; f[i*prim[j]]=f[i]; s[i*prim[j]]=(s[i]+1ll*e[i*prim[j]]*f[i]%mod)%mod; es[i*prim[j]]=1ll*es[i]*sqr(prim[j])%mod; fs[i*prim[j]]=fs[i]; ss[i*prim[j]]=(ss[i]+1ll*es[i*prim[j]]*fs[i]%mod)%mod; break; } e[i*prim[j]]=1ll*prim[j]*(prim[j]-1)%mod; f[i*prim[j]]=s[i]; s[i*prim[j]]=1ll*s[i]*(e[i*prim[j]]+1)%mod; es[i*prim[j]]=sqr(prim[j]); fs[i*prim[j]]=ss[i]; ss[i*prim[j]]=1ll*ss[i]*(es[i*prim[j]]+1)%mod; } } FOR(i,1,n+1)add(s[i],s[i-1]),add(ss[i],ss[i-1]); } int SS(int n){ if(n<=N)return ss[n]; register int idd=id(n),i=1,j,ans=0,pre=0,cur=0; if(dps[idd])return dps[idd]; while(i<=n){ j=n/(n/i); add(ans,1ll*((cur=sm2(j))-pre+mod)*(n/i)%mod); pre=cur; i=j+1; } return dps[idd]=ans; } int S(int n){ if(n<=N)return s[n]; register int idd=id(n),ans=SS(n),i=2,j,cur=0,pre=1; if(dp[idd])return dp[idd]; while(i<=n){ j=n/(n/i); dec(ans,1ll*((cur=sm1(j))-pre+mod)*S(n/i)%mod); i=j+1; pre=cur; } return dp[idd]=ans; } int cal(int n){ if(n<=0)return 0; nn=n; limt=sqrt(nn); register int i =1,j,ans=0,c=nn/limt+limt; FOR(i,1,c+1)dps[i]=dp[i]=0; return 1ll*(S(n)+n)*qpow(2,mod-2,mod)%mod; } int main(){ int a,b; rd(&a,&b); N=pow(b,2.0/3.0); init(N=max(N,10000)); wrn((cal(b)-cal(a-1)+mod)%mod); }