整除和同余相关理论

下面的所有知识均属数论范畴。

---

$Lucas$ 定理

当 $p$ 为质数时有：

$$C_n^m\mathrm{m}\mathrm{o}\mathrm{d}p=C_{n/p}^{m/p}\ast C_{n\mathrm{m}\mathrm{o}\mathrm{d}p}^{m\mathrm{m}\mathrm{o}\mathrm{d}p}\mathrm{m}\mathrm{o}\mathrm{d}p$$

代码就不贴了 （实际上是懒）。

---

如果 $p$ 不为质数呢？

那么我们先用唯一分解定理将 $p$ 分解：
$$p=\sum _{k_i\in Prime}{k}_i^{{\alpha }_i}$$
然后把原式拆分为若干个同余方程，最后用中国剩余定理合并即可。

现在我们讨论如何求 $C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{k}_i^{{\alpha }_i}$，而 $C_n^m=\frac{n!}{m!(n-m)!}$，所以问题本质上是求阶乘的取模。

设我们要求的是 $P!\mathrm{m}\mathrm{o}\mathrm{d}{k}_i^{{\alpha }_i}$，那么我们可以把 $P!$ 包含的所有 $k_i$ 质因子取出，由于剩下的数与模数互质，所以我们就可以直接求逆元了。

所以最终要求的问题是求：

$$\frac{\frac{n!}{k_i^{{\beta }_1}}}{\frac{m!}{k_i^{{\beta }_2}}\cdot \frac{(n-m)!}{k_i^{{\beta }_3}}}\cdot k_i^{{\beta }_1-{\beta }_2-{\beta }_3}\mathrm{m}\mathrm{o}\mathrm{d}{k}_i^{{\alpha }_i}$$

不过我们还是先从研究 $P!\mathrm{m}\mathrm{o}\mathrm{d}{k}_i^{{\alpha }_i}$ 开始。

同 $OI-wiki$，我们用 $22!\mathrm{m}\mathrm{o}\mathrm{d}{3}^2$ 举例。

我们知道：
$$22!=1\times 2\times 3\times 4\times 5\times 6\times 7\times 8\times 9\times 10\times 11\times 12\times 13\times 14\times 15\times 16\times 17\times 18\times 19\times 20\times 21\times 22.$$

我们把 $3$ 的倍数提出来，化简一下有：
$$22!=3^7\times (1\times 2\times 3\times 4\times 5\times 6\times 7)\times (1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19\times 20\times 22)$$

第一个部分是 $k_i$ 的 $\lfloor \frac{P}{k_i}\rfloor$ 次幂，这个部分只需要记录指数，不需要求具体的值；

第二个部分是 $\lfloor \frac{P}{k_i}\rfloor !$，可以递归求解；

第三个部分是 $1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17$，仔细观察可以发现：
$$1\times 2\times 4\times 5\times 7\times 8\equiv 10\times 11\times 13\times 14\times 16\times 17(\mathrm{m}\mathrm{o}\mathrm{d}{k}_i^{{\alpha }_i})$$
形式化地，有：（$t\in Z_{+}$）
$$\prod _{i=1,\gcd (i,k_i)=1}^{k_i^{{\alpha }_i}}i\equiv \prod _{i=1,\gcd (i,k_i)=1}^{k_i^{{\alpha }_i}}(i+t{k}_i^{{\alpha }_i})(\mathrm{m}\mathrm{o}\mathrm{d}{k}_i^{{\alpha }_i})$$
这个循环节一共循环了 $\lfloor \frac{P}{k_i^{{\alpha }_i}}\rfloor$ 次，我们直接求出来之后快速幂即可。

最后是 $19\times 20\times 22$，循环节的余项，表示出来即为：
$$\prod _{i=1,\gcd (i,k_i)=1}^{P\mathrm{m}\mathrm{o}\mathrm{d}{k}_i^{{\alpha }_i}}i$$
分步求解即可。

如果多测并且模数相同，可以对模数的每一个质因子 $k_i$，求出所有 $P\le k_i^{{\alpha }_i}$ 的 ${\prod }_{i=1,\gcd (i,k_i)=1}^{P}i$。

贴上代码：

#include<iostream>
#include<cstdio>
#include<cstdlib>

using LL = long long;
using namespace std;

LL n, m, p;

inline LL Get(LL k, LL p) { if(k >= p) k -= p; return k; }

inline void Exgcd(LL a, LL b, LL &x, LL &y)
{
	if(b == 0LL) return x = 1LL, y = 0LL, void();
	Exgcd(b, a % b, y, x);
	return y -= a / b * x, void();
}

inline LL Pow(LL a, LL b, LL p)
{
	LL res = 1LL;
	while(b)
	{
		if(b & 1) res = res * a % p;
		a = a * a % p;
		b >>= 1;
	}
	return res;
}

inline LL Fac(LL k, LL pi, LL pk)
{
	if(!k) return 1LL;
	LL res = 1LL;
	for(register LL i = 2LL; i <= pk; i++) if(i % pi) res = res * i % pk;
	res = Pow(res, k / pk, pk);
	for(register LL i = 2LL; i <= k % pk; i++) if(i % pi) res = res * i % pk;
	return res * Fac(k / pi, pi, pk) % pk; 
}

inline LL Inv(LL k, LL p)
{
	LL x, y; Exgcd(k, p, x, y);
	return Get(x % p + p, p);
}

inline LL CRT(LL k, LL t) { return k * (p / t) % p * Inv(p / t, t) % p; }

inline LL C(LL n, LL m, LL pi, LL pk)
{
	LL v1 = Fac(n, pi, pk), v2 = Fac(m, pi, pk), v3 = Fac(n - m, pi, pk);
	LL k = 0LL;
	for(register LL i = n; i; i /= pi) k += i / pi;
	for(register LL i = m; i; i /= pi) k -= i / pi;
	for(register LL i = n - m; i; i /= pi) k -= i / pi;
	return v1 * Inv(v2, pk) % pk * Inv(v3, pk) % pk * Pow(pi, k, pk) % pk;
}

inline LL ExLucas(LL n, LL m, LL p)
{
	LL ans = 0LL, t = p, pt = 1;
	for(register int i = 2; i * i <= p; i++) if(t % i == 0)
	{
		while(t % i == 0) pt *= i, t /= i;
		ans = Get(ans + CRT(C(n, m, i, pt), pt), p);
		pt = 1;
	}
	if(t > 1) ans = Get(ans + CRT(C(n, m, t, t), t), p);
	return ans;
}

int main()
{
	#ifdef FILE
		freopen("Input.in", "r", stdin);	
	#endif
	scanf("%lld%lld%lld", &n, &m, &p);
	printf("%lld\n", ExLucas(n, m, p));
	return 0;
}

---

$BSGS$ 算法

当 $a\perp p$ 时有如下算法：

设 $t=\lfloor \sqrt{p}\rfloor$，则方程可以化为：
$$a^{it-k}\equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)$$

进一步有：（如果 $a\perp b$ 不成立，这一步得出的结论不能推至上一步）

$$(a^t{)}^i\equiv b{a}^k(\mathrm{m}\mathrm{o}\mathrm{d}p)$$

因为 $k\in [0,t)$，所以可以把同余式右侧的式子用哈希表存储（如有重复选择最大的 $k$ 即可），然后从小到大枚举同余式左侧的 $i$，如果对于某个 $i$，等式左边的值可以在哈希表中查到，那么此时的 $it-k$ 就是最小的 $x$。

代码不贴 ，因为懒。

---

如果 $a\perp p$ 不成立怎么办？

根据模的分配律，我们设 $g=\gcd (a,p)$，那么原方程可以变成：
$$a^{x-1}\cdot \frac{a}{g}\equiv \frac{b}{g}(\mathrm{m}\mathrm{o}\mathrm{d}\frac{p}{g})$$
所以当 $b\ne 1$ 且 $g\nmid b$ 时，该方程无解。

如果保证一直有解，那么我们将这个方程一直化开，直到 $a\perp \frac{p}{{\sum }_{i=1}^k{g}_i}$。

此时有：（设 $q={\sum }_{i=1}^k\frac{a}{g_i}$）
$$a^{x-k}\equiv \frac{b}{q{\sum }_{i=1}^k{g}_i}(\mathrm{m}\mathrm{o}\mathrm{d}\frac{p}{{\sum }_{i=1}^k{g}_i})$$
这个时候就可以用 $BSGS$ 求解了，最后的答案记得要加上 $k$。

当然，算法过程中还需要特判 $q=\frac{b}{{\sum }_{i=1}^k{g}_i}$ 的情况，此时左侧 $a$ 的次数为 $0$，答案即为 $k$。

贴上代码：

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<cctype>
#include<map>

using namespace std;

int a, p, b;
map<int, int> Hash;

inline int Read()
{
	int v = 0; char c = getchar();
	while(!isdigit(c)) c = getchar();
	while(isdigit(c)) v = (v << 3) + (v << 1) + c - '0', c = getchar();
	return v;
}

inline int gcd(int a, int b) { return b ? gcd(b, a % b) : a; }

inline void Exgcd(int a, int b, int &x, int &y)
{
	if(b == 0) return x = 1, y = 0, void();
	Exgcd(b, a % b, y, x);
	return y -= a / b * x, void();
}

inline int Inv(int a, int b)
{
	int x, y; Exgcd(a, b, x, y);
	return ((x % b) + b) % b;
}

inline int Pow(int a, int b, int p)
{
	int res = 1;
	while(b)
	{
		if(b & 1) res = 1LL * res * a % p;
		a = 1LL * a * a % p;
		b >>= 1;
	}
	return res;
}

inline int BSGS(int a, int b, int p)
{
	Hash.clear(), b %= p;
	int t = (int)(sqrt(p) + 1), w = 1;
	for(register int i = 0; i < t; i++)
	{
		int val = 1LL * b * w % p;
		Hash[val] = i;
		w = 1LL * w * a % p;
	}
	int val = 1; a = Pow(a, t, p), val = 1;
	if(!a) return b ? -1 : 1;
	for(register int i = 1; i <= t; i++)
	{
		val = 1LL * val * a % p;
		int j = Hash.find(val) == Hash.end() ? -1 : Hash[val]; 
		if((~j) && i * t - j >= 0) return i * t - j;
	}
	return -1;
}

inline int ExBSGS(int a, int b, int p)
{
	if(b == 1 || p == 1) return 0;
	int g = gcd(a, p), k = 0, q = 1;
	while(g ^ 1)
	{
		if(b % g) return -1;
		++k; b /= g, p /= g;
		q = 1LL * q * (a / g) % p;
		if(q == b) return k;
		g = gcd(a, p);
	}
	int res = BSGS(a, 1LL * b * Inv(q, p) % p, p);
	return (~res) ? res + k : res;
}

int main()
{
	#ifdef FILE
		freopen("Input.in", "r", stdin);
	#endif
	a = Read(), p = Read(), b = Read();
	while(a || p || b)
	{
		a %= p, b %= p;
		int ans = ExBSGS(a, b, p);
		if(ans == -1) puts("No Solution");
		else printf("%d\n", ans);
		a = Read(), p = Read(), b = Read();
	}
	return 0;
}

---

中国剩余定理

如果保证所有的 $p_i$ 之间两两互质，有如下算法：

设 $M={\prod }_{i=1}^m{p}_i$，$M_i=\frac{M}{p_i}$，$M_i^{-1}$ 为 $M_i$ 在模 $p_i$ 意义下的逆元。

设每个方程的特解为 $x_i$，则方程组特解为 ${\sum }_{i=1}^m{x}_i{M}_i{M}_i^{-1}$。

代码：

#include<iostream>
#include<cstdio>
#include<cstdlib>

#define MAXN 1005

using namespace std;

int n;
long long a[MAXN], b[MAXN], sum = 1, ans;

void exgcd(long long a, long long b, long long &x, long long &y)
{
    if(b == 0)
    {
        x = 1;
        y = 0;
        return;
    }
    exgcd(b, a % b, x, y);
    int z = x;
    x = y;
    y = z - (a / b) * y;
}

int main()
{
    scanf("%d", &n);
    for(int i = 1; i <= n; i++)
    {
        scanf("%lld%lld", &a[i], &b[i]);
        sum *= a[i];
    }
    for(int i = 1; i <= n; i++)
    {
        long long m = sum / a[i], x = 0, y = 0;
        exgcd(m, a[i], x, y);
        ans += b[i] * m * (x < 0 ? x + a[i] : x);
    }
    printf("%lld", ans % sum);
    return 0;
}

---

如果不保证 $p_i$ 两两互质该怎么办？

假设现在我们已经得出了前 $i-1$ 个方程的通解 $x+kM(k\in Z)$，其中 $M$ 显然为 $LC{M}_{j=1}^{i-1}{p}_j$，那么对于第 $i$ 个方程，我们就是要找到一个 $k$ 满足：

$$x+kM\equiv a_i(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i)$$

显然这可以化为一个二元一次方程，用 $exgcd$ 求解即可（同时可以用裴蜀定理判断无解情况）。

代码：（用上了龟速乘避免溢出）

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cctype>

#define MAXN 100005

#define up(i, l, r) for(int i = (l); i <= (r); i = -~i)
#define down(i, l, r) for(int i = (l); i >= (r); i--)

using i64 = long long;
using namespace std;

namespace IO
{
	int f; char c;

	template<typename T> inline void _Read(T &v)
	{
		v = 0; f = 1; c = getchar();
		while(!isdigit(c)) { if(c == '-') f = -1; c = getchar(); }
		while(isdigit(c)) { v = (v << 3) + (v << 1) + (int)(c - '0'); c = getchar(); }
		v *= f;
		return;
	}

	template<typename T> inline void _Write(T k)
	{
		if(k < 0) { putchar('-'); k = -k; }
		if(k > 9) _Write(k / 10);
		putchar((char)(k % 10 + '0'));
		return;
	}

	inline int Read() { int v; _Read(v); return v; }
	inline i64 Read64() { i64 v; _Read(v); return v; }
	inline void Write(int v, char ed = '\n') { _Write(v); putchar(ed); return; }
	inline void Write64(i64 v, char ed = '\n') { _Write(v); putchar(ed); return; }
}

using IO :: Read;
using IO :: Read64;
using IO :: Write;
using IO :: Write64;

int n;
i64 b[MAXN], a[MAXN], mul, ans;

inline void Mul(i64 &res, i64 a, i64 b, i64 p)
{
	res = 0LL;
	while(b)
	{
		if(b & 1LL) res = (res + a) % p;
		a = (a << 1) % p;
		b >>= 1;
	}
	return;
}

inline i64 Exgcd(i64 a, i64 b, i64 &x, i64 &y)
{
	if(b == 0LL) return x = 1LL, y = 0LL, a;
	i64 c = Exgcd(b, a % b, y, x);
	return y -= (a / b) * x, c;
}

int main()
{
	n = Read();
	up(i, 1, n) a[i] = Read64(), b[i] = Read64();
	mul = a[1]; ans = b[1];
	up(i, 2, n)
	{
		i64 A = mul, B = a[i], C = ((b[i] - ans) % B + B) % B, x, y;
		i64 D = Exgcd(A, B, x, y), E = B / D;
		Mul(x, x, C / D, E);
		ans += mul * x; mul *= E;
		ans = (ans % mul + mul) % mul;
	}
	Write64(ans);
	return 0;
}

---

原根

• 定义部分

设 $a，p$ 互质且 $p>1$，则使得同余式 $a^n\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 成立的最小的 $n$ 被称为 $a$ 模 $p$ 的阶，记为 ${\delta }_p(a)$。

由欧拉定理 $a^{\phi (p)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，可知 ${\delta }_p(a)|\phi (p)$。

若存在 $g\perp p$ 使得 ${\delta }_p(g)=\phi (p)$，则称 $g$ 为模 $p$ 的原根。

• 原根存在定理

一个大于 $1$ 的整数 $n$ 存在原根，当且仅当 $n=2,4,p^{\alpha },2p^{\alpha }$（其中 $p$ 为奇素数，$\alpha$ 为正整数）。

（由于笔者太菜，这里的证明省去。事实上，这里的证明的确十分复杂，并不建议深究。）

• 原根个数

若 $n$ 有原根，其个数为 $\phi (\phi (n))$。

如何证明？

设任一原根为 $g$，而阶存在性质 ${\delta }_n(g^k)=\frac{{\delta }_n(g)}{\gcd ({\delta }_n(g),k)}$ （这里也略去证明，感兴趣的可以上网查找），说明如果 $k\perp \phi (n)$，那么 ${\delta }_n(g^k)={\delta }_n(g)=\phi (n)$，即 $g^k$ 也为 $n$ 的原根。

可以发现这样的 $k$ 恰好有 $\phi (\phi (n))$ 个，并且由于原根的定义，这些数互不相同。

而对于这样的一些 $k$，所有的 $g^{kv}\mathrm{m}\mathrm{o}\mathrm{d}n(v\in [1,\phi (n)])$ 也会互不相同，因为 $g^k$ 模 $n$ 的阶为 $\phi (n)$。

所以，这些数确实为 $n$ 的所有原根。

得证。

• 原根求法

我们先枚举求出 $n$ 的最小原根 $g$。

如何判定一个数 $g$ 是否为 $n$ 的原根呢？

如果 $n=2$，原根只能为 $1$，可以直接判定。

否则，$g$ 为模 $n$ 的原根的充要条件是：对于 $\phi (n)$ 的每个质因数 $d$，都不满足 $g^{\frac{\phi (n)}{d}}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}n)$。（再次略去证明）

这样，我们就可以先把 $\phi (n)$ 分解质因数，然后暴力枚举判断了。

为什么暴力枚举不会超时？

• 我国数学家王元于 $1959$ 年证明了若 $n$ 有原根，其最小原根是不超过 $n^{0.25}$ 级别的。所以枚举的时间复杂度在可承受范围内。

然后对于所有与 $\phi (n)$ 互质的 $k\in [1,\phi (n))$， $g^k\mathrm{m}\mathrm{o}\mathrm{d}n$ 均为 $n$ 的原根。

贴上洛谷模板题代码：

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cassert>
#include<algorithm>

#define MAXN 1000001

using namespace std;

int t, n, d, phi[MAXN], len, pr[MAXN], cnt, ans[MAXN];
bool np[MAXN], rt[MAXN];

inline void Sieve()
{
    for(register int i = 2; i < MAXN; i++) phi[i] = i;
    for(register int i = 2; i < MAXN; i++)
    {
        if(!np[i])
        {
            rt[i] = true; if((i << 1) < MAXN) rt[i << 1] = true; phi[i] = i - 1;
            for(register int j = 2; i * j < MAXN; j++) np[i * j] = true, phi[i * j] = phi[i * j] / i * (i - 1);
            if(i & 1 && i <= 1000)
            {
                int tmp = i * i;
                while(tmp < MAXN) { rt[tmp] = true; if((tmp << 1) < MAXN) rt[tmp << 1] = true; tmp *= i; }
            }
        }
    }
    return;
}

inline void Divide(int k)
{
    len = 0;
    for(register int i = 2; i * i <= k; i++)
    {
        if(k % i == 0)
        {
            pr[++len] = i;
            while(k % i == 0) k /= i;
        }
    }
    if(k > 1) pr[++len] = k;
    return;
}

inline int pw(int a, int b, int p)
{
    int res = 1;
    while(b)
    {
        if(b & 1) res = 1LL * res * a % p;
        a = 1LL * a * a % p; 
        b >>= 1;
    }
    return res;
}

inline bool Check(int a, int b, int p)
{
    for(register int i = 1; i <= len; i++) if(pw(a, b / pr[i], p) == 1) return false;
    return true;
}

inline int gcd(int a, int b) { if(b == 0) return a; return gcd(b, a % b); }

int main()
{
    #ifdef FILE
        freopen("Input.in", "r", stdin);
        freopen("Output.out", "w", stdout);
    #endif
    Sieve();
    scanf("%d", &t);
    while(t--)
    {
        cnt = 0;
        scanf("%d%d", &n, &d);
        if(!rt[n]) { puts("0\n"); continue; }
        Divide(phi[n]);
        int g; for(register int i = 1; i < n; i++) if(pw(i, phi[n], n) == 1 && Check(i, phi[n], n)) { g = i; break; }
        int val = 1;
        for(register int i = 1; i <= phi[n]; i++)
        {
            val = 1LL * val * g % n;
            if(gcd(i, phi[n]) == 1) ans[++cnt] = val;
        }
        printf("%d\n", cnt), sort(ans + 1, ans + cnt + 1);
        for(register int i = d; i <= cnt; i += d) printf("%d ", ans[i]);
        putchar('\n');
    }
    return 0;
}

---

二次剩余

下文仅介绍 $p$ 为奇素数时的 $Cipolla$ 算法。

如无特殊情况，所有的定义、运算都是在 $\mathrm{m}\mathrm{o}\mathrm{d}p$ 意义下的。

---

• 定义部分

如果对于非 $0$ 整数 $n$，使得 $x^2\equiv n(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 无解，称 $n$ 为二次非剩余，否则为二次剩余。

否则，我们来讨论一下解的情况。

设有一组不等的解为 $x_1,x_2$，则有 $x_1^2\equiv x_2^2$，移项后有 $(x_1-x_2)(x_1+x_2)\equiv 0$，而前者不等于 $0$，所以 $x_1,x_2$ 必为相反数。

并且模意义下，由于 $p$ 为奇素数，因此它们互相之间是不相等的。

也就是说，方程最多有两互异解，并且互为相反数。这样，我们就可以判断唯一解的情况了：直接判断 $n=0$ 是否成立，成立则 $x=0$ 为唯一解。

进一步地，我们发现：任意一对相反数对应唯一一个二次剩余，所以二次剩余的个数为 $\frac{p-1}{2}$，其它的非 $0$ 整数即为二次非剩余，数量与其相等。

• 欧拉判别准则

若 $x\perp p$ ，当 $x^{\frac{p-1}{2}}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 成立时，$x$ 为模 $p$ 的二次剩余，否则为二次非剩余。

即 “$x$ 是模 $p$ 的二次剩余” 是 “$x^{\frac{p-1}{2}}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 成立” 的充要条件。

充分性：设 $x\equiv k^2(\mathrm{m}\mathrm{o}\mathrm{d}p)$，由费马小定理可知 $x^{\frac{p-1}{2}}\equiv k^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

必要性：设 $g$ 为 $p$ 的一个原根，由于 $x\perp p$，则必定存在一个正整数 $k$ 使得 $x\equiv g^k(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 。又因为 $x^{\frac{p-1}{2}}\equiv g^{\frac{k(p-1)}{2}}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，所以 $\frac{k(p-1)}{2}$ 一定是 $\phi (p)=p-1$ 的倍数，说明 $k$ 是一个偶数。所以令 $m=g^{\frac{k}{2}}$ 就有 $x\equiv m^2(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

综上，我们只要判断 $x^{\frac{p-1}{2}}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 是否成立，就可以求出 $x$ 是否为二次剩余。

• 求法

判定是否为二次剩余后，我们怎样求出方程的解呢？

所以，下面介绍 $Cipolla$ 算法。

对于二次剩余方程 $x^2\equiv n(\mathrm{m}\mathrm{o}\mathrm{d}p)$，我们先随机出一个 $k$ 使得 $k^2-n$ 为二次非剩余，由于其个数与二次剩余相等，因此期望 $2$ 次就能找到。

然后我们类似虚数，定义 $i^2\equiv k^2-n(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

接下来摆定理：$x=(k+i{)}^{\frac{p+1}{2}}$ 为原方程的一个解。

证明：

首先通过二项式定理展开，再用 $Lucas$ 定理化简得 $(k+i{)}^p\equiv k^p+i^p(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

根据欧拉判别法有：
$$\begin{array}{rl}& (i^2{)}^{\frac{p-1}{2}}\equiv i^{p-1}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)\\ & i^p\equiv -i(\mathrm{m}\mathrm{o}\mathrm{d}p)\end{array}$$
由费马小定理可知 $k^p\equiv k(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

所以：
$$\begin{array}{rl}{x}^2& \equiv (k+i{)}^{p+1}\\ & \equiv (k+i{)}^p(k+i)\\ & \equiv (k^p+i^p)(k+i)\\ & \equiv (k-i)(k+i)\\ & =k^2-i^2=n\end{array}$$
得证。

贴上代码：

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cassert>
#include<ctime>

using LL = long long;
using namespace std;

int t;
LL si, n, p;

struct Complex
{
	LL a, b;
	
	Complex operator * (const Complex &t) const
	{
	 	Complex res = (Complex){0LL, 0LL};
		res.a = (a * t.a % p + b * t.b % p * si % p) % p;
		res.b = (a * t.b % p + b * t.a % p) % p;
		return res;
	}
};

inline LL Pow_Real(LL a, LL b)
{
	LL res = 1LL;
	while(b)
	{
		if(b & 1LL) res = res * a % p;
		a = a * a % p;
		b >>= 1;
	}
	return res;
}

inline LL Pow_Complex(Complex a, LL b)
{
	Complex res = (Complex){1LL, 0LL};
	while(b)
	{
		if(b & 1LL) res = res * a;
		a = a * a;
		b >>= 1;
	}
	return res.a;
}

inline LL Solve(LL n, LL p)
{
	if(Pow_Real(n, (p - 1LL) / 2LL) != 1LL) return -1LL;
	LL k;
	while(true)
	{
		k = 1LL * rand() * rand() % p;
		si = ((k * k % p - n) % p + p) % p;
		if(Pow_Real(si, (p - 1LL) / 2LL) != 1LL) break;
	}
	Complex val = (Complex){k, 1};
	return Pow_Complex(val, (p + 1LL) / 2LL);
}

int main()
{
	#ifdef FILE
		freopen("Input.in", "r", stdin);
	#endif
	scanf("%d", &t);
	while(t--)
	{
		scanf("%lld%lld", &n, &p), n %= p;
		if(n == 0LL) { puts("0"); continue; }
		LL ans1 = Solve(n, p), ans2;
		if(ans1 == -1LL) { puts("Hola!"); continue; }
		else
		{
			ans2 = p - ans1;
			if(ans1 > ans2) swap(ans1, ans2);
			printf("%lld %lld\n", ans1, ans2);
			continue;
		}
	}
	return 0;
}

---

$ToBeContinued\cdots \cdots$