((x+y)equiv bpmod p)
((x imes y)equiv cpmod p)
由第一个式子可知：(x+y=b~or~x+y=b+p)
先任选一个代入到第二个式子里得
解二次剩余方程 (q^2equiv apmod p)
因为这个方程去查了很多资料
1. 欧拉准则
对于(x^2equiv apmod p)
证明：
先不考虑a为0的情况。
已知 ((p-x)^2equiv x^2pmod p)， 这是因为(p^2-2xp+x^2equiv x^2 pmod p)
所有有(p-1over 2) 个不同的二次剩余，即(1^2,2^2,cdots,({p-1over 2}) pmod p) （因为前一半和后一半相同了）
又(a^{p-1}equiv 1 pmod p), 可以写为(({a^{p-1over 2}-1})({a^{p-1over 2}} + 1)equiv 0pmod p)
上式中的前后两个因子，必须有一个为0，由(Lagrangs's theorem) 可知 k次多项式最多 k 个解，所以(a^{p-1over 2}-1equiv 0pmod p) 有最多 (p-1over 2)个解。
又(x^2equiv a), 所以(a^{p-1over 2} equiv (x^2)^{p-1over 2} equiv 1pmod p)
所以每一个平方剩余都可以使得第一个因子为0，非0平方剩余最少有(p-1over 2)个，所以这与上面方程的解正好对应起来，也就是说使得(a^{p-1over 2}equiv 1)成立的(a) 都是(p) 的二次剩余。同理可知使得(a^{p-1over 2} equiv -1)成立的(a)都是(p)的非二次剩余。
2. 求解二次剩余
本题比较特殊，(pequiv 3pmod 4) ，那么根据上面推出来的
对于更一般的(p) , 可以参考: http://xbgjxt.swu.edu.cn/jsuns/html/jsuns/2019/1/201901009.htm
综述：
- 根据欧拉准则来判断是否有解
- 有解则求出原题中的x 与 y输出
- 无解则输出-1
标程代码：
int t, p = 1000000007;
//快速幂代码
long long pow(long long x, long long y, long long p) {
}
int main() {
cin >> t;
for (int tt = 0; tt < t; tt++) {
long long b, c;
cin >> b >> c;
long long d = (b * b - 4 * c) % p;
if (d < 0) {
d += p;
}
if (d != 0 && pow(d, (p - 1) / 2, p) != 1) {//如果无解，需要注意d为0的情况
cout << -1 << ' ' << -1 << endl;
} else {
long long r = pow(d, (p + 1) / 4, p);
long long x = (b + r) * pow(2, p - 2, p) % p;
long long y = (b - r) * pow(2, p - 2, p) % p;
if (x < 0)x += p;
if (y < 0)y += p;
if (x > y) swap(x, y);
cout << x << ' ' << y << endl;
}
}
}
自己的代码
const ll mod = 1e9+7;
ll pow_mod(ll a,ll i,ll n){
if(i == 0)return 1 % n;
ll tmp = pow_mod(a, i>> 1,n);
tmp = tmp * tmp % n;
if(i & 1)tmp = tmp * a % n;
return tmp;
}
//红宝书模板代码，求解模p下二次剩余为a的解
ll modsqr(ll a,ll n){
if(a == 0)return 0;
ll b,k,i,x;
if(n == 2)return a % n;
if(pow_mod(a,(n-1)/2,n) == 1){
if(n % 4 == 3)//本题中只会进入下面这个case
x = pow_mod(a,(n+1)/4,n);
else{
for(b = 1;pow_mod(b,(n-1)/2,n) == 1;b++);
i = (n-1)/2;
k = 0;
do{
i/=2;
k/=2;
if((pow_mod(a,i,n) * pow_mod(b,k,n) + 1) % n == 0)
k += (n-1) / 2;
}while(i % 2 == 0);
x = (pow_mod(a,(i+1)/2,n) * pow_mod(b,k/2,n))%n;
}
if(x * 2 > n)x = n-x;
return x;
}
return -1;
}
int main()
{
int T;cin>>T;
while(T--){
ll b,c;
scanf("%lld%lld",&b,&c);
ll q,a;
a = ((b*b - c * 4)%mod + mod) % mod;
q = modsqr(a,mod);
if(q == -1){//无解情况
puts("-1 -1");
continue;
}
//很蠢的分了两种情况....其实乘一个2的逆元即可
ll x = (q + b) / 2;
ll y = mod + b - x;
x = (x % mod + mod) % mod;
y = (y % mod + mod) % mod;
if((x * y) % mod == c){
if(x > y)swap(x,y);
printf("%lld %lld
",x,y);continue;
}
x = (q + b + mod)/2;
y = b - x;
x = (x % mod + mod) % mod;
y = (y % mod + mod) % mod;
if((x * y) % mod == c){
if(x > y)swap(x,y);
printf("%lld %lld
",x,y);continue;
}
puts("-1 -1");
}
return 0;
}
参考资料：