SPOJ:[DIVCNT3]Counting Divisors

 

 

 

题目大意:求1~N的每个数因子数的立方和。

题解:由于N过大,我们不能直接通过线性筛求解。我们可以采用洲阁筛。

洲阁筛的式子可以写成:

 

对于F(1~√n),可以直接线性筛求解。

 

对于,我们进行以下DP:

g[i][j]为1~j中,与前i个质数互质的数的F值之和。

dp过程中,有

如果p[i]>j,则g[i][j]=F(1);

如果p[i]*p[i]>j>=p[i],则g[i][j]=g[i-1][j]-p[i]^k*F(1)=g[d][j]-(p[d+1]^k~p[i]^k)*F(1),其中p[d]为最大的p[d]*p[d]<=j的质数;

如果j>=p[i]*p[i],我们老老实实调用式子。

 

对于,我们用类似的方式DP:

f[i][j]为1~j中,只以第i个到第m个质数为质因子的数的F值之和(m为√N以内质数个数)。

类似的,dp过程中,有

如果p[i]>j,则f[i][j]=F(1);

如果p[i]*p[i]>j>=p[i],则f[i][j]=f[i+1][j]+F(p[i])*F(1)=F(1)+F(p[i]~p[e])*F(1),其中p[e]为最大的p[e]<=j的质数;

如果j>=p[i]*p[i],我们老老实实调用式子。

 

代码:

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 long long block,n,lb[2000005],dp[2000005],zs[2000005],ans[2000005],ans2;
 4 bool bo[2000005],bo2[2000005];
 5 int dp2[2000005],dp3[2000005],m,mm,tot,j,bo3[2000005],tt;
 6 void xxs(int n)
 7 {
 8     ans[1]=1;
 9     for(int i=2;i<=n;i++)
10     {
11         if(bo[i]==0){ mm++; zs[mm]=i; ans[i]=4; bo2[i]=1; bo3[i]=4; }
12         for(int j=1;j<=mm;j++)
13         if(zs[j]*i>n)break;else
14         if(i%zs[j]==0)
15         {
16             if(bo2[i]==1)bo2[i*zs[j]]=1; bo3[i*zs[j]]=bo3[i]+3; ans[i*zs[j]]=ans[i]/bo3[i]*(bo3[i]+3);
17             bo[i*zs[j]]=1; break;
18         }else { ans[i*zs[j]]=ans[i]*4; bo[i*zs[j]]=1; bo3[i*zs[j]]=4; }
19     }
20 }
21 int dy(long long x){ if(x<=block)return x;else return tot-n/x+1; }
22 long long get(int i,int j)
23 {
24     if(dp2[j]==i)return dp[j];else 
25     if(lb[j]<zs[i]){ if(j>0)return 1; return 0; } 
26     return dp[j]-i+dp2[j];
27 }
28 long long get2(int i,int j)
29 {
30     if(dp2[j]==i)return dp[j];else 
31     if(lb[j]<zs[i])return 1; 
32     return (dp3[j]-i+1)*4+1;
33 }
34 int main()
35 {
36     scanf("%d",&tt); xxs(500000);
37     for(int ii=1;ii<=tt;ii++)
38     {
39         scanf("%lld",&n); block=(int)sqrt(n); ans2=0;
40         for(m=1;m<=mm;m++)if(zs[m]>block)break; m--; tot=0;
41         for(int i=1;i<=block;i++){ lb[++tot]=i; if(1ll*i*i<n)lb[++tot]=n/i; }
42         sort(lb+1,lb+tot+1);
43         for(int i=1;i<=tot;i++)dp[i]=lb[i],dp2[i]=0;
44         for(int i=1;i<=m;i++)
45         {
46             for(int j=tot;j>=1;j--)
47             {
48                 if(lb[j]<zs[i]*zs[i])break; dp2[j]=i;
49                 dp[j]=dp[j]-get(i-1,dy(lb[j]/zs[i]));
50             }
51         }
52         for(int i=1;i<=block;i++)
53         ans2=ans2+ans[i]*(get(m,tot-i+1)-1)*4;
54         j=1; 
55         for(int i=1;i<=tot;i++)
56         {
57             dp[i]=1; dp2[i]=m+1;
58             while((j<=m)and(zs[j]<=lb[i]))j++; dp3[i]=j-1;
59         }
60         for(int i=m;i>=1;i--)
61         {
62             for(int j=tot;j>=1;j--)
63             {
64                 if(lb[j]<zs[i]*zs[i])break; dp[j]=get2(i+1,j); dp2[j]=i;
65                 long long t=j,l=1; 
66                 while(lb[t]>=zs[i])
67                 {
68                     t=dy(lb[t]/zs[i]); l+=3;
69                     dp[j]=dp[j]+l*get2(i+1,t);
70                 }
71             }
72         }
73         ans2=ans2+get2(1,tot);
74         printf("%lld\n",ans2);
75     }
76 } 
View Code
posted @ 2018-01-16 19:23  GhoStreach  阅读(319)  评论(0编辑  收藏  举报