夏休みの思い出（２） [yukicoder No.551] - はまやんはまやんはまやん

解説放送

前提知識

平方剰余

解説

https://yukicoder.me/submissions/193014

２次方程式を平方完成する。

Ax^2+Bx+C=0
A(x^2+BA^-1x)+C=0
A(x^2+BA^-1x+B^2(4A^2)^-1-B^2(4A^2)^-1)+C=0
A(x + B(2A)^-1)^2 - B^2(4A)^-1+C=0
A(x + B(2A)^-1)^2 = B^2(4A)^-1-C
(x + B(2A)^-1)^2 = (B^2(4A)^-1-C)/A
　
ここで
a = x + B(2A)^-1
b = (B^2(4A)^-1-C)/A
と置くと

a^2 = b (mod P)

これは平方剰余と呼ばれる問題になる。
この合同方程式が解ければ、x = a - B(2A)^-1 で答えが求まるので、これを計算することを考える。
　
まず、解が存在するかであるが、これはオイラーの基準という式でルジャンドル関数を作って判定する。
ルジャンドル関数が1を返せば解があり（平方剰余）、-1を返せば解が無い（平方非剰余）
これを使って、解が無いならば、"-1"を吐いて終了
　
あとは、実際に問題を解くが、色々なアルゴリズムがある。
自分の解法ではpekempeyさんのCipolla's algorithm実装をほぼそのまま使っている。
出た答えをxに戻して答えると答えとなる。

legendre(int b, int P) := a^2 = b(mod P)についてのルジャンドル関数
modsqrt(int b, int P) := a^2 = b(mod P)を満たすaを返す

int P, R, Q, A, B, C;
//---------------------------------------------------------------------------------------------------
int f() { return mul(sub(mul(B, B, modinv(mul(4, A))), C),modinv(A)); }
int rev(int x) { return sub(x, mul(B, modinv(mul(2, A)))); }
//---------------------------------------------------------------------------------------------------
void _main() {
    cin >> P >> R >> Q;
    mod = P;

    rep(_, 0, Q) {
        cin >> A >> B >> C;

        int b = f();
        if (legendre(b, P) < 0) {
            printf("-1\n");
            continue;
        }

        int a = modsqrt(b, P);
        if (a == 0) printf("%d\n", rev(0));
        else {
            int aa = P - a;
            a = rev(a), aa = rev(aa);
            if (a > aa) swap(a, aa);
            printf("%d %d\n", a, aa);
        }
    }
}