数论入门2——gcd,lcm,exGCD,欧拉定理,乘法逆元,(ex)CRT,(ex)BSGS,(ex)Lucas,原根,Miller-Rabin,Pollard-Rho

数论入门2

另一种类型的数论...

GCD,LCM

定义\(gcd(a,b)\)为a和b的最大公约数，\(lcm(a,b)\)为a和b的最小公倍数，则有：

将a和b分解质因数为\(a=p1^{a1}p2^{a2}p3^{a3}...pn^{an},b=p1^{b1}p2^{b2}p3^{b3}...pn^{bn}\)，那么\(gcd(a,b)=\prod_{i=1}^{n}pi^{min(ai,bi)},lcm(a,b)=\prod_{i=1}^{n}pi^{max(ai,bi)}\)(0和任何数的最大公约数是这个数，最小公倍数是0)

显然成立的吧。

所以我们有\(gcd(a,b)*lcm(a,b)=ab\)。

求\(gcd(a,b)\)的方法：

首先我们有：\(gcd(a,b)=gcd(b,a-b)(a>=b)\)

显然的啊，两个数相减怎么可能影响他们的最大公约数，毕竟是“公”约数，相减肯定还存在这个因子啊。

所以拓展一下我们就有：\(gcd(a,b)=gcd(b,a \% b)(a>=b)\)

取模就是不断相减，所以也是成立的呀。

所以我们可以写出简洁的代码：

int gcd(int a,int b){return !b?a:gcd(b,a%b);}

很显然这个复杂度是log的，每次至少减少一半。

lcm的话，\(a/gcd(a,b)*b\)即可。

然后呢，gcd和lcm具有交换律和结合律：

\(gcd(a,b,c)=gcd(gcd(a,b),c)=gcd(c,b,a)\)

所以求多个数的gcd或lcm的时候就可以递推求啦。

还有一个东西就是\(gcd(ka,kb)=k*gcd(a,b)\)，也是显然的。

exGCD

对于方程\(ax+by=c\)，显然若\(gcd(a,b)|c\)的时候才有解。

因为不整除的话\(ax+by\)肯定有\(c\)不包含的因子，那样的话怎么可能有解呢？

然后呢，我们可以通过这样的方式快速借助求gcd的过程求出\(ax+by=gcd(a,b)\)的一组解。

当\(a=gcd(a,b),b=0\)时，显然\(x=1,y=0\)。

令\(a = b, b = a \% b\)，则有方程\(b *x1 +(a \% b) * y1 = gcd(b, a \% b)\)

又因为\(gcd(a, b) = gcd(b, a \% b)\)，且\(a \% b = a - b * ⌊a / b⌋\)

则\(b * x1 + (a - b * ⌊a / b⌋) * y1 =gcd(a, b)\)

整理得：\(a * y1 +b * (x1 - ⌊a / b⌋ *y1) = gcd(a, b)\)

所以原方程中：\(x = y1, y = x1 - ⌊a / b⌋ *y1\)。于是我们只要递归求出\(x1, y1\)就能求出\(x, y\)。

我们现在已经求得了\(ax +by = gcd(a, b)\)的解,那么对于方程\(ax + by = c (gcd(a, b) | c)\)呢？

因为已经知道\(a *x1 +b * y1 = gcd(a, b)\)的解\(x1, y1\)，左右两边同乘以\(c / gcd(a, b)\) 得：

\(a * x1 * c / gcd(a, b) +b * y1 * c / gcd(a, b) = c\)

则原方程的一组解\(，x2 = x1 * c / gcd(a, b)， y2 = y1 * c / gcd(a, b)\)

由此得出解集\(，{(x, y) | x = x2 + k * b / gcd(a, b)， y = y2 - k * a / gcd(a, b), k ∈ Z}\)

对于线性同余方程\(ax\equiv c\quad(mod\quad b)\)，我们怎么做呢？

显然可以转化成\(ax+by=c\)的形式，然后直接用扩欧解决啦。

贴个代码：

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

欧拉定理

费马小定理：

若p是质数，且\(gcd(a,p)=1\)，则\(a^{p-1}\equiv1(mod\quad p)\)。

欧拉定理：

若\(gcd(a,p)=1\)，则\(a^{\varphi(p)}\equiv1(nod \quad p)\)。

扩展欧拉定理：

\(a^b\equiv \begin{cases} a^{b\%\varphi(p)}~~~~~~~~~~~gcd(a,p)=1\\ a^b~~~~~~~~~~~~~~~~~~gcd(a,p)\neq1,b<\varphi(p)\\ a^{b\%\varphi(p)+\varphi(p)}~~~~gcd(a,p)\neq1,b\geq\varphi(p) \end{cases}~~~~~~~(mod~p)\)

\(\varphi\)为欧拉函数。

证明：欧拉定理和扩展欧拉定理的证明

有个作用就是在算\(a^b\)的时候对b进行取模，使运算更高效，具体见luogu欧拉定理模板题

乘法逆元

若\(gcd(a,b)=1,\)\(ax\equiv 1(mod~p)\)，则称x为a在模p意义下的乘法逆元。

那么求法就很显然了。

欧拉定理：a在模p意义下的逆元为\(a^{\varphi(p)-1}\)。

exGCD：当做一个不定方程求，但是要通过+kp的方式变为一个正数。

此外还有一个线性求逆元的算法。

首先，1在任何模数下的逆元都是1。

然后设\(p=i*k+r,r<i,1<i<p\)。

那么我们有\(i*k+r\equiv 0(mod~p)\)

两边同时乘上\(i^{-1}*r^{-1}\)得：

\(k*r^{-1}+i^{-1}\equiv 0(mod~p)\)

\(i^{-1}\equiv -k* r^{-1}(mod~p)\)

因为\(k=\lfloor p/i\rfloor,r=p\%i\)，所以我们得到了逆元的递推式：

a[i]=-(p/i)*a[p%i]%p;

以及提一句，逆元是完全积性函数，从定义式里面就能看出来。

对于一个前缀积的序列的话，我们可以求出最后一项的逆元，然后\(O(n)\)递推回去，实际上逆元可以理解为模意义下的\(\frac{1}{x}\)，怎么递推的话也很显然了，以阶乘为例：

fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%P;
ifa[n]=ksm(fac[n],P-2);
for(int i=n-1;~i;i--) ifa[i]=ifa[i+1]*(i+1)%P;

还有就是：

\((a/b)\%p=(a\%(bp))/ p\)

证明：

\((a/b)\%p=a/b-⌊(a/b)/p⌋*p\)

\(=a/b-⌊a/(b*p)⌋*p\)

\(=a/b-⌊a/(bp)⌋*b*p/b\)

\(=(a\%(bp) )/p\)

CRT

已知系数全部为1的线性同余方程组：

\(x \equiv a_i(mod ~p_i)\)，其中\(p_i\)两两互质，求x的最小非负整数解。

那么我们考虑这样一个思路：

从1到n考虑每个同余方程组，我们设\(M=\prod_{i=1}^{n}p_i\)，那么当我们考虑到第i个方程的时候，我们设\(T=\frac{M}{p_i}\)，那么我们发现把第i个方程的解\(x\)加上\(kT\)的话，\(x+kT\)带到别的同余方程组的余数不会改变。

假设我们已经求出了前\(i-1\)个方程的通解\(x\)，求出了第\(i\)个方程的一个解\(x_i\)，那么显然，\(Tx_i+x\)是前\(i\)个方程的解。

我们用exGCD求出第i个线性同余方程组\(Tx_i\equiv a_i(mod~p_i)\)的解，一步步合并即可，因为\(p_i\)显然和\(T\)互质，所以实际上是求出逆元之后乘上\(a_i\)。

exCRT

已知系数全部为1的线性同余方程组：

\(x \equiv a_i(mod ~p_i)\)，其中\(p_i\)都是正整数，求x的最小非负整数解。

上面的算法失效了。为什么？都不一定互质了怎么可能还去求逆元？

所以我们要考虑一个新的思路：

假设我们已经求出了前\(i-1\)个方程的通解\(x\)，设\(M\)为前i-1个数的lcm，那么我们发现\(x+kM\)如果能够满足第\(i\)个方程，那么前i个方程就都解出来啦。

我们写个式子：\(x+kM \equiv a_i(mod~p_i)\)，移项得\(kM \equiv a_i-x(mod ~p_i)\)。

那么我们解出来k之后，前i个方程的通解就出来啦，一步步做下去即可。

贴个代码：

ll mul(ll x,ll y,ll p){ll g=0;for(;y;y>>=1,x=(x+x)%p) if(y&1) g=(g+x)%p;return g;}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b){x=1;y=0;return a;}
    ll g=exgcd(b,a%b,y,x);y-=a/b*x;
    return g;
}
for(int i=1;i<=n;i++) scanf("%lld%lld",&b[i],&a[i]);
ans=a[1];M=b[1];
for(int i=2;i<=n;i++)
{
    ll A=M,B=b[i],c=(a[i]%B-ans%B+B)%B,x,y;
    ll gcd=exgcd(A,B,x,y),d=B/gcd;
    if(c%gcd) return 0;x=mul(x,c/gcd,d);
    ans+=x*M;M*=d;ans=(ans%M+M)%M;
}

注意一下，因为是求lcm，所以那个B要除以gcd。

顺便扩展一下，如果\(x\)对应的系数不为1，那么怎么办呢？

其实就是这个题：【NOI2018】屠龙勇士

我们可以用扩欧求出每个方程的解集，这个解集显然是\(x=x_i+k*b\)，写成同余方程的形式就是\(x \equiv xi(mod ~b)\)

然后这不就可以用\(exCRT\)解决了吗。

BSGS

对于这样一个指数型的同余方程组：\(a^x \equiv b~(mod~p),gcd(a,p)=1\)，求x的最小正整数解。

首先我们有：\(a^{\varphi(p)} \equiv 1(mod ~p)\)。

怎么办？\(1-\varphi(p)\)枚举？TLE。

我们设\(x=im-r\)，其中\(m=\lceil \sqrt{\varphi(p)}\rceil\)。

那么原方程可以转化为\(a^{im} \equiv b*a^{r}(mod ~p)\)。

这样的话，\(i,r\)都只有\(m-1\)或\(m\)种取值，我们可以先枚举\(r\)，把\(b*a^r\)存在hash表中，映射值为\(r\)，然后我们枚举\(i\)，递推算\(a^{im}\)，这样的话，我们在hash表里面如果能找到，就可以直接返回\(im-r\)为答案啦。

时间复杂度\(O(\sqrt{\varphi(p)})\)。

exBSGS

对于这样一个指数型的同余方程组：\(a^x \equiv b~(mod~p)\)，求x的最小正整数解。

刚才那个算法失效了，为什么？欧拉定理不成立了啊！

所以我们可以这么想：让\(a,p\)互质，这样就可以用\(BSGS\)做啦。

设\(d=gcd(a,p)\)，那么显然\(d|b\)的时候方程才有解。

我们可以对这个方程整体除以d，设\(c=\frac{a}{d}\)，\(b'=\frac{b}{d},p'=\frac{p}{d}\)得：

\(c*a^{x-1} \equiv b'(mod~p')\)。

我们记录c，然后让b和p不断除以d，等到a和p互质的时候，直接上\(BSGS\)就好啦。

那个c并不影响我们进行\(BSGS\)呢，枚举的时候乘上就好啦。

贴个代码：Luogu模板题

#include <bits/stdc++.h>
using namespace std;
int a,b,p,g;unordered_map<int,int>mmp;
int bsgs(int A,int B,int P)
{
    if(B==1) return 0;//特判B=1
    int cnt=0,C=1,d;
    while((d=__gcd(A,P))!=1)
    {//相除到互质为止，并判断无解的情况
        if(B%d) return -1;
        P/=d;B/=d;C=1ll*C*A/d%P;
        cnt++;if(C==B) return cnt;//可能除的过程中就出解了
    }
    int n=sqrt(P)+1,t=1;mmp.clear();//正常的BSGS
    for(int i=0;i<n;i++) mmp[1ll*t*B%P]=i,t=1ll*t*A%P;//先存B*A^i
    C=1ll*t*C%P;
    for(int i=1;i<=n;i++)//再用C*A^(i*n)在hash表里面寻找解
    {
        if(mmp.find(C)!=mmp.end()) return i*n-mmp[C]+cnt;
        C=1ll*t*C%P;
    }
    return -1;
}
int main()
{
    while(scanf("%d%d%d",&a,&p,&b)==3&&a&&b&&p)
    {
        if(~(g=bsgs(a,b,p))) printf("%d\n",g);
        else puts("No Solution");
    }
}

原根

假设一个数g对于P来说是原根，那么\(g^i mod P\)的结果两两不同,且有 \(1<g<P, 1<i<P,\)那么g可以称为是P的一个原根。

所以呢，原根的一个作用就是可以把1-m-1这些数的某种相乘同余关系转化成原根的幂次相乘关系，进而转化成幂次相加关系，比如说可以神奇地方便卷积呀。（别的作用我太菜了目前还没有找到）

推荐一道题：【SDOI2015】序列统计

然后怎么求原根呢

一个显然但是低效的方法就是枚举原根是多少，然后枚举循环节是多长\((1-\varphi(p))\)，然后判断一下这个循环节的\(a^x\)是否和1在模p意义下同余，如果有说明当前枚举的数字不是原根。

但是这样太低效了，有个结论就是：枚举i只需要枚举φ(p)的因数就好了。

截的这位博主的证明Orz

所以我们把\(\varphi(x)\)分解质因数，用所有质因数判断即可。

贴个代码：

#include <bits/stdc++.h>
using namespace std;
int p;
int phi(int x)
{
    int g=x,X=x;
    for(int i=2;i*i<=X;i++) if(x%i==0)
    {
        g=g/i*(i-1);
        while(x%i==0) x/=i;
    }
    if(x>1) g=g/x*(x-1);
    return g;
}
int ksm(int x,int y,int P){int g=1;for(;y;y>>=1,x=1ll*x*x%P) if(y&1) g=1ll*g*x%P;return g;}
int G(int x)
{
    int ph=phi(x),t=ph,tot=0;static int a[100005];
    for(int i=2;i*i<=ph;i++) if(t%i==0)
    {
        a[++tot]=ph/i;
        while(t%i==0) t/=i;
    }
    if(t!=1) a[++tot]=ph/t;
    for(int g=2,ok;;g++)
    {
        if(ksm(g,ph,x)!=1) continue;ok=1;
        for(int j=1;j<=tot;j++) if(ksm(g,a[j],x)==1){ok=0;break;}
        if(ok) return g;
    }
}
int main()
{
    scanf("%d",&p);
    printf("%d\n",G(p));
}

Lucas

\(Lucas\)定理可以求解组合数取模的问题，就是求\(C_n^m~mod~p\)，p为质数。

上递归式：\(C_n^m=C_{n\%p}^{m\%p}*C_{n/p}^{m/p}~\%p\)。

时间复杂度\(O(log_p(n)*p)\)。

用于n,m很大但是p不是很大的情况下。

贴个代码：Luogu模板

#include <bits/stdc++.h>
using namespace std;
int k,n,m,p;
long long a[200005],b[200005];
long long lucas(int x,int y)
{
    if(x<y) return 0;
    else if(x<p) return b[x]*a[y]*a[x-y]%p;
    else return lucas(x/p,y/p)*lucas(x%p,y%p)%p;
}
int main()
{
    scanf("%d",&k);
    while(k--)
    {
        scanf("%d%d%d",&n,&m,&p);
        a[0]=a[1]=b[0]=b[1]=1;
        for(int i=2;i<=n+m;i++) b[i]=b[i-1]*i%p;
        for(int i=2;i<=n+m;i++) a[i]=(p-p/i)*a[p%i]%p;
        for(int i=2;i<=n+m;i++) a[i]=a[i-1]*a[i]%p;
        printf("%lld\n",lucas(n+m,m));
    }
}

（其实这只是Lucas定理引出的一个组合数取模计算式，并不是Lucas定理）

真正的Lucas定理是这样子的：

把n和m写成p进制数的样子（如果长度不一样把短的补成长的那个的长度）：

\(n=(a_0a_1\dots a_k)_p \\ m=(b_0b_1\dots b_k)_p\)

那么：

\(C_n^m\equiv \prod _{i=0}^k C_{a_i}^{b_i} \mod p\)

证明戳这里

那么看一下上面那个式子，是不是挺显然的呀。

顺便这里我们可以得出组合数奇偶性的规律：

\(C_n^m\)为奇数的充要条件为\(n\&m=\)\(m\)。

至于原因，把n和m对齐写成二进制形式套用Lucas定理就可以很显然看出来。

暂时还没有发现别的应用哎（我太菜了）

exLucas

如果p不是质数该怎么办呢？

考虑这样做：

我们将p分解质因数得到：\(p=\prod _i p_i^{a_i}\)。

这样一个显然的性质就是每一个\(p_i^{a_i}\)都是互质的。

我们如果对于每个单个的\(p_i^{a_i}\)求出了\(C_n^m~ mod~p_i^{a_i}\)的结果，也即\(x \equiv c_i~(mod ~p_i^{a_i})\)，那么用\(CRT\)一合并，不就得出了问题的答案吗？

所以我们现在考虑怎么求\(C_n^m~ mod~p_i^{a_i}\)。

\(p_i^{a_i}\)也不是质数，没法套用\(Lucas\)定理，我们可以这样考虑：把组合数拆成三个阶乘相乘的形式，然后分别算出来，其中两个阶乘还要乘上逆元。

现在考虑怎么求\(n!~mod~p^a\)。

一个著名的例子就是\(n=19,p=3,a=2\)。

\(\begin{aligned} 19!&=1*2*3*\cdots *17*18*19 \\ &=1*2*4*5*7*8*10*11*13*14*16*17*19*(3*6*9*12*15*18) \\ &=(1*2*4*5*7*8)*(10*11*13*14*16*17)*19*3^6*(1*2*3*4*5*6) \end{aligned}\)

我们发现\(10*11*13*14*16*17 \equiv 1*2*4*5*7*8~(mod ~3^2)\)，那么这一部分可以\(O(n/p^a)\)直接暴力算出来然后快速幂一发。

还会剩下半截，这个例子就是一个\(19\)，这个暴力算出来即可。

还有一个\(3^6\)，这个我们在算阶乘的时候不算它，原因后面会提到。

剩下的那个\(6!\)递归即可。

时间复杂度\(O(p^a*log_{p}(n)*log_2(n))\)。

为什么不算那个\(3^6\)呢，我们考虑这样一个问题：算出来的\(m!,(n-m)!\)都要乘上关于\(p_i ^{a_i}\)的逆元，那么就要和\(p_i ^{a_i}\)互质。所以阶乘中我们要把关于\(p_i\)的因子都提取出来，到后面再乘回去。

怎么快速计算\(n!\)中有多少个p呢？一个显然的式子就是：\(\lfloor\frac{n}{p^1}\rfloor+\lfloor\frac{n}{p^2}\rfloor+\lfloor\frac{n}{p^3}\rfloor+...\)。

因为一次方要记一次数，二次方要记两次数，......，然后i次方计数个数就是\(\lfloor\frac{n}{p^i}\rfloor\)。

这个式子写成递推式就是\(f(n)=f(\lfloor\frac{n}{p}\rfloor)+\lfloor\frac{n}{p}\rfloor\)。

所以最后的式子就是\(g_n*inv(g_m)*inv(g_{n-m})*p^{f_n-f_m-f_{n-m}}\)。

至此问题已经解决，看着很麻烦，实际上代码并不长。

时间复杂度\(O(\sum p^a*log_{p}(n)*log_2(n)+nlog_2(p))\)，可以接受。

#include <bits/stdc++.h>
using namespace std;
#define ll long long
ll n,m,p;
ll ksm(ll x,ll y,ll P){ll g=1;for(;y;y>>=1,x=x*x%P) if(y&1) g=g*x%P;return g;}
ll exgcd(ll a,ll b,ll &x,ll &y){if(!b){x=1;y=0;return a;}ll g=exgcd(b,a%b,y,x);y-=a/b*x;return g;}
ll inv(ll x,ll P){ll y,z;exgcd(x,P,y,z);return (y%P+P)%P;}
ll crt(ll a,ll b){return (p/b)*inv(p/b,b)%p*a%p;}
ll fac(ll n,ll pi,ll pk)
{
    if(!n) return 1;ll g=1;
    for(int i=1;i<=pk;i++) if(i%pi) (g*=i)%=pk;
    g=ksm(g,n/pk,pk);
    for(int i=1;i<=n%pk;i++) if(i%pi) (g*=i)%=pk;
    return g*fac(n/pi,pi,pk)%pk;
}
ll C(ll n,ll m,ll pi,ll pk)
{
    ll pa=fac(n,pi,pk),pb=fac(m,pi,pk),pc=fac(n-m,pi,pk),x=0;
    for(ll i=n;i;i/=pi) x+=i/pi;
    for(ll i=m;i;i/=pi) x-=i/pi;
    for(ll i=n-m;i;i/=pi) x-=i/pi;
    return pa*inv(pb,pk)%pk*inv(pc,pk)%pk*ksm(pi,x,pk)%pk;
}
ll exlucas(ll n,ll m,ll P)
{
    ll g=0,x=P,t;
    for(ll i=2;i*i<=P;i++) if(x%i==0)
    {
        t=1;while(x%i==0) t*=i,x/=i;
        (g+=crt(C(n,m,i,t),t))%=P;
    }
    if(x>1) (g+=crt(C(n,m,x,x),x))%=P;
    return g;
}
int main()
{
    scanf("%lld%lld%lld",&n,&m,&p);
    printf("%lld\n",exlucas(n,m,p));
}

Miller-Rabin

这是一种可以在与数字大小几乎无关的复杂度下判断一个数是否为质数的算法。

首先我们有费马小定理：\(a^{p-1}\equiv 1\pmod{p}(a<p,p\ is\ a\ prime)\)

那么我们是不是可以多选几个a，保证a与p互质，然后判断一个数是不是质数？

这样就太low了。因为有些数字无论如何都能通过这样的测试，比如说\(561=51*11\)就可以通过所有这样的测试。

所以有一个准确率更高的方法：二次探测法。

先上结论：如果p为质数,\(x<p\)，\(x^2 \equiv1(mod~ p)\)，那么$x=1 _orp-1 $。

证明：\(x^2 \equiv1(mod~ p)~~->~~(x+1)(x-1) \equiv0(mod~p)\)。

所以\(p|(x+1)(x-1)\)，因为p是大于x的质数，所以\(p=x+1~~or~~p \equiv x-1(mod~ p)\)。

所以\(x=1~or~p-1\)。

有了这个结论，我们可以随机几个a，a选为质数最好，然后将\(p-1\)不断除以2直到\(p-1\)是个奇数为止，记\(p-1=2^s*q\)，那么我们先令\(x=a^q\)，如果这个时候\(x^2~mod ~p=1\)，那么我们可以看\(x\)的值，如果x不是\(1\)或者\(p-1\)那么就可以返回合数了。最后我们还可以用上面那个费马小定理的方法顺便提高一下正确率。

时间复杂度\(k*log_2(p)\)。其中k是测试用的a的个数。

正确率不会证明，但是好像是\(1-(\frac{1}{4})^k\)的。

代码并不难实现：

bool MR(ll n)
{
    if(n<2) return 0;//三种特判
    if(n==2||n==3||n==5||n==7||n==11||n==13||n==17||n==19) return 1;
    if(!(n%2)||!(n%3)||!(n%5)||!(n%7)||!(n%11)||!(n%13)||!(n%17)||!(n%19)) return 0;
    ll s=0,x=n-1;while(!(x&1)) x>>=1,s++;//去掉n-1含有的的质因子2
    for(int i=1;i<=8;i++)//枚举底数
    {
        if(pr[i]>n) return 1;//大于n了可以直接判为质数，因为小于n的所有质数都通过了测试
        ll r=ksm(pr[i],x,n),l=r;
        for(int j=1;j<=s;j++)
        {
            r=msc(r,r,n);
            if(r==1&&l!=1&&l!=n-1) return 0;//二次探测法
            l=r;
        }
        if(r!=1) return 0;//用费马小定理验证
    }
    return 1;
}

Pollard-Rho

这是一种可以期望在\(O(n^{\frac{1}{4}})\)的复杂度内找出n的一个非平凡因子的算法。

把它用到质因数分解上，效率肯定会高于直接暴力根号分解。

假设我们已经实现了Pollard-Rho算法，那么我们可以写出如下的solve函数去分解质因数：

void solve(ll n)
{
    if(n==1) return;
    if(MR(n)){a[++tot]=n;return;}
    ll t=n;while(t==n) t=PR(t);//非平凡
    while(n%t==0) n/=t;//去掉这个因子
    solve(t);solve(n);//分解质因数
}

然后呢我们考虑一个递推序列\(x_i=(x_{i-1}^2+c)~mod~p\)，其中\(c\)是一个随机常数，这样一直生成下去的话一定会出现一个环（太菜了不会证为什么），我们考虑该怎么寻找这个环。

假设有两个人在一个特定的轨道上面赛跑，如果赛道没有环，那么跑得快的人一定一直在最前面，如果赛道有环，那么跑得快的人一定和跑得慢的人相遇。我们考虑相对运动，可以让一个人站着不动，然后跑的快的人一直跑，并且设置一个步长k，如果跑得快的人k步之后还没有遇到环，那么就让跑得慢的人站在跑得快的人现在的位置，跑得快的人继续跑。这样最终一定可以找到一个环。

考虑怎么和寻找因子的过程结合，设跑得慢的人现在站在x上，跑得快的人现在站在y上，那么我们判断\(|y-x|\)和\(n\)的最大公约数是否大于1，如果大于1，就说明找到了一个不是1的因子，就可以返回结果了。

但是怎么判断找没找到环呢？如果\(x=y\)，那么返回的值是n，我们结合上面的solve函数，如果找到的值是n就重新进行算法，直到找到的因子不是n为止。

一个玄学的优化就是每次步长k可以\(<<=1\)，我也不知道为什么反正就是可以变得很快。

贴个代码：

ll PR(ll n)
{
    ll x=1,y=1,d=1,t=1,c=1ll*rand()*rand()*rand()*rand()%(n-1)+1;
    for(ll k=1;;k*=2,t=1,y=x)//brent找环
    {
        for(int i=1;i<=k;i++)
        {
            x=msc(x,x,n);x=(x+c)%n;//生成序列
            t=msc(t,abs(x-y),n);//累乘（卡常用）
            if(i%127==0){d=__gcd(t,n);if(d>1) break;}//卡常，判断是否找到因子了
        }
        if(d>1||(d=__gcd(t,n))>1) break;//细节，有可能是上面if跳出来的也有可能是没有达到128但是也满足条件了，都要判
    }
    return d;
}

题目

填个坑，待填