1. 程式人生 > 實用技巧 >淺談二次剩餘

淺談二次剩餘

淺談二次剩餘

定義:對於正整數\(p,n\),若存在\(x\)使得\(x^2 \equiv n(\bmod\ p)\).則稱\(n\)是模\(p\)的二次剩餘。(在本文中我們只考慮\(p\)為奇質數的情況)

勒讓德括號,尤拉判別準則

下面引入勒讓德括號,它可以簡化討論:

\[\left(\frac{n}{p}\right)=\left\{\begin{array}{ll}1 & (p \nmid n \text { 且 } n \text { 是模 } p \text { 的二次剩餘 }) \\ -1 & (p \nmid n \text { 且 } n \text { 是模 } p \text { 的二次非剩餘 }) \\ 0 & (p \mid n)\end{array}\right. \]

定理(尤拉判別準則): 當\(p\)為奇素數時,\(\left( \frac{n}{p}\right) \equiv n^{\frac{p-1}{2}} (\bmod\ p)\)

證明: \(p|n\)的情況是顯然的,我們考慮\(p \nmid n\)時的情況。顯然\(n^{\frac{p-1}{2}}=\plusmn 1\).若\(n\)是二次剩餘,則\(n^{\frac{p-1}{2}}=(x^2)^{\frac{p-1}{2}}=x^{p-1}=1\). 若\(n^{\frac{p-1}{2}}=1\),把\(n\)用模\(p\)下的原根\(g\)表示為\(n=g^k\),那麼\(g^{k\frac{p-1}{2}}=1\)

.又因為\(g^{p-1}=1\),由原根的性質(見此處定理3.1)得到\(p-1|k\frac{p-1}{2}\),那麼\(k\)必定為偶數,令\(x=g^{\frac{k}{2}}\)即可,說明\(n\)是二次剩餘。也就是說 \(n^{\frac{p-1}{2}}=1\)\(n\)是二次剩餘的充要條件。反之-1的情況成立,上面的\(-1\)在模意義下其實就是\(p-1\)

從這個定理我們已經得到二次剩餘的一種求法,求出原根\(g\)後用BSGS求出\(g^k=n\)的解,那麼\(g^{\frac{k}{2}}\)就是一個解。複雜度\(O(\sqrt{p})\)

Cipolla演算法

上面這個\(O(\sqrt{p})\)

的演算法還不夠優秀,我們用Cipolla演算法來求解\(x^2 \equiv n(\bmod\ p)\)。演算法流程如下:

  1. 不斷隨機\(a\),直到找到一個\(a\)使得\(a^2-n\)是非二次剩餘(用尤拉判別準則判斷)(我們之後會證明期望隨機次數為2).
  2. 定義\(\omega^2=a^2-n\).但是我們知道\(a^2-n\)是非二次剩餘,那麼我們可以類比複數域對無法開根的數\(-1\)定義\(\mathrm{i}^2=-1\),我們把所有數表示為\(p+q\omega\)的形式,運算規則類似複數。那麼\(x=\plusmn (a+\omega)^{\frac{p+1}{2}}\).

引理1:\(x^2 \equiv n(\bmod\ p)\)有且僅有兩組解,且滿足\(x_1+x_2 \equiv n(\bmod\ 0)\)

證明:假設對於\(x_1 \neq x_2\)\(x_1^2 \equiv x_2^2\),移項得\((x_1-x_2)(x_1+x_2)\equiv 0\)。 這很類似實數的開根。又因為每個\(n\)都對應兩個不同的\(x_1,x_2\),一共只有\(p-1\)個數,那麼二次剩餘的數量恰為\(\frac{p-1}{2}\).那麼隨機到非二次剩餘的概率為\(\frac{p-1}{2p} \approx \frac{1}{2}\).期望次數為\(\sum_{i=0}\frac{i}{2^i}=2-\frac{N+2}{2^N}\),趨近於2. 這樣演算法的複雜度就是\(O(\log p)\)

引理2:$ \omega^{p} \equiv-\omega$
證明: \(\omega^{p} \equiv \omega\left(\omega^{2}\right)^{\frac{p-1}{2}} \equiv \omega\left(a^{2}-n\right)^{\frac{p-1}{2}} \equiv-\omega\)

引理3:\((A+B)^{p} \equiv A^{p}+B^{p}\)
證明:根據二項式定理,由於\(p\)是質數,除了\(C_{p}^{0}, C_{p}^{p}\) 外的組合數分子上的階乘沒法消掉\(p\),模\(p\)都為\(0,\) 只有\(C_{p}^{0} A^{0} B^{p}+C_{p}^{p} A^{p} B^{0}\)

根據引理2,3:

\[(a+i)^{p+1} \equiv\left(a^{p}+i^{p}\right)(a+i)\equiv\left(a\cdot a^{p-1}+i^{p}\right)(a+i) \equiv(a-i)(a+i) \equiv a^{2}-i^{2} \equiv n \]

那麼\(x=\plusmn (a+i)^{\frac{p+1}{2}}\)(可以用反證法證明虛部為0,這裡不再贅述)

struct com { //cipolla用的類似複數的東西
    ll real;
    ll imag;//a+b*sqrt(w)
    com() {

    }
    com(ll _real,ll _imag) {
        real=_real;
        imag=_imag;
    }
};
inline com mul(com a,com b,ll w) {//過載乘法
    return com((a.real*b.real%mod+a.imag*b.imag%mod*w%mod)%mod,(a.real*b.imag%mod+a.imag*b.real%mod)%mod);
}
inline com fpow(com x,ll k,ll w) {
    com ans=com(1,0);
    while(k) {
        if(k&1) ans=mul(ans,x,w);
        x=mul(x,x,w);
        k>>=1;
    }
    return ans;
}
inline ll cipolla(ll x) {
    if(fast_pow(x,(mod-1)/2)==mod-1) return -1;//x不存在二次剩餘
    while(1) {
        ll a=((rand()<<15)|rand())%mod;
        ll w=(a*a%mod-x+mod)%mod;
        if(fast_pow(w,(mod-1)/2)==mod-1) { //x-a^2不存在二次剩餘
            return fpow(com(a,1),(mod+1)/2,w).real%mod;
        }
    }
    return -1;
}

完整程式碼:LuoguP5491

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
using namespace std;
typedef long long ll;

int T;
ll n,mod;
inline ll fast_pow(ll x,ll k) {
    ll ans=1;
    while(k) {
        if(k&1) ans=ans*x%mod;
        x=x*x%mod;
        k>>=1;
    }
    return ans;
}
inline ll inv(ll x) {
    return fast_pow(x,mod-2);
}

struct com { //cipolla用的類似複數的東西
    ll real;
    ll imag;//a+b*sqrt(w)
    com() {

    }
    com(ll _real,ll _imag) {
        real=_real;
        imag=_imag;
    }
};
inline com mul(com a,com b,ll w) {
    return com((a.real*b.real%mod+a.imag*b.imag%mod*w%mod)%mod,(a.real*b.imag%mod+a.imag*b.real%mod)%mod);
}
inline com fpow(com x,ll k,ll w) {
    com ans=com(1,0);
    while(k) {
        if(k&1) ans=mul(ans,x,w);
        x=mul(x,x,w);
        k>>=1;
    }
    return ans;
}
inline ll cipolla(ll x) {
    if(fast_pow(x,(mod-1)/2)==mod-1) return -1;//x不存在二次剩餘
    while(1) {
        ll a=((rand()<<15)|rand())%mod;
        ll w=(a*a%mod-x+mod)%mod;
        if(fast_pow(w,(mod-1)/2)==mod-1) { //x-a^2不存在二次剩餘
            return fpow(com(a,1),(mod+1)/2,w).real%mod;
        }
    }
    return -1;
}

int main() {
    scanf("%d",&T);
    while(T--){
        scanf("%lld %lld",&n,&mod);
        if(n==0){
            puts("0");
            continue;
        }
        ll ans=cipolla(n);
        if(ans==-1) puts("Hola!");
        else{
            ans=(ans%mod+mod)%mod;
            printf("%lld %lld\n",min(ans,mod-ans),max(ans,mod-ans));
        }
    }
}