P2155 [SDOI2008] 沙拉公主的困惑

\(\text{Description}\)

• 給定 \(t,r\)，有 \(t\) 次詢問，每次給定 \(n,m\)，請求出 \([1,n!]\) 中與 \(m!\) 互質的數的數量，答案對 \(r\) 取模.

• 對於 \(100\%\) 的資料，\(1\le m\le n\le 10^7,1\le t\le 10^4,2\le r\le 10^9+10\) 且 \(r\) 為質數.

\(\text{Solution}\)

首先因為 \(m\le n\)，所以 \(m!\mid n!\)，則與 \(m!\) 互質的有 \(\dfrac{n!}{m!}\cdot\varphi(m!)\)

個.

我們來整理一下：

\[\begin{aligned} \dfrac{n!}{m!}\cdot\varphi(m!) &=\dfrac{n!}{m!}\cdot(m!\cdot\prod\limits_{質數p\mid m!}\dfrac{p-1}{p})\\ &=n!\cdot\prod\limits_{質數p\mid m!}\dfrac{p-1}{p} \end{aligned} \]

因為 \(m!=1\times 2\times \cdots\times m\)，所以質數 \(p\) 要整除 \(m!\) 實際上就是要滿足質數 \(p\le m\).

所以就是

\[n!\cdot\prod\limits_{質數p\le m}\dfrac{p-1}{p}\\ =\dfrac{n!\cdot\prod\limits_{質數p\le m}(p-1)}{\prod\limits_{質數p\le m}p} \]

看起來預處理階乘，字首積和逆元，然後詢問時三個一乘就完事了.

求連續 \(n\) 個數的逆元見 【數學】逆元.

具體地，用 \(jc_x\) 表示 \(x!\)，\(phi_x\) 表示 \(\prod\limits_{前x個質數p}(p-1)\)（好吧我知道這不是 \(phi\) 但我們暫且叫它 \(phi\)），\(inv_x\) 表示 \(inv\left(\prod\limits_{前x個質數p}p\right)\).

查詢時用 \(\operatorname{upper\_bound}\) 找到質數中第一個大於 \(m\) 的，則在這之前的所有質數必然都是小於等於 \(m\) 的，記它的位置是 \(y\)，\(x=y-1\)，然後答案就是 \((jc_n\cdot phi_x\cdot inv_x)\bmod r\)

.

scanf("%d%d", &n, &m);
int x = upper_bound(p + 1, p + p[0] + 1, m) - (p + 1);
printf("%d\n", 1ll * jc[n] * phi[x] % r * inv[x] % r);

然而兔隊提供了 \(\text{hack}\) 資料，\(\text{hack}\) 的原理是，當分子與分母中含有相同個因數 \(r\) 時，正常計算會約分掉，但程式中如果把三個乘起來的話這三個中有 \(0\)，算出來就是 \(0\)，這是錯的.

所以重點在於考慮如何處理這種情況.

因為 \(r\) 為質數，所以 \(\prod\limits_{質數p\le m}(p-1)\) 不可能含有因數 \(r\)，所以當 \(n!\) 中因數 \(r\) 的個數比 \(\prod\limits_{質數p\le m}p\) 多的話，約分後還是 \(r\) 的倍數，輸出 \(0\).

否則說明它們的因數 \(r\) 的個數一定相等！那麼所有 \(r\) 都會被約掉.

在求階乘和逆元時遇到 \(r\) 就不乘.

這樣算還得特判 \(m=1\) 的情況，答案為 \(jc_n\).

\(\text{Code}\)

#include <iostream>
#include <cstdio>
#include <algorithm>
using namespace std;

const int MAXN = 1e7 + 5;

int p[MAXN];
bool vis[MAXN];

void get_prime()
{
	for (int i = 2; i < MAXN; i++)
	{
		if (!vis[i])
		{
			p[++p[0]] = i;
		}
		for (int j = 1; j <= p[0] && i * p[j] < MAXN; j++)
		{
			vis[i * p[j]] = true;
			if (i % p[j] == 0)
			{
				break;
			}
		}
	}
}

int t, r, n, m;

int ksm(int a, int b)
{
	int base = a, ans = 1;
	while (b)
	{
		if (b & 1)
		{
			ans = 1ll * ans * base % r;
		}
		base = 1ll * base * base % r;
		b >>= 1;
	}
	return ans;
}

int jc[MAXN];

void get_jc()
{
	jc[0] = 1;
	for (int i = 1; i < MAXN; i++)
	{
		if (i != r)
		{
			jc[i] = 1ll * jc[i - 1] * i % r;
		}
		else
		{
			jc[i] = jc[i - 1];
		}
	}
}

int phi[MAXN];

void get_phi()
{
	phi[0] = 1;
	for (int i = 1; i <= p[0]; i++)
	{
		phi[i] = 1ll * phi[i - 1] * (p[i] - 1) % r;
	}
}

int pro[MAXN], inv[MAXN];

void get_inv()
{
	pro[0] = 1;
	for (int i = 1; i <= p[0]; i++)
	{
		if (p[i] != r)
		{
			pro[i] = 1ll * pro[i - 1] * p[i] % r;
		}
		else
		{
			pro[i] = pro[i - 1];
		}
	}
	inv[p[0]] = ksm(pro[p[0]], r - 2);
	for (int i = p[0] - 1; i >= 1; i--)
	{
		if (p[i + 1] != r)
		{
			inv[i] = 1ll * inv[i + 1] * p[i + 1] % r;
		}
		else
		{
			inv[i] = inv[i + 1];
		}
	}
}

int main()
{
	scanf("%d%d", &t, &r);
	get_prime();
	get_jc();
	get_phi();
	get_inv();
	while (t--)
	{
		scanf("%d%d", &n, &m);
		if (n / r > m / r)
		{
			puts("0");
			continue;
		}
		else if (m == 1)
		{
			printf("%d\n", jc[n]);
			continue;
		}
		int x = upper_bound(p + 1, p + p[0] + 1, m) - (p + 1);
		printf("%d\n", 1ll * jc[n] * phi[x] % r * inv[x] % r);
	}
	return 0;
}