二次剩余

即模意义下的开根

定义

对于 \(a,p\)，若存在 \(x\) 使得 \(x^2\equiv a\pmod p\)，则 \(a\) 为模 \(p\) 下的二次剩余，否则为非二次剩余

数量

假定 \(p\) 为奇质数

则二次剩余有 \(\frac {p+1}2\) 个，非二次剩余有 \(\frac{p-1}2\) 个

勒让德记号

定义勒让德记号 \(\left(\frac ap\right)\)，当 \(p\mid a\) 时值为 \(0\)，若 \(a\) 为模 \(p\) 的二次剩余则值为 \(1\)，否则为 \(-1\)

根据欧拉判别准则，有 \(\left(\frac ap\right)\equiv a^{\frac{p-1}2}\bmod p\)

Cipolla 算法

求出一个 \(x\) 使得 \(x^2\equiv a\pmod p\)

首先判断是否是二次剩余

然后找到一个 \(w\) 满足 \(w^2-a\) 为模 \(p\) 下的非二次剩余（每次随机选一个，约 \(\frac 12\) 的概率非二次剩余，期望选 \(O(1)\) 次）

定义虚数单位 \(i\) 为满足 \(i^2\equiv w^2-a\pmod p\) 的数（类比负数开根得到的虚数），则 \((w+i)^{\frac{p-1}2}\bmod p\) 即为 \(x^2\equiv a\pmod p\) 的一个解（计算过程类似复数运算，实部和虚部分别对 \(p\) 取模，可证结果的虚部一定为 \(0\)）

模板

namespace modint_nm {

    namespace modint_help {
        template <u32 mod> inline constexpr bool is_quad_residue(modint<mod> a){return (a ^ ((mod - 1) >> 1)) == 1;}
        template <u32 mod>
            struct mod_complex {
                modint<mod> real, imag;
                inline constexpr mod_complex(u32 R = 0, u32 I = 0) noexcept : real(R), imag(I) {}
                inline constexpr mod_complex(modint<mod> R, modint<mod> I = 0) noexcept : real(R), imag(I) {}
                inline constexpr mod_complex& operator = (mod_complex o) noexcept { real = o.real;imag = o.imag;return *this; }
            };
        template <u32 mod>
            inline constexpr mod_complex<mod> mult(mod_complex<mod> a, mod_complex<mod> b, modint<mod> i_2) noexcept
                { return mod_complex<mod>(fmam(a.imag * b.imag, i_2, a.real, b.real), fmam(a.real, b.imag, a.imag, b.real)); }
        template <u32 mod>
            inline constexpr mod_complex<mod> pow(mod_complex<mod> a, u64 b, modint<mod> i_2) noexcept
                { mod_complex<mod> ret(1);for (; b; b >>= 1, a = mult(a, a, i_2))if (b & 1)ret = mult(ret, a, i_2);return ret; }
    }

    template <u32 mod> inline constexpr modint<mod> sqrt(modint<mod> a){
        if (a == modint<mod>(0))return 0;
        if (!modint_help::is_quad_residue(a))return modint<mod>::nan;
        // if constexpr ((mod & 3) == 3)return a ^ ((mod + 1) >> 2);
        // else if constexpr ((mod & 7) == 5){//Atkin
        //     modint<mod> b = (a + a) ^ (mod >> 3);
        //     return a * b * (2 * a * b * b - 1);
        // }
        //Cipolla
        modint<mod> r = 0, i_2 = 0; u32 seed = 1u;
        do {r = (seed = u64(seed) * 16807UL % 2147483647UL);/*minstd_rand0*/} while (modint_help::is_quad_residue(i_2 = r * r - a));
        modint<mod> x = modint_help::pow(modint_help::mod_complex<mod>(r, modint<mod>(1)), (mod + 1) >> 1, i_2).real;
        return x.val < mod - x.val? x : -x;
    }

}

参考

1. \(\text{12.11 math1.pdf \; \;by Tx\_Lcy}\)