二次同余方程的解

算法

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

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

证明：由费马小定理，(n^{p-1}-1equiv (n^frac{p-1}2-1)(n^frac{p-1}2+1)equiv 0 (mod p) Rightarrow n^frac{p-1}2 equiv pm 1 (mod 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 (mod p))的解。

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

[(a+sqrt{omega})^{p+1}\ equiv(a+sqrtomega)^p(a+sqrtomega)\ equiv(a-sqrtomega)(a+sqrtomega)\ equiv a^2-omegaequiv n (mod 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
",ans1,ans2);
		else printf("%lld
",ans1);
	}
	return 0;
}