Luckyleaves's Blog

stay hungry,stay foolish greedy and lucky

Miller Robin

前言

玄学的概率测素数, 这是一种随机的算法, 出错率大约在$($1/4$)^n$, 所以当n > 50以后正确率就接近于1了(反正比计算机本身的正确率要高)

题目传送门

前置芝士

费马小定理:

当a于p互质时,$(a^p) \equiv a (\bmod p)$

二次探测定理

内容:如果 φ(p)=p−1,p>1,p>X,且$X^2 ≡ 1(\bmod p)$,那么 $X = 1 || p - 1$

证明:

∵$X^2≡1 \pmod p$

∴ $p|X^2 - 1$

∴$p|(X - 1)(X+ 1)$

∵p是大于X的质数

∴ $p = X + 1 || p \equiv X-1\pmod p$ 即$X = 1 || p - 1$

证毕

原理

好, 那么接下来让我们走进它的原理。

首先,如果当满足$(a^p-1) \equiv a \pmod p$时p都为质数那这个问题是不是就迎刃而解了呢?(就可以和愚蠢的二次探测定理说再见了)

但愿望是美好的,现实是骨感的,当你自信地将这个代码交上去以后,你会发现:你Wa了。

我们来看一下这个数据:$2^{340} \equiv 1 \pmod {341}$,然而 341=31∗11

……

实现

  • 1.将$p - 1$提出所有2的因数, 我们假设有t个。然后将剩余的部分定义为res(用于二次探测定理)
  • 2.枚举一个数a,并定义一个数$x = a^{res} \pmod p $
  • 3.如果 $\forall y = x^2 \pmod p, y != p - 1$ 那么p就不是一个质数

  • 4.当我们的底数已经足够多了时,就可以跳出了

    主函数及其他流程

    牛B的大佬就可以去刷题了(别来找茬了

  • 我们先判断当前数是否为素数

  • 如果不是素数的话,就找它的因子

  • 递归对该因子和约出来的另一个因子进行分解(直到为质数)

    因子从何而来

    如果不嫌弃的话,我们可以一个一个试

    咳,我们来说正解。

  • 我们假设要找的因子为a

  • 我们运用随机的艺术find一个$x,y$,并不停随机$x$, 具体的法子一般是$x = x^2 + c$(c就是随机艺术的产物)

  • 使$a = gcd(y - x, n) && a \in(1, n)$则我们就找到了一个因子(至于为什么……)

  • 那如果$x = y$出现了呢,这就说明出现了循环,所以我们就要判环,运用倍增的思想,让$y$记住$x$,当$x = {x_0}^2$时,y = x,所以当x跑到y时,已经跑完一个圈。

  • 定义一个$i = 1,j = 2$,当执行$gcd$时$i ++$,如果i == j, 则$y = x, j <<= 1$

    代码

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    #include<iostream>
    #include<cstring>
    #include<cstdio>
    #include<cstdlib>
    #include<cmath>
    #include<algorithm>
    #include<ctime>
    #define LL long long
    using namespace std;
    int T;
    LL n, ans;
    LL a[30] = {2,3,5,7,11,13,17,41,61,24251};
    LL Gcd(LL x, LL y)
    {
    return y ? Gcd(y, x % y) : x;
    }
    LL mul(LL x, LL k, LL mod)
    {
    LL res = 0;
    while(k)
    {
    if(k & 1) res += x, res %= mod;
    x += x;
    x %= mod;
    k >>= 1;
    }
    return res % mod;
    }
    LL qpow(LL x, LL k, LL mod)
    {
    LL res = 1;
    while(k)
    {
    if(k & 1) res = mul(res, x, mod);
    x = mul(x, x, mod);
    k >>= 1;
    }
    return res % mod;
    }
    bool Robin(LL p)
    {
    if(p < 2) return false;
    if(p == 2) return true;
    if(p % 2 == 0) return false;
    LL d = p - 1;
    int s = 0;
    while(!(d & 1)) d >>= 1, s ++;
    for(int i = 0;i <= 9 && a[i] < p;i ++)
    {
    LL x = qpow(a[i], d, p), y = 0;
    for(int j = 1;j <= s;j ++)
    {
    y = mul(x, x, p);
    if(y == 1 && x != 1 && x != (p - 1)) return false;
    x = y;
    }
    if(y != 1) return false;
    }
    return true;
    }
    inline LL pollard(LL p)
    {
    while(1)
    {
    LL x = rand() % (p - 2) + 1;
    LL y = rand() % (p - 2) + 1;
    LL c = rand() % (p - 2) + 1;
    LL i = 0, j = 1;
    while(++ i)
    {
    x = (mul(x, x, p) + c) % p;
    if(x == y) break;
    LL d = Gcd(abs(y - x), p);
    if(d > 1) return d;
    if(i == j)
    {
    y = x;
    j <<= 1;
    }
    }
    }
    }
    void Find(LL p)
    {
    if(p == 1) return ;
    if(Robin(p))
    {
    ans = max(ans, p);
    return ;
    }
    LL k = p;
    while(p <= k) k = pollard(k);
    Find(k), Find(p / k);
    }
    void print(LL x)
    {
    if(x < 0) x = -x, putchar('-');
    if(x > 9) print(x / 10);
    putchar(x % 10 + '0');
    }
    int main()
    {
    scanf("%d", &T);
    srand(time(0));
    while(T --)
    {
    ans = 0;
    scanf("%lld", &n);
    Find(n);
    if(ans == n) puts("Prime");
    else print(ans),puts("");
    }
    return 0;
    }

我们会发现:咦怎么T了

这时就要卡常(打表)了,经过被巨佬一顿嘲讽,我们终于问到了卡常的办法

快速乘

1
2
3
4
5
6
7
8
typedef long double ld;
typedef long long ll;
typedef unsigned long long ull;
inline ull qmul(ull a,
ull b,const ull mod){
ll c=(ll)(a)*b-(ll)((ull)((ld)(a)*b/mod)*mod);
return c<0? c+mod:(c<mod? c:c-mod);
}

连简单的$\log(n)$的龟速乘都不放过……

然后我们就可以得到优秀的93分

最后的优化

93分了,冲鸭

向哪冲呢?

还记得我们的$Gcd$吗,当然,$Gcd$函数是不能再优化了,但我们调用$Gcd$的次数是可以再优化的。

  • 我们都知道, 龟速乘的模数都是质数, 但我们的模数可能并不是质数
  • 所以我们可以根据取模的性质:如果模数和被模的数都含有一个公约数$k$,那么这次模运算的结果必然也会是这个公约数$k$的倍数。所以如果我们将若干个$(y - x)$ 相乘,因为模数是$n$ ,所以如果若干个$(y - x)$中有一个与$n$有公约数,最后的结果定然也会含有这个公约数。
  • 所以我们可以多算几次$(y - x)$, 再来求$Gcd$(63次擦不多吧)
  • 记得一边欧拉$(y - x)$,一边倍增,判环哟

    AC代码

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    #include<iostream>
    #include<cstring>
    #include<cstdio>
    #include<cstdlib>
    #include<cmath>
    #include<algorithm>
    #include<ctime>
    #define LL long long
    using namespace std;
    int T;
    LL n, ans;
    LL a[30] = {2,3,5,7,11,13,17,41,61,24251};
    LL Gcd(LL x, LL y)
    {
    return y ? Gcd(y, x % y) : x;
    }
    typedef long double ld;
    typedef long long ll;
    typedef unsigned long long ull;
    inline ull mul(ull a,ull b, ull mod){
    ll c=(ll)(a)*b-(ll)((ull)((ld)(a)*b/mod)*mod);
    return c<0? c+mod:(c<mod? c:c-mod);
    }
    LL qpow(LL x, LL k, LL mod)
    {
    LL res = 1;
    while(k)
    {
    if(k & 1) res = mul(res, x, mod);
    x = mul(x, x, mod);
    k >>= 1;
    }
    return res % mod;
    }
    bool Robin(LL p)
    {
    if(p < 2) return false;
    if(p == 2) return true;
    if(p % 2 == 0) return false;
    LL d = p - 1;
    int s = 0;
    while(!(d & 1)) d >>= 1, s ++;
    for(int i = 0;i <= 9 && a[i] < p;i ++)
    {
    LL x = qpow(a[i], d, p), y = 0;
    for(int j = 1;j <= s;j ++)
    {
    y = mul(x, x, p);
    if(y == 1 && x != 1 && x != (p - 1)) return false;
    x = y;
    }
    if(y != 1) return false;
    }
    return true;
    }
    inline LL pollard(LL p)
    {
    while(1)
    {
    LL x = rand() % (p - 2) + 1;
    LL y = rand() % (p - 2) + 1;
    LL c = rand() % (p - 2) + 1;
    LL i = 0, j = 1, b = 1;
    while(++ i)
    {
    x = (mul(x * 1ull, x * 1ull, p * 1ull) + c) % p;
    b = mul(b * 1ull, abs(y - x) * 1ull, p * 1ull);
    if(!b || x == y) break;
    if(i == j || !(i % 63))
    {
    LL d = Gcd(b, p);
    if(d > 1) return d;
    if(i == j) y = x, j <<= 1;
    }
    }
    }
    }
    void Find(LL p)
    {
    if(p <= ans || p == 1) return ;//小优化
    if(Robin(p))
    {
    ans = max(ans, p);
    return ;
    }
    LL k = pollard(p);
    while(p % k == 0) p /= k;
    Find(k), Find(p);
    }
    void print(LL x)
    {
    if(x < 0) x = -x, putchar('-');
    if(x > 9) print(x / 10);
    putchar(x % 10 + '0');
    }
    int main()
    {
    scanf("%d", &T);
    srand(time(0));
    while(T --)
    {
    ans = 0;
    scanf("%lld", &n);
    Find(n);
    if(ans == n) puts("Prime");
    else print(ans),puts("");
    }
    return 0;
    }