高次同余方程，二次同余方程学习笔记

写在前面

文章作者实力有限，本文可能有个别错误，如有错误请友好地指出。
高次同余方程就是(x^aequiv b(mod p))
二次同余方程就是(x^2 equiv b(mod p))
我们接下来讨论解这两种方程的方法。
那么有一个问题。既然知道了高次同余方程的解法，就可以直接用解高次同余的方法解二次剩余方程。为什么要单独学二次同余方程呢。
因为我区间加区间修改用的是线段树不是树套树。即问题特殊化之后可以使用一些特殊的方法，这种方法可能会比一般方法高效，简便。

正文

高次同余方程

首先需要知道阶和原根

阶

定义

满足(a^xequiv 1(mod p))的最小的正整数x是a关于模p的阶，接下来(a)的阶表示为(<a>)。

条件

(gcd(a,p)=1)，这是显然的。

性质

1：(<a>mid varphi(p))。
证明：
因为(a^{varphi(p)}equiv 1(mod p))，(a^{<a>}equiv 1(mod p))
所以(a^{varphi(p)-<a>}equiv 1(mod p))
进而得出(a^{varphi(p) mod <a>}equiv 1(mod p))
假设(<a>)不是(varphi(p))的约数。
(varphi(p) mod<a> eq 0)
((varphi(p) mod<a>) < (<a>))与(<a>)的最小性矛盾，与假设矛盾原命题成立。
证毕。

2：设a关于模p的阶为(<a>)则(a^0)，(a^1)，(a^2)，...，(a^{<a>-1})两两不同。
证明：
假设(a^xequiv1(mod p))，(a^yequiv 1(mod p)),且(a^xequiv a^y(mod p))，(x<y<<a>)
则(a^{y-x}equiv 1(mod p))。(y-x<<a>)与阶的最小性矛盾，假设不成立，原命题成立。
证毕。

3：设a关于模p的阶为(<a>)，则(a^xequiv a^y(mod p))的充要条件是(xequiv y(mod <a>))。
两边一直除(a^{<a>})因为(a^{<a>}equiv 1(mod p))所以结论显然。

原根

（因为原根的定义涉及到阶，所以默认互质，又因为原根的性质3，p为质数）

定义

满足g关于模p的阶等于(varphi(p))的g是p的一个原根。

性质

1：设g是模p的一个原根，那么对于任意一个小于p的正整数x，都存在唯一一个值r使得r是以g为底的x对模p的指数。
其实就是说(g^requiv x(mod p))的r是唯一的，反过来也一样。其中g为p的一个原根，(x<p)，(r<varphi (p))其实就是阶的第二个性质。不再赘述。
2：对于一个质数，所拥有的原根个数是(varphi(varphi(p)))。
3：一个模数存在原根的从要条件是，这个模数可以表示为(1,2,4,p,2p,p^n)其中(p)是奇质数，n是任意正整数
（这两个性质的证明可能会用到抽象代数的知识，我不会555...）

---

那么如何找原根？

考虑到最小的原根一般比较小，所以我们枚举，原根然后判断。
假设我们当前枚举到(i)，判断的时候只需要枚举(varphi(p))的质因子(p1,p2,p3...pk)。然后判断(i^{frac{varphi(p)}{pj}})是不是全部都不是1，如果全部都不是1，(i)就是(p)的一个原根。
为什么？
我们这个过程其实是枚举(varphi(p))的约数，因为阶的性质1，所以i的阶是(varphi(p))的一个约数，我们只需要看有没有一个(varphi(p))的约数(x)使得(i^x equiv 1(mod p))就能判断i是不是原根了。
那为什么不直接枚举约数，用这种枚举方法对吗？
(varphi(p))的约数至少比(varphi(p))缺少一个质因子，因为如果(x^aequiv 1)那么(x^{ab}equiv 1)，所以少多个质因子的情况被少一个质因子的情况包含，所以枚举所有少掉一个质因子的情况即(frac{varphi(p)}{pj})就行。

---

解高次同余方程的步骤

首先要会BSGS和exgcd
(x^aequiv b(mod p))
1：先求出(p)的一个原根(g)。
2：求出以(g)为底的b关于模p的指数r，即(g^requiv b(mod p))（离散对数，用BSGS求解）
3：令(xequiv g^y(mod p)),带回原式就是求(g^{ay} equiv g^r(mod p))
4：根据阶的性质3就是求(ayequiv r(mod p))（一次同余方程用exgcd求解）
5：在0~(varphi(p)-1)中求得y的解，带回(xequiv g^y(mod p))解完方程。

有些值得注意的点：

因为原根的性质1求出的答案没有重复。
为什么要用原根代换b，别的数不行吗？我们需要用原根表示b。见原根的性质1。
因为要求原根且p为质数所以试用的p的范围很小，见原根的性质3。
只要求求出通解的复杂度为(O(sqrt p))
下面是输出所有解的代码。。。

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=201010;
unordered_map<int,int> hsh;
int prime[N],tot,T,phi,ans[N],num,a,b,p,g,x,y;
int read(){
	int sum=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){sum=sum*10+ch-'0';ch=getchar();}
	return sum*f;
}
void work(int x){
	int tmp=x;
	tot=0;
	for(int i=2;i*i<=x;i++){
		if(tmp%i==0){
			prime[++tot]=i;
			while(tmp%i==0)tmp/=i;
		}
	}
	if(tmp>1)prime[++tot]=tmp;
}
int ksm(int x,int b){
	int tmp=1;
	while(b){
		if(b&1)tmp=tmp*x%p;
		b>>=1;
		x=x*x%p;
	}
	return tmp;
}
int BSGS(int a,int b){
	hsh.clear();
	int block=sqrt(p)+1;
	int tmp=b;
	for(int i=0;i<block;i++,tmp=tmp*a%p)hsh[tmp]=i;
	a=ksm(a,block);
	tmp=1;
	for(int i=0;i<=block;i++,tmp=tmp*a%p){
		if(hsh.count(tmp)&&i*block-hsh[tmp]>=0)return i*block-hsh[tmp];
	}
}
int exgcd(int &x,int &y,int a,int b){
	if(b==0){
		x=1;y=0;
		return a;
	}
	int gcd=exgcd(x,y,b,a%b);
	int z=x;
	x=y;y=z-(a/b)*y;
	return gcd;
}
signed main(){
	T=read();
	while(T--){
		p=read(),a=read(),b=read();
		b%=p;
		phi=p-1;
		work(phi);
		for(int i=1;i<=p;i++){
			bool flag=false;
			for(int j=1;j<=tot;j++)if(ksm(i,phi/prime[j])==1){flag=true;break;}
			if(flag==false){g=i;break;}
		}
		int r=BSGS(g,b);
		int gcd=exgcd(x,y,a,phi);
		if(r%gcd!=0){printf("No Solution
");continue;}
		x=x*r/gcd;
		int k=phi/gcd;
		x=(x%k+k)%k;
		num=0;
		while(x<phi){ans[++num]=ksm(g,x),x+=k;}
		sort(ans+1,ans+1+num);
		for(int i=1;i<=num;i++)printf("%lld ",ans[i]);
		printf("
");
	}
	return 0;
}

二次同余方程

（p为素质数至于为什么要有素，接着看就知道了）
(x^2equiv A(mod p))

一些引理

引理1：x有正整数解的充要条件是(A^{frac{p-1}{2}} equiv 1(mod p))
证明：
由费马小定理得

[A^{p-1} equiv 1(mod p) ]

[A^{p-1}-1 equiv 0(mod p) ]

[(A^{frac{p-1}{2}}+1)+(A^{frac{p-1}{2}}-1) equiv 0 Rightarrow A^{frac{p-1}{2}}equiv pm1(mod p) ]

设(x^2 equiv A(mod p))则(A^{frac{p-1}{2}} equiv x^{p-1}(mod p))
由费马小定理

[x^{p-1} equiv 1(mod p) ]

[A^{frac{p-1}{2}}equiv 1(mod p) ]

证毕。
推论：x无正整数解的时(A^{frac{p-1}{2}} equiv -1(mod p))
引理2：一共只有((p-1)/2)（不考虑0）个数存在模(p)意义下的二次方根，对于其中的每个数都存在两个互为相反数的解。
证明：
假设有两个数u，v，(u^2-v^2 equiv 0(mod p))那么((u+v)*(u-v) equiv 0 (mod p))，因为u，v(subset [1,p-1])（不考虑0），所以(-p+1leq u-vleq p-1)无法被(p)整除，所以(u+v equiv 0(mod p))。
所以u，v互为相反数且对应唯一一个A，反过来也一样。因为不考虑0且p为奇数，这样的数有(frac{p-1}{2})对，且对应着(frac{p-1}{2})个A。
证毕。
此时我们就可以有一个解法。

解法一

1：我们求出(p)的一个原根(g)。然后求出以(g)为底的A关于模p的指数r，即(g^requiv A(mod p))（离散对数，用BSGS求解）
2：(A)满足(A^{frac{p-1}{2}}equiv1(mod p))带入得(g^{frac{r(p-1)}{2}}equiv 1(mod p))由原根的定义得(p-1 mid frac{r(p-1)}{2} Rightarrow r mid 2)。然后(g^{frac{r}{2}})就是一个解，因为(rmid2)所以可以直接快速幂，另一个解就是相反数。

---

又一些引理

定义：如果(x^2 equiv n (mod p))有解就说n是模p意义下的一个二次剩余，无解就说n是模p意义下的一个非二次剩余。
然后就是一个非常神的推导。
我们定义设a满足(w=a^2-n)是一个非二次剩余。
引理3：((a+sqrt{w})^p equiv a-sqrt{w}(mod p))这实际上拓展了数域，可以用复数去理解。
证明：
把((a+sqrt{w})^p)用二项式定理展开为(sum_{j=0}^p C_p^ja^jw^{p-j})。因为p是一个质数所以当(j!=0且j!=p)时(C_p^j)都为(0)否则为1，因为(C_p^j)分子中有(p)消不掉最后(mod p)得0。那么((a+sqrt{w})^p equiv a^p+sqrt{w}^pequiv a^{p-1}*a+w^{frac{p-1}{2}}*sqrt{w}(mod p))由费马小定理和引理1的推论(别忘了w是一个非二次剩余）最后化简为(a-sqrt{w})。

推导

[(a+sqrt{w})^{p+1}=(a+sqrt{w})^{p}*(a+sqrt{w}) ]

[=(a-sqrt{w})(a+sqrt{w}) ]

[=a^2-w ]

[=a^2-(a^2-n) ]

[=n ]

然后我们随机一个(a)满足条件，(xequiv (a+sqrt{w})^{frac{p+1}{2}}(mod p))根据拉格朗日定理可得答案中w的系数必然为0（然而我并不知道拉格朗日定理，所以模板拓展数域的乘法是背的），就可以解得(x)了,又因为引理2所以期望随机次数为2。复杂度可以保证。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstdlib>
#include<ctime>
using namespace std;
#define random(a,b) (rand()%(b-a+1)+a)
#define int long long
int n,p,w;
bool flag;
struct complex{
	int x,y;
	complex (int xx=0,int yy=0){
		x=xx;y=yy;
	}
};
complex operator *(complex a,complex b){
	return complex(((a.x*b.x%p+w*a.y%p*b.y%p)%p+p)%p,((a.x*b.y%p+a.y*b.x%p)%p+p)%p);
}
int read(){
	int sum=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){sum=sum*10+ch-'0';ch=getchar();}
	return sum*f;
}
complex ksm(complex x,int b){
	complex tmp(1,0);
	while(b){
		if(b&1)tmp=tmp*x;
		b>>=1;
		x=x*x;
	}
	return tmp;
}
int ksm(int x,int b){
	int tmp=1;
	while(b){
		if(b&1)tmp=tmp*x%p;
		b>>=1;
		x=x*x%p;
	}
	return tmp;
}
int work(int n){
	if(p==2)return n;
	if(ksm(n,(p-1)/2)+1==p)flag=true;
	int a;
	while(233){
		a=random(0,p-1);
		w=((a*a-n)%p+p)%p;
		if(ksm(w,(p-1)/2)+1==p)break;
	}
	complex res(a,1);
	complex ans(0,0);
	ans=ksm(res,(p+1)/2);
	return ans.x;
}
signed main(){
	srand((unsigned)time(NULL));
	p=read();n=read();
	n%=p;
	int ans1=work(n);
	int ans2=p-ans1;
	if(flag){printf("No Solution
");return 0;}
	if(ans1==ans2)printf("%lld
",ans1);
	else printf("%lld %lld",min(ans1,ans2),max(ans1,ans2));
	return 0;
}