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

嗯这道题拖了好几天了...

设dp[i][j]表示枚举到第i个数，这i个数乘积为j的方案数。

发现j的取值范围太大，数组根本开不下。考虑缩小j的范围。

发现题目中有关i个数乘积的部分中，有用的只有乘积与k的gcd。

因此修改dp[i][j]为枚举到第i个数，这i个数的乘积与k的gcd是k的第几个约数。

考虑转移方程，发现方程即为dp[i][j]=dp[i-1][k]*dp[1][a[j]/a[k]是k的第几个约数]。

因此只需预处理出dp[1][sth]，也就是在1~m中有多少个数与k的gcd为a[sth]。

枚举的话肯定行不通，要用更快捷的方法。

考虑容斥原理。

由于博主太蒟蒻，无法叙述清楚容斥系数，烦请各位大佬自行查看代码...

#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()
{
    freopen("orbits.in","r",stdin);
    freopen("orbits.out","w",stdout);
    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++)
        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;
}