各种友(e)善(xin)数论总集,从入门到绝望2

目录

因为原来的那篇已经很多了,所以在此写上第二篇。

这一章可以说是紧紧围绕的素数的主旨展开的。

@

前置芝士

二进制的GCD

二进制的GCD用的是更相减损法。

首先,我们有两个数字\(x,y\)。

  1. 当\(x\%2==0,y\%2==0,gcd(x,y)=gcd(x>>1,y>>1)*2\)
  2. 当\(x\%2==0,y\%2==1,gcd(x,y)=gcd(x>>1,y)\)
  3. 当\(x\%2==1,y\%2==0,gcd(x,y)=gcd(x,y>>1)\)
  4. 当\(x\%2==1,y\%2==1,x>=y,gcd(x,y)=gcd((x-y)>>1,y)\),其实这里的\(x-y\)不用除\(2\)也可以,只是这样写更快,反正都是偶数。
inline  LL  gcd(LL  x,LL  y)
{
    int  ans=0;
    while(x  &&  y)
    {
        if(x&1  &&  y&1)
        {
            y>x?x^=y^=x^=y:0;
            x=(x-y)>>1;
        }
        else  if(x&1)y>>=1;
        else  if(y&1)x>>=1;
        else  x>>=1,y>>=1,ans++;
    }
    return  (x+y)<<ans;
}

常数好像比GCD更小,复杂度都是\(log\),但是貌似跑起来差不多,写写吧,反正都差不多了,还更稳一点,你说是吧。

快速乘

博主博主,平常\(O(1)\)都已经如此之快,难道可以\(O(0)\)?

不不不,都一样,只不过算的是\(x*y\%z\),因为有时候\(x*y\)溢出了long long,但是结果并没有,所以发明了快速乘。

O(1)

\(O(1)\)版的非常简单。\(x*y\%z=x*y-(x*y/z)*z\)(在C++环境下),但是这个有什么特殊的吗?

就是用溢出对待溢出,我们先用long double(16位)得出(x*y/z)(你在比赛的时候,可以先用\(sizeof(long\) \(double)\)得出你电脑的long double是几位的,如果不是\(16\)位的话,那评测机应该也是,为了稳妥还是用\(log\)的吧)。

然后我们乘一下,两边可能会溢出,但是我们还是能减出正确的结果。

另外还有一坨精度问题,模数大的话还是少用吧。

inline  LL  ksc(LL  x,LL  y,LL  z)
{
    LL  c=(LD)x*y/z+0.5;
    LL  ans=x*y-c*z;
    return  ans<0?ans+z:ans;
}

O(log)

有没有什么稳得一批又好用的快速乘?

当然后,假设又有个\(x,y,z\)。

我们把\(x\)拆成二进制:\(c_{n}*2^{n}+c_{n-1}*2^{n-1}+....\),而\(c\)的取值只能为\(0,1\)
然后\(x*y\)在乘法分配率一下:\(y*c_{0}*2^{0}+y*c_{1}*2^{1}+...\),那么我们只要边乘边模不就好起来了吗。

inline ll ksc(ll x,ll y,ll p){//计算x乘y的积
    ll res=0;//加法初始化
    while(y){
        if(y&1)res=(res+x)%p;//模仿二进制
        x=(x<<1)%p; y>>=1;//将x不断乘2达到二进制
    }return res;
}
// ll 表示 long long
//这里的代码用的是https://www.cnblogs.com/812-xiao-wen/p/10543023.html

floyd提出的判环法

判环

你以为是floyd,不是,是这样的,假设一个链表有环,怎么判环,我们这样想:\(y\)以\(x\)两倍的速度奔跑,那么当\(x,y\)相遇时,\(y\)刚好跑完几圈了,就退出。

找环

其实还有个扩展,如何找到环的起始位置。

我们设链表头走到环开始的地方步数为\(m\),从链表头走到相遇地点的步数是\(m+k\),然后环的长度为\(n\),\(x,y\)分别走了\(X,Y\)圈。

那么\(S_{x}=m+k+Xn,S_{y}=m+k+Yn\),然后又因为\(S_{y}=2S_{x}\),所以\(S_{x}=(Y-X)n\)。
也就是说两人走的距离肯定是\(n\)的倍数。

然后我们再把\(x\)提到了链表开始的地方,让两个人继续开始走,当走了\(m\)步(\(x\)在环开始的地方),那么因为说了\(S_{y}=2S_{x}=2(Y-X)n\),也就是说\(y\)走的步数应该是\(n\)的倍数,也就是说当第一次相遇的时候,他应该离走完这个环到环开始的地方为\(m\),所以\(y\)也到了环开始的地方,所以\(x,y\)将会在环开始的地方相遇。

可惜这里不用。

生日悖论

首先我们来看看,现在我们来选数字,在1-100之间,如果我们选一个数字,那么是1的概率则是\(\frac{1}{100}\),但是如果我们选两个数字,然后取差的绝对值,会怎样?我们选一个数字,然后选到他周围的两个数字的概率就变成了\(\frac{1}{50}\)了!(忽略第一个数字\(1\)和\(100\)的情况,那还是\(\frac{1}{100}\)),难道多元能增加概率!

没错,这就是生日悖论的内容。

生日悖论的重要思想是什么,\(1-x\)的范围,如果有\(\sqrt{x}\)个数字的话,重复的概率就会高达\(50\%\),恐不恐怖。

至于证明,在这里贴上大佬的证明。
各种友(e)善(xin)数论总集,从入门到绝望2

Miller_rabin

相信在第一章里面,你们已经学会了费马小定律了,那就不讲了QMQ。

前言

我们都知道判断一个数字是不是素数,有一种方法就是试除法,直接从\(2\)枚举到\(\sqrt{p}\),但是有没有一种方法,能比\(O(\sqrt{p})\)还快有准确无误呢?

答案是并没有,但是如果你要求的是很大概率的话,打我可以告诉你的是,Miller_rabin就是这么一种算法,基本上准确无误,就连强伪素数都能跑过去,是什么呢?

二次探测

我们都知道,选取一个\(p\)数字,然后用费马小定理判断一下,如果不是\(1\),那就不是素数,但是存在这么一种数字,能满足费马小定理但是不是素数的一类数字,我们又要怎么判断呢?

这里就要引入一个定理了,这个定理可以很大概率的判断是不是素数,加上费马小定理。

如果\(p\)是质数且\(a^2≡1(\mod p)(a<p)\),那么\(a=1,p-1\)。
我们可以来证明一下:
\[a^2≡1(\mod p)\]
\[a^2-1≡(\mod p)\]
\[(a+1)(a-1)≡0(\mod p)\]
那么,因为\(a<p\)且\(p\)是质数,所以\(a=1,p-1\)。

那么我们就可以把\(p-1\)分成\(2^{t}*k\),然后随机选取一个值\(x\),然后计算\(x^{k}\),然后继续不断取平方:\(x^{2^{i}*k}\),然后不断的用二次探测来检测,更重要的是我们最后还可以用用费马小定理,当然,\(x\)我们可以手动取几个素数来多判几次,不知道为什么,素数成功概率大一点QMQ。

很明显是\(log\)的。

代码

inline  LL  ksc(LL  x,LL  y,LL  z)
{
    LL  c=(LD)x*y/z+0.5;
    LL  ans=x*y-c*z;
    return  ans<0?ans+z:ans;
}
inline  LL  ksm(LL  x,LL  m,LL  mod)
{
    if(m==0)return  1%mod;
    LL  ans=1;
    while(m>1)
    {
        m&1?ans=ksc(ans,x,mod):0;
        x=ksc(x,x,mod);m>>=1;
    }
    return  ksc(ans,x,mod);
}
inline  int  log2(LL  &x)
{
    int  ans=0;
    while(x%2==0)ans++,x>>=1;
    return  ans;
}
int  su[]={2,3,5,7,11,23,29,61};
inline  bool  pd(LL  x)//判断一个素数 
{
    for(int  i=0;i<=7;i++)
    {
        if(x<=su[i])return  1;
        LL  y=x-1;int  tt=log2(y);
        y=ksm(su[i],y,x);
        while(tt--)
        {
            LL  z=ksc(y,y,x);
            if(z==1  &&  y!=1  &&  y!=x-1)return  0;
            y=z;
        }
        if(y!=1)return  0;
    }
    return  1;
}

Pollard-Rho

前言

你是否想快速的分解一个素数?

想吗?少年。

---来自SaDiao博主的一席话。

例题

都看到了,就是想咯QMQ。

例题

假设我们要分的数字是\(p\)

构建随机数列

我们构建一个随机函数\(x_{i}=x_{i-1}*x_{i-1}+C(\mod p)\)(\(C\)为我们自己定的常数),这么强?

同时\(x_{1}=2\),我们发现这个序列每个都跟前面的数字有关,那么不就是类似链表了吗?而且因为是模了后序列,所以会有环(根据生日悖论,出现相同的数字概率是\(O(\sqrt{p})\)),就可以用\(floyd\)判环法了。

我们再设一个函数:\(y_{i}=x_{i-1}*x_{i-1}+C(\mod q)\)(\(q\)为\(p\)的一个质因数),\(y_{1}=2\),那么这个序列也会出现循环节的,我们再在两个循环节上找到两个位置\(i,j\),使得\(i<j,l(i)=l(j)\),然后我们就会发现\(|x_{i}-x_{j}|\)含有\(q\),也就是\(gcd(|x_{i}-x_{j}|,p)≠1\),那么我们就可以找到一个素数了。

但是我们并不知道\(q\),我们又怎么找到\(i,j\)呢,我们会发现其实\(i,j\)就是\(y\)数列循环节上的对应位置,而上面也只是提供了一种可行性,也就是说我们可以用\(floyd\)找环法来在\(x\)数列中找,如果\(gcd\)为\(1\),那么继续找环,如果\(gcd\)为\(p\)(差为0就会这样),说明我们找到了模数为\(q\)的环,可惜也是模数为\(p\)的环,那么我们就退出,然后\(C++\),如果两个都不是,那么我们不就找到了一个质因数了吗?

再看看概率有多大,原本我们找两个数字的差找到质因数的概率应该比较小,但是如果我们是\(gcd≠1\)的话,那概率不就大了吗?而且期望的循环节大小为\(\sqrt{p}\),不就好起来了吗?

而且判断一个数字是不是素数就靠Miller_rabin了。

优(ka)化(chang)

我们要发现一个事情:\(x\%z\),可以等同于\(z*((x/gcd(x,z))\%(z/gcd(x,z)))\)。

那么这有什么用呢?我们可以把几个差乘起来,然后模一下(如果出现了质因子的话不会因为模了而消失掉,上面写了),至于乘几次,我选择的是\(127\),当然,\(pow(p,0.1)\)也有人用,不要太大就可以了,我们后面则叫乘了\(step\)次。

那么我们也是乘完后GCD。

  1. 如果为\(1\)继续。
  2. 如果为\(n\)的话,说不定中间有质因子呢,我们也回到\(step\)次一起,不乘起来,一次次慢慢来。
  3. 如果两个都不是,恭喜,喜提质因子一枚。

然后自我感觉良好,优化了不少的常数。

注意事项

我们会发现,打完代码有事后还是会卡住的。

为什么,因为\(4\)是个神奇的数字,我们可以把\(C\)从\(1\)到\(4\)枚举一遍,会发现差统统不会涉及到\(2\),而其他\(2\)的次方,比如\(2^{i}\),在\(C=2^{i}-4\)的时候,肯定跳一次就能得到结果,\(x=0,y=C\),然后就可以把\(2\)筛出来,为什么\(4\)不行,因为\(C=0\)了,所以我们就只会一遍遍的得到\(4\)然后重来。

那我们特判\(4\)不就行了?不不不,特判要讲究艺术。

你想想,我们要是特判\(\%2==0\)不是一样的吗,而且还造福了其他的数字,尤其是\(2\)的次方,不加这个很有可能就老是到\(2^{i}-4\)才跳出来。

时间复杂度证明

一次的Pollard-Rho的复杂度是多少?(这道题目得用几次Pollard-Rho)

设\(N=A*B(A<=B)\),那么\(A<=\sqrt(N)\),我们是在\(x\)找\(A\)的循环节,期望复杂度为\(O(\sqrt{A})\),那么不就是\(O(N^{\frac{1}{4}})\)吗?

但是其实不然,因为加上GCD什么乱起八糟的,正宗的应该是:
\(O(\frac{N^{\frac{1}{4}}logN}{step}+logN)\),但原本就是玄学算法你加这么多干嘛,而且在long long范围内我们的\(step\)肯定大于\(logN\),毕竟我们的\(step\)原本就是\(127\)吗。

所以我们还是写\(O(N^{\frac{1}{4}})\),好看又好写。

代码

开了O2在luogu跑了600+ms,快的飞起,不开也有1.44s了,这不就快的飞起了吗。

可以试试这个跑不跑得出结果46856248255981。

强伪素数呀,都跑过去了。

#include<cstdio>
#include<cstring>
#include<cmath>
using  namespace  std;
typedef  long  double  LD;
typedef  long  long  LL;
inline  LL  zabs(LL  x){return  x<0?-x:x;}
inline  LL  mymin(LL  x,LL  y){return  x<y?x:y;}
inline  LL  mymax(LL  x,LL  y){return  x>y?x:y;}
inline  LL  ksc(LL  x,LL  y,LL  z)
{
    LL  c=(LD)x*y/z+0.5;
    LL  ans=x*y-c*z;
    return  ans<0?ans+z:ans;
}
inline  LL  ksm(LL  x,LL  m,LL  mod)
{
    if(m==0)return  1%mod;
    LL  ans=1;
    while(m>1)
    {
        m&1?ans=ksc(ans,x,mod):0;
        x=ksc(x,x,mod);m>>=1;
    }
    return  ksc(ans,x,mod);
}
inline  int  log2(LL  &x)
{
    int  ans=0;
    while(x%2==0)ans++,x>>=1;
    return  ans;
}
int  su[]={2,3,5,7,11,23,29,61};
inline  bool  pd(LL  x)//判断一个素数 
{
    for(int  i=0;i<=7;i++)
    {
        if(x<=su[i])return  1;
        LL  y=x-1;int  tt=log2(y);
        y=ksm(su[i],y,x);
        while(tt--)
        {
            LL  z=ksc(y,y,x);
            if(z==1  &&  y!=1  &&  y!=x-1)return  0;
            y=z;
        }
        if(y!=1)return  0;
    }
    return  1;
}
inline  LL  gcd(LL  x,LL  y)//实测二进制版GCD只比原来的快了20+ms,估计是因为优化减少了GCD的调用次数,凸显不出优势。
{
    int  ans=0;
    while(x  &&  y)
    {
        if(x&1  &&  y&1)
        {
            y>x?x^=y^=x^=y:0;
            x=(x-y)>>1;
        }
        else  if(x&1)y>>=1;
        else  if(y&1)x>>=1;
        else  x>>=1,y>>=1,ans++;
    }
    return  (x+y)<<ans;
}
inline  LL  Pol(LL  now,LL  step,LL  add)
{
    if(now%2==0)return  2;//防止毒瘤的4的情况 
    LL  x=2,y=2,d=1;
    while(1)
    {
        LL  tx=x,ty=y;
        for(int  i=1;i<=step;i++)
        {
            x=ksc(x,x,now)+add;x>=now?x-=now:0;
            y=ksc(y,y,now)+add;y>=now?y-=now:0;
            y=ksc(y,y,now)+add;y>=now?y-=now:0;
            d=ksc(d,zabs(x-y),now);
        }
        d=gcd(d,now);
        if(d==1)continue;
        else  if(d!=now)return  d;
        x=tx;y=ty;
        for(int  i=1;i<=step;i++)
        {
            x=ksc(x,x,now)+add;x>=now?x-=now:0;
            y=ksc(y,y,now)+add;y>=now?y-=now:0;
            y=ksc(y,y,now)+add;y>=now?y-=now:0;
            d=gcd(zabs(x-y),now);
            if(d!=1)return  d%now;
        }
    }
}
inline  LL  work(LL  n)
{
    if(pd(n)  ||  n==1)return  n;
    LL  tmp=0,step=127/*玄学步数*/,add=1;
    while(!tmp)tmp=Pol(n,step,add++);
    //
    if(n/tmp<tmp)tmp=n/tmp;//使得n/tmp>=tmp
    LL  ans=work(n/tmp);
    if(ans>=tmp)return  ans;
    return  mymax(ans,work(tmp));
    //实现上的一个优化,优化空间小,但是能优化,而且不会耗多少空间,基本正优化
}
LL  n;
int  main()
{
    int  T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld",&n);
        LL  ans=work(n);
        if(ans==n)printf("Prime\n");
        else  printf("%lld\n",ans);
    } 
    return  0;
}
上一篇:CQNK信心赛(2019.8.20)


下一篇:xin-小太阳