提高B组模拟轨道

wjh15101051's space

卧薪尝胆，厚积薄发。

提高B组模拟轨道

Date: Sun Aug 12 16:05:32 CST 2018

In Category: NoCategory

Description：

$v=\frac{a_1\times a_2\times \cdots\times a_n}{k}$ 已知 $n$ 和 $k$ ，求这 $n$ 个正整数都不大于 $m$ 的情况下有多少种选择方式使得 $v$ 为与 $k$ 互质正整数。

$1\le n\le 3000$ $1\le m\le 10^9$ $1\le k\le 10^7 $

Solution：

$DP$ ，设状态为 $f[i][j]$ 表示选了 $i$ 个数乘积与 $k$ 的最大公约数为 $k$ 的第 $j$ 个约数且乘积除以公约数与 $k$ 互质的方案数，设 $a[i]$ 为 $k$ 的第 $i$ 个约数， $b[i]$ 为 $i$ 是 $k$ 的第几个约数，转移方程为 $\begin{align}f[i][j]=\sum_{l=1}^j f[i-1][l]\times f[1][b[a[j]/a[l]]]\end{align}$ ，预处理 $a[j]$ 的约数转移。含义是选了 $i-1$ 个数与 $k$ 的最大公约数是 $k$ 的第 $l$ 个约数且乘积除以 $a[l]$ 与 $k$ 互质，则后面再选一个数使得与 $k$ 的最大公约数是 $k$ 的第 $j$ 个约数且乘积除以公约数与 $k$ 互质，如果不考虑乘积除以公约数与 $k$ 互质，则所有使得含有一个约数 $a[p]$ 且 $lcm(a[p],a[l])=a[j]$ 的数字合法，但是由于要求互质，所以乘积除以 $k$ 之后不能再有 $k$ 的约数，所以只能是 $a[j]/a[l]$ ，就能保证最后的最大公约数是 $a[j]$ 且除掉后与 $k$ 互质。现在的问题就变成了如何求 $f[1][l]$ 即与 $k$ 的最大公约数是 $k$ 的第 $l$ 个约数且除掉 $a[l]$ 后与 $k$ 互质的数的个数。即 $gcd(w,k)=a[l]$ ，且 $gcd(\lfloor\frac w{a[l]}\rfloor,k)=1$ ，发现满足第二个条件一定满足第一个条件，问题就变成了求 $1$ ~ $\lfloor\frac m{a[l]}\rfloor$ 中与 $k$ 互质的数的个数，容斥原理做就可以了。

Code：


#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n,m,k;

#define MAXN 3010

#define MOD 10007

int a[MAXN],tot = 0;

int prime[MAXN],cnt = 0;

int f[MAXN][MAXN];

int res;

int b[10000010];

void dfs(int cur,int curlcm,int border,int type)

{

	if(curlcm > border)return;

	if(cur == 0)

	{

		res += type * (border / curlcm);

		return;

	}

	if(curlcm * prime[cur] <= border)dfs(cur - 1,curlcm * prime[cur],border,-type);

	dfs(cur - 1,curlcm,border,type);

	return;

}

bool cmp(int x,int y){return x > y;}

inline void add(int &x,int y){x = x + y;x =  (x >= MOD ? x - MOD : x);return;}

struct num

{

	int v,nxt;

}p[1000010];

int pnum = 0;

int lin[MAXN] = {0};

inline void ins(int a,int b)

{

	++pnum;p[pnum].v = b;p[pnum].nxt = lin[a];lin[a] = pnum;

	return;

}

inline int read()

{

	register int res = 0;

	register char c = getchar();

	while(!isdigit(c))c = getchar();

	while(isdigit(c))

	{

		res = (res << 1) + (res << 3) + (c ^ 48);

		c = getchar();

	}

	return res;

}

int main()

{

	freopen("orbits.in","r",stdin);

	freopen("orbits.out","w",stdout);

	n = read();m = read();k = read();

	for(register int i = 1;i * i <= k;++i)if(k % i == 0)

	{

		a[++tot] = i;

		if(i * i != k)a[++tot] = k / i;

	}

	register int s = k;

	for(register int i = 2;i * i <= k;++i)

	{

		if(s % i == 0)prime[++cnt] = i;

		while(s % i == 0)s /= i;

	}

	if(s != 1)prime[++cnt] = s;

	sort(a + 1,a + 1 + tot);

	for(register int i = 1;i <= tot;++i)b[a[i]] = i;

	sort(prime + 1,prime + 1 + cnt,cmp);

	for(register int i = 1;i <= tot;++i)

	{

		res = 0;

		dfs(cnt,1,m / a[i],1);

		f[1][i] = res % MOD;

	}

	for(register int j = 1;j <= tot;++j)

		for(register int l = 1;a[l] * a[l] <= a[j];++l)

			if(a[j] % a[l] == 0)

			{

				ins(j,l);

				register int pp = a[j] / a[l];

				if(pp != a[l])ins(j,b[pp]);

			}

	for(register int i = 2;i <= n;++i)

		for(register int j = 1;j <= tot;++j)

			for(register int l = lin[j];l != 0;l = p[l].nxt)

				add(f[i][j],f[i - 1][p[l].v] * f[1][b[a[j] / a[p[l].v]]] % MOD);

	printf("%d\n",f[n][tot]);

	fclose(stdin);

	fclose(stdout);

	return 0;

}

In tag: 数学-容斥原理 DP-DP

wjh15101051