数学（二）

BSGS

BSGS用于求解形如$A^x\equiv B\pmod{P}$的高次剩余问题。
在这里我们先引入一些基础的数论知识。

欧拉函数及欧拉定理

欧拉函数$\phi(x)$表示x内与x互质的数的个数，比如$\phi(6)=2$，$\phi(8)=4$
欧拉函数的计算公式：
设$n=p_1^{k_1}p_2^{k_2}…p_m^{k_m}$，则有

$$phi(n)=n\times\dfrac{p_1-1}{p_1}\times\dfrac{p_2-1}{p_2}\times...\times\dfrac{p_m-1}{p_m}$$

欧拉定理：$若a,n互质，则a^{\phi(n)}\equiv 1\pmod{n}$

欧拉降幂

首先我们有广义欧拉定理

$$a^b=\begin{cases} a^{b\bmod \phi(p)} & \gcd(a,p)=1\\ a^b & \gcd(a,p)\not=1,b<\phi(p)\\ a^{b\bmod\phi(p)+\phi(p)} & \gcd(a,p)\not=1,b\ge\phi(p) \end{cases}\pmod p$$

因此对于指数非常大的情况我们可以通过欧拉降幂使得指数变小。
对于形容$w_1^{w_2^{w_3^{…}}}$这样的指数形式，实际上不断求$\phi$经过$O(\log p)$次后就会求到1，也就是说模意义下有效的指数层数不超过$O(\log p)$。
实际实现的时候，其实并不需要考虑gcd，因为指数加个$\phi(p)$并不影响答案。
下面给一个模板，长度为n的序列，每次查询$[l,r]$的$w_l^{w_{l+1}^{w_{l+2}^{ {…}^{w_r} } } }$在模p意义下的值

ll p, w[N], phi[N];
int n, tot;
ll ksm(ll a, ll b, ll p) {
    ll ans = 1;
    for (; b; b >>= 1, (a *= a) %= p)
        if (b & 1)(ans *= a) %= p;
    return ans;
}
void getphi(ll x, int step) {
    ll ans = 1;
    tot = step;
    for (ll i = 2; i * i <= x; ++i) {
        if (x % i == 0) {
            ans *= i - 1;
            x /= i;
            while (x % i == 0)
                ans *= i, x /= i;
        }
    }
    if (x != 1)ans *= x - 1;
    phi[step] = ans;
    if (ans == 1)return;
    getphi(ans, step + 1);
}
bool judge(ll a, ll b, ll p) {
    ll ans = 1;
    for (; b; b >>= 1, a *= a) {
        if (b & 1)ans *= a;
        if (a >= p || ans >= p)return true;
    }
    return false;
}

int main() {
    n = fast_IO::read();
    p = fast_IO::read();
    phi[0] = p;
    getphi(p, 1);
    for (int i = 1; i <= n; ++i)w[i] = fast_IO::read();
    int q = fast_IO::read();
    while (q--) {
        int l = fast_IO::read(), r = fast_IO::read();
        if (p == 1) {
            puts("0");
            continue;
        }
        ll ans = w[min(l + tot, r)] % phi[min(tot, r - l)] + (w[min(l + tot, r)] >= phi[min(tot, r - l)] ? phi[min(tot, r - l)] : 0);
        for (int i = min(l + tot, r) - 1, now = i - l; i >= l; --i, --now)
            ans = ksm(w[i] % phi[now], ans, phi[now]) + (judge(w[i], ans, phi[now]) ? phi[now] : 0);
        printf("%lld\n", ans % p);
    }
    return 0;
}

整数的阶

设a与n是互质的整数，其中$a\not=0,n>0$，使得$a^x\equiv 1\pmod{n}$成立的最小正整数x称为a模n的阶，符号为${ \rm { ord } }_na$
根据欧拉定理，$\phi(n)$是$a^x\equiv 1\pmod{n}$的一个解，但未必是最小，所以$\phi(n)$未必是阶。
下面给出阶相关的几个定理：

• $若\gcd(a,n)=1，则x_0是方程a^x\equiv 1\pmod{n}的一个解的充要条件是{ \rm{ ord } }_na|x_0$
• $若\gcd(a,n)=1，则{ \rm {ord} }_na|\phi(n)$
• $若\gcd(a,n)=1，则a^i\equiv a^j\pmod{n}的充要条件是i\equiv j\pmod{ { \rm {ord} }_na }$
• $若\;{ \rm {ord} }_na=t，且u为正整数，则有{ \rm {ord} }_n(a^u)=\dfrac{t}{\gcd(t, u)}$

原根

当$ { \rm {ord} }_na=\phi(n)$时，称a为n的原根。
咕咕咕

BSGS算法

懂了这些基础知识后，我们来看看BSGS算法是怎样求解的
首先使用BSGS的前提条件是P为素数，于是我们知道${ \rm {ord} }_PA\le \phi(P)$，所以我们就考虑$0\sim\phi(P)$的解。
我们知道，小于P的数可以表示成$x=a\lceil\sqrt{P}\rceil+b,0\le a,b < \lceil\sqrt{P}\rceil$
那么原式子就可以变形为$A^{a\lceil\sqrt{P}\rceil}\equiv B\times A^{-b}\pmod P$
那么我们只需要枚举所有的b将$B\times A^{-b}$插入哈希表，然后枚举a在哈希表中查找$A^{a\lceil\sqrt{P}\rceil}$就行了。
复杂度为$O(\sqrt{P})$
当然，如果我们将x表示成$x=a\lceil\sqrt{P}\rceil-b,1\le a \le \lceil\sqrt{P}\rceil,0\le b< \lceil\sqrt{P}\rceil$，就不用求A的逆元了
此时原式子变形为$A^{a\lceil\sqrt{P}\rceil}\equiv B\times A^b\pmod P$

ll ksm(ll a, ll b, ll p) {
	ll ans = 1;
	for (; b; (a *= a) %= p, b >>= 1)
		if (b & 1)(ans *= a) %= p;
	return ans;
}
ll BSGS(ll A, ll B, ll P) {
	unordered_map<ll, ll>mp;
	unordered_map<ll, bool>mmp;
	B %= P;
	ll s = ceil(sqrt(P));
	for (ll i = 0; i < s; ++i)
		mp[B * ksm(A, i, P) % P] = i,
		mmp[B * ksm(A, i, P) % P] = true;
	A = ksm(A, s, P);
	for (ll i = 1; i <= s; ++i) {
		ll now = ksm(A, i, P);
		if (mmp[now] && i * s - mp[now] >= 0)return i * s - mp[now];
	}
	return -1;
}

exBSGS

当P不为素数时，直接BSGS就有问题了，因为可能会出现$A^k\equiv 0\pmod{P}$
当$\gcd(A,P)\not|B$时，仅当$B=1$时有解（显然至少有一个解为0），其余均无解，而当$\gcd(A,P)=1$时直接进行BSGS，否则

$$A^x\equiv B\pmod{P}\\ A^{x-1}\times\dfrac{A}{\gcd(A,P)}\equiv \dfrac{B}{\gcd(A,P)}\pmod{\dfrac{P}{\gcd(A,P)}}\\ A^{x-1}\equiv \dfrac{B}{\gcd(A,P)}\times(\dfrac{A}{\gcd(A,P)})^{-1}\pmod{\dfrac{P}{\gcd(A,P)}}\\$$

而这个式子可看作

$$A'\equiv B'\pmod{P'}$$

反复进行直到$\gcd(A’,P’)=1$，再进行BSGS就可以了。
注意，最终我们BSGS求出的是$x-t$，我们在前面的过程中记录t的大小，最后加上BSGS的答案即可

ll exBSGS(ll A, ll B, ll P) {
	ll x, y, gcd, add = 0;
	if (B == 1)return 0;
	for (; (gcd = exGCD(A, P, x, y)) != 1; ++add) {
		if (B % gcd)return -1;
		B /= gcd, P /= gcd;
		exGCD(A / gcd, P, x, y);
		x = (x % P + P) % P;
		B = B * x % P;
	}
	ll ans = BSGS(A, B, P);
	if (~ans)return ans + add;
	return -1;
}

Lucas

卢卡斯定理用于求解当$p$为素数时$x\equiv C_n^m\pmod{p}$。
定理：

$$C_n^m\equiv C_{\lfloor n/p\rfloor}^{\lfloor m/p\rfloor}\times C_{n\bmod p}^{m\bmod p}\pmod{p}$$

证明：
首先我们需要证明$C_n^m=\dfrac{n}{m}\times C_{n-1}^{m-1}$对$1\le m\le n - 1$恒成立。

$$C_n^m=\dfrac{n!}{m!\times(n-m)!}=\dfrac{n}{m}\times\dfrac{(n-1)!}{(m-1)!\times(n-m)!}=\dfrac{n}{m}\times C_{n-1}^{m-1}$$

再结合$(1+x)^n$的二项式展开，有

$$(1+x)^p=1+C_p^1 x+C_p^2 x^2+...+C_p^{p-1} x^{p-1}+x^p=1+\dfrac{p}{1}C_{p-1}^0 x+\dfrac{p}{2}C_{p-1}^1 x^2+...+\dfrac{p}{p-1}C_{p-1}^{p-2} x^{p-1}+x^p\equiv 1+x^p\pmod{p}$$

对于$(1+x)^n$，有

$$(1+x)^n=((1+x)^p)^{\lfloor n/p\rfloor}\cdot(1+x)^{n\bmod p}\equiv (1+x^p)^{\lfloor n/p\rfloor}\cdot(1+x)^{n\bmod p}\pmod{p}$$

然后我们考虑上面得到的式子以及$(1+x)^n$中$x^m$的系数是什么，很容易得到

$$C_n^mx^m\equiv C_{\lfloor n/p\rfloor}^{\lfloor m/p\rfloor}x^{\lfloor m/p\rfloor \times p}\cdot C_{n\bmod p}^{m\bmod p}x^{m\bmod p}\pmod{p}\\ 即C_n^m\equiv C_{\lfloor n/p\rfloor}^{\lfloor m/p\rfloor}\times C_{n\bmod p}^{m\bmod p}\pmod{p}$$

所以我们只需要递归求解即可，注意边界条件

ll dfs(ll n, ll m, ll p) {
	if (n < m)return 0;
	if (n < p)return mul[n] * inv[m] % p * inv[n - m] % p;
	return dfs(n / p, m / p, p) * dfs(n % p, m % p, p) % p;
}
ll Lucas(ll n, ll m, ll p) {
	mul[0] = mul[1] = 1;
	inv[0] = inv[1] = 1;
	for (int i = 2; i <= n; ++i)
		mul[i] = mul[i - 1] * i % p;
	//线性递推求逆元
	for (int i = 2; i <= n; ++i)
		inv[i] = (p - p / i) * inv[p % i] % p;
	for (int i = 2; i <= n; ++i)
		(inv[i] *= inv[i - 1]) %= p;
	return dfs(n, m, p);
}

exLucas

当p不为素数的时候，我们就不能用Lucas求了，于是我们引入扩展卢卡斯定理。
首先声明一点，exLucas和Lucas没有本质联系，只是求解的问题类型相同。
我们可以将p进行质因数分解，得到$p=p_1^{k_1}\cdot p_2^{k_2}\cdot…\cdot p_t^{k_t}$
那么对于下面的模方程

$$\begin{cases} x\equiv C_n^m\pmod{p_1^{k_1}}\\ x\equiv C_n^m\pmod{p_2^{k_2}}\\ ...\\ x\equiv C_n^m\pmod{p_t^{k_t}} \end{cases}$$

我们可以用CRT进行求解，最终得到的x就是$C_n^m\pmod{p}$
因此扩展卢卡斯定理的重点在于如何求解

$$\begin{cases} C_n^m\pmod{p_1^{k_1}}\\ C_n^m\pmod{p_2^{k_2}}\\ ...\\ C_n^m\pmod{p_t^{k_t}} \end{cases}$$

假设我们当前要求$C_n^m\pmod{p^{k}}$即$\dfrac{n!}{m!(n-m)!}\pmod{p^k}$
那么这东西可以直接求解吗？显然是不能的，因为a对p的逆元存在的前提条件是gcd(a,p)=1
所以我们对这个式子变形，改成如下形式

$$C_n^m\equiv\dfrac{\dfrac{n!}{p^x}}{\dfrac{m!}{p^y}\times\dfrac{(n-m)!}{p^z}}\times p^{x-y-z}\pmod{p}$$

其中$\dfrac{x!}{p^y}$表示x!除掉其中所有的质因子p，且质因子p的个数为y。
不难证明，$x\ge y+z$，因此我们提出来的p的幂用快速幂求解即可。
现在的问题只剩下如何求解$\dfrac{x!}{p^y}\pmod{p^k}$了
而x!中p的倍数显然有$\lfloor\dfrac{x}{p}\rfloor$个，于是我们就得到

$$\dfrac{x!}{p^{\lfloor\frac{x}{p}\rfloor}}=(\lfloor\dfrac{x}{p}\rfloor)!\prod\limits_{i=1,i\not\equiv 0\pmod{p}}^{n}i$$

而最终我们要的结果是在模$p^k$意义下的，所以对于最右边求积的部分，我们可以将其改写成

$$\prod\limits_{i=1,i\not\equiv 0\pmod{p}}^{n}i=(\prod\limits_{i=1,i\not\equiv0\pmod{p}}^{p^k}i)^{\lfloor \frac{n}{p^k}\rfloor}(\prod\limits_{i=p^k\times\lfloor \frac{n}{p^k}\rfloor,i\not\equiv0\pmod{p}}^{n}i)$$

对后面的部分求解完后，看看前面的，$(\lfloor\dfrac{x}{p}\rfloor)!$中可能还有p这个质因子，因此我们递归下去继续求解。
到此，整个求解过程就结束了。

下面给出模板

ll inv(ll n, ll p) {
	//exgcd求逆元
	ll x, y;
	exGCD(n, p, x, y);
	return (x % p + p) % p;
}
ll fac(ll n, ll p, ll mul) {
	//求%mul意义下的阶乘
	if (!n)return 1;
	ll ans = 1;
	for (ll i = 2; i <= mul; ++i)
		if (i % p)(ans *= i) %= mul;
	//循环节
	ans = ksm(ans, n / mul, mul);
	//余项
	for (ll i = 2; i <= n % mul; ++i)
		if (i % p)(ans *= i) %= mul;
	return ans * fac(n / p, p, mul) % mul;//递归求解
}
ll C(ll n, ll m, ll p, ll mul) {
	//求%mul意义下的组合数C(n,m)
	ll up = fac(n, p, mul), d1 = fac(m, p, mul), d2 = fac(n - m, p, mul);
	ll k = 0;
	for (ll i = n; i; i /= p)k += i / p;//x
	for (ll i = m; i; i /= p)k -= i / p;//y
	for (ll i = n - m; i; i /= p)k -= i / p;//z
	return up * inv(d1, mul) % mul * inv(d2, mul) % mul * ksm(p, k, mul) % mul;
}
ll exLucas(ll n, ll m, ll mod) {
	ll lim = sqrt(mod), tmp = mod;
	*p = 0;
	for (int i = 2; i <= lim; ++i) {
		if (tmp % i == 0) {
			p[++*p] = i; mulp[*p] = 1;
			while (tmp % i == 0)mulp[*p] *= i, tmp /= i;
		}
	}
	if (tmp != 1)p[++*p] = tmp, mulp[*p] = tmp;
	for (int i = 1; i <= *p; ++i)
		b[i] = mulp[i], a[i] = (C(n, m, p[i], mulp[i]) % b[i] + b[i]) % b[i];
	return CRT(*p);//构成的模方程组为x=a(mod b)
}