Mr素性测试+Prho分解小记

· · 个人记录

一些玄学的东西……

Mr素性测试·前置芝士

(下面的p都指质数)

费马小定理:a^{p-1}=1(mod\ \ p)

那么反过来,a^{p-1}≠1(mod\ \ p)的数一定不是素数。

二次探测定理:

理解的话,大概是$1$在膜质数意义下开方一定得到$1$或者$p-1$(即$-1$) **爆LL的乘法**: - 可以使用所谓龟速乘,把乘法分治成加法,复杂度$O(logn)$,**跑的很慢**,主要是多个$log$让复杂度很不美观。 - `__int128`带走,实测常数极大,而且有一些OJ/比赛不支持。 - 使用下面的魔法: (~~好像被大佬Hack了,慎用~~) ```cpp inline ll guisuMul(ll a,ll b,ll x){ return (a*b-(ll)((long double)a*b/x+0.1)*x)%x; } ``` 具体原理百度“骆可强 快速乘”即可,这里就不再赘述了。 ## Mr素性测试·干货 Mr素性测试**效果** : 不用提前筛,$O(log^2n)$就能测出是否是素数。 设待测的数为$n$。 (先暴力试除掉一些因子,可以避免一些边界情况) 首先,你可以找一个数$a$,然后计算$a^{n-1}$,假如不等于1的话,则判定为非素数。 可以想见:多随机几个$a$,如果都`hack`失败,那么基本可以确定$n$就是质数。 这里有一个有趣的结论 : > 对于合数$m$,满足$a^{m-1}=1(mod\ \ m)$的$a$称为坏的,因为它不能帮助我们`hack`这个$m

首先,在模m意义下,所有数的m-1次方形成一个乘法群Z_*^m。封闭性 : a^{m-1}b^{m-1}=(ab)^{m-1}

那么,所有不好的数的m-1次方也形成一个乘法群。封闭性 : a^{m-1}b^{m-1}=1,仍然在群里。

那么,如果我们存在一个好的数,就说明不是所有数都是坏的,那么不好的群就是Z_*^m的真子群,其大小是|Z_*^m|的约数(除去其本身)

那么坏的数至多占|Z_*^m|这个群的\frac{1}{2}。

这就是Fermat素性测试。

问题是,存在一些(较多)神奇的合数,能通过任意的a的测试。

为了保证正确率,我们考虑利用二次探测定理:

令(m-1)=k*2^t,我们把m-1的所有因子2都提取出来(相当于进行一系列开方)

首先计算a^k\mod m,假如得到1或n-1,hack失败。

否则,这只可能是n-1开方得来的,不断平方看看能否变为n-1即可。

这玩意利用二次剩余理论证明一下也可以得到每次测试概率大概是\frac{1}{2},同样也存在大量反例。

但是,把两个定理相配合,相当于把两个反例集合交集,基本变成了空集……

这样子错误率就大大减小了,对于一个合数,可以证明一次测试成功hack的概率约为0.75。

那么试个O(logn)次就基本万无一失了。

实战中,选择a=\{2,3,5,7,13,19,61,233\}基本上10^{18}次方以内都没问题了。

Code:

并没有使用快速乘。

#include<algorithm>
#include<cstdio>
#define ll long long
using namespace std;
ll powM(ll a,ll t,ll m)
{
  ll ans=1;
  while(t){
      if (t&1)ans=ans*a%m;
      a=a*a%m;
      t>>=1;
  }return ans; 
}
bool mr(ll n,ll a)
{
  ll t=n-1;
  while(!(t&1))t>>=1;
  ll buf=powM(a,t,n);
  if (buf==1||buf==n-1)return 0;
  while((t<<=1)<n-1){
    buf=buf*buf%n;
    if (buf==n-1)return 0;
  }return 1;
}
const int testp[8]={2,3,5,7,13,19,61,233};
bool ptest(ll n)
{
  if (n<2)return 0;
  for (int i=2;i*i<=min(n,10000ll);i++)
    if (n%i==0)return 0;
  if (n<=10000)return 1;
  for (int i=0;i<8;i++)
    if (mr(n,testp[i]))return 0;
  return 1;
}
int n;
int main()
{
  scanf("%d%d",&n,&n);
  for (int i=1,a;i<=n;i++){
      scanf("%d",&a);
    if (ptest(a))puts("Yes");
      else puts("No");
  }return 0;
}

Prho分解·问题引入

P4718 【模板】Pollard-Rho算法

大概就是给你一个long long范围内的数,要求你分解。

暴力试除法O(\sqrt{n})肯定没戏的……

Pro分解效果 : O(n^{0.25})把一个数质因数分解。

Prho分解·前置芝士

Mr素性测试(两者关系并不十分紧密)

生日悖论(怪论):

随便找23个人,他们之中有任意两人生日相同的概率已经大于50\%。

也就是 : 随机写[1,N]之内的整数,期望写到第\sqrt{N}个时出现重复。

Prho分解·干货

我们先考虑通过某种方法得到n的某个非平凡因数(不为1和它本身)。

拿Mr特判掉n是素数的情况。

我们假设n有一个因子q,那么,我们如果能找到两个数:

$gcd(|x1-x2|,n)$就能得到$n$的一个非平凡因子。 问题在于怎么得到这个东西…… 考虑一个数列$f_0=1;\ \ f_i=F(f_{i-1}) \mod n$,这里值域为$[0,n-1]$的整数。 其中$F$是某种单变量映射,结果和只和参数相关(比如$F(x)=x^2+t$)。 如果$F$的映射方式比较随机,那么根据生日悖论,这个数列在$O(\sqrt{n})$步内就会出现**重复**(环)。 我们定义$g_i=f_i\mod q$,当然,我们**不能直接求出$g$来**,不过这可以帮助我们分析。 由于$f$近似随机,可以想见$g$也是近似随机的。 所以根据生日悖论,$g$在$O(\sqrt{q})$步内,就会有**重复**,设两个$g$相同的位置为$i,j$: 给个图让大家看看 $\large\rho$(rho) 究竟是什么意思 : ![](https://cdn.luogu.com.cn/upload/image_hosting/ce0jnf1f.png) 也就对应着$f_i=f_j(mod\ \ q)$,如果此时恰巧$f_i≠f_j$,那么我们利用$|f_i-f_j|$就找到了一个因子。 如果很不幸$f_i=f_j$,那么分解失败。 事实上,这种情况概率极小,因为$g$的环期望大小远小于$n$的,两个$g_i$相同时,往往在$(mod\ \ n)$下还没有绕到环上。 分解失败的概率趋近于$n/q=O(\frac{1}{\sqrt{n}})$,从hash的角度容易理解,此时只需要略微改动$F(x)$(换个种子)然后重新分解。这部分复杂度可不计。 问题在于我们并不能显式地求出$g$,但是如果我们已经获得了答案(如上图左侧),我们能够通过$gcd$来判定。 实战中我们令$F(x)=x^2+t$,即平方之后加上常数。这个东西生成出来取模之后,可以认为近似随机。 现在要找环,我们有几种方法: - **追及法** 找两个人在$f$上面跳,第一个人速度为1,第二个人速度为2,那么如果$g$有环,前者一定会被后者追 上。 形式化的说,我们每次取$f_i$和$f_{2i}$就好了。 我们如何判定$f_i=f_{2i}(mod\ \ q)$(也即追上?) 如果$f_i=f_{2i}=0$则表示失败了,否则,计算$gcd(|f_i-f_{2i}|,n)$,如果**大于**$1$,则成功。 根据上文,一次测试只需要消耗期望$\sqrt{q}$次递推。 - **倍增法** 每次取$\large f_{2^k}$,然后让$\large f_{(2^k,2^{k+1}]}$与其一一判定。 同样只需要消耗期望$\sqrt{q}$次递推。 如果很不幸分解失败,我们只需要重置上文$F(x)$中的$t$就好了。 由$q\leq \sqrt{n}\Rightarrow \sqrt{q}\leq n^{0.25}$,所以总的复杂度是$O(n^{0.25}logn)

每步都要求gcd代价有点高,这也导致复杂度多个log。

注意到,我们求出的gcd一旦不为1就可以返回了。

又注意到:gcd(a,n)>1||gcd(b,n)>1\Leftrightarrow gcd(ab\mod n,n)>1

我们可以把连续lim个要测试的值乘起来膜n,然后求与n的gcd。

假如大于1,那么这一段内就找到了解,回头一个一个跳就好了。

取lim=O(logn)(实战中lim=128)就可以优化掉一个log了,好开心呢!

通过上述算法,我们就能得到n的一个非平凡因子d。

然后我们继续分治分解(n/d),d,直到产生素数为止。

通常先暴力试除一些素数来减小常数,避开边界。

Code:

追及法:

#include<algorithm>
#include<cstdio>
#define ll long long
using namespace std;
inline ll mul(ll a,ll b,ll m){
  ll r=((long double)a/m*b+0.5);
  r=a*b-r*m;
  return r<0?r+m:r;
}
ll powM(ll a,ll t,ll m)
{
  ll ans=1;
  while(t){
      if (t&1)ans=mul(ans,a,m);
      a=mul(a,a,m);
      t>>=1;
  }return ans; 
}
bool mr(ll n,ll a)
{
  ll t=n-1;
  while(!(t&1))t>>=1;
  ll buf=powM(a,t,n);
  if (buf==1||buf==n-1)return 0;
  while((t<<=1)<n-1){
    buf=mul(buf,buf,n);
    if (buf==n-1)return 0;
  }return 1;
}
const int testp[8]={2,3,5,7,13,19,61,233};
bool ptest(ll n)
{
  if (n<2)return 0;
  for (int i=2;i*i<=min(n,10000ll);i++)
    if (n%i==0)return 0;
  if (n<=10000)return 1;
  for (int i=0;i<8;i++)
    if (mr(n,testp[i]))return 0;
  return 1;
}
inline ll gcd(ll a,ll b){
  if(a==0)return b;
  if(a<0)a=-a;
  ll t;
  while(b){t=a%b;a=b;b=t;}
  return a;
}
ll sav[150];int tot;
const int lim=128;
ll prho(ll n,ll c)
{
  ll x1=(c+1)%n,x2=(mul(x1,x1,n)+c)%n,buf;
  tot=0;
  while(1){
    buf=mul(x1-x2,buf,n);
    sav[++tot]=x1-x2;
    if (tot==lim){
      if (gcd(buf,n)>1)break;
      tot=0;
    }
    x1=mul(x1,x1,n)+c;
    x2=mul(x2,x2,n)+c;
    x2=mul(x2,x2,n)+c;
  }
  for (int i=1;i<=tot;i++){
    buf=gcd(sav[i],n);
    if (buf>1)return buf;
  }return n;
}
long long ans;
void solve(ll n)
{
  if (n<=ans)return ;
  if (ptest(n)){ans=max(ans,n);return ;}
  ll sav=prho(n,rand()%(n-1)+1);
  while(sav==n)sav=prho(n,rand()%(n-1)+1);
  solve(sav);solve(n/sav);
}
int main()
{
  srand(233);
  int T;scanf("%d",&T);
  for (int qt=1;qt<=T;qt++){
    ll a,sav;
      scanf("%lld",&a);
      if (ptest(a)){
      puts("Prime");
      continue;
    }sav=a;ans=0;
      for (int i=2;i*i<min(10000ll,a);i++)
        if (a%i==0){
          ans=i;
          while(a%i==0)a/=i;
      }
      solve(a);
    printf("%lld\n",ans);
  }return 0;
}

倍增法:

ll prho(ll n,ll c)
{
  tot=0;
  ll x,y=0,mt=1;
  for (int st=1;;st<<=1){
    x=y;
    for (int i=0;i<st;i++){
      y=mul(y,y,n)+c;
      mt=mul(y-x,mt,n);
      sav[++tot]=y-x;
      if (tot==lim){
        if (gcd(mt,n)>1)break;
        tot=0;
      }
    }if (tot==lim)break;
  }
  for (int i=1;i<=tot;i++){
    mt=gcd(sav[i],n);
    if (mt>1)return mt;
  }return n;
}