timus 1132 Square Root (二次剩余)

求$n$在模$p$($p$是素数)下的开方。
若$w=a^2-n$不是p的二次剩余，则$( a + \sqrt{w} ) ^ { (p+1) / 2 }$是$x^2=n(mod \ p)$的一个解

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <vector>
#include <queue>
#include <map>
#include <set>
#include <algorithm>
#include <ctime>
#include <cstdlib>
#include <functional>
#pragma comment(linker,"/STACK:102400000,102400000")
using namespace std;


#define eps 1e-10
#define N 1600020
#define B 20
#define M 10000020
#define inf 0x3f3f3f3f
#define LL long long
#define pii pair<int, int>
#define MP make_pair
#define fi first
#define se second
#define ls (i << 1)
#define rs (ls | 1)
#define md (ll + rr >> 1)
#define lson ll, md, ls
#define rson md + 1, rr, rs
int w;


int qpow(int x, int k, int p) {
    int ret = 1;
    while(k) {
        if(k & 1) ret = 1LL * ret * x % p;
        k >>= 1;
        x = 1LL * x * x % p;
    }
    return ret;
}
pii mul(pii x, pii y, int p) {

    pii ret;
    ret.fi = 1LL * x.fi * y.fi % p + 1LL * x.se * y.se % p * w % p;
    if(ret.fi >= p) ret.fi -= p;

    ret.se = 1LL * x.se * y.fi % p + 1LL * y.se * x.fi % p;
    if(ret.se >= p) ret.se -= p;
    return ret;
}
pii qpow(pii x, int k, int p) {
    pii ret(1, 0);
    while(k) {
        if(k & 1) ret = mul(ret, x, p);
        k >>= 1;
        x = mul(x, x, p);
    }
    return ret;
}
int legendre(int a, int n) {
    return qpow(a, (n - 1) / 2, n);
}
int main() {
    int cas;
    scanf("%d", &cas);
    while(cas--) {
        int a, n;
        scanf("%d%d", &a, &n);
        a %= n;
        if(n == 2) {
            puts("1");
            continue;
        }

        if(legendre(a, n) == n - 1) {
            puts("No root");
            continue;
        }
        int ans;
        while(1) {
            int x = rand() % (n - 1) + 1;
            w = 1LL * x * x % n - a;
            if(w == 0) {
                ans = x;
                break;
            }
            if(w < 0) w += n;
            if(legendre(w, n) == 1) continue;
            pii t = MP(x, 1);
            t = qpow(t, (n + 1) / 2, n);
            ans = t.fi;
            break;
        }
        int x = ans;
        int y = n - ans;
        if(x > y) swap(x, y);
        printf("%d %d\n", x, y);
    }
    return 0;
}