数论算法 剩余系相关 学习笔记 (基础回顾,(ex)CRT,(ex)lucas,(ex)BSGS,原根与指标入门,高次剩余,Miller_Rabin+Pollard_Rho)

注：转载本文须标明出处。

原文链接https://www.cnblogs.com/zhouzhendong/p/Number-theory.html

数论算法 剩余系相关 学习笔记 (基础回顾,(ex)CRT,(ex)lucas,(ex)BSGS,原根与指标入门,高次剩余,Miller_Robin+Pollard_Rho)

本文概要

　　1.　基础回顾

　　2.　中国剩余定理 (CRT) 及其扩展

　　3.　卢卡斯定理 (lucas) 及其扩展

　　4.　大步小步算法 (BSGS) 及其扩展

　　5.　原根与指标入门

　　6.　从简到难解决高次剩余问题

　　7.　大素数判定与大数分解

　　8.　相关推荐

　　9.　鸣谢

基础回顾

欧几里得算法与扩展欧几里得算法

int gcd(int x,int y){
    return y?gcd(y,x%y):x;
}

int exgcd(int a,ini b,int &x,int &y){
    if (!b){
        x=1,y=0;
        return a;
    }
    int res=exgcd(b,a%b,y,x);
    y-=(a/b)*x;
    return res;
}

快速幂

int Pow(int x,int y,int mod){
	int ans=1;
	for (;y;y>>=1,x=1LL*x*x%mod)
		if (y&1)
			ans=1LL*ans*x%mod;
	return ans;
}

欧拉函数

int phi(int n){
	int res=n;
	for (int i=2;i*i<=n;i++)
		if (n%i==0){
			res=res/i*(i-1);
			while (n%i==0)
				n/=i;
		}
	if (n>1)
		res=res/n*(n-1);
	return res;
}

欧拉定理与费马小定理

　　欧拉定理 ：如果 $a,p$ 互质，则 $a^{\varphi(p)} \equiv 1 \pmod p$ 。

　　费马小定理 ：当 $p$ 为质数的时候，则 $a^{p-1} \equiv 1 \pmod p$ 。

逆元

　　给定 $a,p$ ，求 $a^{-1} \pmod p$ 。

　　求法：如果 $a,p$ 不互质，显然无解。否则可以通过 扩展欧几里得 或 费马小定理/欧拉定理 来求。

　　当 $p$ 为质数的时候，运用费马小定理可以得到：

$$a^{-1}\equiv a^{p-2} \pmod p$$

中国剩余定理

问题模型

　　给定同余方程组

$$\begin{cases}x&\equiv&x_1&\pmod {p_1}\\x&\equiv&x_2&\pmod {p_2}\\ &&\vdots\\x&\equiv&x_n&\pmod {p_n}\end{cases}$$

　　求解 $x$ 在对 ${\rm lcm}(p_1,p_2,\cdots ,p_n)$ 取模的意义下的值，或判断无解。

模数两两互质的情况 (CRT)

　　我们考虑一种构造方法。

　　首先，设 $P=p_1p_2\cdots p_n={\rm lcm}(p_1,p_2,\cdots ,p_n)$ ，$A_i=\cfrac{P}{p_i}$ ，$B_i\equiv A_i^{-1} \pmod {p_i}$ 。

　　则，

$$x\equiv \sum_{i=1}^{n} x_iA_iB_i \pmod P$$

　　这就是著名的孙子定理，即中国剩余定理。

　　接下来我们来证明这个构造的正确性。

　　首先，由于 $p_i$ 与 $A_i$ 互质，所以必然存在逆元 $B_i$ 。

　　然后，由于 $p_j|A_i (j\neq i)$ ，所以 $(\sum_{i=1}^{n}x_iA_iB_i)\equiv x_jA_jB_j \pmod {p_j}$ 。

　　又由于 $B_j\equiv A_j^{-1} \pmod {p_j}$ ，所以 $x_jA_jB_j\equiv x_j \pmod {p_j}$ 。

　　所以构造成立。

不保证模数互质的情况 (exCRT)

　　接下来我们来解决扩展中国剩余定理。（下述证明摘自博主另一篇博文：https://www.cnblogs.com/zhouzhendong/p/exCRT.html）

　　我们只需要知道如何合并两个方程，就可以将原方程组逐个合并了，所以下面只介绍合并两个方程的方法。

$$\begin{eqnarray*}x&\equiv& a_1&\pmod {p_1}\\x&\equiv& a_2&\pmod {p_2}\end{eqnarray*}$$

　　则 $x=a_1+k_1p_1=a_2+k_2p_2$ 。

　　令 $p={\rm lcm}(p_1,p_2)$，需要求形同 $x\equiv a\pmod p$的一个式子。

$$\begin{eqnarray*}a_1+k_1p_1&\equiv&a_2+k_2p_2 &\pmod {p}\\a_2-a_1&\equiv&k_1p_1-k_2p_2 & \pmod {p}\end{eqnarray*}$$

　　我们先用 exgcd 求解方程 $p_1x+p_2y=\gcd(p_1,p_2)$。设 $g=\gcd(p_1,p_2),z=a_2-a_1$。

　　由于 $g|p$ ，所以，如果 $g|z$ 不成立，那么显然 $g+kp|z$ 也不成立，故无解。

　　否则有解。记 $t=\frac zg$ ，则 $k_1=tx_1,k_2=-tx_2$，于是，$a\equiv a_1+k_1p_1\equiv a_2+k_2p_2\pmod p$ ，单次合并完成。

扩展中国剩余定理代码

bool CRT(LL w1,LL p1,LL w2,LL p2,LL &w,LL &p){
    LL x,y,z=w2-w1,g=ex_gcd(p1,p2,x,y);
    if (z%g)
        return 0;
    LL t=z/g;
    x=Mul(x,t,p2/g);
    p=p1/g*p2;
    w=((w1+Mul(x,p1,p))%p+p)%p;
    return 1;
}

例题传送门

POJ2891 - 模板题

BZOJ2010

NOI2018Day2T1 - 经典题

51Nod1362

　

卢卡斯定理

定义

　　设 $n,m$ 为正整数， $p$ 为素数，若

$$n=a_kp^k+a_{k-1}p^{k-1}+\cdots a_1p+a_0$$

$$m=b_kp^k+b_{k-1}p^{k-1}+\cdots b_1p+b_0$$

　　所有的 $a_i,b_i$ 都为非负整数，则：

$$\binom{n}{m} \equiv \prod _{i=0}^{k} \binom {a_i}{b_i}\pmod p$$

　　或者用递归表示：

$$\binom{n}{m} \equiv \binom {\left\lfloor\frac{n}{p}\right\rfloor}{\left\lfloor\frac{m}{p}\right\rfloor}\binom{n \bmod p}{m \bmod  p}$$

证明

　　下面首先摘录一段“冯志刚版《初等数论》第37页”的内容来证明。貌似也可以从组合意义方向yy。

　　由 $p$ 为素数，可知对 $1\leq j\leq p-1 $ ，都有 $\binom{p}{j}=\cfrac{p}{j}\binom{p-1}{j-1}\equiv 0 \pmod p $，于是，

$$\begin{eqnarray*}(1+x)^p&=&1+\binom{p}{1}x+\cdots + \binom{p}{p-1}x^{p-1}+x^p \\ &\equiv & 1 + x^p \pmod p \end{eqnarray*}$$

　　利用上述结果，可知

$$\begin{eqnarray*}(1+x)^a&=&(1+x)^{a_0}\times((1+x)^p)^{a_1}\times \cdots \times ((1+x)^{p^k})^{a_k} \\ &\equiv& (1+x)^{a_0} \times (1+x^p)^{a_1} \times \cdots \times (1+x^{pk})^{a_k} \pmod p \end{eqnarray*}$$

　　对比两边 $x^m$ 项的系数（用二项式定理及 $p$ 进制数的性质）可得

$$\binom{n}{m} \equiv \prod _{i=0}^{k} \binom {a_i}{b_i}\pmod p$$

代码

LL Lucas(LL n,LL m){
	if (n<mod&&m<mod)
		return C(n,m);
	return C(n%mod,m%mod)*Lucas(n/mod,m/mod);
}

扩展卢卡斯

　　考虑模数不是质数的情况。

　　首先对模数分解质因数，令 $P=\prod p_i^{k_i}$ ，其中 $p_i$ 都是指数。我们只需要对于每一个 $p_i^{k_1}$ 为模数的时候求出答案，最后 CRT 合并就可以了。

　　接下来考虑对于 $p^k$ 形式的模数取模的做法。其中 $p$ 为质数。

　　由于 $ap^x$ 形式的数在 对于 $p^k$ 取模意义下没有逆元，所以我们要对于这一类数特殊处理。扩展卢卡斯的关键在于除去和统计阶乘中含 $p$ 因子的个数。

　　

　　首先，我们考虑计算除去含 $p$ 因子后的阶乘在对于 $p^k$ 取模意义下的值。

　　考虑一个数列 $1$~$n$ ：

$$1,2,\cdots p-1,p,p+1,\cdots ,2p,\cdots ,kp, \cdots ,\left\lfloor\frac{n}{p} \right \rfloor p, \cdots ,n$$

　　其中，所有不为 $p$ 倍数的数都直接对于我们要算的东西有直接贡献，我们直接把它累乘到答案中；一般来说 $p^k$ 都不会很大，所以我们可以枚举 $1$ 到 $p^k$ 中与 $p$ 互质的数累乘，超过 $p^k$ 之后，就相当于从 $0$ 开始循环了。

　　对于 $p$ 的倍数的数，我们将它们都除掉一个 $p$ ，然后递归调用上述过程，直到没有这种数了为止。时间复杂度 $O(p^k)$ 。

　　然后，我们统计组合数中包含的 $p$ 因子的个数。这个比较容易。

　　考虑 $1$~$n$ 的自然数中，至少含有一个因子 $p$ 的数的个数为 $\left\lfloor \frac {n}{p} \right\rfloor $ 。然后除去这些数的因子 $p$ ，不断求 $1$~$\left\lfloor \frac {n}{p} \right\rfloor $ 之间的至少含有一个因子 $p$ 的数的个数即可。

代码

　　代码貌似还不短啊。但是比起后面的高次剩余算不了什么了。

LL Pow(LL x,LL y,LL mod){
    if (y==0)
        return 1LL;
    LL xx=Pow(x,y/2,mod);
    xx=xx*xx%mod;
    if (y&1LL)
        xx=xx*x%mod;
    return xx;
}
void ex_gcd(LL a,LL b,LL &x,LL &y){
    if (!b)
        x=1,y=0;
    else
        ex_gcd(b,a%b,y,x),y-=a/b*x;
}
LL Inv(LL X,LL mod){
    if (!X)
        return 0;
    LL a=X,b=mod,x,y;
    ex_gcd(a,b,x,y);
    x=(x%b+b)%b;
    return x;
}
LL ex_lucas(LL n,LL pi,LL pk){
    if (!n)
        return 1LL;
    LL ans=1;
    for (LL i=2;i<=pk;i++)
        if (i%pi)
            ans=ans*i%pk;
    ans=Pow(ans,n/pk,pk);
    for (LL i=2;i<=n%pk;i++)
        if (i%pi)
            ans=ans*i%pk;
    return ans*ex_lucas(n/pi,pi,pk)%pk;
}
LL C(LL n,LL m,LL pi,LL pk){
    if (m>n)
        return 0;
    LL a=ex_lucas(n,pi,pk),b=ex_lucas(m,pi,pk),c=ex_lucas(n-m,pi,pk);
    LL k=0,ans;
    for (LL i=n;i;i/=pi,k+=i);
    for (LL i=m;i;i/=pi,k-=i);
    for (LL i=n-m;i;i/=pi,k-=i);
    ans=a*Inv(b,pk)%pk*Inv(c,pk)%pk*Pow(pi,k,pk)%pk;
    return ans*(P/pk)%P*Inv(P/pk,pk)%P;
}
LL C(LL n,LL m){
    LL ans=0;
    for (int i=1;i<=cnt;i++)
        ans=(ans+C(n,m,px[i],Pow(px[i],py[i],P+1)))%P;
    return ans;
}

例题

BZOJ2142

大步小步算法 (BSGS)

问题模型

　　求解 $x$ ，满足 $A^x\equiv B \pmod p$ 。$1\leq A,B,p\leq 10^9$

BSGS

　　BSGS 算法可以解决的是 $A,p$ 互质的情况。

　　首先确定一个阀值 $M$ （一般取 $\sqrt p$ 的近似值），则最终的解都可以表示成 $aM-b$ 的形式。令

$$A^{aM-b}\equiv B \pmod p$$

　　则可以得到

$$A^{aM}\equiv B\times A^b \pmod p$$

　　于是，我们只需要先枚举 $b$ 的值，把对应的 $B\times A^b$ 扔到一个 $Map$ 里面，然后再枚举左侧的 $a$ ，查询满足上式条件的 $b$ 的值即可。

　　时间复杂度 $O(\sqrt p \log p)$ 。用无序Map (unordered_map) 可以省掉 $\log$ ，但是常数较大。推荐的方法是手写 Hash 表，时间复杂度 $O(\sqrt p)$ 。

代码

struct hash_map{
    static const int Ti=233,mod=1<<16;
    int cnt,k[mod+1],v[mod+1],nxt[mod+1],fst[mod+1];
    int Hash(int x){
        int v=x&(mod-1);
        return v==0?mod:v;
    }
    void clear(){
        cnt=0;
        memset(fst,0,sizeof fst);
    }
    void update(int x,int a){
        int y=Hash(x);
        for (int p=fst[y];p;p=nxt[p])
            if (k[p]==x){
                v[p]=a;
                return;
            }
        k[++cnt]=x,nxt[cnt]=fst[y],fst[y]=cnt,v[cnt]=a;
        return;
    }
    int find(int x){
        int y=Hash(x);
        for (int p=fst[y];p;p=nxt[p])
            if (k[p]==x)
                return v[p];
        return 0;
    }
    int &operator [] (int x){
        int y=Hash(x);
        for (int p=fst[y];p;p=nxt[p])
            if (k[p]==x)
                return v[p];
        k[++cnt]=x,nxt[cnt]=fst[y],fst[y]=cnt;
        return v[cnt]=0;
    }
}Map;
int BSGS(int A,int B,int P){
    int M=max((int)(0.8*sqrt(1.0*P)),1),AM=Pow(A,M,P);
    Map.clear();
    for (int b=0,pw=B;b<M;b++,pw=1LL*pw*A%P)
        Map.update(pw,b+1);
    for (int a=M,pw=AM;a-M<P;a+=M,pw=1LL*pw*AM%P){
        int v=Map.find(pw);
        if (v)
            return a-(v-1);
    }
    return -1;
}

把 Hash 表改成 unordered_map 就只有这么一点点代码了：

unordered_map <int,int> Map;
int BSGS(int A,int B,int P){
    int M=max((int)(0.8*sqrt(1.0*P)),1),AM=Pow(A,M,P);
    Map.clear();
    for (int b=0,pw=B;b<M;b++,pw=1LL*pw*A%P)
        Map[pw]=b+1;
    for (int a=M,pw=AM;a-M<P;a+=M,pw=1LL*pw*AM%P)
        if (Map[pw])
            return a-(Map[pw]-1);
    return -1;
}

扩展BSGS

　　不保证 $A,p$ 互质。

　　令 $g=\gcd(A,p)$ ，则：

$$\frac{A}{g} A^{x-1} \equiv \frac {B}{g} \pmod {\frac{p}{g}}$$

　　当 $B$ 不是 $g$ 的倍数的时候，显然无解。

　　于是又得到了一个新的问题，性质与原问题一样。于是可以不断递归求解，直到 $g=1$ 为止。

　　最终会得到一个形如

$$\frac {A^k} {d} A^{x-k} \equiv \frac{B}{d} \pmod {\frac {p}{d}}$$

　　的方程。

　　其中 $d$ 为各个 $g$ 的乘积。

　　于是问题转化成了求解方程

$$ a A^x \equiv b \pmod {p^\prime}$$

　　于是 BSGS 就可以解决的问题了，最终把 $k$ 加回去就可以了。

　　注意，上述做法不能解决 $x<k$ 的情况，对于这种情况，我们只需要特殊处理一下就可以了。

代码

　　此处略去 Hash 表。

int ExBSGS(int A,int B,int P){
    A%=P,B%=P;
    int k=0,v=1;
    while (1){
        int g=gcd(A,P);
        if (g==1)
            break;
        if (B%g)
            return -1;
        k++,B/=g,P/=g,v=1LL*v*(A/g)%P;
        if (v==B)
            return k;
    }
    if (P==1)
        return k;
    int M=max((int)sqrt(1.0*P),1),AM=Pow(A,M,P);
    Map.clear();
    for (int b=0,pw=B;b<M;b+=1,pw=1LL*pw*A%P)
        Map.update(pw,b+1);
    for (int a=M,pw=1LL*v*AM%P;a-M<P;a+=M,pw=1LL*pw*AM%P){
        int v=Map.find(pw);
        if (v)
            return a-(v-1)+k;
    }
    return -1;
}

例题

BZOJ2480 - 扩展bsgs模板题

原根与指标入门

　　声明：这里摘取了部分 PoPoQQQ 的课件《原根与指标 v1.0.1 .ppt》的内容，感谢 PoPoQQQ。读者如果想要看这个课件，请自行移步 UOJ 群文件中下载。

阶

定义

　　对于一个正整数 $a$ ，求满足 $a^x \equiv 1 \pmod m$ 的最小正整数 $x$ 。保证 $a,m$ 互质。这个最小正整数 $x$ 称作 $a$ 在对 $m$ 取模意义下的阶，记做 ${\rm ord}_m(a)$ ，在 $m$ 的值十分明确的时候，可以记做 ${\rm ord}(a)$ 。

性质

　　1.　$a^n\equiv 1 \pmod m$ 的充要条件是 ${\rm ord}_m(a) | n$ 。推论： ${\rm ord}_m(a) | \varphi (m) $ 。

　　2.　若 $a\equiv b \pmod m$ ，则 ${\rm ord} _m(a) = {\rm ord}_m(b)$ 。

　　3.　若 $a^n\equiv a^i \pmod m $，则 $n\equiv i \pmod {{\rm ord}_m(a)}$ 。

　　4.　令 $n= {\rm ord} _m (a) $ ，则 $a^0,a^1,\cdots ,a^{n-1}$ 对 $m$ 取模两两不同。

　　由于上述性质的证明比较简单，所以这里就不写证明了。

原根

定义

　　若 $g$ 是模 $p$ 意义下的原根，则 $g$ 满足 ${\rm ord}_p(g) = \varphi(p)$ 。

性质

　　1.　模 $p$ 意义下存在原根，当且仅当 $p$ 是如下形式的数： $2,4,x^a,2x^a$ 。（$x$ 为奇素数， $a$ 为正整数）

　　2.　当 $p$ 为奇素数时，模 $p$ 意义下的原根个数为 $\varphi(\varphi(p))$ 。

　　3.　对于质数 $p$ ， $\varphi (p)=p-1$ ，将 $p-1$ 分解质因数，得到 $p-1=\prod {p_i^{q_i}}$ ，则正整数 $g$ 是模 $p$ 意义下的原根的充分必要条件是：对于所有 $i$ ，$g^{\frac{p-1}{p_i}}\not\equiv 1 \pmod p$ 。　　证明： 充分性很显然。必要性：首先考虑阶的第 1 点性质，可以得知 ${\rm ord}_p(g)| p-1$ ，那么，如果这个值比 $p-1$ 小，必然可以找到一个 $i$ ，使得 ${\rm ord}_p(g) | \frac{p-1}{p_i} $ ，那么 $g^{\frac{p-1}{p_i}}\equiv 1 \pmod p$ ，故 $g$ 不是原根，否则，说明 ${\rm ord }_p(m) = p-1 =\varphi(p)$ ，$g$ 是原根。

如何求原根

　　由于原根一般都很小，所以我们直接暴搜原根！

　　具体多小？根据 PoPoQQQ 的课件，在 $10^9$ 范围内，最小的原根大小大约不超过 $n^{0.25}$ 。

代码

bool Get_g_Check(int P,int x){
    for (int i=1;i<=Fac_tot;i++)
        if (Pow(x,(P-1)/Fac_p[i],P)==1)
            return 0;
    return 1;
}
int Get_g(int P){
    Fac_tot=0;
    int v=P-1;
    for (int i=1;prime[i]*prime[i]<=v&&i<=pcnt;i++)
        if (v%prime[i]==0){
            Fac_p[++Fac_tot]=prime[i];
            while (v%prime[i]==0)
                v/=prime[i];
        }
    if (v>1)
        Fac_p[++Fac_tot]=v;
    for (int i=2;;i++)
        if (Get_g_Check(P,i))
            return i;
    return -1;
}

　　

指标

定义

　　设 $g$ 是模 $p$ 意义下的一个原根，若 $a,p$ 互质，则存在一个唯一的正整数 $x$ 满足 $g^x\equiv a \pmod m $ 。这个 $x$ 就叫做 “以 $g$ 为底的 $a$ 对模的一个指标"，记作 $x={\rm ind}_g(a) $ 。

求法

　　直接 BSGS 。

高次剩余

问题模型

　　求解方程

$$x^A \equiv B \pmod p$$

　　$1\leq A,B,p \leq 10^9$

当 $p$ 为质数时

　　传送门 - 51Nod1038题意与代码

　　首先，我们可以找到模 $p$ 意义下的原根 $g$ 。

　　令 $t={\rm ind}_g(B)$ ，设 $x=g^y$ ，则

$$g^{Ay} \equiv g^t \pmod p$$

　　即

$$Ay \equiv t \pmod {\varphi (p)}$$

　　把 $y$ 解出来即可得到 $x$ 。

　　代码见传送门。

当 $p$ 为奇数时

　　传送门 - BZOJ2219题意与代码

　　注：这里这道 bzoj 题只要求输出解的个数。

　　

　　考虑对 $p$ 分解质因数。

　　假设 $p=\prod p_i^{q_i}$ ，$p_i$ 为奇素数。那么我们显然可以对于每一个 $p_i$ 分开求解，最后可以通过中国剩余定理得到：在原模数意义下的解的总数，为所有 模 $p_i^{q_i}$ 意义下的解的个数 的乘积。

　　考虑求解形如 $x^A \equiv B \pmod {p^k}$ 的方程。

　　如果 $\gcd(B,p^k) = p^k$

　　　　那么，只需要保证 $p^{\left\lceil\frac kA\right\rceil}|x$ 即可。

　　　　令 $base = p^{\left\lceil\frac kA\right\rceil}$ ，则解的形式为 $\alpha\cdot base$ ，解的个数为 $p^{k-\left\lceil\frac kA\right\rceil} $ 。

　　如果 $\gcd(B,p^k) = p^q \ \ \ (0\leq q<k)$

　　　　则显然有 $A | q$ 。

　　　　令 $B^\prime = \frac B {p^q}$ ，则可以把原方程转化成求

$$p^q y^A \equiv p^q B^\prime \pmod {p^k}$$

　　　　即

$$y^A \equiv B^\prime \pmod {p^{k-q}}$$

　　　　由于 $\gcd(B^\prime , p^{k-q})=1$ ，那么求解 $y$ 显然可以用之前介绍的方法来得到。

　　　　最终， $x \equiv y \times p^{\frac A q} \pmod {p^k}$ 。于是，如果 $y$ 在模 $p^{k-q}$ 意义下的解的个数为 $ans$ ，那么 $x$ 在模 $p^k$ 意义下的解的个数为 $p^{q-\frac Aq}\times ans$ 。

　　代码见传送门。

当 $p$ 没有限制时

　　传送门 - 51Nod1123题意与代码

　　其他都一样，在求解 $x^A \equiv B \pmod {p^k}$ 这一类问题时，当 $p=2$ 时，问题会比较特殊。

　　考虑到 $x^A \equiv B \pmod {2^k}$ 的充要条件是 $x^A \equiv B \pmod {2^{k-1}}$ 。

　　反过来，如果 $x^A \equiv B \pmod {2^{i-1}}$ ，那么只可能导致 $x^A \equiv B \pmod {2^i}$ 或者 $(x+2^{i-1})^A \equiv B \pmod {2^i}$ 。只需要写个 dfs 就好了。

　　但是由于这题要求出具体的解，所以代码有点长。

　　代码见传送门。

小结

　　高次剩余的求解中，要用到之前提到过的几乎所有算法。数论算法在我们眼中突然变成了码农题，240行高次剩余，瞬间虐爆数据结构。

　　”告辞“！ 剩余！

如果您还想来一道模板题的话

51Nod1039

大素数判定和大整数分解算法 - Miller_Rabin & Pollard_Rho

　　懒得写了。

　　自行百度。

　　网上的博客写的还是挺详细的。

　　我来这里只是贴个板子。

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
typedef long double LD;
namespace Pollard_Rho{
    int prime[9]={2,3,5,7,11,13,17,19,23};
    ULL RR;
    int Pcnt;
    LL p[70];
    vector <LL> res;
    LL R(LL mod){
        return (RR+=4179340454199820289LL)%mod;
    }
    LL Mul(LL x,LL y,LL mod){
        LL d=(LL)floor((LD)x*y/mod+0.5);
        LL res=x*y-d*mod;
        if (res<0)
            res+=mod;
        return res;
    }
    LL Pow(LL x,LL y,LL mod){
        LL ans=1%mod;
        for (;y;y>>=1,x=Mul(x,x,mod))
            if (y&1)
                ans=Mul(ans,x,mod);
        return ans;
    }
    bool Miller_Rabin(LL n){
        if (n<=1)
            return 0;
        for (int i=0;i<9;i++)
            if (n==prime[i])
                return 1;
        LL d=n-1;
        int tmp=0;
        while (!(d&1))
            d>>=1,tmp++;
        for (int i=0;i<9;i++){
            LL x=Pow(prime[i],d,n),p=x;
            for (int j=1;j<=tmp;j++){
                x=Mul(x,x,n);
                if (x==1&&p!=1&&p!=n-1)
                    return 0;
                p=x;
            }
            if (x!=1)
                return 0;
        }
        return 1;
    }
    LL f(LL x,LL c,LL mod){
        return (Mul(x,x,mod)+c)%mod;
    }
    LL gcd(LL x,LL y){
        return y?gcd(y,x%y):x;
    }
    LL Get_Factor(LL c,LL n){
        LL x=R(n),y=f(x,c,n),p=n;
        while (x!=y&&(p==n||p==1)){
            p=gcd(n,max(x-y,y-x));
            x=f(x,c,n);
            y=f(f(y,c,n),c,n);
        }
        return p;
    }
    void Pollard_Rho(LL n){
        if (n<=1)
            return;
        if (Miller_Rabin(n)){
            res.push_back(n);
            return;
        }
        while (1){
            LL v=Get_Factor(R(n-1)+1,n);
            if (v!=n&&v!=1){
                Pollard_Rho(v);
                Pollard_Rho(n/v);
                return;
            }
        }
    }
    void work(LL n){
        res.clear();
        Pollard_Rho(n);
    }
}

相关推荐

　　基础数论函数相关 - 莫比乌斯反演 - 杜教筛 - https://www.cnblogs.com/zhouzhendong/p/Number-theory-Function.html

　　快速傅里叶变换/快速数论变换 入门与经典例题 - https://www.cnblogs.com/zhouzhendong/p/Fast-Fourier-Transform.html

鸣谢

　　冯志刚版《初等数论》

　　yyb - BSGS - http://www.cnblogs.com/cjyyb/p/8810050.html

　　PoPoQQQ 原根与指标 v1.0.1 .ppt

　　感谢 51Nod 和 BZOJ 提供大量例题来源。