數論彙總

最前面：\(\LaTeX\)可能須要加載一會，請耐心等待o~

前言

　　數論在\(\text{OI}\)中十分重要。

　　什麼是數論？

\[ 研究整數的理論 ——zzq \]

　　本文包含全部側邊目錄中呈現的內容。

　　下面直奔主題。

整除

　　若\(a\)是\(b\)的因數，或\(b\)是\(a\)的倍數，則\(a\)整除\(b\)，記做\(a\mid b\)。

　　關於整除，有如下幾點：

　　一、若\(a\mid b\)，\(b\mid c\)，則\(a\mid c\)。

　　二、若\(a\mid b\)，\(a\mid c\)，則\(a\mid b+c\)。

　　三、若\(a\mid b\)，對於任意的\(k \in \Z\)，都有\(a\mid kb\)。

同餘

　　（小學生都懂）

　　給定整數\(a\)，\(b\)和模數\(p\)，若\(a \bmod p = b \bmod p\)，則稱\(a\)和\(b\)在模\(p\)意義下同餘，簡寫成\(a\equiv b \pmod p\)。

　　同餘有如下性質：

　　一、\((a \bmod p) \plusmn (b \bmod p) = (a \plusmn b) \bmod p\)

　　二、\((a \bmod p) \times (b \bmod p) = (a \times b) \bmod p\)

　　三、\(\frac{a \bmod p}{b \bmod p} = \frac{a}{b} \bmod p\)（若是你直接算的話，數學老師會揍你，正確姿式放後面）

　　四、\(ax \equiv bx \pmod p \Leftrightarrow a \equiv b \pmod { \frac{ p }{ \gcd(x,p) } }\)

　　\(Proof\)：

\[\because ax \equiv bx \pmod p\]

\[\therefore (a-b)x \equiv 0 \pmod p\]

\[\therefore p\mid (a-b)x\]

\[\therefore \frac{p}{\gcd(x,p)}\mid \frac{x}{\gcd(x,p)}(a-b)\]

\[\because \gcd(\frac{p}{\gcd(x,p)},\frac{x}{\gcd(x,p)})=1，也就是互質\]

\[\therefore \frac{p}{\gcd(x,p)} \mid (a-b)\]

\[\therefore a-b \equiv 0 \pmod{ \frac{p}{\gcd(x,p)} }\]

\[因而a\equiv b \pmod{ \frac{p}{\gcd(x,p)} }得證\]

裴蜀定理

　　先定義線性組合：形如\(ax+by\)的式子，其中\(a,b,x,y\in \Z\)。

　　裴蜀定理就是：若是\(a\)、\(b\)是不爲\(0\)的任意整數，則\(\gcd(a,b)\)是線性組合\(ax+by\)的最小正整數。

　　\(Proof\):

\[由\gcd(a,b)\mid a和\gcd(a,b)\mid b得\]

\[\gcd(a,b)\mid ax+by\]

\[設d=\gcd(a,b)，ax+by的最小正值爲s，由d\mid ax+by恆成立獲得d\mid s\]

\[令q=\lfloor \frac{a}{s} \rfloor，r=a \bmod s=a-q(ax+by)=a(1-qx)+b(-qy)，發現r也爲a,b的線性組合\]

\[又\because 0 \leq r < s\]

\[\therefore r = 0（最小正整數爲s，而r<s，說明r=0)\]

\[\therefore q=\frac{a}{s}，即s\mid a，相似的可獲得s\mid b\]

\[\therefore s\mid ax+by，而d\mid ax+by\]

\[\therefore s = d\]

\[ax+by的最小正值爲\gcd(a,b)得證\]

　　這必定理十分經常使用。

展轉相除法 && 擴展歐幾里得算法（Ex-GCD）

　　推論：\(\gcd(a,b)=\gcd(b, a \bmod b)\)。

　　\(Proof\):

\[由裴蜀定理得ax+by的最小正整數爲\gcd(a,b) \tag{*}\]

\[同時獲得bx+(a \bmod b)y的最小正整數爲\gcd(b, a\bmod b) \tag{**}\]

\[而bx+(a \bmod b)y = bx + (a - \lfloor \frac{a}{b} \rfloor \times b)y=bx+ay-\lfloor \frac{a}{b}\rfloor by = ay+b(x-\lfloor \frac{a}{b}\rfloor y)\]

\[該式又是關於a,b的線性組合，最小正值爲\gcd(a,b)，與(*)，(**)同樣。\]

　　同時，這個跟\(Ex-GCD\)的推導大同小異：爲了求解\(ax+by=\gcd(a,b)\)方程，利用遞歸求解，即新的\(x'=y,y'=x-\lfloor \frac{a}{b} \rfloor y\)，當\(b=0\)時是邊界，顯然\(x=1,y=0\)，此時知足\(ax+by=\gcd(a,b)=a\)。

　　相似的，也能夠獲得：

　　推論\(1\):\(\forall a,b\in \Z : \gcd(a,b)=\gcd(a-b,b)\)；

　　推論\(2\):\(\forall a,b爲偶數:\gcd(a,b)=2 \times \gcd(a/2,b/2)\)；

　　推論\(3\):\(\forall a爲偶數,b爲奇數:\gcd(a,b)=\gcd(a/2,b)\)。以上式子能夠進一步推廣成：若\(a=2^i \times a'\)，\(b=2^j \times b'\)，則\(\gcd(a,b)=2^{\min(i,j)} \times \gcd(a',b')\)。

　　證實略

羣

　　\((G,·)\)表示\(G\)對定義在集合\(G\)上的運算\(·\)構成羣，它須要知足：

　　（1）封閉性：若\(a,b \in G\)，則存在惟一肯定的\(c \in G\)，使得\(a·b=c\)。

　　（2）結合律：對於\(a,b,c \in G : (a·b)·c = a·(b·c)\)

　　（3）存在單位元：存在\(e \in G\)，使得\(\forall a \in G : a·e=e·a\)，\(e\)被稱做單位元。

　　（4）存在逆元：\(\forall a \in G, \exist b \in G : a·b=b·a=e\)，\(b\)又寫做\(a^{-1}\)。

逆元

定義

　　根據羣中的理論，在模\(p\)意義下，對於一個數\(a\)，若是存在一個數\(b\)，使得\(a \times b \equiv 1 \pmod p\)，則\(b\)是\(a\)的逆元，記做\(a^{-1}\)。同時，\(a\)和\(a^{-1}\)互爲逆元。

　　逆元存在的條件：\((a,p)=1\)。

　　\(Proof\):

\[上式可當作ax+py=1，只不過y應爲負數\]

\[由裴蜀定理得：(a,p)\mid 1，因此(a,p)=1，\text{Q.E.D}\]

　　上面的\(\frac{a}{b}\)能夠寫成\(a\times b^{-1}\)，假設逆元存在，因而前面式\(3\)就能夠成了。

性質

　　一、存在惟一性

　　模\(p\)意義下，對於\(a\)而言，若是其存在逆元，則逆元只有一個\(a'=a^{-1}\)。

　　證實很顯然：假設存在另外一個逆元\(a''\times a \equiv 1\)，由\(a'\ \times a \equiv 1 \pmod p\)得：\((a''-a') \times a \equiv 0 \pmod p\)，有\(a \neq 0\)，因此有\(a''-a'\equiv 0 \pmod p\)，因而\(a' \equiv a'' \pmod p\)。

　　二、徹底積性函數

　　也就是\((ab)^{-1} \equiv a^{-1} \times b^{-1} \pmod p\)。

　　證實也很簡單：考慮有\(a \times a^{-1} \equiv 1 \pmod p\)和\(b \times b^{-1} \equiv 1 \pmod p\)，乘到一塊兒就是\((ab)(a^{-1}\times b^{-1}) \equiv 1 \pmod p\)，因此\((ab)^{-1} \equiv a^{-1}\times b^{-1} \pmod p\)。

　　如何求解逆元呢？

求解1：暴力

　　一個一個試唄（很少說）

　　複雜度：\(\text{O}(p)\)。

求解2：費馬小定理

　　定理：當\(p\)爲素數時，\(a^{p-1}\equiv 1 \pmod p\)

　　\(Proof\)：

\[構造一個縮系A=\{x\mid (x,p)=1\}，則\mid A\mid =\varphi(p)=p-1，\varphi(p)表示與p互質的數個數\]

\[對於a \neq 0構造集合B=\{ax\mid x \in A\}，易證在元素模p意義下：A=B\]

\[先證\forall x \in B，x' \in A：x由惟一的x'·a獲得\]

\[與證惟一性同樣，假設有x''·a\equiv x，則有(x'-x'')a \equiv 0，有a\neq 0，故有x\equiv x''，因此惟一\]

\[\therefore 在模p意義下，\prod_{x \in A}x \equiv \prod_{y \in B}y \equiv \prod_{x \in A}ax \equiv a^{p-1}\prod_{x \in A}x \pmod p\]

\[即證a^{p-1}\equiv 1 \pmod p\]

　　因此要求\(a^{-1}\)，只需求\(a^{p-2}\)便可。

　　複雜度：\(\text{O}(\log p)\)。

　　這裏能夠推廣成\(a^{\varphi(p)}\equiv 1 \pmod p\)。（縮系大小爲\(\varphi(p)\)）

求解3：擴展歐幾里得

　　寫個式子你就明白了：\(ax+py=1\)，\(x\)就是\(a^{-1}\)。

　　複雜度：\(\text{O}(\log p)\)。

求解4：歐拉篩

　　由徹底積性的性質獲得的方法。這個方法能夠批量求解，可是對於\(p\)是素數的狀況下仍然須要前面的方法，複雜度\(\text{O}(n+素數\times \log n)\)

求解5：線性遞推

　　假如咱們要求\(a\)在模\(p\)意義下的逆元，構造式子\(p=qa+r\)，發現\(q=\lfloor \frac{p}{a} \rfloor\)，\(r=p \bmod a\)。

　　而後咱們有：

\[r\times r^{-1} \equiv 1 \pmod p\]

　　用\(r=p-qa\)替換得：

\[(p-qa)\times r^{-1} \equiv 1 \pmod p\]

　　而後：

\[-qa\times r^{-1} \equiv 1 \pmod p\]

　　因此咱們就有了：

\[a^{-1} \equiv -q \times r^{-1} \pmod p\]

　　用最上面的東西替換就獲得了：

\[a^{-1} \equiv -\lfloor \frac{p}{a} \rfloor \times (p \bmod a)^{-1} \pmod p\]

　　因而咱們就能夠線性求解了。

n = read(), p = read(); // 咱們要求1~n在模p意義下的逆元

inv[1] = 1; // 首先1的逆元爲1
rep(i, 2, n)
    inv[i] = (ll)(p-p/i) * inv[p % i] % p; // 日後線性推出來

　　複雜度：\(\text{O}(p)\)。

　　可是這裏有缺陷：就是逆元必須從前日後推，若是隨機性地從中挑些數出來詢問，而\(p\)又很大，就會\(\text{T}\)。因此根據題目須要進行選擇（可能還須要巧妙的方法）。

解法6：另外一種方法

　　前提是獲得了模數\(p\)和要詢問的數組\(A\)，要求\(A\)的逆元。

　　令\(S=\prod A_i\)，\(A_i\)的逆元就是\(\frac{S/A_i}{S}\)。咱們不但願用費馬小定理求屢次逆元，這樣複雜度很大。咱們只要求出\(S^{-1}\)便可，這一步的複雜度爲\(\text{O}(\log p)\)。

　　對於\(S/A_i\)，咱們維護一個前綴\(P_i=\prod_{x<i}A_x\)，維護一個後綴\(Q_i=\prod_{x>i}A_x\)，則\(S/A_i = P_i\times Q_i\)（正好少乘了一個\(A_i\)），因而問題得解。這一步的複雜度只與查詢數多少有關，爲\(\text{O}(n)\)。

　　總複雜度：\(\text{O}(n+\log p)\)。正好彌補瞭解法\(5\)的缺陷。

解同餘方程

　　即求出\(ax\equiv b \pmod p\)的解。

　　首先根據裴蜀定理必需要有\((a,p)\mid b\)，要否則無解。若是\((a,p)=1\)，顯然\(a\)存在逆元，則\(x \equiv ba^{-1} \pmod p\)；

　　若是\((a,p)\neq 1\)，令\(d=(a,p)\)，

\[\because d\mid a, d\mid p, d\mid b\]

\[\therefore 由同餘性質獲得:\frac{a}{d} x \equiv \frac{b}{d}\pmod{ \frac{p}{d} }\]

　　令\(a'=\frac{a}{d}\)，\(b'=\frac{b}{d}\)，\(p'=\frac{p}{d}\)，此時有\((a', p')=1\)，\(a'\)存在逆元，因此\(x \equiv b'a'^{-1} \pmod {p'}\)。若是逆元不能用費馬定理求解（好比說\(p'\)不是質數），就用\(\text{Ex-GCD}\)求解。

歐拉函數

　　前面提過：與\(n\)互質的數的個數，記做\(\varphi(n)\)。

　　這個函數的求解：\(\varphi(n)=n\prod_{p_i}(1-\frac{1}{p_i})\)，其中\(p_i\)是\(n\)的素因子。

　　\(Proof\)：

\[咱們用容斥\]

\[\varphi(n)=n-\frac{n}{p_1}-\frac{n}{p_2}-\dotsb+\frac{n}{p_1p_2}+\frac{n}{p_1p_3}-\dotsb\frac{n}{p_1p_2p_3}-\dotsb\]

\[（這裏全部的除法均能整除）\]

\[提出其中一個素因子：\]

\[ \begin{array}{ll} \varphi(n)&=n[(1-\frac{1}{p_1}-\frac{1}{p_2}\dotsb-\frac{1}{p_{s-1}}+\frac{1}{p_1p_2}+\frac{1}{p_1p_3}+\dotsb)-\frac{1}{p_s}(1-\frac{1}{p_1}-\frac{1}{p_2}\dotsb-\frac{1}{p_{s-1}}+\frac{1}{p_1p_2}+\frac{1}{p_1p_3}+\dotsb)]\\ \\ &=n(1-\frac{1}{p_s})(1-\frac{1}{p_1}-\frac{1}{p_2}\dotsb-\frac{1}{p_{s-1}}+\frac{1}{p_1p_2}+\frac{1}{p_1p_3}+\dotsb)\\ \\ &=\dotsb \\ \\ &=n\prod_{p_i}(1-\frac{1}{p_i}) \end{array} \]

　　解釋下：就是每次在容斥的式子裏面提出一個素因子，剩下的因式與沒有該素因子的式子同樣，能夠合併，因而命題得證。

　　由此，咱們看出\(\varphi(n)\)是積性函數，能夠用線性篩求解。

const int N = 10000000;

int p[N], check[N], phi[N];
void init() {
    check[1] = 1;
    phi[1] = 1;
    rep(i, 2, N) {
        if (!check[i]) phi[p[++p[0]] = i] = i-1;
        rep(j, 1, p[0]) {
            if (i*p[j] > N) break;
            check[i*p[j]] = 1;
            if (!(i % p[j])) { phi[i*p[j]] = phi[i] * p[j]; break; } // 質因數p[j]在i中已經出現，互質個數增長p[j]倍
            else phi[i*p[j]] = phi[i] * (p[j] - 1); // 不然乘上p[j]-1
        }
    }
}

　　 當\(n>2\)時，互質的數是成對出現的，若\((n,m)=1\)，則\((n,n-m)=1\)，顯然\(m\neq n-m\)，這一有用的性質能夠證實：全部與\(n\)互質的數的和爲\(\frac{n\varphi(n)}{2}\)。\(n\leq 2\)時特殊代入一下便可。

歐拉定理

　　上面提過：當\((a,n)=1\)時，有\(a^{\varphi(n)}\equiv 1 \pmod n\)，證實略了。

　　還有拓展歐拉定理：任意狀況下，\(a^{b+\varphi(n)} \equiv a^{b \bmod \varphi(n) + \varphi(n)} \pmod n\)。咱們採用分解化歸再合併的方法證實該命題。

　　\(Proof\)：

\[對於(a,n)\neq 1，取a的質因子p，將n分解使得n=s\times p^r，且(a,s)=1\]

\[由積性函數性質得:\varphi(s)\mid \varphi(n)\]

\[\therefore 有p^{\varphi(n)} \equiv 1 \pmod s\]

\[\forall k \in \Z:p^{k\varphi(n)} \equiv 1 \pmod s\]

\[由同餘的性質4:p^{k\varphi(n)} \times p^r \equiv p^r \pmod {s\times p^r}\]

\[\therefore p^{k\varphi(n)+r} \equiv p^r \pmod n\]

\[\therefore \forall c \in \Z_+:p^{k\varphi(n)+r+c} \equiv p^{r+c} \pmod n\]

\[也就是說:\forall c\in \Z, c\geq r:p^c \equiv p^{k\varphi(n)+c} \pmod n\]

\[對於a=p_1^{\alpha_1}p_2^{\alpha_2}\dotsb，先考慮其中任意一個因子p^k，p是素數，假設r\leq \varphi(n) \leq c，有:\]

\[(p^k)^c \equiv p^{kc} \equiv p^{kt\varphi(n)+kc} \equiv (p^k)^{t\varphi(n)+c} \pmod n\]

\[綜上獲得:(p^k)^c \equiv (p^k)^{c \bmod{\varphi(n)} + \varphi(n)}，加上一個\varphi(n)是爲了避免讓指數<c。\]

\[將全部質因子乘起來:(\prod_{p_i}p_i^{\alpha_i})^c \equiv \prod_{p_i}(p_i^{\alpha_i})^c \equiv \prod_{p_i}(p_i^{\alpha_i})^{c \bmod \varphi(n) + \varphi(n)} \equiv (\prod_{p_i}p_i^{\alpha_i})^{c \bmod \varphi(n) + \varphi(n)} \pmod n\]

\[用b+\varphi(n)替換c（由於有c\geq \varphi(n)），即證a^{b+\varphi(n)} \equiv a^{b \bmod \varphi(n) + \varphi(n)}\]

威爾遜定理

　　若\(p\)是質數，恆有\((p-1)! \equiv -1 \pmod p\)成立。

　　\(Proof\)：

\[發現a和a^{-1}互爲逆元，獲得逆元的成對出現\]

\[但若是a\equiv a^{-1}怎麼辦，咱們能夠求出來\]

\[也就是說aa^{-1}\equiv a^2 \equiv 1 \pmod p\]

\[咱們變個形:a^2-1 \equiv 0 \pmod p\]

\[也就是:p \mid a^2-1 \Rightarrow p \mid (a-1)(a+1)\]

\[\therefore a-1 \equiv 0 \pmod p 或者 a+1 \equiv 0 \pmod p，不可能存在第三者由於p是質數\]

\[獲得a \equiv 1 \pmod p或a\equiv -1 \pmod p\]

\[對於p=2，兩式等價，(p-1)!\equiv 1 \equiv p-1 \pmod p，成立\]

\[對於p\neq 2，有兩解爲\pm 1（1或p-1）\]

\[對於其它數必定能與逆元配對，乘積爲1\]

\[\therefore (p-1)! \equiv 1^{\frac{p-3}{2}} \times 1 \times -1 \equiv -1 \pmod p，這裏\frac{p-3}{2}指逆元配對數\]

\[\text{Q.E.D}\]

中國剩餘定理

　　用來求解：

\[ \begin{cases} x \equiv a_1 \pmod {p_1} \\ x \equiv a_2 \pmod {p_2} \\ \dotsb \\ x \equiv a_n \pmod {p_n} \\ \end{cases} \]

　　的最小正整數解（或者通解）。

解法一：構造

　　讓\(M=\prod p_i\)，\(m_i = \frac{M}{p_i}\)，令\(m_i^{-1}\)爲\(m_i\)在模\(p_i\)意義下的逆元，則有\(m_im_i^{-1} \equiv 1 \pmod {p_i}\)，即\(a_im_im_i^{-1} \equiv a_i \pmod {p_i}\)，將\(a_im_im_i^{-1}\)記做\(L_i\)。這裏能夠看出\(m_i^{-1}\)的存在必須知足：\((m_i,p_i)=1\)，也就是\(p_i\)兩兩互質才行。

　　構造\(x=\sum L_i\)，\(x \bmod M\)就是同餘方程組的最小正整數解。

　　\(\forall i,j:L_j \bmod p_i \equiv \begin{cases} 0,j \neq i（m_i中有因子p_i）\\ a_i,j=i（由上面定義可知）\end{cases}\)

N = read();

M = 1;
rep(i, 1, N) {
    p[i] = read(), a[i] = read(); // p爲模數，a爲餘數
    M *= p[i]; 
}

ans = 0;
rep(i, 1, N) {
    ll m = M / p[i], invm, y, d;
    exgcd(m, p[i], invm, y, d); // 求出mi和mi^-1
    ans = (ans + (ll)a[i] * m % M * (invm+p[i])) % M; // 加起來
}

printf("%lld", ans);

　　構造的侷限性也很顯然：模數要兩兩互質才行，不然沒法構造。

　　複雜度：\(\text{O}(n\log (\prod p_i))\)，由於要求逆。

解法二：合併

　　構造能夠應用於各類情形（不互質也ok）。思想就是將兩個同餘式合併成一個同餘式，合併\(n-1\)次只剩下一個同餘式，就是咱們要的答案。

　　令\(k_i=\text{LCM}_{x\leq i}\ b_x\)，假設咱們合併了前面的\(i-1\)個方程，通解爲\(x\)，也就是說這個\(x\)代入前\(i-1\)個方程恆成立，爲了讓第\(i\)個方程成立，咱們有：

\[x+t\times k_{i-1}\equiv a_i \pmod {p_i}\]

　　記\(x'=x+t\times k_{i-1}\)，發現\(x'\)對於前\(i\)個式子恆成立，並且\(k_{i-1}\)與前\(i-1\)個模數取餘爲\(0\)。解出來\(t\)，而後一直算下去。

　　可是不互質時可能會無解，當且僅當\((k_{i-1}, p_i)\nmid a_i-x\)。

int N;
ll a[maxn], p[maxn];

void inc(ll &a, ll b, ll p) { // 同餘加
    a += b;
    if (a >= p) a -= p;
}

ll dec(ll v, ll p) { return v < 0 ? v+p : v; } // 同餘減後的處理

ll Mult(ll a, ll b, ll p) { // 快(gui)速乘
    ll res = 0;
    for (ll k = a; b; inc(k, k, p), b >>= 1)
        if (b & 1) inc(res, k, p);
    return res;
}

void exgcd(ll a, ll b, ll &x, ll &y, ll &d) {
    b ? (exgcd(b, a % b, y, x, d), y -= a/b*x) : (x = 1, y = 0, d = a); // 求ax+by=gcd(a,b)的一組解
}

void excrt() {
    ll x = a[1], k = p[1]; // 初始值
    rep(i, 2, N) {
        ll t, y, d, c, _k;
        exgcd(k, p[i], t, y, d); // 求解方程t*k+p_i*y=(k,p_i)的解
        if ((c = dec(a[i]-x%p[i], p[i])) % d) return; // 若是(k,p_i)不整除c，同餘方程無解
        _k = k, k = k/d*p[i]; // 用_k記錄原來的k
        t = Mult((t+p[i])%p[i], c/d, p[i]); // 更新t的值（按比例放大）
        inc(x, Mult(t, _k, k), k); // 用x'代替x
    }
    printf("%lld", x);
}

原根

　　定義原根以前，咱們先定義個東西。

　　整數的階：給定\(a\)，求使得\(a^x \equiv 1 \pmod p\)成立的最小的正整數\(x\)，這個\(x\)成爲\(a\)模\(p\)的階，同時要知足\((a,p)=1\)，記做\(\text{ord}_pa\)。由歐拉定理易知它確定存在，且能推出\(\text{ord}_pa \mid \varphi(p)\)。

　　階的性質：

　　一、\(\text{ord}_pa \mid \varphi(p)\)。

　　二、\(a\)，\(a^2\)，\(a^3\dotsb a^{\text{ord}_pa}\)在模\(p\)意義下互不相等。顯然若是相等會違背階的定義。

　　三、\(a^r \equiv a^{r'} \pmod p \Leftrightarrow r \equiv r' \pmod{\text{ord}_pa}\)。一樣很顯然。

　　四、\(\text{ord}_ma^k=\frac{\text{ord}_ma}{(\text{ord}_ma,k)}\)。證實用一下裴蜀定理便可。

　　五、\(\text{ord}_m(ab) = \text{ord}_ma \times \text{ord}_mb \Leftrightarrow (\text{ord}_ma, \text{ord}_mb) = 1\)。用整除和階的性質證一證便可。

　　如何求\(\text{ord}_pa\)？枚舉\(\varphi(p)\)的質因子試除判斷便可。

// 事先咱們求出了phi[N]，pri表示分解phi[N]獲得的質因數
int ord(int n) {
    int ret = phi[N];
    for (register int i = 0; i < (int)pri.size(); ++i) {
        if (qpow(n, ret/pri[i], N) == 1) ret /= pri[i]; // 試除成功
    }
    return ret;
}

　　原根：若\(\text{ord}_pa=\varphi(p)\)時，則稱\(a\)是\(p\)的原根。

　　有原根的數爲\(1\)，\(2\)，\(4\)，\(t\)，\(2t\)，\(t^k\)，\(2t^k\)，其中\(t\)爲奇素數，\(k \in \Z_+\)。

　　由此咱們能夠看出：原根\(g\)的\(t\)次方在模\(p\)意義下會遍歷全部\(\varphi(p)\)個與\(p\)互質的數。

　　若是在模數\(p\)下有原根，則原根的個數爲\(\varphi(\varphi(p))\)。

　　\(Proof\)：

\[咱們找一個任意的原根g，這個g的任意次方構成p的縮系，大小爲\mid \varphi(p) \mid\]

\[且根據階的定義，會獲得\{g,g^2,g^3,\dotsb g^{\varphi(p)}\}，兩兩不一樣\]

\[由於有g^{\varphi(p)}\equiv 1，因此g^x \equiv g^{x \bmod \varphi(p)}，構成一個環\]

\[咱們考慮其中哪些元素是原根\]

\[設原根爲g'，易知其指數\leq \varphi(p)時，會取遍縮系中元素\]

\[由於g'在縮系中，因此g'能夠表示成g^k，其中1\leq k \leq \varphi(p)\]

\[g'是原根當且僅當g^k,g^{2k},\dotsb,g^{k\varphi(p)}兩兩不一樣\]

\[等價於k,2k,\dotsb,k\varphi(p)在模\varphi(p)意義下兩兩不一樣\]

\[命題的成立必須知足(k,\varphi(p))=1\]

\[能夠理解成從環上起點開始每次走k步走回到起點的最小次數爲\varphi(p)\]

\[也就是解kx \equiv 0 \pmod{\varphi(p)}，獲得最小正數解x=\frac{\varphi(p)}{(\varphi(p),k)}\]

\[而使得x=\varphi(p)，當且僅當(\varphi(p),k)=1\]

\[這樣的k有\varphi(\varphi(p))個，即證原根個數只有\varphi(\varphi(p))個\]

　　經過這個還能夠證在縮系中，若是\(d \mid \varphi(p)\)，則恰有\(\varphi(d)\)個數的階爲\(d\)，這裏不證了（同上）。

　　判斷\(a\)是否爲原根？一樣試除。因爲\(\text{ord}_pa \mid \varphi(p)\)，咱們枚舉\(\varphi(p)\)的質因子\(p_i\)試除而後看\(a^{\frac{\text{ord}_pa}{p_i}}\)是否同餘爲\(1\)，若是對於全部質因子除完後都不成立，則\(a\)是\(p\)的原根。證實略。

　　如何求原根？原根比較多，且相對分佈均勻，因此直接暴力找。

BSGS

　　\(\text{BSGS}\)（\(\text{Baby Step Giant Step}\)）是用來解決\(a^x \equiv b \pmod p\)的問題。

(a,p)=1

　　令\(m=\lceil \sqrt{p} \rceil\)，則\(x\)能夠表示成\(mq+r\)，其中\(0 \leq q \leq m\)，\(0 \leq r < m\)。

　　那麼\(a^x\equiv b \pmod p\)能夠寫成\(a^{mq+r} \equiv b \pmod p\)。

　　移項得：\(a^{mq} \equiv ba^{-r} \pmod p\)。

　　如今問題變成了：找出一組\(q\)，\(r\)，使得上式成立，咱們枚舉\(q\)，將左邊獲得的式子的全部種取值哈希起來，而後再枚舉\(r\)，判斷右邊的取值是否在左邊的取值中出現，若是出現就找到了解。

　　枚舉部分的複雜度\(\text{O}(\sqrt p)\)，求逆元的複雜度爲\(\text{O}(\log p)\)，因此總複雜度爲\(\text{O}(\sqrt p + \log p)\)。

　　若是你不肯意求逆也能夠，只要設\(x=mq-r\)，而後\(0 < q \leq m\)，\(0 \leq r < m\)，而後移項過程當中\(-r\)變成\(+r\)，因而就不用逆元了。

// 前面定義了Hash類，支持不重複插入和查詢
int BSGS(int a, int b) { // 解a^x = b
    if (a == 0) return b ? -1 : 1;
    Hash::clear(); // 清空
    int m = (int)ceil(sqrt(P+0.5)); // 儘量使得m大
    int am = qpow(a, m), t = 1; // 求a^m
    rep(i, 0, m) {
        Hash::insert(t, i); // 一個一個插入
        t = (ll)t * am % P;
    }
    t = b;
    int inva = inv(a), ans = -1, q; // 求出a的逆元，標記答案
    rep(r, 0, m-1) {
        if (~(q = Hash::exist(t))) { // 哈希表存在t
            int x = m*q + r; // 求解指數
            if (ans == -1 || x < ans) ans = x; // 更新答案（可能不是最小的指數）
        }
        t = (ll)t * inva % P;
    }
    return ans; // 若是ans=-1代表無解；反之有解且爲最小解。
}

(a,p)≠1

　　這種狀況下，\(a\)的逆元不存在，因此咱們要約分兩邊使得\(a\)的逆元出現。

　　對於\(a^x \equiv b \pmod p\)，提出一個\(d=(a,p)\)並約分獲得\(a^{x-1}\frac{a}{d} \equiv \frac{b}{d} \pmod {\frac{p}{d}}\)。固然要有\(d \mid b\)，不然無解。

　　繼續，若是\(d'=(a,\frac{p}{d}) \neq 1\)，再拿出一個\(a\)提取因子，獲得\(a^{x-2}\frac{a}{d}\frac{a}{d'} \equiv \frac{b}{dd'} \pmod {\frac{p}{dd'}}\)

　　以此類推，直到其互質，此時有\(a^{x-t}\frac{a^t}{\prod d} \equiv \frac{b}{\prod d} \pmod {\frac{p}{\prod d}}\)

　　顯然此時\(\frac{a^t}{\prod d}\)存在逆元，因而能夠移過去獲得\(a^{x-t} \equiv \frac{b}{\prod d}\frac{a^{-t}}{\prod d} \pmod {\frac{p}{\prod d}}\)，而後就成了上面的問題。

　　這裏約分時因爲約分和求公約數，多了一個\(\text{O}(\log^2 p)\)的複雜度，其餘不變。

int exBSGS(int a, int b, int p) {
    int t = 0, d = gcd(a, p), l = 1; // 如上所述，l表示約分後左邊的因子
    while (d > 1) {
        if (b % d) return -1; // 若是不能整除說明必定無解
        t++; p /= d, b /= d; // 約分，並記錄指數
        l = (ll)l * (a / d) % p; //
        d = gcd(a, p);
    }
    if (!l) return t; // 若是左邊的因子爲0，顯然有a=0，b=0（同餘意義下），而後解就是t，這個特判掉
    a %= p, b = (ll)b * inv(l, p) % p; // 移項
    int ans = BSGS(a, b, p); // 轉化成了普通的bsgs問題
    return ~ans ? ans+t : -1; // 若是無解返回-1，有解返回ans+t
}

剩餘問題

　　解決形如\(x^k \equiv a \pmod p\)的問題。暴力找彷佛代價過高了。咱們先從簡單的入手。

二次剩餘

　　若\(x^2 \equiv a \pmod p\)有解，則稱\(a\)是\(p\)的二次剩餘。不然稱其爲非二次剩餘。

　　定義\(\text{Legendre}\)符號：

\[ \left( \frac{a}{p} \right)=\begin{cases} 1 & a是p的二次剩餘 \\ -1 & a不是p的二次剩餘 \\ 0 & p \mid a \end{cases} \]

　　（長相奇怪，其實只是記號）

　　先討論\(p\)是奇素數的狀況。此時有：

\[\left( \frac{a}{p} \right) \equiv a^\frac{p-1}{2} \pmod p\]

　　\(Proof\)：

\[當p \mid a時命題顯然成立\]

\[當p \nmid a時，若是a是p的二次剩餘，則有x^2\equiv a \pmod p\]

\[\therefore (x^2)^{\frac{p-1}{2}} \equiv a^{\frac{p-1}{2}} \pmod p\]

\[\therefore x^{p-1} \equiv a^{\frac{p-1}{2}} \pmod p\]

\[由費馬小定理易知左邊同餘1，因此右邊同餘1\]

\[若是a不是p的二次剩餘，則不存在x使得x^2 \equiv a \pmod p成立\]

\[\because p是奇質數\]

\[\therefore 必定能找到x \neq y，使得xy \equiv a \pmod p\]

\[而共有\varphi(p)=p-1個數，其能夠兩兩配對使得上式成立（由逆元的惟一性獲得）\]

\[\therefore (p-1)! \equiv a^{\frac{p-1}{2}} \pmod p\]

\[又由威爾遜定理得:(p-1)! \equiv -1 \pmod p\]

\[\therefore a^{\frac{p-1}{2}} \equiv -1 \pmod p\]

　　咱們證實了全部的二次剩餘獲得的結果都是\(1\)，全部的非二次剩餘獲得的結果都是\(0\)，而又不存在第三者，因此命題能成立。

　　而後有\(p\)是奇質數下，二次剩餘和非二次剩餘的數量同樣多。這個很顯然，讓\(1\)到\(p-1\)的平方在模\(p\)意義下構成集合，而咱們有：

\[x^2 \equiv (p-x)^2 \pmod p\]

　　這個東西展開下就證完了，因此咱們只要考慮\(1\)到\(\frac{p-1}{2}\)之間的平方便可。

　　如有\(x^2 \equiv y^2 \pmod p\)，\(x \neq y\)，則\(x^2-y^2 \equiv 0 \pmod p \Rightarrow p \mid (x-y)(x+y)\)，因爲\(p\)是質數，因此有\(p \mid x-y\)或\(p \mid x+y\)而獲得\(x \equiv y\)或\(-y\)，第一種狀況不存在，第二種狀況發現兩數必定不在同一個部分（即在\(\frac{p-1}{2}\)兩邊），不合假設。因此獲得\(1\)到\(\frac{p-1}{2}\)之間，兩兩的平方不一樣，而這樣獲得的二次剩餘只有\(\frac{p-1}{2}\)個，非二次剩餘也只有\(\frac{p-1}{2}\)個。因此得證。

Cipolla

　　該算法是用來解決二次剩餘問題的，思想基於構造。注意下文解決的是\(x^2 \equiv n \pmod p\)，右邊的是\(n\)而不是\(a\)了。

　　若是\(n\)是二次非剩餘，顯然無解；

　　若是\(n\)是二次剩餘，咱們隨機一個\(a\)，使得\(\left( \frac{a^2-n}{2} \right) = -1\)。由非二次剩餘和二次剩餘數量同樣多，因此有\(\frac{1}{2}\)的機率取到，而指望下\(2\)次即能獲得。證實不重要，放在後面。定義\(\omega = \sqrt{a^2-n}\)，事實上\(\omega\)在整數域中沒有數存在，因此相似於定義虛數單位\(i\)同樣，存數也先存成\(a+b\omega\)這種形式，以後會消掉這個東西。

　　咱們要\(x^2 \equiv a \pmod p\)，先要知道這幾個引理。

　　引理1：\((x+\omega)^p \equiv x^p + \omega^p \pmod p\)。

　　\(Proof\)：

\[由二項式定理得：(x+\omega)^p = \sum_{i=0}^p {p \choose i} x^i\omega^{p-i} \equiv x^p+\omega^p \pmod p\]

\[中間若干項項因爲組合數係數被模掉了，這樣就\text{Q.E.D}了\]

　　引理2：\(\omega^p \equiv -\omega\)。

　　\(Proof\)：

\[注意\omega可能不存在於整數域中，因此咱們不能直接運用費馬小定理\]

\[\omega^p = \omega^{p-1}\omega = (\omega^2)^\frac{p-1}{2}\omega = (a^2-n)^\frac{p-1}{2}\omega \equiv -\omega \pmod p\]

\[這些咱們運用的是上面的定義推導出來的\]

　　結論：\(x=(a+\omega)^\frac{p+1}{2}\)

　　\(Proof\)：

\[(a+\omega)^{p+1} = (a+\omega)^p(a+\omega) = (a^p+\omega^p)(a+\omega)=(a-\omega)(a+\omega)=a^2-\omega^2=a^2-(a^2-n)=n \equiv x^2 \pmod p\]

\[\Rightarrow (a+\omega)^{p+1} \equiv x^2 \pmod p\]

\[\Rightarrow x \equiv (a+\omega)^\frac{p+1}{2} \pmod p\]

\[\text{Q.E.D}\]

　　那這個式子咋求？直接上覆數。根據拉格朗日定理，最後虛部必定爲\(0\)。不會拉格朗日定理？其實我也不會那咱們反證下吧。

　　\(Proof\)：

\[假設存在x=A+B\omega，B \neq 0，使得x^2 = (A+B\omega)^2 \equiv n \pmod p\]

\[那麼必定有A^2 + 2AB\omega + B^2\omega^2 = A^2 + 2AB\omega + B^2(a^2-n) \equiv n \pmod p\]

\[\Rightarrow A^2+B^2(a^2-n)-n \equiv -2AB\omega \pmod p\]

\[\because式子左邊的虛部爲0，代表式子的右邊虛部也爲0\]

\[\therefore必須有:-2AB \equiv 0 \pmod p\]

\[\because B \not\equiv 0\]

\[\therefore A \equiv 0\]

\[\therefore B^2\omega^2 \equiv n \pmod p\]

\[\therefore \omega^2 \equiv nB^{-2} \pmod p\]

\[n是二次剩餘，B^2是二次剩餘，因此B^{-2}也應爲二次剩餘\]

\[根據\text{Legendre}符號能夠推出二次剩餘與二次剩餘的積爲二次剩餘，而\omega^2卻不是二次剩餘\]

\[矛盾，因此B=0\]

　　還有一種狀況就是\(p=2\)，此時它是質數但不是奇質數，因此特判掉便可。

　　若是\(x\)有解，其實存在兩組解，另外一組解爲\(p-x\)。

int T, N, P;

namespace Cipo {
#define rd() (1ll * rand() * RAND_MAX + rand()) // 隨機數可能範圍不夠
    int n, P; // x^2=n (%p)
    int a, t; // t = w^2 = a^2-n
    struct comp { // 自定義複數類
        int x, y;
        comp operator * (const comp &rhs) const { // (x1,y1)(x2,y2)=(x1x2+y1y2t,x1y2+x2y1)，這是複數乘，注意取模
            return (comp){(1ll*x*rhs.x + 1ll*y*rhs.y%P*t)%P, (1ll*x*rhs.y+1ll*y*rhs.x)%P};
        }
    };
    comp qpow(comp a, int b) { // 複數類快速乘
        comp res = (comp){1, 0};
        for (comp k = a; b; k = k*k, b >>= 1)
            if (b & 1) res = res*k;
        return res;
    }
    int qpow(int a, int b) { // 普通快速乘
        int res = 1;
        for (register int k = a; b; k = (ll)k*k%P, b >>= 1)
            if (b & 1) res = (ll)res*k%P;
        return res;
    }
    int sqrt(int n, int P) { // 開根號
        Cipo::n = n, Cipo::P = P;
        if (qpow(n, (P-1)>>1) == P-1) return -1; // n不是二次剩餘
        while (a = rd() % P, qpow(t = (1ll*a*a-n+P)%P, (P-1)>>1) == 1); // 判斷a^2-n獲得的是否是二次剩餘，若是是從新隨機a
        return qpow((comp){a, 1}, (P+1)>>1).x; // 返回結果
    }
}

int main() {
    srand(time(0));
    T = read();
    while (T--) {
        N = read(), P = read(); N %= P;
        if (!N) { printf("0\n"); continue; } // 若是n=0代表x=0(%p)
        if (P == 2) { printf("1\n"); continue; } // p=2特判
        int ans = Cipo::sqrt(N, P); // 開根
        if (ans == -1) printf("Hola!\n"); // 無解
        else printf("%d %d\n", min(ans, P-ans), max(ans, P-ans)); // 不然有兩解
    }
    return 0;
}

　　該算法的複雜度：\(\text{O}(\log p)\)。

　　附指望是\(2\)的證實：

\[ \sum_{i=1}^{\infty}\frac{1}{2^i}\times i （這個由指望的定義來：機率×次數）\\ =\sum_{i=0}^{\infty}\frac{1}{2^i} + \frac{1}{2}\sum_{i=0}^{\infty}\frac{1}{2^i}+\dotsb - (1+\frac{1}{2}+\frac{1}{4}+\dotsb) \\ =2+1+\frac{1}{2}+\dotsb-2 \\ =2\times 2-2=2 \]

K次剩餘

　　定義同上，只不過把\(2\)換成了\(k\)，但沒法套用上面的解法求解。這裏仍然假定\(p\)是奇質數。

　　發現構造走不通了！不要緊，咱們學習了原根！

　　衆所周知原根遍歷了全部與\(p\)互質的數，因此咱們用原根替換它們。

　　好比說\(x^k \equiv a \pmod p\)，設\(g\)爲原根，\(a=g^s\)，\(x=g^m\)，則原式變成了\((g^m)^k \equiv g^s \pmod p\)，也就是\(g^{km} \equiv g^s \pmod p\)，\(s\)能夠經過\(\text{BSGS}\)解得；\(m\)未知。而後由階的性質獲得\(km \equiv s \pmod{\varphi(p)}\)，也就是\(km \equiv s \pmod{p-1}\)，解\(m\)便可。

　　想要全部解？剛剛的同餘方程的解在模\(p-1\)的意義下的全部取值即爲全部解。事實上題目通常要求找出任意解，正常解同餘方程便可。若是上面的任意一步都出了問題（好比說\(s\)無解，或者同餘方程沒有解），那就沒有解了。

// 碼量驚人，博主寫了140+行。這裏只放最核心部分
// 這裏用K次剩餘解決二次剩餘問題
int solve(int k, int a, int P) { // x^k=a (%p)
    int s = BSGS(g, a, P); // 找出s使得g^s=a
    int d = gcd(k, P-1);
    if (s == -1 || s % d) return -1; // BSGS無解或同餘方程無解，說明該問題無解
    return qpow(g, 1ll*(s/d)*inv(k/d, (P-1)/d)%P, P); // 解同餘方程km=s(%phi(p))
}

int main() {
    init(); // 初始化質數表

    T = read();
    while (T--) {
        N = read(), P = read(); N %= P;
        if (!N) { printf("0\n"); continue; }
        if (P == 2) { printf("1\n"); continue; }
        getG(P); // 獲得模P的原根
        int ans = solve(2, N, P); // 解決K次剩餘就寫成solve(K, N, P)
        if (ans == -1) printf("Hola!\n");
        else printf("%d %d\n", min(ans, P-ans), max(ans, P-ans));
    }

    return 0;
}

Miller Rabin

　　若是咱們想要快速判斷一個數是否是質數，你可能會想到\(\text{O}(\sqrt n)\)複雜度的算法，但這個算法太劣了；或者想到運用費馬小定理，由於當\(p\)爲質數時有\(a^{p-1} \equiv 1 \pmod p\)，因此猜想逆命題也成立。咱們經過隨機\(a\)判斷式子是否成立來判斷\(p\)是否是質數。若是不許確，一次中了不表明必定是質數，那麼咱們多測幾回，能夠將其爲合數的機率降到很低。但有一類特殊的合數能讓式子恆成立，不管你怎麼選擇\(a\)，該同餘式子一直成立，這類合數叫作\(\text{Carmichael}\)數，例如\(561\)。

　　那怎麼辦呢？先介紹一下二次探測原理。

　　若\(p\)是質數，且有\(a^2 \equiv 1 \pmod p\)，那麼\(a \equiv \pm 1 \pmod p\)。前文有相似的證實，也就是移項化成整除式，這裏就不證了。

　　咱們經過這個方法來驗證。

　　把\(p-1\)分解成\(2^k \times t\)。而後隨機一個\(a\)而後讓其自乘成\(a^2\)，結合二次探測原理來判斷，若是在某一時刻自乘出\(1\)，則前一次必須是\(1\)或\(p-1\)，不然違背了原理，這個必定是合數。測試\(t\)輪後就獲得了\(a^{p-1}=a^{2^k \times t}\)。若是最後獲得的數不是\(1\)，那麼就違背了費馬小定理，\(p\)是合數。

　　聽說測試\(n\)次的正確率爲\(1-\left( \frac{1}{4} \right)^n\)。

int test[10] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}; // 聽說用這幾個數在OI的ll範圍內出錯率爲0%

int MR(ll P) {
    if (P == 1) return 0; // 毒瘤的1:他不是質數，但會卡後面
    ll t = P-1; int k = 0;
    while (!(t & 1)) k++, t >>= 1; // 分解出2
    rep(i, 0, 9) {
        if (P == test[i]) return 1; // test中的數都是質數。若是這個數出現過那它必定是質數
        else if (test[i] > P) return 0; // 測試的數比模數大，說明測試失敗了
        ll a = pow(test[i], t, P), nxt = a; // 初始化a=test[i]^t%P
        rep(j, 1, k) {
            nxt = Mult(a, a, P); // 計算nxt=a^2%P。因爲模數可能很大，這裏的Mult用的是gui速乘
            if (nxt == 1 && a != P-1 && a != 1) return 0; // 若是違背了二次探測定理，測試失敗
            a = nxt;
        }
        if (a != 1) return 0; // 違背了費馬小定理，測試失敗
    }
    return 1; // 經過了全部的測試
}

　　忽略快（gui）速乘的複雜度，該複雜度爲\(\text{O}(\log P)\)。

Pollard Rho

　　\(\text{Pollard-Rho}\)是一個用來快速分解質因數的算法，它巧妙的運用了「生日悖論」這一現象將\(N\)分解，配合着\(\text{Miller-Rabin}\)能夠將複雜度降爲\(\text{O}(N^{\frac{1}{4}})\)。

　　描述這個算法，咱們先從「生日悖論」提及。

　　班級裏有\(50\)名學生，保證生日是隨機的，不存在暗箱操做。假定一年爲\(365\)天，問出現兩我的生日相同的機率是多少？你可能脫口而出：\(\frac{50}{365}\)。

　　但事實上機率很是之高以致於幾乎確定有兩我的的生日在同一天，這個機率\(P=1-\frac{365 \choose 50}{365^{50}} \approx 97.142\%\)。爲啥直覺跟計算出來的差別如此之大？緣由在於要保證兩兩生日都不等的機率佔樣本空間的比重很是之小。換個問題：若是你走進一個\(50\)人的班級裏，你與班級裏同窗生日相等的機率又是多少呢？這時就是\(\frac{50}{365}\)了。

　　而後有一個結論：出現兩我的生日在同一天的指望人數是\(\text{O}(\sqrt N)\)。這個證實博主不會，直接跳過。

　　進入到\(\text{Pollard-Rho}\)前，還要介紹一個算法中的一個隨機函數。

　　設隨機函數\(f(x)\)，在模\(P\)的意義下，知足\(f(x)=x^2+c\)，\(c\)爲隨機參數，那麼構成了\([0,P)\rightarrow[0, P)\)下的映射，且有\(f(x)=f(x+P)\)。藉助初始參數\(x_0\)構造隨機數列\(\{ x_i\}\)，知足：\(\forall i>0 : x_i=f(x_{i-1})\)。顯然其在模\(N\)意義下一定會走向一個循環，由於定義域和值域大小有限。這個函數生成的數列仍是很隨機的，而後根據根據生日悖論有指望環大小爲\(\text{O}(\sqrt P)\)，緣由是若是出現了以前出現過的數字，顯然就會走入循環，根據上面的結論，出現兩個相同數字的指望次數就爲\(\text{O}(\sqrt P)\)。把\(x\)的一步一步映射畫出來圖，發現長得十分像\(\rho\)，算法由此得名。

　　對於將要被分解的數\(N\)，咱們選取兩個數字\(s\)和\(t\)，若是有\((\mid s-t \mid, N) \neq 1,N\)，很顯然咱們找到了一個因子，能夠遞歸分解兩部分了。可是畢竟數字太多，如何合理地隨機\(s\)和\(t\)呢？就運用隨機函數。

　　若是咱們隨機選取一些數字兩兩判斷髮現時間複雜度和空間複雜度又會彈回去。因此咱們採用「\(\text{Floyd}\)判圈法」：有兩我的在起點\(x_0\)處，在環上一我的一次走一步，另外一個一次走兩步。當他們再次相遇的時候代表一我的至少走完了一圈。這時一我的的數字記爲\(s\)，另外一個記爲\(t\)，再繼續如上判斷。運用此方法能更好判環。

　　指望下咱們須要\(\text{O}(N^{\frac{1}{4}})\)找出來一個因子。假設\(P\)可被分解，有\(P=a_1a_2,a_1 \leq a_2\)，在模\(P\)意義下，生成出來的序列出現環的指望爲\(\text{O}(\sqrt N)\)。對於環上的任意兩個數\(x'\)和\(x''\)，\((\mid x'-x'' \mid, P)=a_1\)，當且僅當：\(x' \equiv x'' \pmod {a_1}\)成立。而生成的數列\(\{x_i\}\)在模\(N\)下是隨機的，在模\(a_1\)下也是隨機的，再次運用結論獲得出現\(x'\equiv x''\)的指望爲\(\text{O}(\sqrt{a_1})\)，又由於有\(a_1 \leq \sqrt P\)，因此獲得指望爲\(\text{O}(N^{\frac{1}{4}})\)。

　　可是還有種狀況：\(x'=x''\)，會致使\((\mid x' - x'' \mid, P)=P\)。這代表兩我的相遇，選取的是環上的同一個數。遇到這種狀況代表咱們隨機數選取的很差，要更換隨機參數\(c\)從新操做。至此最原始的\(\text{Pollard-Rho}\)就出來了。

　　分析複雜度時咱們只關心最大的數分解的複雜度，分解後的數字繼續分解對總複雜度的貢獻不大。但要不斷求\(\gcd\)，總複雜度爲\(\text{O}(N^{\frac{1}{4}}\log N)\)。有\(\log\)的複雜度很差看，考慮削掉。

　　爲了減小\(\gcd\)的次數，咱們將若干次獲得的\(\mid s-t \mid\)的結果乘起來。\((a,b)=d\)時，\((ac \bmod b,b)\geq d\)必定成立。記乘積模\(N\)的值爲\(prod\)，若是\((prod,N)=1\)，顯然每一次\(\mid s-t \mid\)都與\(N\)互質；若是\(=N\)，代表其中一定有因子或\(0\)；不然就是\(N\)的因子。此時咱們退回來一步一步判斷，找出因子或者是發現相遇了。調整適當的次數連乘，均攤下來複雜度就變成了\(\text{O}(N^{\frac{1}{4}})\)。

// 下列程序段實現求N的最大質因子
inline ll f(ll x, ll c, ll p) { return inc(Mult(x, x, p), c, p); } // 隨機函數

ll Floyd(ll N) {
    ll c = rand()&15, s = 1, t = f(s, c, N); // 隨機c，s=1表示x0=1，t表示從起點跳一步
    for (;;) {
        ll _s = s, _t = t, prod = 1;
        for (int i = 1; i <= 128; i++)
            prod = Mult(prod, abs(_s - _t), N), _s = f(_s, c, N), _t = f(f(_t, c, N), c, N); // 連續處理128步
        ll d = gcd(prod, N); // d=(prod, N)
        if (d == 1) { s = _s, t = _t; continue; } // 若是d=1，表示
        for (int i = 1; i <= 128; i++) {
            if ((d = gcd(abs(s - t), N)) != 1) return d; // !=1返回
            s = f(s, c, N), t = f(f(t, c, N), c, N); // 不然繼續
        }
    }
}

ll PR(ll N) { // Pollard-Rho主過程
    if (MR(N)) return N; // 若是N是質數，直接返回
    ll ans;
    while ((ans = Floyd(N)) == N); // 找到的是同一個點，從新來
    return max(PR(ans), PR(N / ans)); // 遞歸，返回最大值
}

int main() {
    srand(time(0));
    scanf("%d", &T);
    while (T--) {
        scanf("%lld", &N);
        ll ans = PR(N); // 找最大的質因子
        if (ans == N) printf("Prime\n"); else printf("%lld\n", ans);
    }
    return 0;
}

möbius反演

　　博主曾經寫過一篇，可能不夠完整，在這裏就放個連接了：Mobius反演學習

矩陣

機率指望

多項式

數列