算法
问题是解方程(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
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;
}