Luckyleaves's Blog

stay hungry，stay foolish greedy and lucky

Miller Robin

Posted on Views: 17

前言

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

题目传送门

前置芝士

费马小定理：

当a于p互质时，(ap)≡a(modp)

二次探测定理

内容：如果 φ(p)=p−1,p>1,p>X，且X2≡1(modp)，那么 X=1||p−1

证明：

∵X2≡1(modp)

∴ p|X2−1

∴p|(X−1)(X+1)

∵p是大于X的质数

∴ p=X+1||p≡X−1(modp) 即X=1||p−1

证毕

原理

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

首先，如果当满足(ap−1)≡a(modp)时p都为质数那这个问题是不是就迎刃而解了呢？（就可以和愚蠢的二次探测定理说再见了）

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

我们来看一下这个数据：2340≡1(mod341)，然而 341=31∗11

……

实现

  • 1.将p−1提出所有2的因数， 我们假设有t个。然后将剩余的部分定义为res（用于二次探测定理）
  • 2.枚举一个数a，并定义一个数x=ares(modp)

  • 3.如果 ∀y=x2(modp),y!=p−1 那么p就不是一个质数

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

    主函数及其他流程

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

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

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

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

    因子从何而来

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

    咳，我们来说正解。

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

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

  • 使Misplaced &则我们就找到了一个因子（至于为什么……）

  • 那如果x=y出现了呢，这就说明出现了循环，所以我们就要判环，运用倍增的思想，让y记住x，当x=x02时，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;
    }