锐客网

  1. 首页 > it技术 > >

数论|Codeforces 235E Number Challenge 莫比乌斯反演+数论

数论莫比乌斯反演

题意 【数论|Codeforces 235E Number Challenge 莫比乌斯反演+数论】设d(s)表示s的约数个数，给出a,b,c，求
∑i=1a∑j=1b∑k=1cd(ijk)∑ i = 1 a ∑ j = 1 b ∑ k = 1 c d ( i j k )
a,b,c<=2000
分析 题解貌似是一个很鬼畜的做法。。。
约数个数函数 σ0(d)σ 0 ( d ) 有一个小结论，就是
σ0(ij)=∑p|i∑q|j[(p,q)=1]σ 0 ( i j ) = ∑ p | i ∑ q | j [ ( p , q ) = 1 ]
这个东西推广到三个数的成绩同样也是成立的，那么我们要求的就是 ∑(i,j)=1,(j,k)=1,(i,k)=1?ai??bj??ck?∑ ( i , j ) = 1 , ( j , k ) = 1 , ( i , k ) = 1 ? a i ? ? b j ? ? c k ?
枚举i和j，式子变成 ∑i=1a∑j=1b[(i,j)=1]?ai??bj?∑k=1c[(ij,k)=1]?ck?∑ i = 1 a ∑ j = 1 b [ ( i , j ) = 1 ] ? a i ? ? b j ? ∑ k = 1 c [ ( i j , k ) = 1 ] ? c k ?
不妨设 f(s)=∑k=1c[(s,k)=1]?ck?f ( s ) = ∑ k = 1 c [ ( s , k ) = 1 ] ? c k ?
我们只要通过反演在 O(n2logn)O ( n 2 l o g n )把所有 f(s)f ( s )求出来，然后就可以在 O(n2logn)O ( n 2 l o g n )的复杂度内求答案了。
代码

#include #include #include #include #include using namespace std; typedef long long LL; const int N=2005; const int MOD=1073741824; int a,b,c,n,tot,prime[N*N],mu[N*N],f[N],s[N*N]; bool not_prime[N*N]; void mod(int &x) { x-=x>=MOD?MOD:0; }int gcd(int x,int y) { if (!y) return x; else return gcd(y,x%y); }void get_prime(int n) { mu[1]=1; for (int i=2; i<=n; i++) { if (!not_prime[i]) prime[++tot]=i,mu[i]=-1; for (int j=1; j<=tot&&i*prime[j]<=n; j++) { not_prime[i*prime[j]]=1; if (i%prime[j]==0) break; mu[i*prime[j]]=-mu[i]; } } }int main() { scanf("%d%d%d",&a,&b,&c); get_prime(a*b); for (int d=1; d<=c; d++) { int n=c/d; for (int i=1,last; i<=n; i=last+1) { last=n/(n/i); mod(f[d]+=(n/i)*(last-i+1)); } } for (int i=1; i<=c; i++) { if (!mu[i]) continue; int w=mu[i]*f[i]; w+=w<0?MOD:0; for (int j=i; j<=a*b; j+=i) mod(s[j]+=w); } int ans=0; for (int i=1; i<=a; i++) for (int j=1; j<=b; j++) if (gcd(i,j)==1) mod(ans+=(LL)(a/i)*(b/j)*s[i*j]&(MOD-1)); printf("%d",ans); return 0; }

    推荐阅读


    上一篇：树形dp|Codeforces 917D Stranger Trees 树形dp+容斥原理

    下一篇：GDOI2018游记