ZigZagK的博客
[Miller-Rabin+Pollard-Rho]Codeforces1025B【Weakened Common Divisor】题解
2018年8月21日 15:09
查看标签

题目概述

有 $n$ 组数对 $(a_i,b_i)$ ,求一个数使得 $\forall i,d|a_i\lor d|b_i$ 。

解题报告

因为随便求所以找共有的素因子就好了,那么先求出 $GCD$ 表示所有 $lcm(a_i,b_i)$ 的 $gcd$ ,然后找 $GCD$ 的素因子。但 $GCD$ 有 $10^{18}$ 级别,上MR&PR就好啦!更简单的方法就是用 $GCD$ 去匹配 $a_i$ 或 $b_i$ ,如果 $(GCD,a_i)>1$ 就 $GCD=(GCD,a_i)$ 否则 $GCD=(GCD,b_i)$ 。

Miller-Rabin

Miller-Rabin(MR)素数判断是利用费马小定理 $a^{p-1}\equiv1\ (mod\ p)$ 的逆定理:如果 $a^{p-1}\equiv1\ (mod\ p)$ 则 $p$ 是素数。

然而这货是伪命题,但有很大概率正确。所以我们可以随机几个 $a$ 来验证,如果 $a^{n-1}\not\equiv1\ (mod\ n)$ 则 $n$ 不是素数。但这样并不保险……有个二次探测定理:如果 $p$ 是素数,且 $x^2\equiv 1\ (mod\ p)$ 那么 $x=1\lor x=p-1$ ,这个移下项就出来了:$p|(x+1)(x-1)$ 。假如待测数是 $n$ ,我们把 $n-1$ 表示为 $2^kt$ ,先求出 $x=a^t$ ,如果 $x=1$ 那么验证成功,否则检查 $x^2\ mod\ n$ 是不是 $1$ ,如果是 $1$ 那么检查 $x$ 是不是 $n-1$ ,如果不是就验证失败,否则验证成功,以此类推做 $k$ 次。

inline bool MR(LL n){
    for (int i=0;i<8;i++) if (!(n%prime[i])) return n==prime[i]; //可以先判掉几个素数
    for (int t=1;t<=10;t++){
        int k=0;LL x,s;for (x=n-1;!(x&1);x>>=1,k++);x=Pow(rand()%(n-2)+2,x,n);
        if (x==1) continue;
        for (;k;k--,x=s) {s=Mul(x,x,n);if (s==1) {if (x<n-1) return false;break;}}
        if (s>1) return false;
    }return true;
}

Pollard-Rho

泼辣的肉Pollard-Rho(PR)可以期望 $O(n^{1\over 4})$ 找出一个 $n$ 的因子 $d$ (感觉好像贼快啊),为了判断 $d$ 是不是素数需要利用MR。大概是这样的:如果 $(x.n)>1$ 那么 $(x,n)$ 就是 $n$ 的一个因子。

怎么找 $x$ ?随机大法好!随机出 $a$ 和 $r$ ,每次 $a'=(a^2+r)\ mod\ n$ (随机数生成器,其实可以随便搞),如果 $(|a-a'|,n)>1$ 说明找到了。但是 $a$ 可能会出现循环,所以需要神奇的Floyd判环。

干嘛这么做,复杂度又是什么?令生成的随机数列为 $\{a_i\}$ ,另一个数列 $\{b_i=a_i\ mod\ d\}$ ( $d$ 是找到的因子),如果 $b_i=b_j\Leftrightarrow d|(a_j-a_i)$ 就说明找到了一组解。跟据生日悖论,$a$ 期望 $\sqrt n$ 出现循环,而 $b$ 期望 $\sqrt d$ 出现循环,所以 $b$ 更容易出现循环,只要在 $a$ 出现循环之前 $b$ 出现循环就意味着找到了一个因子。当然没找到就多做几次啊2333。

inline LL FindD(LL n){
    LL x=rand()%(n-1)+1,y=x,r=rand()%(n-1)+1;
    for (int i=1,k=1;;i++){
        AMOD(x=Mul(x,x,n),r,n-1);x++;if (x==y) return -1;LL GCD=gcd(Abs(x-y),n);
        if (GCD>1) return GCD;if (i==k) k<<=1,y=x;
    }
}
void PR(LL n){
    if (n==1) return;if (MR(n)) {p[++P]=n;return;} //特判
    LL D;for (D=FindD(n);D<0;D=FindD(n));PR(D);PR(n/D);
}

示例程序

#include<cstdio>
#include<algorithm>
#include<map>
using namespace std;
typedef long long LL;typedef long double DB;
const int maxn=150000,maxp=10000,prime[]={2,3,5,7,11,13,17,19};

int n,a[maxn+5],b[maxn+5];LL LCM[maxn+5];int P;LL p[maxp+5];

inline LL Mul(LL x,LL y,LL MOD) {return ((x*y-(LL)((DB)x/MOD*y)*MOD)%MOD+MOD)%MOD;}
inline LL Pow(LL w,LL b,LL MOD) {LL s=1;for (;b;b>>=1,w=Mul(w,w,MOD)) if (b&1) s=Mul(s,w,MOD);return s;}
inline bool MR(LL n){
    for (int i=0;i<8;i++) if (!(n%prime[i])) return n==prime[i];
    for (int t=1;t<=10;t++){
        int k=0;LL x,s;for (x=n-1;!(x&1);x>>=1,k++);x=Pow(rand()%(n-2)+2,x,n);
        if (x==1) continue;
        for (;k;k--,x=s) {s=Mul(x,x,n);if (s==1) {if (x<n-1) return false;break;}}
        if (s>1) return false;
    }return true;
}
#define Abs(x) ((x)<0?-(x):(x))
LL gcd(LL a,LL b) {return b?gcd(b,a%b):a;}
inline LL lcm(LL a,LL b) {return a/gcd(a,b)*b;}
inline void AMOD(LL &x,LL tem,LL MOD) {if ((x+=tem)>=MOD) x-=MOD;}
inline LL FindD(LL n){
    LL x=rand()%(n-1)+1,y=x,r=rand()%(n-1)+1;
    for (int i=1,k=1;;i++){
        AMOD(x=Mul(x,x,n),r,n-1);x++;if (x==y) return -1;LL GCD=gcd(Abs(x-y),n);
        if (GCD>1) return GCD;if (i==k) k<<=1,y=x;
    }
}
void PR(LL n){
    if (n==1) return;if (MR(n)) {p[++P]=n;return;}
    LL D;for (D=FindD(n);D<0;D=FindD(n));PR(D);PR(n/D);
}
int main(){
    freopen("program.in","r",stdin);freopen("program.out","w",stdout);
    scanf("%d",&n);for (int i=1;i<=n;i++) scanf("%d%d",&a[i],&b[i]),LCM[i]=lcm(a[i],b[i]);
    LL GCD=LCM[1];for (int i=2;i<=n;i++) GCD=gcd(GCD,LCM[i]);if (GCD==1) return puts("-1"),0;
    PR(GCD);return printf("%lld\n",p[1]),0;
}
版权声明:本博客所有文章除特别声明外,均采用 CC BY 4.0 CN协议 许可协议。转载请注明出处!