Cipolla 模板

求解二次剩余问题。

定义勒让德符号 $(\frac n p) = n ^ {\frac {p - 1} 2}$ 。有

$$
(\frac n p) = \begin {cases}
0 & p \mid n \\
1 \mod p & \text {n is a quadratic nonresidue modulo p.}\\
-1 \mod p & \text {n is not a quadratic nonresidue modulo p.}\\
\end {cases}
$$

如果 $n$ 是 $p$ 的二次剩余，考虑求解。随机生成一个数 $a$ 满足 $a ^ 2 - n$ 不是 $p$ 的二次剩余。类似复数，记 $\omega = \sqrt {a ^ 2 - n}$ ，可以证明 $
\pm (a + \omega) ^ {\frac {p + 1} 2} $ 即为答案。

查看代码

#include <cstdio>
#include <ctime>
#include <cstdlib>
#include <algorithm>
using namespace std;
template <class Type>
void read (Type &x)
{
    char c;
    bool flag = false;
    while ((c = getchar()) < '0' || c > '9')
        c == '-' && (flag = true);
    x = c - '0';
    while ((c = getchar()) >= '0' && c <= '9')
        x = (x << 1) + (x << 3) + c - '0';
    flag && (x = ~x + 1);
}
template <class Type, class ...rest>
void read (Type &x, rest &...y) { read(x), read(y...); }
template <class Type>
void write (Type x)
{
    x < 0 && (putchar('-'), x = ~x + 1);
    x > 9 && (write(x / 10), 0);
    putchar('0' + x % 10);
}
typedef long long LL;
int w, n, p;
int binpow (int b, int k = p - 2)
{
    int res = 1;
    for (; k; k >>= 1, b = (LL)b * b % p)
        if (k & 1) res = (LL)res * b % p;
    return res;
}
struct Complex
{
    int x, y;
    Complex (int _x = 1, int _y = 0) { x = _x, y = _y; }
    Complex operator * (const Complex &_) const
        { return Complex(((LL)x * _.x + (LL)y * _.y % p * w) % p, ((LL)x * _.y + (LL)y * _.x) % p); }; 
    friend Complex operator ^ (Complex b, int k)
    {
        Complex res;
        for (; k ; k >>= 1, b = b * b)
            if (k & 1) res = res * b;
        return res;
    }
};
int main ()
{
    srand(time(NULL));
    int T; read(T);
    while (T--)
    {
        read(n, p);
        if (n % p == 0) puts("0");
        else if (binpow(n, p - 1 >> 1) ^ 1) puts("-1");
        else
        {
            int a = rand() % p;
            while (binpow(w = ((LL)a * a - n + p) % p, p - 1 >> 1) ^ p - 1)
                a = rand() % p;
            int x1 = (Complex(a, 1) ^ p + 1 >> 1).x, x2 = p - x1;
            if (x1 > x2) swap(x1, x2);
            write(x1), putchar(' '), write(x2), puts("");
        }
    }
    return 0;
}