【笔寄】n 次同余 - 洛谷保存站

为了赶时间有亿些证明没写，需要的话查看下面的链接

<https://blog.csdn.net/zxyoi_dreamer/article/details/85195819>

<https://www.luogu.com.cn/blog/wengweijie/solution-p5668>

## 目录

1. 0 次同余
2. 1 次同余
3. 2 次同余 / 剩余
   1. $ax^2+bx+c\equiv0\pmod{p^\alpha}$
   2. $y^2\equiv A\pmod{p^\alpha}$
   3. $x^2\equiv a\pmod m,(a,m)=1$
      1. 奇素数的平方（非）剩余
         1. Legendre
         2. Jacobi
         3. Cipolla 算法
      2. 任意模数的平方（非）剩余
4. n 次同余 / 剩余
   1. 指数、原根和指标
      1. $a^x\equiv b\pmod m,(a,m)=1$ / BSGS
      2. $a^x\equiv b\pmod m$ / exBSGS

[TOC]

## 1. 0 次同余

$a\equiv0\pmod p$

啊首先 0 次同余要说吗？

```cpp
#include <cstdio>
int main() {
	int a, p;
	scanf("%d %d", &a, &p);
	puts(a ? "NO" : "YES");
	return 0;
}
    
```

（

## 2. 1 次同余

$ax+b\equiv0\pmod p$

这是经典的 exgcd，[P1516 青蛙的约会](https://www.luogu.com.cn/problem/P1516)，由于不在这篇文章的讨论范围内，略（

## 3. 2 次同余

### 3.1. $ax^2+bx+c\equiv0\pmod p$

令 $f(x)=ax^2+bx+c$

首先应该判断是否有解  
设 $p$ 可以拆分为 $\prod\limits_{i=1}^kp_{i}$，其中 $p_i$ 两两互素，则 $f(x)\equiv0\pmod p$ 与  
$\left\{\begin{aligned}f(x)&\equiv0\pmod {p_1}\\f(x)&\equiv0\pmod{p_2}\\\cdots\end{aligned}\right.$  
等价。
因此我们考虑 $f(x)\equiv0\pmod{p^\alpha}$ 的解情况。

$\tiny\Downarrow\text{一条分割线}$

---

分类讨论：

1. 3.1.1. $p^\alpha\mid(a,b,c)$

易得，则任意整数都满足条件

> 定义 $p^\gamma\mid\mid k$ 为 $p^\gamma\mid k$ 且 $p^{\gamma+1}\nmid k$，称 $p^\gamma$ **恰整除** $k$。

2. 3.1.2. $p^\gamma\mid\mid(a,b,c),\gamma<\alpha$

此时 $ax^2+bx+c\equiv0\pmod{p^\alpha}$ 的解情况与 $\frac a{p^\gamma}x^2+\frac b{p^\gamma}x+\frac c{p^\gamma}\equiv0\pmod{p^{\alpha-\gamma}}$ 一样。

此时 $p\nmid(a,b,c)$

3. 3.1.3. $p\nmid(a,b,c)$

1. 3.1.3.1. $p\mid a,p\mid b$

此时 $p\nmid c$，$ax^2+bx+c\equiv0\pmod{p}$ 无解。

因此 $ax^2+bx+c\equiv0\pmod{p^\alpha}$ 也无解。

2. 3.1.3.2. $p\mid a,p\nmid b$

此时 $f'(x)\equiv0\pmod{p}$ 无解，$p\nmid f'(x)$。

所以 $f(x)\equiv0\pmod{p^\alpha}$ 与 $f(x)\equiv0\pmod p$ 等价

又因为 $p\mid a$ 所以又等价于 $bx+c\equiv0\pmod p$

又因为 $p\nmid b$，所以 $(p,b)=1$，由一次同余，肯定有解

3. 3.1.3.3. $p\nmid a,p>2$

此时 $(p^\alpha,4a)=1$，左右乘上 $4a$：

$y^2-A\equiv0\pmod{p^\alpha}$，其中 $A=b^2-4ac,y=2ax+b$。

此时 $f(x)\equiv0\pmod{p^\alpha}$ 有解说明 $y^2-A\equiv0\pmod{p^\alpha}$

证明一下充分性：

> 设 $y^2-A\equiv0\pmod{p^\alpha}$ 有一解 $y_0$
      >
      > 因为 $(2a,p^\alpha)=1$，$y=2ax+b\equiv y_0\pmod{p^\alpha}$
      >
      > 所以 $(2ax+b)^2-A\equiv0\pmod{p^\alpha}$ 有解
      >
      > 所以 $f(x)\equiv0\pmod{p^\alpha}$ 有解
      >
      > $\tiny\text{这里有点绕，你品}$

4. 3.1.3.4. $p\nmid a,p=2$

1. 3.1.3.4.1. $2\nmid b$

则 $f'(x)=2ax+b\equiv0\pmod2$ 无解，$2\nmid f'(x)$

所以 $f(x)\equiv0\pmod{2^\alpha}$ 与 $f(x)\equiv0\pmod2$ 等价。

由费马小定理，$x^2\equiv x\pmod2$，所以又等价于 $(a+b)x+c\equiv0\pmod2$

此时 $f(x)\equiv0\pmod{2^\alpha}$ 有解的充要条件为 $2\mid c$

2. 3.1.3.4.2. $2\mid b$

我们设 $b=2b_1$

> 下一行的 $\alpha$ 和 $a$ 是不一样的，请仔细区分！

因为 $(2^\alpha,a)=1$，原式与 $a\cdot ax^2+a\cdot 2b_1x+a\cdot c\equiv0\pmod{2^\alpha}$ 等价。

转换一下可得 $y^2-A\equiv0\pmod{2^\alpha}$，其中 $y=ax+b_1,A=b_1^2-ac$。

类比 3.1.3.3 的必要性证明可以证明这个必要性

---

$\tiny\Uparrow\text{一条分割线}$

所以刚刚我们干了什么？

我们将 $f(x)=ax^2+bx+c\equiv0\pmod{p^\alpha}$ 的判定转换为了 $y^2-A\equiv0\pmod{p^\alpha}$ 的判定！！1

### 3.2. $y^2-A\equiv0\pmod{p^\alpha}$

分类讨论：

1. 3.2.1. $p^\alpha\mid A$

原式等价于 $y^2\equiv0\pmod{p^\alpha}$，易求出解

2. 3.2.2. $p^\alpha\nmid A$

设 $p^\beta\mid\mid A,A=p^\beta A_1,\alpha>\beta\geq0$

当 $\beta\geq1$，$p\mid y$，设 $p^r\mid\mid y,y=p^rt$，则代入可得：

$p^{2r}t^2-p^\beta A_1\equiv0\pmod{p^\alpha},p\nmid t,p\nmid A_1$

因此 $(p^\beta A_1,p^\alpha)=(p^{2r}t^2,p^\alpha)$

又因为 $(p^\beta A_1,p^\alpha)=p^\beta$，所以 $\beta=\min(2r,a)$

又因为 $\alpha>\beta\geq0$，所以 $\beta=2r$，所以 $\beta$ 为偶数才有解。

所以原式是否有解与 $t^2-A_1\equiv0\pmod{p^{\alpha-\beta}}$ 有解等价。

其中 $(A_1,p^{\alpha-\beta})$

以上，我们成功的将普通的二次同余方程转换为了 $x^2\equiv a\pmod m,(a,m)=1$

> 定义，若 $x^2\equiv a\pmod m,(a,m)=1$ 有解，则称 $a$ 为模 $m$ 下的**平方剩余**，否则称为**平方非剩余**
>
> $\tiny\text{其实就是是否能在模意义下开方}$

### 3.3. $x^2\equiv a\pmod m,(a,m)=1$

#### 3.3.1. 奇素数的平方（非）剩余

奇素数，即 $p\in \text{prime},p>2$

**欧拉判别条件**：

1. $a$ 为模 $m$ 下的平方剩余的充要条件为 $a^\frac{p-1}2\equiv1\pmod p$。
2. $a$ 为模 $m$ 下的平方非剩余的充要条件为 $a^\frac{p-1}2\equiv-1\pmod p$。
3. 当 $a$ 为模 $m$ 下的平方剩余时，原式有且仅有**两个解。**

证明：~~我也想写证明，可是这里太小写不下~~以后再补。

模 $p$ 的剩余系中，平方剩余个数与平方非剩余个数都是 $\frac{p-1}2$，且这 $\frac{p-1}2$ 个平方剩余**分别**与序列 $1^2,2^2,\cdots,(\frac{p-1}2)^2$ 中恰好一个数同余。

证明：咕！

现在可以开始讲求解了，但在这之前，我们先了解以下 Legendre 与 Jacobi 符号。

为了什么？手算，所以你不看也没事。

#### 3.3.1.1. Legendre

定义 $\left(\frac ap\right)=\begin{cases}1,&a\text{ 是模 }p\text{ 下的平方剩余}\\-1,&a\text{ 是模 }p\text{ 下的平方非剩余}\\0,&p\mid a\\\end{cases}$

咕！

#### 3.3.1.2. Jaboci

定义 $\left(\dfrac am\right)=\prod\limits_{i=1}^r\left(\dfrac a{p_i}\right),m=\prod\limits_{i=1}^rp_i$，$p_i$ 是素数， $\left(\dfrac a{p_i}\right)$ 是 $a$ 对 $p_i$ 的 Legendre 符号。

咕！

#### 3.3.1.3. Cipolla 算法

对于任何的整数 $b$，令 $i^2=b^2-a$，若  $i^2$ 平方非剩余，则同余方程 $x^2\equiv a\pmod p,(a,p)=1$ 的解为 $x_0=(b+i)^\frac{p+1}2,x_1=p-x_0$。

证明：

> $$\begin{array}{rll}&x_0&\\\equiv&a^\frac12&\\\equiv&(b^2-(b^2-a))^\frac12&\\\equiv&(b^2-i^2)^\frac12&\\\equiv&((b-i)\cdot(b+i))^\frac12&\\\equiv&((b^p-i)\cdot(b+i))^\frac12&\scriptsize\text{欧拉定理}\\\equiv&((b^p+((i^2)^\frac{p-1}2\cdot i))\cdot(b+i))^\frac12&\scriptsize\text{欧拉判别公式}\\\equiv&((b^p+i^p)\cdot(b+i))^\frac12&\\\equiv&((b+i)^p\cdot(b+i))^\frac12&\scriptsize\text{二项式展开}\\\equiv&(b+i)^\frac{p+1}2&\pmod p\\\end{array}$$
>
> 因为 $(p-x)^2\equiv x^2\pmod p$，所以 $x_1=p-x_0$

实现的时候，由于你只知道 $i^2$ 的值而不知道 $i$ 的值，因此你需要像复数一样创建一个 `struct` 来计算快速幂，可以保证最后的结果不包含 $i$ 的项。

代码实现（[【模板】二次剩余](https://www.luogu.com.cn/problem/P5491)）：

```cpp
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <algorithm>
using namespace std;
long long p;
inline long long qmod(const long long a, const long long p) {
	// [-p, 2p-1]
	if (a >= p) return a - p;
	if (a < 0) return a + p;
	else return a;
}
long long ii;
struct complex {
	long long r, i;
	complex(long long R = 0, long long I = 0) {r = R; i = I;}
	complex operator + (complex sec) const {
		return complex(qmod(r + sec.r, p), qmod(i + sec.i, p));
	}
	complex operator * (complex sec) const {
		// (a+bi)*(c+di)=ac+i^2bd+(ad+bc)i
		return complex(qmod(r * sec.r % p + ii * i % p * sec.i % p, p), qmod(r * sec.i % p + i * sec.r % p, p));
	}
};
complex qkpow(complex a, long long b) {
	complex ans = 1;
	while (b) {
		if (b & 1) ans = ans * a;
		a = a * a; b >>= 1;
	}
	return ans;
}
long long qkpow(long long a, long long b) {
	long long ans = 1;
	while (b) {
		if (b & 1) ans = ans * a % p;
		a = a * a % p; b >>= 1;
	}
	return ans;
}
void Cipolla(long long a) {
	if (a == 0) {
		puts("0"); return;
	}
	if (qkpow(a, p >> 1) == p - 1) {
		puts("Hola!"); return;
	}
	long long b;
	do {
		b = rand() % p;
		ii = b * b - a;
		ii = (ii % p + p) % p;
	} while (qkpow(ii, p >> 1) != p - 1);
	complex c(b, 1);
	c = qkpow(c, p + 1 >> 1);
	long long ans[2] = {c.r, p - c.r};
	sort(ans, ans + 2);
	if (ans[0] == ans[1]) {
		printf("%lld\n", ans[0]);
	}
	else {
		printf("%lld %lld\n", ans[0], ans[1]); 
	}
}
int main() {
	int t;
	scanf("%d", &t);
	srand(time(NULL));
	while (t--) {
		long long a;
		scanf("%lld %lld", &a, &p);
		Cipolla(a);
	}
	return 0;
}
```

#### 3.3.1.4. Tonelli-Shanks 算法

要推主线！！1

要咕！！1

#### 3.3.2. 任意模数的平方（非）剩余

考虑拆分模数 $P$。

分类讨论：

1. 当 $P$ 为奇素数

详见 3.3.1。

2. 当 $P$ 为奇素数的幂

这时候设 $P=p^\alpha$。

这时有两种解决办法：

1. Tonelli-Shanks 算法改版：

咕！

2. ？算法：

分类讨论

1. $p\mid a$：

设 $a=p^cm,p\nmid m$，此时 $2\mid c$。

设 $x=p^tn,p\nmid n$

那么 $p^{2t}n2\equiv p^cn\pmod{p^\alpha}$

代码（[U124195 任意模数二次剩余](https://www.luogu.com.cn/problem/U124195)）：

时间复杂度为 $O(\log_2a+\frac{\sqrt m}{\ln m}\cdot\log_2\log_2m)$

```cpp
#include <cstdio>
#include <algorithm>
using namespace std;
long long a, p;
const int P = 100005;
int prime[P], pcnt;
bool mark[P];
void primeinit() {
	mark[1] = 1;
	for(int i = 2; i <= 100000; ++i) {
		if (!mark[i]) prime[++pcnt] = i;
		for (int j = 1; j <= pcnt && i * prime[j] <= 100000; ++j) {
			mark[i * prime[j]] = 1;
			if (i % prime[j] == 0) break;
		}
	}
}
long long mul(long long a, long long b, long long mod) {
	long long res = a * b - (long long)((long double)a / mod * b) * mod;
	return (res + mod) % mod;
//	   (a*b)%m
//	-> (a%m*b)
//	-> (a-a/m*m)*b
//	-> a*b-(a/m*b)*m
}
long long ii;
struct complex {
    long long r, i;
    complex(long long R = 0, long long I = 0) {r = R; i = I;}
    complex operator * (pair<complex, long long> sec) const {
        return complex((mul(r, sec.first.r, p) + mul(mul(ii, i, p), sec.first.i, p)) % p, (mul(r, sec.first.i, p) + mul(i, sec.first.r, p)) % p);
    }
};
complex qkpow(complex a, long long b, long long p) {
    complex ans = 1;
    while (b) {
        if (b & 1) ans = ans * make_pair(a, p);
        a = a * make_pair(a, p); b >>= 1;
    }
    return ans;
}
long long qkpow(long long a, long long b, long long p) {
    long long ans = 1;
    while (b) {
        if (b & 1) ans = mul(ans, a, p);
        a = mul(a, a, p); b >>= 1;
    }
    return ans;
}
long long exgcd(long long a, long long b, long long &x, long long &y) {
	if (b) {
		long long g = exgcd(b, a % b, y, x);
		y -= (a / b) * x; return g;
	}
	y = 0; x = 1; return a;
}
long long inv(long long a, long long mod) {
	long long x, y;
	exgcd(a, mod, x, y); // ax+ymod=1 -> ax=1(mod mod)
	return x;
}
long long Cipolla(long long a, long long p) {
	// solution for x^2=a(mod p) p>2 p=prime
	a %= p; long long b;
	if (qkpow(a, p >> 1, p) == p - 1) return -1;
	do {
		b = rand() % p;
		ii = (mul(b, b, p) - a + p) % p;
	} while (qkpow(ii, p >> 1, p) != p - 1);
	long long ans = qkpow(complex(b, 1), p + 1 >> 1, p).r;
	return min(ans, p - ans);
}
long long solutionPK(long long a, int k, long long p) {
	// solution for x^2=a(mod p^k) p>2 p=prime
	int pk = qkpow(p, k, 114514114514114514ll);
	if (a % pk == 0) return 0;
	long long mmod = 1; int times = 0;
	while (a % p == 0) a /= p, times++, mmod *= (times & 1) ? p : 1;
	if (times & 1) return -1;
	k -= times; pk /= mmod * mmod;
	long long res = Cipolla(a, p); if (res == -1) return -1;
	ii = a; complex pows = qkpow(complex(res, 1), k, pk);
	return mul(pows.r, inv(pows.i, pk), pk) * mmod;
}
long long solution2K(long long a, int k) {
	// solution for x^2=a(mod 2^k) 
	if ((a & ((1 << k) - 1)) == 0) return 0;
	a &= ((1 << k) - 1);
	int t = __builtin_ctz(a); long long res = 1; a >>= t;
	if ((a & 7) != 1) return -1;
	if (t & 1) return -1;
	k -= t;
	for (int i = 4; i <= k; i++) {
		res = (res + (a % (1 << i) - res * res) / 2) % (1 << k);
	}
	if (res < 0) res += 1 << k;
	return res << (t >> 1);
}
long long re[20], mods[20];
int cnts;
long long mainsolution(long long a, long long p) {
	a %= p; cnts = 0; long long pp = p;
	for (int i = 1; i <= pcnt && prime[i] * prime[i] <= p; i++) {
		long long k = 0;
		while (p % prime[i] == 0) p /= prime[i], k++;
		if (k == 0) continue;
		mods[++cnts] = qkpow(prime[i], k, 114514114514114514ll);
		if (i == 1) re[cnts] = solution2K(a, k);
		else re[cnts] = solutionPK(a, k, prime[i]);
		if (re[cnts] == -1) return -1;
	}
	if (p != 1) {
		mods[++cnts] = p;
		if (p == 2) re[cnts] = solution2K(a, 1);
		else re[cnts] = solutionPK(a, 1, p);
		if (re[cnts] == -1) return -1;
	}
	// CRT:
	long long ans = 0; p = pp;
	for (int i = 1; i <= cnts; i++) {
		ans = (ans + mul(mul(p / mods[i], inv(p / mods[i], mods[i]), p), re[i], p)) % p;
	}
	return ans;
}
int main() {
	primeinit();
	int t;
	scanf("%d", &t);
	while (t--) {
		scanf("%lld %lld", &a, &p);
		printf("%lld\n", mainsolution(a, p));
	}
	return 0;
}
```

## 4. n 次同余

### 4.1. 指数、原根和指标

由欧拉定理，若 $(a,m)=1,m>1$，则 $a^{\varphi(m)}\equiv1\pmod m$

因此，若 $(a,m)=1,m>1$，则存在正整数 $\gamma$ 使得 $a^{\gamma}\equiv1\pmod m$

所以存在最小的正整数满足上述要求。

定义：若 $m>1,(a,m)=1$，则是同余式 $a^{\gamma}\equiv1\pmod m$ 成立的最小正整数 $\gamma$ 称为 $a$ 对模 $m$ 的指数，又称 $a$ 对模 $m$ 的阶，记作 $\delta_m(a)$ 或 $\text{ord}_m(a)$。

定义：若 $\delta_m(a)=\varphi(m)$，则 $a$ 是模 $m$ 的一个原根（primitive root）。

一些关于指数与原根的定理：

/\* 日后再补 \*/

讲一下原根的判定：

> 对于任意模数 $m$，原根不一定存在。
>
> 可以证明，当且仅当 $m$ 为 $2,4,p^\alpha,2p^\alpha$（$p$ 是奇素数）之一时，才存在它的原根。

证明：

> 补

原根判定定理：

> 设 $m>1$，$\varphi(m)$ 的所有不同素因数是 $q_1,q_2,\cdots,q_k$，如果 $(g,m)=1$，那么 $g$ 是模 $m$ 的一个的充要条件是：
>
> $g^\frac{\varphi(m)}{q_i}\not\equiv1\pmod m$

证明：

> 补

因此代码实现直接模拟即可

代码（[SP3713 PROOT - Primitive Root](https://www.luogu.com.cn/problem/SP3713)）：

~~这个代码是 wa 的，明天（5 月 3 日）调。~~ ac 了。

时间复杂度 $O(T\cdot \log_2n)$

```cpp
#include <cstdio>
const int MAXHP = 46666;
int prime[MAXHP], pcnt;
bool flag[MAXHP];
void primeinit() {
	for (int i = 2; i <= 46340; i++) {
		if (flag[i] == 0) prime[++pcnt] = i;
		for (int j = 1; j <= pcnt && prime[j] * i <= 46340; j++) {
			flag[i * prime[j]] = 1;
			if (i % prime[j] == 0) break;
		}
	}
}
int fraclist[10], fcnt;
int qkpow(int a, int b, int p) {
	int ans = 1;
	while (b) {
		if (b & 1) ans = 1ll * ans * a % p;
		a = 1ll * a * a % p; b >>= 1;
	}
	return ans;
}
int gcd(int a, int b) {
	return b ? gcd(b, a % b) : a;
}
int main() {
	primeinit();
	int n, m;
	while (scanf("%d %d", &n, &m), n || m) {
		fcnt = 0;
		int phin = n - 1, pphin = n - 1;
		for (int i = 1; i <= pcnt && prime[i] * prime[i] <= pphin; i++) {
			if (pphin % prime[i] == 0) fraclist[++fcnt] = prime[i];
			while (pphin % prime[i] == 0) pphin /= prime[i];
		}
		if (pphin > 1) fraclist[++fcnt] = pphin;
		while (m--) {
			int m;
			scanf("%d", &m);
			if (gcd(m, n) != 1) {
				puts("NO"); continue;
			}
			bool flg = 1;
			for (int i = 1; i <= fcnt; i++) {
				if (qkpow(m, phin / fraclist[i], n) == 1) {
					flg = 0; break;
				}
			}
			puts(flg ? "YES" : "NO");
		}
	}
	return 0;
}
```

刚刚我们已经学会了如何判定原根

如果想找一个最小原根的话其实可以按照上述方法暴力找的

由证明，一个数 $m$ 如果拥有原根，则其最小原根 $\leq m^\frac14$。

因此我们可以按照上面的方法大胆暴力

至此，时间复杂度为 $O(m^\frac14\cdot \log_2^2m)$。

我们获得了一个特解之后可以推出通解：

> 假设我们知道了的模 $m$ 下的最小原根 $g$，则所有的原根即为 $g^x\bmod m$，其中 $(x,\varphi(m))=1$。

通过这个，我们可以判断 $\varphi(m)$ 次互素即可。

时间复杂度为 $O(\varphi(m)\log_2m)$，考虑优化。

由于之前求最小原根的时候我们已经将 $\varphi(m)$ 素因数拆分过了，所以我们可以类似于埃式筛一样晒出与 $\varphi(m)$ 互素的数。时间复杂度 $O(\varphi(m)\log_2\log_2\varphi(m))$，

线性筛可以做到 $O(\varphi(m))$ 的时间复杂度，但是由于上述算法优越的常数，线性筛实际运行效果较慢。

总共有 $\varphi(\varphi(m))$ 个数。

为了使复杂度不爆表，排序换成基数排序。

代码（加优化）（[P6091 【模板】原根](https://www.luogu.com.cn/problem/P6091)）：

自带大常数（

```cpp
#include <cstdio>
#include <cstring>
using namespace std;
const int MAXHP = 1005;
int prime[MAXHP], pcnt;
bool flag[MAXHP];
void primeinit() {
	for (int i = 2; i <= 1000; i++) {
		if (flag[i] == 0) prime[++pcnt] = i;
		for (int j = 1; j <= pcnt && prime[j] * i <= 1000; j++) {
			flag[i * prime[j]] = 1;
			if (i % prime[j] == 0) break;
		}
	}
}
int getphi(int n) {
	int ans = n;
	for (int i = 1; i <= pcnt && prime[i] * prime[i] <= n; i++) {
		if (n % prime[i] == 0) {
			ans = ans / prime[i] * (prime[i] - 1);
			while (n % prime[i] == 0) n /= prime[i];
		}
	}
	if (n > 1) ans = ans / n * (n - 1);
	return ans;
}
int qkpow(int a, int b, int p) {
	int ans = 1;
	while (b) {
		if (b & 1) ans = 1ll * ans * a % p;
		a = 1ll * a * a % p; b >>= 1;
	}
	return ans;
}
int frac[20], top, tm[20];
int f[20], top2, t[20];
int gcd(int a, int b) {
	return b ? gcd(b, a % b) : a;
}
bool check(int n, int phin, int m) {
	if (gcd(m, n) != 1) return 0;
	bool flg = 1;
	for (int i = 1; i <= top; i++) {
		if (qkpow(m, phin / frac[i], n) == 1) {
			flg = 0; break;
		}
	}
	return flg;
}
bool shaizi[1000005];
void shai(int phin) {
	for (int i = 1; i <= top; i++) {
		for (int j = frac[i]; j <= phin; j += frac[i]) {
			shaizi[j] = 1;
		}
	}
}
int _b_[1000005];
void sort(int *a, int n) {
        int i,cnt0[256]={0},cnt1[256]={0},cnt2[256]={0},cnt3[256]={0}; 
        for(i=0;i<n;i++)
                cnt1[(a[i]>>8)&255]++,cnt2[(a[i]>>16)&255]++,cnt3[(a[i]>>24)&255]++,cnt0[a[i]&255]++;
	for(i=1;i<256;i++) cnt0[i]+=cnt0[i-1],cnt1[i]+=cnt1[i-1],cnt2[i]+=cnt2[i-1],cnt3[i]+=cnt3[i-1];
	for(i=n-1;~i;i--) _b_[--cnt0[a[i]&255]]=a[i];
	for(i=n-1;~i;i--) a[--cnt1[(_b_[i]>>8)&255]]=_b_[i];
	for(i=n-1;~i;i--) _b_[--cnt2[(a[i]>>16)&255]]=a[i];
	for(i=n-1;~i;i--) a[--cnt3[(_b_[i]>>24)&255]]=_b_[i];
}
int rs[1000005], ttp;
long long qkkpow[1000005];
int main() {
	primeinit();
	int T; scanf("%d", &T);
	while (T--) {
		top = top2 = 0;
		int n, d; scanf("%d %d", &n, &d);
		if (n == 2) {
			puts("1");
			if (d == 1) puts("1"); else puts("");
		}
		else if (n == 4) {
			puts("1");
			if (d == 1) puts("3"); else puts("");
		}
		else {
			int phin = getphi(n), phinn = phin;
			for (int i = 1; i <= pcnt && prime[i] * prime[i] <= phinn; i++) {
				if (phinn % prime[i] == 0) {
					frac[++top] = prime[i];
					tm[top] = 0;
					while (phinn % prime[i] == 0) {
						phinn /= prime[i];
						tm[top]++;
					}
				}
			}
			if (phinn > 1) frac[++top] = phinn, tm[top] = 1;
			int nn = n;
			for (int i = 1; i <= pcnt && prime[i] * prime[i] <= nn; i++) {
				if (nn % prime[i] == 0) {
					f[++top2] = prime[i];
					t[top2] = 0;
					while (nn % prime[i] == 0) {
						nn /= prime[i];
						t[top2]++;
					}
				}
			}
			if (nn > 1) f[++top2] = nn, t[top2] = 1;
			if (top2 > 2) puts("0"), puts("");
			else if (top2 == 1 && f[1] == 2) puts("0"), puts("");
			else {
				if (top2 == 2 && (f[1] != 2 || t[1] != 1)) puts("0"), puts("");
				else {
					int ming = 1;
					while (check(n, phin, ming) == 0) {
						ming++;
					}
					ttp = 1; rs[1] = ming;
					shai(phin);
					qkkpow[1] = ming;
					for (int i = 2; i <= phin; i++) {
						qkkpow[i] = qkkpow[i - 1] * ming % n;
						if (shaizi[i] == 0) {
							rs[++ttp] = qkkpow[i];
						}
						else shaizi[i] = 0;
					}
					printf("%d\n", ttp);
					sort(rs + 1, ttp);
					for (int i = d; i <= ttp; i += d) {
						printf("%d ", rs[i]);
					}
					puts("");
				}
			}
		}
	}
	return 0;
}
/*
10
999959 114514
999959 114514
999959 114514
999959 114514
999959 114514
999959 114514
999959 114514
999959 114514
999959 114514
999959 114514
*/
```

指标定义自己去翻（怒）

#### 4.1.1. $a^x\equiv b\pmod m,(a,m)=1$ / BSGS

如果按照指标的定义去求，恭喜你，时间复杂度 $O(m)$。

所以我们考虑将此式子平衡规划一下。

因为 $(a,m)=1$，我们将 $x$ 拆成 $A\lceil\sqrt m\rceil-B$，其中 $A,B\in[0,\lceil\sqrt m\rceil]$。

由于逆元存在， 原式等价于 $a^{A\lceil\sqrt m\rceil}\equiv ba^B\pmod m$

预处理 $ba^B$，然后枚举 $A$，求出 $a^{A\lceil\sqrt m\rceil}$，按照合法的 $A,B$ 推回 $x$。

如果使用 map，时间复杂度 $O(\sqrt m\log_2m)$。

如果手写 hash 表，时间复杂度 $O(\sqrt m)$，非常优秀（？）。

实现的时候可以先存一下 $a^{\sqrt m}$，然后每次累乘 $a^{\sqrt m}$ 即可 $O(\sqrt m)$ 得到 $a^{A\lceil\sqrt m\rceil}$。

代码（[P3846 [TJOI2007] 可爱的质数 /【模板】BSGS](https://www.luogu.com.cn/problem/P3846)）：

```cpp
#include <cstdio>
#include <cmath>
using namespace std;
int p = 1000000007;
int sqrtq = 31623;
namespace hash {
	const int p1 = 1919839, p2 = 1919843, p3 = 2947331;
	int tong[p1], tong2[p2], tong3[p3];
	int mem[100000], mem2[100000], mem3[100000], tp;
	void clear() {
		while (tp--) {
			tong[mem[tp]] = tong2[mem2[tp]] = tong3[mem3[tp]] = 0;
		}
	}
	void insert(int x, int v) {
		mem[tp] = x % p1;
		mem2[tp] = x % p2;
		mem3[tp++] = x % p3;
		tong[x % p1] = v;
		tong2[x % p2] = v;
		tong3[x % p3] = v;
	}
	int get(int x) {
		if (tong[x % p1] == tong2[x % p2] && tong[x % p1] == tong3[x % p3]) {
			return tong[x % p1];
		}
		else return 0;
	}
}
int main() {
	int t;
//	scanf("%d", &t);
	t = 1;
	while (t--) {
		hash::clear();
		scanf("%d", &p);
		sqrtq = sqrt(p) + 0.99;
		int a, b;
		scanf("%d %d", &a, &b);
		long long pt = b; hash::insert(pt, 0);
		for (int i = 1; i <= sqrtq; i++) {
			pt = pt * a % p;
			hash::insert(pt, i);
		}
		long long ppt = 1;
		for (int i = 1; i <= sqrtq; i++) {
			ppt = ppt * a % p;
		}
		pt = ppt; bool flg = 0; int B;
		for (int i = 1; i <= sqrtq; i++) {
			if (B = hash::get(pt)) {
				flg = 1; printf("%d\n", i * sqrtq - B); break;
			}
			pt = pt * ppt % p;
		}
		if (flg == 0) puts("no solution");
	}
	return 0;
}
/*
5
6 2
2 3
2 6
3 4
7 9
*/
```

#### 4.1.2 $a^x\equiv b\pmod m,(a,m)\neq1$ / exBSGS

我们令 $d_1=(a,m)\neq1$，则根据裴蜀定理，当 $d_1\nmid b$ 时无解。

当 $d_1\mid b$，有 $\frac a{d_1}\cdot a^{x-1}\equiv \frac b{d_1}\pmod{\frac m{d_1}}$。

因为 $(\frac a{d_1},\frac m{d_1}) = 1$，逆元存在，所以 $a^{x-1}\equiv \frac b{d_1}\cdot(\frac a{d_1})^{-1}\pmod{\frac m{d_1}}$。

当 $(a,\frac m{d_1})=1$ 时可以使用 BSGS 求解。

否则，继续令 $d_2=(a,\frac m{d_1})$……

最终得到：
$$
D = \prod\limits_{i=1}^kd_i
$$

$$
a^{x-k}\equiv\frac bD\cdot(\frac{a^k}D)^{-1}\pmod{\frac mD},(a,\frac mD)=1
$$

> 注意上式中存在 $a^{x-k}$，而 $k$ 是根据 $a,m$ 而定的。
>
> 因此可能存在 $x<k$ 的情况
>
> 如果你的 BSGS 可以解决解为 $0$ 的情况那么只需要枚举 $0\sim k-1$ 的情况
>
> 否则枚举 $0\sim k$ 的情况
>
> 为什么时间复杂度正确？$k\leq\log_2m$

使用 BSGS 即可解决。

时间复杂度：由于每次需要除 gcd，因此时间复杂度为 $O(\log_2^2m)$

总时间复杂度： $O(\log_2^2m+\sqrt m)$

注意有一些题目要求 $0^0\equiv0\pmod1$，就极其离谱（

代码（[P4195 【模板】扩展 BSGS/exBSGS](https://www.luogu.com.cn/problem/P4195)）：

```cpp
#include <cstdio>
#include <cmath>
using namespace std;
int gcd(int a, int b) {
	return b ? gcd(b, a % b) : a;
}
int qkpow(long long a, int b, int m) {
	long long ans = 1;
	while (b) {
		if (b & 1) ans = ans * a % m;
		a = a * a % m; b >> 1;
	}
	return ans;
}
namespace hash {
	const int p1 = 1919839, p2 = 1919843, p3 = 2947331;
	int tong[p1], tong2[p2], tong3[p3];
	int mem[100000], mem2[100000], mem3[100000], tp;
	void clear() {
		while (tp) {
			tp--;
			tong[mem[tp]] = tong2[mem2[tp]] = tong3[mem3[tp]] = 0;
		}
	}
	void insert(int x, int v) {
		mem[tp] = x % p1;
		mem2[tp] = x % p2;
		mem3[tp++] = x % p3;
		tong[x % p1] = v;
		tong2[x % p2] = v;
		tong3[x % p3] = v;
	}
	int get(int x) {
		if (tong[x % p1] == tong2[x % p2] || tong[x % p1] == tong3[x % p3] || tong2[x % p2] == tong3[x % p3]) {
			if (tong[x % p1] == tong2[x % p2]) return tong[x % p1];
			else return tong3[x % p3];
		}
		else return 0;
	}
}
int BSGS(int a, int b, int p, int ext = 1) {
	hash::clear();
	int sqrtq = sqrt(p) + 0.99;
	a %= p; b %= p;
	if (a == 0) {
		if (b == 0) return 1;
		else return -1;
	}
	if (b == ext) {
		return 0;
	}
	long long pt = b; hash::insert(pt, 0 + 1);
	for (int i = 1; i < sqrtq; i++) {
		pt = pt * a % p;
		hash::insert(pt, i + 1);
	}
	long long ppt = 1;
	for (int i = 1; i <= sqrtq; i++) {
		ppt = ppt * a % p;
	}
	pt = ppt * ext % p; bool flg = 0; int B;
	for (int i = 1; i <= sqrtq; i++) {
		if (B = hash::get(pt)) {
			flg = 1; return i * sqrtq - B + 1;
		}
		pt = pt * ppt % p;
	}
	if (flg == 0) return -1;
	else return -2;
}
int exBSGS(int a, int b, int m) {
	a %= m; b %= m;
	if (b == 1) return 0;
	if (m == 1) return 0;
	int d = gcd(a, m), k = 0;
	long long inva = 1;
	while (d != 1) {
		if (b % d != 0) return -1;
		b /= d, m /= d;
		k++; inva = (inva * (a / d)) % m;
		if (inva == b) return k;
		d = gcd(a, m);
	}
	int ans = BSGS(a, b, m, inva);
	if (ans == -1) return -1;
	else return ans + k;
}
int main() {
	int a, p, b;
	while (scanf("%d %d %d", &a, &p, &b), a || p || b) {
		int ans = exBSGS(a, b, p);
		if (~ans) printf("%d\n", ans);
		else puts("No Solution");
	}
	return 0;
}
```

正在渲染内容...