二次剩余

二次剩余

参考： 二次剩余Cipolla算法学习笔记

#include <bits/stdc++.h>
using namespace std;

const int mod=1e9+9;
namespace TwoRemain {
    template <typename A, typename B> inline int add(A x, B y) {if(x + y < 0) return x + y + mod; return x + y >= mod ? x + y - mod : x + y;}
    template <typename A, typename B> inline void add2(A &x, B y) {if(x + y < 0) x = x + y + mod; else x = (x + y >= mod ? x + y - mod : x + y);}
    template <typename A, typename B> inline int mul(A x, B y) {return 1ll * x * y % mod;}
    template <typename A, typename B> inline void mul2(A &x, B y) {x = (1ll * x * y % mod + mod) % mod;}
    int fmul(int a, int p, int Mod = mod) {
        int base = 0;
        while(p) {
            if(p & 1) base = (base + a) % Mod;
            a = (a + a) % Mod; p >>= 1;
        }
        return base;
    }
    int fp(int a, int p, int Mod = mod) {
        int base = 1;
        while(p) {
            if(p & 1) base = fmul(base, a, Mod);
            p >>= 1; a = fmul(a, a, Mod);
        }
        return base;
    }
    int f(int x) {
        return fp(x, (mod - 1) >> 1);
    }
    struct MyComplex {
        int a, b;
         int cn;
        MyComplex operator * (const MyComplex &rhs)  {
            return {
                add(fmul(a, rhs.a), fmul(cn, fmul(b, rhs.b, mod))),
                add(fmul(a, rhs.b), fmul(b, rhs.a)),
                cn
            };
        }
    };
    MyComplex fp(MyComplex a, int p) {
        MyComplex base = {1, 0, a.cn};
        while(p) {
            if(p & 1) base = base * a;
            a = a * a; p >>= 1;
        }
        return base;
    }
    int TwoSqrt(int n) {
        if(f(n) == mod - 1) return -1;
        if(f(n) ==  0) return  0;
        int a = -1, val = -1;
        while(val == -1) {
            a = rand() << 15 | rand();
            val = add(mul(a, a), -n);
            if(f(val) != mod - 1) val = -1;
        }
        return fp({a, 1, val}, (mod + 1) / 2).a;
    }
}
using namespace TwoRemain;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout<<TwoSqrt(5)<<endl;
    return 0;
}