c++|51nod 1227 平均最小公倍数 杜教筛

Description 求
∑i=ab1n∑j=1ilcm(i,j)∑ i = a b 1 n ∑ j = 1 i l c m ( i , j )
Solution 明天回家,现在有点浑浑噩噩
终于推出了与题解一致的柿子,so moved
照例推柿子,区间可以变成前缀和之差,枚举gcd

ans=∑d=1b∑i=1?bd?∑j=1ij?[gcd(i,j)=1]a n s = ∑ d = 1 b ∑ i = 1 ? b d ? ∑ j = 1 i j ? [ gcd ( i , j ) = 1 ]
然后后面这一坨实际上就是 ∑ni=1i?[gcd(n,i)=1]=φ(n)?n2∑ i = 1 n i ? [ g c d ( n , i ) = 1 ] = φ ( n ) ? n 2的应用,带进去就得到 ans=∑d=1b∑i=1?bd?φ(i)?i2a n s = ∑ d = 1 b ∑ i = 1 ? b d ? φ ( i ) ? i 2
需要注意的是当i=1的时候后面的贡献应当是1但是柿子体现不出来,需要我们单独加上去。这个和除二一起先不管
令 S(n)=sumni=1φ(i)?iS ( n ) = s u m i = 1 n φ ( i ) ? i,那么 ans=b+12∑d=1bS(?bd?)a n s = b + 1 2 ∑ d = 1 b S ( ? b d ? )
然后问题就变成了怎么求 S(n)S ( n ),这个是非常经典的杜教筛问题我就不写了 【c++|51nod 1227 平均最小公倍数 杜教筛】
Code

#include #include #include #define rep(i,st,ed) for (int i=st; i<=ed; ++i)typedef long long LL; const int MOD=1000000007; const int N=10000005; std:: map map; int prime[N/10]; LL f[N+5],ny6,ny2; bool not_prime[N+5]; void pre_work(int n) { f[1]=1; rep(i,2,n) { if (!not_prime[i]) { prime[++prime[0]]=i; f[i]=(i-1); } for (int j=1; i*prime[j]<=n&&j<=prime[0]; j++) { not_prime[i*prime[j]]=1; if (i%prime[j]==0) { f[i*prime[j]]=f[i]*prime[j]%MOD; break; } f[i*prime[j]]=f[i]*(prime[j]-1)%MOD; } }// f[0]=1; rep(i,1,n) f[i]=(f[i-1]+f[i]*(LL)i%MOD)%MOD; }LL ksm(LL x,LL dep) { LL ret=1; while (dep) { if (dep&1) ret=ret*x%MOD; x=x*x%MOD; dep/=2; } return ret; }LL get_f(LL n) { if (n<=N) return f[n]; if (map[n]) return map[n]; LL ret=n*(n+1)%MOD*(2*n%MOD+1)%MOD*ny6%MOD; for (LL i=2,j; i<=n; i=j+1) { j=n/(n/i); LL tmp=(j*(j+1)%MOD-(i-1)*i%MOD+MOD)%MOD*ny2%MOD; ret=(ret+MOD-tmp*get_f(n/i)%MOD)%MOD; } ret=(ret%MOD+MOD)%MOD; return map[n]=ret; }LL solve(LL n) { if (n==0) return 0; LL ret=0; for (LL i=1,j; i<=n; i=j+1) { j=n/(n/i); ret=(ret+(j-i+1)*get_f(n/i)%MOD)%MOD; } ret=(ret+n)%MOD; ret=ret*ny2%MOD; return ret; }int main(void) { ny2=ksm(2,MOD-2); ny6=ksm(6,MOD-2); pre_work(N); LL a,b; scanf("%lld%lld",&a,&b); LL ans=solve(b); ans=ans-solve(a-1); ans=(ans%MOD+MOD)%MOD; printf("%lld\n", ans); return 0; }

    推荐阅读