锐客网

  1. 首页 > it技术 > >

莫比乌斯反演|51nod 1227 平均最小公倍数 莫比乌斯反演+杜教筛

莫比乌斯反演杜教筛

题意 Lcm(a,b)表示a和b的最小公倍数，A(n)表示Lcm(n,i)的平均数(1 <= i <= n)，
例如：A(4) = (Lcm(1,4) + Lcm(2,4) + Lcm(3,4) + Lcm(4,4)) / 4 = (4 + 4 + 12 + 4) / 4 = 6。
F(a, b) = A(a) + A(a + 1) + …… A(b)。(F(a,b) = ∑A(k), a <= k <= b)
例如：F(2, 4) = A(2) + A(3) + A(4) = 2 + 4 + 6 = 12。
给出a,b，计算F(a, b)，由于结果可能很大，输出F(a, b) % 1000000007的结果即可。
1 <= a <= b <= 10^9
分析 我们要求的实际就是
∑i=1n∑j=1ijgcd(i,j)
=∑d=1n∑i=1?nd?∑j=1ij?[gcd(i,j)=1]
到这里为止就有两种方法，一种是转化为用 φ(d)来做，这样会好做很多，但我没想到，于是用了反演的做法。
按套路反演一波很容易得到 ans=12?∑k=1nμ(k)?k?∑d=1?nk?∑i=1?nkd?(i2+i)
设 s(d)=μ(d)?d,f(d)=d?(d+1)?(2d+1)6+n?(n+1)2
那么 ans=12?∑k=1ns(k)?∑d=1?nk?f(?nkd?)
设 T=kd可以得到 ans=12?∑T=1nf(?nT?)?∑d|Ts(d)
现在问题就在于如何求 g(T)=∑d|Tμ(d)?d的前缀和。
画一画不难发现 ∑d=1ng(d)=∑d=1nS(?nd?)
其中 S(d)=∑ni=1s(i)。
注意到 s(d)和 g(d)都是积性函数，可以用线性筛预处理一部分，然后其余的用杜教筛即可。 【莫比乌斯反演|51nod 1227 平均最小公倍数 莫比乌斯反演+杜教筛】
代码

#include #include #include #include #include #include using namespace std; typedef long long LL; const int MOD=1000000007; const int ny2=500000004; const int ny6=166666668; const int N=10000005; int l,r,prime[N],tot,s[N],g[N]; bool not_prime[N]; map w,wg; void get_prime(int n) { s[1]=g[1]=1; for (int i=2; i<=n; i++) { if (!not_prime[i]) prime[++tot]=i,s[i]=-i,g[i]=s[i]+s[1]; for (int j=1; j<=tot&&i*prime[j]<=n; j++) { not_prime[i*prime[j]]=1; if (i%prime[j]==0) { g[i*prime[j]]=g[i]; break; } s[i*prime[j]]=s[i]*s[prime[j]]; g[i*prime[j]]=(LL)g[i]*g[prime[j]]%MOD; } } for (int i=2; i<=n; i++) s[i]+=s[i]<0?MOD:0,s[i]+=s[i-1],s[i]-=s[i]>=MOD?MOD:0; for (int i=2; i<=n; i++) g[i]+=g[i]<0?MOD:0,g[i]+=g[i-1],g[i]-=g[i]>=MOD?MOD:0; }int get_s(int n) { if (n<=10000000) return s[n]; if (w[n]) return w[n]; int ans=1; for (int i=2,last; i<=n; i=last+1) { last=n/(n/i); (ans-=(LL)((LL)last*(last+1)/2%MOD-(LL)i*(i-1)/2%MOD)*get_s(n/i)%MOD)%=MOD; } ans+=ans<0?MOD:0; return w[n]=ans; }int get_g(int n) { if (n<=10000000) return g[n]; if (wg[n]) return wg[n]; int ans=0; for (int i=1,last; i<=n; i=last+1) { last=n/(n/i); ans+=(LL)(last-i+1)*get_s(n/i)%MOD; ans-=ans>=MOD?MOD:0; } return wg[n]=ans; }int get_f(int n) { int ans=(LL)n*(n+1)%MOD*(n*2+1)%MOD*ny6%MOD+(LL)n*(n+1)/2%MOD; ans-=ans>=MOD?MOD:0; return ans; }int solve(int n) { int ans=0; for (int i=1,last; i<=n; i=last+1) { last=n/(n/i); ans+=(LL)(get_g(last)+MOD-get_g(i-1))*get_f(n/i)%MOD; ans-=ans>=MOD?MOD:0; } ans+=ans<0?MOD:0; return (LL)ans*ny2%MOD; }int main() { get_prime(10000000); scanf("%d%d",&l,&r); printf("%d",(solve(r)-solve(l-1)+MOD)%MOD); return 0; }

    推荐阅读


    上一篇：bootstrap添加progressbar及提示

    下一篇：# 2019 GDUT Rating Contest #I E. Convention