JZOJ5787轨道
Description
2018年1月31日,152年一遇的超级大月全食在中国高空出现(没看到的朋友真是可惜),小B看到月食,便对月球的轨道产生了兴趣。他上网查重力加速度的公式,公式如下:
就在这个时候,他想到了一个跟这个差不多的问题,那就是对于以下公式:
已知n和k,求这n个正整数在都不大于m的情况下有多少种选择方式,使得v为与k互质正整数?
![](https://i-blog.csdnimg.cn/blog_migrate/cba393e4c3644343b60661d19b66b05e.png)
就在这个时候,他想到了一个跟这个差不多的问题,那就是对于以下公式:
![](https://i-blog.csdnimg.cn/blog_migrate/4bb66cac4018585a1bbbb99722607d47.png)
已知n和k,求这n个正整数在都不大于m的情况下有多少种选择方式,使得v为与k互质正整数?
Input
一行三个正整数n,m,k(意义见题目描述)。
Output
输出一个答案,代表方案数。因为答案可能会很大,所以输出方案数mod 10007的值。
数据范围
对于20%的数据 1<=n,m<=8 k<=100
对于40%的数据 1<=n<=50 1<=m<=10^6 1<=k<=10^4
对于70%的数据 1<=n<=100 1<=m<=10^9 1<=k<=10^7
对于100%的数据 1<=n<=3000 1<=m<=10^9 1<=k<=10^7
对于20%的数据 1<=n,m<=8 k<=100
对于40%的数据 1<=n<=50 1<=m<=10^6 1<=k<=10^4
对于70%的数据 1<=n<=100 1<=m<=10^9 1<=k<=10^7
对于100%的数据 1<=n<=3000 1<=m<=10^9 1<=k<=10^7
题解
模拟赛的题,然后我就GG了
设dp[i][j]为前i个数的乘积与k的gcd是k的第j个约数(且乘积除以公约数与k互质)的方案数。
然后转移方程是dp[i][j]=sigema dp[i-1][k]*dp[1][第j个约数/第k个约数是第几个约数]这里的k是一个枚举的变量。第j个约数可以整除第k个约数。
然后考虑初始化dp[1][...];
dp[1][j]代表的是和1~m中与k的gcd为k的第j个约数的数的数量。
但是枚举绝对会T。
我们其实要求的是gcd(x,k)=第j个约数(1<=x<=m)的方案数。
把公式化一下化为gcd(x.k)=1(1<=x<=m/第j个约数(向下取整)(以后称为m1))的方案数。
然后使用容斥。
怎么用容斥呢,举个例子。
设k的质因数为A,B,C
方案数为m1/1
- m1/A
- m1/B
- m1/C
+ m1/(A*B)
+ m1/(B*C)
+ m1/(A*C)
- m1/(A*B*C)
很简单的容斥,仔细想想就能明白。具体实现还是看代码吧。
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<cstdlib> 5 #include<cmath> 6 #include<algorithm> 7 #include<vector> 8 #define MOD 10007 9 using namespace std; 10 int n,m,k,fac[4001],kpri[4001],fsf[4001][4001],mp[10000001]; 11 int m1,sum; 12 int dp[3001][4001]; 13 void dfs(int cnt,int p_m,int assemble) 14 { 15 if(cnt>kpri[0]) {sum+=m1/assemble*p_m;return;} 16 dfs(cnt+1,p_m,assemble); 17 dfs(cnt+1,-p_m,assemble*kpri[cnt]); 18 } 19 int main() 20 { 21 scanf("%d%d%d",&n,&m,&k); 22 int sqtk=sqrt(k); 23 for(int i=1;i<=sqtk;i++) 24 if(k%i==0) 25 { 26 fac[++fac[0]]=i; 27 if(k/i>sqtk) fac[++fac[0]]=k/i; 28 } 29 sort(fac+1,fac+1+fac[0]); 30 int tmp=k; 31 for(int i=2;i<=sqtk;i++) 32 { 33 if(tmp==1) break; 34 if(tmp%i==0) 35 { 36 kpri[++kpri[0]]=i; 37 while(tmp%i==0) tmp/=i; 38 } 39 } 40 if(tmp!=1) kpri[++kpri[0]]=tmp; 41 sort(kpri+1,kpri+1+kpri[0]); 42 for(int i=1;i<=fac[0];i++) 43 { 44 mp[fac[i]]=i; 45 sum=0,m1=m/fac[i]; 46 dfs(1,1,1); 47 dp[1][i]=sum%MOD; 48 } 49 for(int i=1;i<=fac[0];i++){ 50 cout<<dp[1][i]<<" "; 51 } 52 cout<<endl; 53 for(int i=1;i<=fac[0];i++) 54 for(int j=1;j<=i;j++) 55 if(fac[i]%fac[j]==0) fsf[i][++fsf[i][0]]=j; 56 for(int i=2;i<=n;i++) 57 for(int j=1;j<=fac[0];j++) 58 { 59 if(fsf[j][0]==0) continue; 60 for(int k=1;k<=fsf[j][0];k++) 61 (dp[i][j]+=dp[i-1][fsf[j][k]]*dp[1][mp[fac[j]/fac[fsf[j][k]]]])%=MOD; 62 } 63 printf("%d\n",dp[n][fac[0]]); 64 return 0; 65 }