JZOJ5787轨道（容斥+DP）

JZOJ5787轨道

Description

2018年1月31日，152年一遇的超级大月全食在中国高空出现（没看到的朋友真是可惜），小B看到月食，便对月球的轨道产生了兴趣。他上网查重力加速度的公式，公式如下：


就在这个时候，他想到了一个跟这个差不多的问题，那就是对于以下公式：

已知n和k，求这n个正整数在都不大于m的情况下有多少种选择方式，使得v为与k互质正整数？
 

Input

一行三个正整数n,m,k（意义见题目描述）。

Output

输出一个答案，代表方案数。因为答案可能会很大，所以输出方案数mod 10007的值。
Data Constraint
数据范围
对于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 }

 

posted @ 2018-09-05 16:20  Xu-daxia  阅读(390)  评论(0)  编辑  收藏  举报