二次同余方程的解

算法

问题是解方程\(x^2 \equiv n \ (\bmod p)\)，其中\(p\)是奇质数。

引理：\(n^{\frac{p-1}2}\equiv \pm 1\ (\bmod p)\)

证明：由费马小定理，\(n^{p-1}-1\equiv (n^\frac{p-1}2-1)(n^\frac{p-1}2+1)\equiv 0\ (\bmod p) \Rightarrow n^\frac{p-1}2 \equiv \pm 1\ (\bmod p)\)

引理：方程有解当且仅当\(n^\frac{p-1}2 \equiv 1\)

定理：设\(a\)满足\(\omega=a^2-n\)不是模\(p\)的二次剩余，那么\(x=(a+\sqrt{\omega})^\frac{p+1}2\)是二次剩余方程\(x^2 \equiv n \ (\bmod p)\)的解。

证明：由\((a+\sqrt{\omega})^p=a^p+\omega^\frac{p-1}2\sqrt{w}=a-\sqrt{\omega}\)，前面的等号用二项式定理和\(\binom pi=0\ (\bmod p),i \neq 0,p\)得到，后面的等号用费马小定理和\(\omega\)是模\(p\)的二次非剩余。然后

\[(a+\sqrt{\omega})^{p+1}\\ \equiv(a+\sqrt\omega)^p(a+\sqrt\omega)\\ \equiv(a-\sqrt\omega)(a+\sqrt\omega)\\ \equiv a^2-\omega\equiv n\ (\bmod p) \]

在算法实现的时候，对\(a​\)的选择可以随机，因为大约有一半数是模\(p​\)的二次非剩余。然后快速幂即可。

Timus Online Judge 例题

1132. Square Root

Time limit: 1.0 second
Memory limit: 64 MB

The number x is called a square root of a modulo n (root(a,n)) if x*x = a (mod n). Write the program to find the square root of number a by given modulo n.

Input

One number K in the first line is an amount of tests (K ≤ 100000). Each next line represents separate test, which contains integers a and n (1 ≤ a, n ≤ 32767, n is prime, a and n are relatively prime).

Output

For each input test the program must evaluate all possible values root(a,n) in the range from 1 to n − 1 and output them in increasing order in one separate line using spaces. If there is no square root for current test, the program must print in separate line: ‘No root’.

Sample

input	output
5 4 17 3 7 2 7 14 31 10007 20011	2 15 No root 3 4 13 18 5382 14629

Problem Author: Michael Medvedev

p=2时特判输出1，不然会WA。龟速乘会TLE。

#include<bits/stdc++.h>
#define rg register
#define il inline
#define co const
template<class T>il T read(){
    rg T data=0,w=1;
    rg char ch=getchar();
    while(!isdigit(ch)){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(isdigit(ch))
        data=data*10+ch-'0',ch=getchar();
    return data*w;
}
template<class T>il T read(rg T&x){
    return x=read<T>();
}
typedef long long ll;

ll n,mod;
ll add(ll x,ll y){
	return (x+y)%mod;
}
//ll mul(ll x,ll y){
//	x%=mod,y%=mod;
//	ll re=0;
//	for(;y;y>>=1,x=add(x,x))
//		if(y&1) re=add(re,x);
//	return re;
//}
ll mul(ll x,ll y){
	return x*y%mod;
}
ll pow(ll x,ll k){
	x%=mod,k%=(mod-1);
	ll re=1;
	for(;k;k>>=1,x=mul(x,x))
		if(k&1) re=mul(re,x);
	return re;
}
ll omega;
struct complex{ll a,b;};
complex operator*(co complex&x,co complex&y){
	return (complex){add(mul(x.a,y.a),mul(x.b,mul(y.b,omega))),add(mul(x.a,y.b),mul(x.b,y.a))};
}
complex operator^(complex x,ll k){
	complex re=(complex){1,0};
	for(;k;k>>=1,x=x*x)
		if(k&1) re=re*x;
	return re;
}
ll sqrt(ll n){
	if(mod==2) return 1;
	n%=mod;
	if(!n) return 0;
	ll a=rand()%mod;
	while(pow(add(mul(a,a),mod-n),(mod-1)/2)!=mod-1) a=rand()%mod;
	omega=add(mul(a,a),mod-n);
	return ((complex){a,1}^(mod+1)/2).a;
}
int main(){
//	freopen(".in","r",stdin),freopen(".out","w",stdout);
	int kase=read<int>();
	while(kase--){
		read(n),read(mod);
		if(pow(n,(mod-1)/2)!=1) {puts("No root");continue;}
		ll ans1=sqrt(n),ans2=mod-ans1;
		if(ans1>ans2) std::swap(ans1,ans2);
		if(ans2!=ans1) printf("%lld %lld\n",ans1,ans2);
		else printf("%lld\n",ans1);
	}
	return 0;
}