SPOJDIVCNT2: Counting Divisors(莫比乌斯反演)
http://acm.tzc.edu.cn/acmhome/vProblemList.do?method=problemdetail&oj=SPOJ&pid=DIVCNT2
给出n求
其中
是除数函数,0代表0次方.
1 #include<algorithm> 2 #include<cstdio> 3 #include<cmath> 4 #include<cstring> 5 #include<iostream> 6 #define ll long long 7 ll n,sig[100000005],N,a[200005]; 8 int p[10000005]; 9 bool mark[100000005]; 10 char mul[100000005]; 11 int smu[100000005],L; 12 int read(){ 13 int t=0,f=1;char ch=getchar(); 14 while (ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();} 15 while ('0'<=ch&&ch<='9'){t=t*10+ch-'0';ch=getchar();} 16 return t*f; 17 } 18 ll Read(){ 19 ll t=0,f=1;char ch=getchar(); 20 while (ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();} 21 while ('0'<=ch&&ch<='9'){t=t*10+ch-'0';ch=getchar();} 22 return t*f; 23 } 24 void Init(int n){ 25 sig[1]=1;smu[1]=1;mul[1]=1; 26 for (int i=2;i<=n;i++){ 27 if (!mark[i]){ 28 p[++p[0]]=i; 29 mul[i]=-1; 30 smu[i]=1; 31 sig[i]=2; 32 } 33 for (int j=1;j<=p[0]&&p[j]*i<=n;j++){ 34 mark[p[j]*i]=1; 35 if (i%p[j]==0){ 36 sig[p[j]*i]=sig[i]/(smu[i]+1)*(smu[i]+2); 37 smu[p[j]*i]=smu[i]+1; 38 mul[p[j]*i]=0; 39 break; 40 } 41 sig[p[j]*i]=sig[i]*2; 42 smu[p[j]*i]=1; 43 mul[p[j]*i]=-mul[i]; 44 } 45 } 46 for (int i=1;i<=n;i++) sig[i]+=sig[i-1],smu[i]=smu[i-1]+std::abs((int)mul[i]); 47 } 48 49 ll R(ll x){ 50 if (x<=L) return sig[x]; 51 ll res=0; 52 for (register ll i=1,j;i<=x;i=j+1){ 53 j=(x/(x/i)); 54 res+=(j-i+1)*(x/i); 55 } 56 return res; 57 } 58 ll Smu(ll x){ 59 if (x<=L) return smu[x]; 60 ll res=0; 61 for (register ll i=1;i*i<=x;i++) 62 if (mul[i]) res+=mul[i]*(x/(i*i)); 63 return res; 64 } 65 void solve(ll n){ 66 int m=sqrt(n); 67 ll ans=0; 68 ll pre=smu[m],tt; 69 for (int i=1;i<=m;i++) if (mul[i]) ans+=R(n/i); 70 for (ll i=m+1,j;i<=n;i=j+1){ 71 j=(n/(n/i)); 72 tt=Smu(j); 73 ans+=(tt-pre)*(R(n/i)); 74 pre=tt; 75 } 76 printf("%lld\n",ans); 77 } 78 int main(){ 79 int T=read(); 80 ll mx=0; 81 N=1000000000000; 82 for (int i=1;i<=T;i++) 83 a[i]=Read(),mx=std::max(mx,a[i]); 84 if (mx<=10000) L=10000;else L=pow(N,2.0/3.0); 85 Init(L); 86 for (int i=1;i<=T;i++){ 87 solve(a[i]); 88 } 89 }
