【算法总结】积性函数相关

【线性筛】

〖模板代码〗

[线性筛质数]

 1 int main()
 2 { 
 3     n=read();
 4     for(int i=2;i<=n;i++)
 5     {
 6         if(!f[i])pri[++cnt]=i;
 7         for(int j=1;j<=cnt;j++)
 8         {
 9             if(i*pri[j]>n)break;
10             f[i*pri[j]]=1;
11             if(i%pri[j]==0)break;
12         }
13     }
14 }

View Code

[线性求约数和]

 1 void pre()
 2 {
 3     miu[1]=1;f[1].first=1;
 4     for(int i=2;i<=mx;i++)
 5     {
 6         if(!np[i]){pri[++cnt]=i;d[i]=i;f[i].first=i+1;miu[i]=-1;}
 7         for(int j=1;i*pri[j]<=mx;j++)
 8         {
 9             np[i*pri[j]]=true;
10             if(i%pri[j]==0)
11             {
12                 miu[i*pri[j]]=0;
13                 d[i*pri[j]]=d[i]*pri[j];
14                 f[i*pri[j]].first=f[i].first+f[i/d[i]].first*d[i*pri[j]];
15                 break;
16             }
17             miu[i*pri[j]]=-miu[i];
18             d[i*pri[j]]=pri[j];
19             f[i*pri[j]].first=f[i].first*f[pri[j]].first;
20         }
21     }
22     for(int i=1;i<=mx;i++)f[i].second=i;
23 }

View Code

〖相关题目〗

1.【Technocup 2018 - Elimination Round 2】F. Paths

题意：给定数字n，建立一个无向图。对于所有1~n之间的数字，当数字gcd(u,v)≠1时将u、v连一条边，边权为1。d(u, v)表示u到v的最短路，求所有d(u, v)的和，其中1 ≤ u < v ≤ n。

分析：Zsnuoの博客

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #define LL long long
 5 using namespace std;
 6 const int N=1e7+5;
 7 int n,m,tot,now,pri[N],p[N],phi[N],num[N],sum[N]; 
 8 LL one,two,three;
 9 int main()
10 {
11     scanf("%d",&n);
12     phi[1]=1;
13     for(int i=2;i<=n;i++)
14     {
15         if(!p[i]){p[i]=pri[++tot]=i;phi[i]=i-1;}
16         for(int j=1;j<=tot;j++)
17         {
18             if(i*pri[j]>n)break;
19             p[i*pri[j]]=pri[j];
20             if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
21             else phi[i*pri[j]]=phi[i]*(pri[j]-1);
22         }
23     } 
24     for(int i=2;i<=n;i++)one+=i-1-phi[i];
25     for(int i=2;i<=n;i++)num[p[i]]++;
26     for(int i=2;i<=n;i++)sum[i]=sum[i-1]+num[i];
27     for(int i=2;i<=n;i++)two+=1ll*num[i]*sum[n/i];
28     for(int i=2;i<=n;i++)if(1ll*p[i]*p[i]<=n)two--;
29     two=two/2-one;m=n-1;
30     for(int i=tot;i>=1;i--)
31         if(pri[i]*2>n)m--;
32         else break;
33     three=1ll*m*(m-1)/2-one-two;
34     printf("%I64d
",one+two*2+three*3);
35     return 0;
36 }

View Code

2.【51nod1584】加权约数和

题意：见原题

分析：KsClaの博客

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=1e6;
 8 const int mod=1e9+7;
 9 int T,n,tot,now,pri[N+5],miu[N+5];
10 int mn[N+5],h[N+5],p[N+5],g[N+5],f[N+5];
11 bool np[N+5];
12 int read()
13 {
14     int x=0,f=1;char c=getchar();
15     while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
16     while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
17     return x*f;
18 }
19 int main()
20 {
21     h[1]=p[1]=miu[1]=g[1]=1;
22     for(int i=2;i<=N;i++)
23     {
24         if(!np[i])
25         {
26             pri[++tot]=i;
27             miu[i]=-1;mn[i]=i;h[i]=i+1;
28             p[i]=(1ll*i*i%mod+i+1)%mod;
29         }
30         for(int j=1;i*pri[j]<=N;j++)
31         {
32             now=i*pri[j];np[now]=true;
33             if(i%pri[j]==0)
34             {
35                 miu[now]=0;mn[now]=mn[i]*pri[j];
36                 h[now]=(h[i]+1ll*h[i/mn[i]]*mn[now]%mod)%mod;
37                 p[now]=(1ll*(mn[now]+mn[i])%mod*mn[now]%mod*p[i/mn[i]]%mod+p[i])%mod;
38                 break;
39             }
40             miu[now]=-miu[i];mn[now]=pri[j];
41             h[now]=1ll*h[i]*h[pri[j]]%mod;
42             p[now]=1ll*p[i]*p[pri[j]]%mod;
43         }
44     }
45     for(int i=2;i<=N;i++)p[i]=1ll*p[i]*i%mod;
46     for(int i=2;i<=N;i++)p[i]=(p[i]+p[i-1])%mod;
47     for(int i=2;i<=N;i++)miu[i]=1ll*(miu[i]+mod)%mod*i%mod*i%mod;
48     for(int i=2;i<=N;i++)g[i]=(g[i-1]+h[i])%mod;
49     for(int i=2;i<=N;i++)g[i]=1ll*i*g[i]%mod*h[i]%mod;
50     for(int i=1;i<=N;i++)
51         for(int j=1;i*j<=N;j++)
52             f[i*j]=(f[i*j]+1ll*miu[i]*g[j]%mod)%mod;
53     for(int i=2;i<=N;i++)f[i]=(f[i-1]+f[i])%mod;
54     T=read();
55     for(int i=1;i<=T;i++)
56     {
57         n=read();
58         n=(2*f[n]%mod-p[n]+mod)%mod;
59         printf("Case #%d: %d
",i,n);
60     }
61     return 0;
62 }

View Code

【杜教筛】

〖相关资料〗

《浅谈一类积性函数的前缀和》《论逗逼的自我修养之寒假颓废记》

〖相关题目〗

1.【51nod1239】欧拉函数之和

题意：求φ的前缀和。

分析：ONION_CYCの博客

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=5e6;
 8 const int mod=1e9+7;
 9 const int inv=(mod+1)/2;
10 int s,tot,pri[N+5],phi[N+5],a[N+5],b[N+5];
11 bool isp[N+5];
12 LL n;
13 int memory(LL x)
14 {
15     if(x<=s)return a[x];
16     else return b[n/x];
17 }
18 void add(LL x,int ans)
19 {
20     if(x<=s)a[x]=ans;
21     else b[n/x]=ans;
22 }
23 int solve(LL n)
24 {
25     if(n<=N)return phi[n];
26     if(~(memory(n)))return memory(n);
27     int ans=1ll*(n%mod)*((n+1)%mod)%mod*inv%mod;
28     LL pos;
29     for(LL i=2;i<=n;i=pos+1)
30     {
31         pos=n/(n/i);
32         ans=(ans-1ll*(pos-i+1)%mod*solve(n/i)%mod+mod)%mod;
33     }
34     add(n,ans);
35     return ans;
36 }
37 int main()
38 {
39     memset(a,-1,sizeof(a));
40     memset(b,-1,sizeof(b));
41     scanf("%lld",&n);
42     s=sqrt(n);phi[1]=1;
43     for(int i=2;i<=N;i++)
44     {
45         if(!isp[i]){pri[++tot]=i;phi[i]=i-1;}
46         for(int j=1;j<=tot;j++)
47         {
48             if(i*pri[j]>N)break;
49             isp[i*pri[j]]=true;
50             if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
51             else phi[i*pri[j]]=phi[i]*(pri[j]-1);
52         }
53         phi[i]=(phi[i]+phi[i-1])%mod;
54     }
55     if(n<=N)printf("%d",phi[n]);
56     else printf("%d",solve(n));
57     return 0;
58 }

View Code

2.【51nod1244】莫比乌斯函数之和

题意：求μ的前缀和。

分析：无

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=5e6;
 8 int s,tot,pri[N+5],miu[N+5],a[N+5];
 9 bool np[N+5];
10 LL nn,L,R;
11 LL solve(LL n)
12 {
13     if(n<=N)return miu[n];
14     if(~a[nn/n])return a[nn/n];
15     LL ans=1,pos;
16     for(LL i=2;i<=n;i=pos+1)
17     {
18         pos=n/(n/i);
19         ans=ans-(pos-i+1)*solve(n/i);
20     }
21     return a[nn/n]=ans;
22 }
23 int main()
24 {
25     miu[1]=1;
26     for(int i=2;i<=N;i++)
27     {
28         if(!np[i]){miu[i]=-1;pri[++tot]=i;}
29         for(int j=1;i*pri[j]<=N;j++)
30         {
31             np[i*pri[j]]=true;
32             if(i%pri[j]==0){miu[i*pri[j]]=0;break;}
33             miu[i*pri[j]]=-miu[i];
34         }
35         miu[i]+=miu[i-1];
36     }
37     scanf("%lld%lld",&L,&R);
38     nn=L-1;memset(a,-1,sizeof(a));L=solve(L-1);
39     nn=R;memset(a,-1,sizeof(a));R=solve(R);
40     printf("%lld",R-L);
41     return 0;
42 }

View Code

3.【bzoj3944】Sum

题意：求φ和μ的前缀和。

分析：无

 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=2e6;
 8 LL T,n,nn,tot,pri[N+5];
 9 LL miu[N+5],phi[N+5],mmiu[N+5],mphi[N+5];
10 bool np[N+5];
11 struct node{LL phi,miu;}ans,tmp;
12 node solve(LL n)
13 {
14     if(n<=N)return (node){phi[n],miu[n]};
15     if(~mphi[nn/n])return (node){mphi[nn/n],mmiu[nn/n]};
16     LL ansp=1ll*n*(n+1)/2,ansm=1,pos;
17     for(LL i=2;i<=n;i=pos+1)
18     {
19         pos=n/(n/i);tmp=solve(n/i);
20         ansp=ansp-(pos-i+1)*tmp.phi;
21         ansm=ansm-(pos-i+1)*tmp.miu;
22     }
23     mphi[nn/n]=ansp;mmiu[nn/n]=ansm;
24     return (node){ansp,ansm};
25 }
26 void work()
27 {
28     scanf("%lld",&n);nn=n;
29     memset(mphi,-1,sizeof(mphi));
30     ans=solve(n);
31     printf("%lld %lld
",ans.phi,ans.miu);
32 }
33 int main()
34 {
35     phi[1]=miu[1]=1;
36     for(int i=2;i<=N;i++)
37     {
38         if(!np[i]){pri[++tot]=i;miu[i]=-1;phi[i]=i-1;}
39         for(int j=1;i*pri[j]<=N;j++)
40         {
41             np[i*pri[j]]=true;
42             if(i%pri[j]==0)
43             {
44                 miu[i*pri[j]]=0;
45                 phi[i*pri[j]]=phi[i]*pri[j];
46                 break;
47             }
48             miu[i*pri[j]]=-miu[i];
49             phi[i*pri[j]]=phi[i]*(pri[j]-1);
50         }
51         miu[i]+=miu[i-1];
52         phi[i]+=phi[i-1];
53     }
54     scanf("%lld",&T);
55     while(T--)work();
56     return 0;
57 }

View Code

4.【51nod1238】最小公倍数之和V3

题意：输出小于等于N的所有数，两两之间的最小公倍数之和。

分析：jiry_2の博客

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=5e6;
 8 const int mod=1e9+7;
 9 const int inv2=(mod+1)/2;
10 const int inv6=(mod+1)/6;
11 int tot,now,pri[N+5],mn[N+5];
12 bool np[N+5];
13 LL n,nn,sum,phi[N+5],a[N+5];
14 LL calc1(LL a)
15 {
16     LL suma=a%mod;
17     return suma*(suma+1)%mod*inv2%mod;
18 }
19 LL calc2(LL a,LL b)
20 {
21     LL suma=(a-1)%mod;
22     suma=suma*(suma+1)%mod*((2*suma+1)%mod)%mod*inv6%mod;
23     LL sumb=b%mod;
24     sumb=sumb*(sumb+1)%mod*((2*sumb+1)%mod)%mod*inv6%mod;
25     return (sumb-suma+mod)%mod;
26 }
27 LL calc3(LL a,LL b)
28 {
29     LL suma=(a-1)%mod;
30     suma=suma*(suma+1)%mod*inv2%mod;
31     suma=suma*suma%mod;
32     LL sumb=b%mod;
33     sumb=sumb*(sumb+1)%mod*inv2%mod;
34     sumb=sumb*sumb%mod;
35     return (sumb-suma+mod)%mod;
36 }
37 LL solve(LL n)
38 {
39     if(n<=N)return phi[n];
40     if(~a[nn/n])return a[nn/n];
41     LL ans=0,pos;
42     for(LL i=1;i<=n;i=pos+1)
43     {
44         pos=n/(n/i);
45         ans=(ans+calc1(n/i)*calc3(i,pos))%mod;
46     }
47     for(LL i=2;i<=n;i=pos+1)
48     {
49         pos=n/(n/i);
50         ans=(ans-solve(n/i)*calc2(i,pos)%mod+mod)%mod;
51     }
52     return a[nn/n]=ans;
53 }
54 int main()
55 {
56     scanf("%lld",&n);
57     nn=n;phi[1]=1;
58     for(int i=2;i<=N;i++)
59     {
60         if(!np[i])
61         {
62             pri[++tot]=i;
63             phi[i]=(1ll*i*(i-1)%mod+1)%mod;
64             mn[i]=i;
65         }
66         for(int j=1;i*pri[j]<=N;j++)
67         {
68             now=i*pri[j];
69             np[now]=true;
70             if(i%pri[j]==0)
71             {
72                 mn[now]=mn[i]*pri[j];
73                 phi[now]=(phi[i]+phi[i/mn[i]]*mn[now]%mod*(mn[i]*(pri[j]-1)%mod))%mod;
74                 break;
75             }
76             phi[now]=phi[i]*(1ll*pri[j]*(pri[j]-1)%mod+1)%mod;
77             mn[now]=pri[j];
78         }
79     }
80     for(int i=1;i<=N;i++)phi[i]=phi[i]*i%mod;
81     for(int i=2;i<=N;i++)phi[i]=(phi[i]+phi[i-1])%mod;
82     memset(a,-1,sizeof(a));
83     printf("%lld
",solve(n));
84     return 0;
85 }

View Code

 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=5e6;
 8 const int mod=1e9+7;
 9 const int inv2=(mod+1)/2;
10 const int inv6=(mod+1)/6;
11 int tot,pri[N+5];
12 LL n,nn;
13 LL phi[N+5],a[N+5];
14 bool np[N+5];
15 LL calc(LL a,LL b)
16 {
17     LL suma=(a-1)%mod;
18     suma=suma*(suma+1)%mod*((2*suma+1)%mod)%mod*inv6%mod;
19     LL sumb=b%mod;
20     sumb=sumb*(sumb+1)%mod*((2*sumb+1)%mod)%mod*inv6%mod;
21     return (sumb-suma+mod)%mod;
22 }
23 LL solve(LL n)
24 {
25     if(n<=N)return phi[n];
26     if(~a[nn/n])return a[nn/n];
27     LL ans=(n%mod)*((n+1)%mod)%mod*inv2%mod,pos;
28     ans=ans*ans%mod;
29     for(LL i=2;i<=n;i=pos+1)
30     {
31         pos=n/(n/i);
32         ans=(ans-solve(n/i)*calc(i,pos)%mod+mod)%mod;
33     }
34     return a[nn/n]=ans;
35 }
36 int main()
37 {
38     scanf("%lld",&n);
39     nn=n;phi[1]=1;
40     for(int i=2;i<=N;i++)
41     {
42         if(!np[i]){pri[++tot]=i;phi[i]=i-1;}
43         for(int j=1;i*pri[j]<=N;j++)
44         {
45             np[i*pri[j]]=true;
46             if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
47             phi[i*pri[j]]=phi[i]*(pri[j]-1);
48         }
49         phi[i]=phi[i]*i%mod*i%mod;
50         phi[i]=(phi[i]+phi[i-1])%mod;
51     }
52     memset(a,-1,sizeof(a));
53     LL pos,ans=0,sum;
54     for(LL i=1;i<=n;i=pos+1)
55     {
56         pos=n/(n/i);
57         sum=(solve(pos)-solve(i-1)+mod)%mod;
58 //        printf("   %lld %lld
",pos,solve(pos));
59         sum=sum*(((n/i)%mod)*((n/i+1)%mod)%mod*inv2%mod)%mod;
60         ans=(ans+sum)%mod;
61     }
62     printf("%lld
",ans);
63 //    for(int i=1;i<=n;i++)printf("%lld
",phi[i]);
64     return 0;
65 }

View Code

5.【51nod1237】最大公约数之和

题意：输出小于等于N的所有数，两两之间的最大公约数之和。

分析：无

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=5e6;
 8 const int mod=1e9+7;
 9 const int inv=(mod+1)/2;
10 int tot,pri[N+5];
11 LL n,nn,ans,pos,phi[N+5],a[N+5];
12 bool np[N+5];
13 LL calc(LL a,LL b)
14 {
15     LL suma=(a-1)%mod;
16     suma=suma*(suma+1)%mod*inv%mod;
17     LL sumb=b%mod;
18     sumb=sumb*(sumb+1)%mod*inv%mod;
19     return (sumb-suma+mod)%mod;
20 }
21 LL solve(LL n)
22 {
23     if(n<=N)return phi[n];
24     if(~a[nn/n])return a[nn/n];
25     LL ans=(n%mod)*((n+1)%mod)%mod*inv%mod,pos;
26     for(LL i=2;i<=n;i=pos+1)
27     {
28         pos=n/(n/i);
29         ans=(ans-(pos-i+1)%mod*solve(n/i)%mod+mod)%mod;
30     }
31     return a[nn/n]=ans;
32 }
33 int main()
34 {
35     phi[1]=1;
36     for(int i=2;i<=N;i++)
37     {
38         if(!np[i]){pri[++tot]=i;phi[i]=i-1;}
39         for(int j=1;i*pri[j]<=N;j++)
40         {
41             np[i*pri[j]]=true;
42             if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
43             phi[i*pri[j]]=phi[i]*(pri[j]-1);
44         }
45         phi[i]=(phi[i]+phi[i-1])%mod;
46     }
47     scanf("%lld",&n);nn=n;
48     memset(a,-1,sizeof(a));
49     for(LL i=1;i<=n;i=pos+1)
50     {
51         pos=n/(n/i);
52         ans=(ans+calc(i,pos)*solve(n/i)%mod)%mod;
53     }
54     n%=mod;n=n*(n+1)%mod*inv%mod;
55     ans=(ans*2%mod-n+mod)%mod;
56     printf("%lld",ans);
57     return 0;
58 }

View Code

6.【51nod1227】平均最小公倍数

题意：见原题

分析：无

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=1e6;
 8 const int mod=1e9+7;
 9 const int inv2=(mod+1)/2;
10 const int inv6=(mod+1)/6;
11 int tot,pri[N+5];
12 LL a,b,nn,pos,sum,phi[N+5],m[N+5];
13 bool np[N+5];
14 LL calc(LL a,LL b)
15 {
16     LL suma=(a-1)%mod;
17     suma=suma*(suma+1)%mod*inv2%mod;
18     LL sumb=b%mod;
19     sumb=sumb*(sumb+1)%mod*inv2%mod;
20     return (sumb-suma+mod)%mod;
21 }
22 LL solve(LL n)
23 {
24     if(n<=N)return phi[n];
25     if(~m[nn/n])return m[nn/n];
26     LL ans=n%mod,pos;
27     ans=ans*(ans+1)%mod*((ans*2+1)%mod)%mod*inv6%mod;
28     for(LL i=2;i<=n;i=pos+1)
29     {
30         pos=n/(n/i);
31         sum=calc(i,pos)*solve(n/i)%mod;
32         ans=(ans-sum+mod)%mod;
33     }
34     return m[nn/n]=ans;
35 }
36 LL work(LL n)
37 {
38     LL ans=0;
39     for(LL i=1;i<=n;i=pos+1)
40     {
41         pos=n/(n/i);
42         sum=((pos-i+1)%mod)*solve(n/i)%mod;
43         ans=(ans+sum)%mod;
44 //        printf("%lld
",ans);
45     }
46     return (ans+n)%mod*inv2%mod;
47 }
48 int main()
49 {
50     phi[1]=1;
51     for(int i=2;i<=N;i++)
52     {
53         if(!np[i]){pri[++tot]=i;phi[i]=i-1;}
54         for(int j=1;i*pri[j]<=N;j++)
55         {
56             np[i*pri[j]]=true;
57             if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
58             phi[i*pri[j]]=phi[i]*(pri[j]-1);
59         }
60         phi[i]=phi[i]*i%mod;
61         phi[i]=(phi[i]+phi[i-1])%mod;
62     }
63     scanf("%lld%lld",&a,&b);
64     memset(m,-1,sizeof(m));nn=a;a=work(a-1);
65     memset(m,-1,sizeof(m));nn=b;b=work(b);
66     printf("%lld",(b-a+mod)%mod);
67     return 0;
68 }

View Code

7.【bzoj4176】Lucas的数论

题意：见原题

分析：PoPoQQQの博客

 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=6e6;
 8 const int mod=1e9+7;
 9 int S,tot,pri[N+5],miu[N+5],a[N+5];
10 int n,nn,pos,n2,pos2,sum,ans;
11 bool np[N+5];
12 int findmiu(int n)
13 {
14     if(n<=S)return miu[n];
15     if(~a[nn/n])return a[nn/n];
16     int ans=1,pos;
17     for(int i=2;i<=n;i=pos+1)
18     {
19         pos=n/(n/i);
20         ans=(ans-(LL)(pos-i+1)*findmiu(n/i)%mod+mod)%mod;
21     }
22     return a[nn/n]=ans;
23 }
24 int solve(int a,int b)
25 {
26     a=findmiu(a-1);b=findmiu(b);
27     return (b-a+mod)%mod;
28 }
29 int main()
30 {
31     scanf("%d",&n);
32     S=ceil(pow(n,0.75)-1e-7)+1e-7;
33     miu[1]=1;
34     for(int i=2;i<=S;i++)
35     {
36         if(!np[i]){pri[++tot]=i;miu[i]=-1;}
37         for(int j=1;i*pri[j]<=S;j++)
38         {
39             np[i*pri[j]]=true;
40             if(i%pri[j]==0){miu[i*pri[j]]=0;break;}
41             miu[i*pri[j]]=-miu[i];
42         }
43         miu[i]+=miu[i-1];
44     }
45     nn=n;memset(a,-1,sizeof(a));
46     for(int i=1;i<=n;i=pos+1)
47     {
48         pos=n/(n/i);n2=n/i;sum=0;
49         for(int j=1;j<=n2;j=pos2+1)
50         {
51             pos2=n2/(n2/j);
52             sum=(sum+(LL)(pos2-j+1)*(n2/j)%mod)%mod;
53         }
54         sum=(LL)sum*sum%mod;
55         ans=(ans+(LL)solve(i,pos)*sum%mod)%mod;
56     }
57     printf("%d",ans);
58     return 0;
59 }

View Code

8.【51nod1220】约数之和

题意：见原题

分析：无。

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=6e6;
 8 const int mod=1e9+7;
 9 int n,S,nn,tot,pos,ans;
10 int pri[N+5],miu[N+5],a[N+5];
11 bool np[N+5];
12 void pre()
13 {
14     S=ceil(pow(n,0.75)-1e-7)+1e-7;
15     miu[1]=1;
16     for(int i=2;i<=S;i++)
17     {
18         if(!np[i]){pri[++tot]=i;miu[i]=-1;}
19         for(int j=1;i*pri[j]<=S;j++)
20         {
21             np[i*pri[j]]=true;
22             if(i%pri[j]==0){miu[i*pri[j]]=0;break;}
23             miu[i*pri[j]]=-miu[i];
24         }
25         miu[i]=1ll*miu[i]*i%mod;
26         miu[i]=(miu[i-1]+miu[i])%mod;
27         miu[i]=(miu[i]+mod)%mod;
28     }
29 }
30 int calc(int a,int b)
31 {
32     a=1ll*a*(a+1)/2%mod;
33     b=1ll*b*(b+1)/2%mod;
34     return (b-a+mod)%mod;
35 }
36 int find(int n)
37 {
38     if(n<=S)return miu[n];
39     if(~a[nn/n])return a[nn/n];
40     int ans=1,pos,sum;
41     for(int i=2;i<=n;i=pos+1)
42     {
43         pos=n/(n/i);
44         sum=1ll*calc(i-1,pos)*find(n/i)%mod;
45         ans=(ans-sum+mod)%mod;
46     }
47     return a[nn/n]=ans;
48 }
49 int solve(int a,int b)
50 {
51     a=find(a-1);b=find(b);
52     return (b-a+mod)%mod;
53 }
54 int work(int n)
55 {
56     int ans=0,pos;
57     for(int i=1;i<=n;i=pos+1)
58     {
59         pos=n/(n/i);
60         ans=(ans+1ll*calc(i-1,pos)*(n/i)%mod)%mod;
61     }
62     return 1ll*ans*ans%mod;
63 }
64 int main()
65 {
66     scanf("%d",&n);
67     pre();nn=n;memset(a,-1,sizeof(a));
68     for(int i=1;i<=n;i=pos+1)
69     {
70         pos=n/(n/i);
71         ans=(ans+1ll*solve(i,pos)*work(n/i)%mod)%mod;
72     }
73     printf("%d",ans);
74     return 0;
75 }

View Code

9.【bzoj3512】DZY Loves Math IV

题意：见原题

分析：permuiの博客

 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<cmath>
 5 #include<map>
 6 #define LL long long
 7 using namespace std;
 8 const int N=1e6;
 9 const int mod=1e9+7;
10 const int inv=(mod+1)/2;
11 int n,m,ans,tot,now;
12 int pri[N+5],phi[N+5],mn[N+5],a[N+5];
13 bool np[N+5];
14 map<int,int>M[N+5];
15 int power(int a,int b)
16 {
17     int ans=1;
18     while(b)
19     {
20         if(b&1)ans=1ll*ans*a%mod;
21         a=1ll*a*a%mod;b>>=1;
22     }
23     return ans;
24 }
25 int getsum(int n)
26 {
27     if(n<=N)return phi[n];
28     if(~a[m/n])return a[m/n];
29     int ans=1ll*n*(n+1)%mod*inv%mod,pos;
30     for(int i=2;i<=n;i=pos+1)
31     {
32         pos=n/(n/i);
33         ans=(ans-1ll*(pos-i+1)*getsum(n/i)%mod+mod)%mod;
34     }
35     return a[m/n]=ans;
36 }
37 int Phi(int n)
38 {
39     int ans=1,num;
40     if(n<=N){ans=(phi[n]-phi[n-1]+mod)%mod;return ans;}
41     int m=sqrt(n)+0.5;
42     for(int i=1;pri[i]<=m;i++)
43     {
44         if(n%pri[i])continue;
45         num=0;while(n%pri[i]==0)num++,n/=pri[i];
46         ans*=power(pri[i],num-1)*(pri[i]-1);
47     }
48     if(n!=1)ans*=(n-1);
49     return ans;
50 }
51 int S(int n,int m)
52 {
53     if(m==0)return 0;
54     if(n==1)return getsum(m);
55     if(M[n][m])return M[n][m];
56     int ans=0;
57     for(int i=1;i*i<=n;i++)
58     {
59         if(n%i)continue;
60         int j=n/i;ans=(ans+1ll*Phi(j)*S(i,m/i)%mod)%mod;
61         if(i!=j)ans=(ans+1ll*Phi(i)*S(j,m/j)%mod)%mod;
62     }
63     return M[n][m]=ans;
64 }
65 int main()
66 {
67     phi[1]=mn[1]=1;
68     for(int i=2;i<=N;i++)
69     {
70         if(!np[i])pri[++tot]=i,phi[i]=i-1,mn[i]=i;
71         for(int j=1;i*pri[j]<=N;j++)
72         {
73             now=i*pri[j];np[now]=true;
74             if(i%pri[j]==0)
75             {
76                 phi[now]=phi[i]*pri[j];
77                 mn[now]=mn[i];break;
78             }
79             phi[now]=phi[i]*(pri[j]-1);
80             mn[now]=mn[i]*pri[j];
81         }
82         phi[i]=(phi[i-1]+phi[i])%mod;
83     }
84     memset(a,-1,sizeof(a));
85     scanf("%d%d",&n,&m);
86     for(int i=1;i<=n;i++)
87         ans=(ans+1ll*i/mn[i]*S(mn[i],m)%mod)%mod;
88     printf("%d",ans);
89     return 0;
90 }

View Code

【洲阁筛】

〖相关资料〗

《洲阁筛学习》

〖相关题目〗

1.【spoj】DIVCNT3 - Counting Divisors (cube)

题意：见原题

分析：Monster_Yiの博客

  1 #include<cstdio>
  2 #include<cstring>
  3 #include<cmath>
  4 #include<cstdlib>
  5 #include<algorithm>
  6 #define LL long long
  7 using namespace std;
  8 const int N=320000;
  9 const int mx=316240;
 10 const int p=1000003;
 11 int tot,id,pri[N],num[N],ci[N],mn[N];
 12 int cnt,first[p],w[p],nex[p],D[p];
 13 LL T,n,ans1,f[N],f2[p],sum[N],to[p],d[p],g[p];
 14 LL read()
 15 {
 16     LL x=0,f=1;char c=getchar();
 17     while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
 18     while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
 19     return x*f;
 20 }
 21 void ins(int x,LL y)
 22 {
 23     int id=y%p;to[++cnt]=y;w[cnt]=x;
 24     nex[cnt]=first[id];first[id]=cnt;
 25 }
 26 int idx(LL x)
 27 {
 28     for(int i=first[x%p];i;i=nex[i])
 29         if(to[i]==x)return w[i];
 30     return 0;
 31 }
 32 void work1()
 33 {
 34     cnt=id=0;int k;
 35     for(LL i=1;i<=n;i=n/(n/i)+1)first[n/i%p]=0;
 36     for(LL i=1;i<=n;i=n/(n/i)+1)
 37         ins(++id,n/i),d[id]=g[id]=n/i,D[id]=0;
 38     for(int i=1;i<=tot;i++)
 39         for(int j=1;j<=id&&1ll*pri[i]*pri[i]<=d[j];j++)
 40             k=idx(d[j]/pri[i]),g[j]-=g[k]-(i-1-D[k]),D[j]=i;
 41 }
 42 void work2()
 43 {
 44     for(int i=tot;i>=1;i--)
 45         for(int j=1;j<=id&&1ll*pri[i]*pri[i]<=d[j];j++)
 46         {
 47             if(pri[i+1]>d[j])f2[j]=1;
 48             else if(1ll*pri[i+1]*pri[i+1]>d[j])
 49                 f2[j]=(num[min(1ll*mx,d[j])]-num[pri[i+1]-1])*4+1;
 50             for(LL pi=pri[i],c=1;pi<=d[j];pi*=pri[i],c++)
 51             {
 52                 LL k=d[j]/pi,k2;
 53                 if(pri[i+1]>k)k2=1;
 54                 else if(1ll*pri[i+1]*pri[i+1]>k)
 55                     k2=(num[min(1ll*mx,k)]-num[pri[i+1]-1])*4+1;
 56                 else k2=f2[idx(k)];
 57                 f2[j]+=k2*(3*c+1);
 58             }
 59         }
 60 }
 61 LL solve(LL n)
 62 {
 63     if(n<=mx)return sum[n];
 64     ans1=0;work1();work2();
 65     for(int i=1,pos;i<=mx;i=pos+1)
 66     {
 67          int j=idx(n/i);LL k;
 68          if(pri[tot+1]>n/i)k=0;
 69          else k=g[j]-(tot-D[j]+1);
 70          ans1+=(sum[pos=min(1ll*mx,n/(n/i))]-sum[i-1])*k*4;
 71     }
 72     return ans1+f2[1];
 73 }
 74 int main()
 75 {
 76     f[1]=sum[1]=1;
 77     for(int i=2;i<=mx;i++)
 78     {
 79         num[i]=num[i-1]+(!f[i]);
 80         if(!f[i])pri[++tot]=i,f[i]=4,mn[i]=i,ci[i]=1;
 81         for(int j=1,now;(now=i*pri[j])<=mx;j++)
 82         {
 83             if(i%pri[j]==0)
 84             {
 85                 ci[now]=ci[i]+1;
 86                 f[now]=f[i/mn[i]]*(ci[now]*3+1);
 87                 mn[now]=mn[i]*pri[j];
 88                 break;
 89             }
 90             f[now]=f[i]*4;mn[now]=pri[j];ci[now]=1;
 91         }
 92         sum[i]=(sum[i-1]+f[i]);
 93     }
 94     pri[tot+1]=316241;
 95     T=read();
 96     while(T--)
 97     {
 98         n=read();
 99         printf("%lld
",solve(n));
100     }
101     return 0;
102 }

View Code

2.【51nod 1184】第N个质数

题意：给出一个数N，求第N个质数。

分析：WerKeyTom_FTDの博客

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long 
 6 using namespace std;
 7 const int N=12500000;
 8 int n,tot,f[N+5],pri[N+5];
 9 LL l,r,mid,ans,tmp;
10 bool np[N+5];
11 LL g(LL n,int m)
12 {
13     if(!m)return n;
14     if(m==1)return n-n/2;
15     if(n<=N)
16     {
17         if(m>=f[n])return 1;
18         if(m>=f[(int)sqrt(n)])return f[n]-m+1;
19     }
20     return g(n,m-1)-g(n/pri[m],m-1);
21 }
22 LL calc(LL n)
23 {
24     if(n<=N)return f[n];
25     tmp=ceil(sqrt(n));
26     return f[tmp]+g(n,f[tmp])-1;
27 }
28 int main()
29 {
30     for(int i=2;i<=N;i++)
31     {
32         if(!np[i])pri[++tot]=i;
33         for(int j=1;i*pri[j]<=N;j++)
34         {
35             np[i*pri[j]]=true;
36             if(i%pri[j]==0)break;
37         }
38         f[i]=f[i-1]+(!np[i]);
39     }
40     scanf("%d",&n);
41     l=1;r=22801763489;
42     if(n<=500000000)r=11037271757;
43     else l=11037271758;
44     if(n<=750000000)r=16875026921;
45     else l=16875026922;
46     if(n<=850000000)r=19236701629;
47     else l=19236701630;
48     if(n<=900000000)r=20422213579;
49     else l=20422213580;
50     if(n<=940000000)r=22801763489;
51     else l=22801763490;
52     if(n<=970000000)r=22086742277;
53     else l=22086742278;
54     if(n<=990000000)r=22563323957;
55     else l=22563323958;
56     while(l<=r)
57     {
58         mid=(l+r)>>1;
59         if(calc(mid)>=n)ans=mid,r=mid-1;
60         else l=mid+1;
61     }
62     printf("%lld",ans);
63     return 0;
64 }

View Code