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)

很简单的容斥，仔细想想就能明白。具体实现还是看代码吧。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<vector>
#define MOD 10007
using namespace std;
int n,m,k,fac[4001],kpri[4001],fsf[4001][4001],mp[10000001];
int m1,sum;
int dp[3001][4001];
void dfs(int cnt,int p_m,int assemble)
{
    if(cnt>kpri[0]) {sum+=m1/assemble*p_m;return;}
    dfs(cnt+1,p_m,assemble);
    dfs(cnt+1,-p_m,assemble*kpri[cnt]);
}
int main()
{
    scanf("%d%d%d",&n,&m,&k);
    int sqtk=sqrt(k);
    for(int i=1;i<=sqtk;i++)
        if(k%i==0)
        {
            fac[++fac[0]]=i;
            if(k/i>sqtk) fac[++fac[0]]=k/i;
        }
    sort(fac+1,fac+1+fac[0]);
    int tmp=k;
    for(int i=2;i<=sqtk;i++)
    {
        if(tmp==1) break;
        if(tmp%i==0)
        {
            kpri[++kpri[0]]=i;
            while(tmp%i==0) tmp/=i;
        }
    }
    if(tmp!=1) kpri[++kpri[0]]=tmp;
    sort(kpri+1,kpri+1+kpri[0]);
    for(int i=1;i<=fac[0];i++)
    {
        mp[fac[i]]=i;
        sum=0,m1=m/fac[i];
        dfs(1,1,1);
        dp[1][i]=sum%MOD;
    }
    for(int i=1;i<=fac[0];i++){
        cout<<dp[1][i]<<" "; 
    }
    cout<<endl; 
    for(int i=1;i<=fac[0];i++)
        for(int j=1;j<=i;j++)
            if(fac[i]%fac[j]==0) fsf[i][++fsf[i][0]]=j;
    for(int i=2;i<=n;i++)
        for(int j=1;j<=fac[0];j++)
        {
            if(fsf[j][0]==0) continue;
            for(int k=1;k<=fsf[j][0];k++)
                (dp[i][j]+=dp[i-1][fsf[j][k]]*dp[1][mp[fac[j]/fac[fsf[j][k]]]])%=MOD;
        }
    printf("%d
",dp[n][fac[0]]);
    return 0;
}