文档章节

Miller-Rabin素性测试|Pollard's Rho算法

o
 osc_1ee7cxmx
发布于 2018/08/06 15:46
字数 1703
阅读 10
收藏 0

精选30+云产品,助力企业轻松上云!>>>

Miller-Rabin 素性测试

Miller-Rabin 素数测试

一本通上的M-R不透彻啊~

Miller-Rabin是利用随机化算法判断一个数是合数还是素数。

首先,如果一个数N是素数,那么他一定满足费马小定理。

$a^{N-1}\equiv1\pmod N$

我们可以任取数字a,计算这个式子的值来判断N是否为素数。

但是这么做不靠谱啊,有很多合数会被卡~

我们介绍一个相关的引理。

当p是素数且p大于2时,$1\bmod p$的平方根只有1和-1。

证明:

假设x是$1\bmod p$的平方根,则有$x^2\equiv1\pmod p$

则有$(x+1)(x-1)\equiv 0\pmod p$

也就是说$(x+1)(x-1)$能够整除素数p。根据欧几里得引理,(x+1)或(x-1)整除p,也就是$x\equiv1\pmod p$和$x\equiv-1\pmod p$

我们假设n是一个素数(n>2),则n-1是一个偶数,且能被表示成$2^s*d$的形式,d是奇数。

那么对于任意的a和$0\le r\le s-1$,有$a^d\equiv1\pmod n$和$a^{2^rd}\equiv -1\pmod n$。

为什么呢?因为费马小定理。对于一个素数n,有$a^{n-1}\equiv1\pmod n$,由于上面的引理,如果我们不断对$a^{n-1}$取平方根,我们总会得到1或-1。如果我们得到了-1,则$a^{2^rd}\equiv -1\pmod n$成立,判定p是素数。如果我们从未得到-1,通过这个过程我们取遍了2的所有幂次,则$a^d\equiv1\pmod n$成立。

Miller-Rabin素性测试就在于上述原理的逆否,如果我们找到了一个a,使得对于任意的$0\le r\le s-1$,满足以下两个式子:$a^d\not\equiv1\pmod n$,$a^{2^rd}\not\equiv-1\pmod n$,n就不是一个素数。

如果一个数a通过了上述测试,就称a是n是合数的凭证。否则称a是n是素数的强伪证,也就是他可能出锅。(毕竟是随机化算法嘛)我们就需要多找几个a,使得被误判的概率更小。

Pollard's Rho算法

给出一个数N,快速提取N的**[一个]**因数。

生日悖论

在一群学生中,随机选一个学生,则他的生日为4 月 1 号的概率为$\frac 1 {365}$。

转化为古典概型,在[1,365]中随机选取一个数,该数为90的概率为$\frac 1 {365}$

然后我们现在需要随机选取2个人,他们两个人生日相同的概率也是$\frac 1{365}$

然而我们现在随机选取3个人,他们中任意两个人生日相同的概率是:

我们利用单步容斥,任意两人生日相同的概率=1-没有两人生日相同的概率。

没有两人生日相同的概率为$\frac{365364363}{365365365}$,约为0.9917958341152186(使用Windows7的计算器算出),则任意两人的生日相同概率为0.0082041658847814(同上)。

我们现在随机选取k个人,则容易得到的是,任意两个人生日相同的概率为:

$1-\frac{365!}{(365-k)!*365^k}$(k<=365)。

当k>365,P=1

当k=23时,P>50%,大概是$\sqrt N$

当k=42时,P>90%

当k=100是,P=99.99996928%

这个很大了~

因数分解

我们选取k个数,询问是否存在$x_i-x_j$能够整除N

当$k=\sqrt N$,可能性提高到了50%

我们将可能性从$\frac 1 N$提高到$\sqrt{\frac1 N}$。

对于一个$10^{10}$的大整数,我们只需要做$10^5$次比较,除法。这还是炸了。。。

我们选取k个数,询问是否存在$gcd(x_i-x_y,N)>1$。看着是棒了好多了,

下面我们介绍P-R

Pollard's Rho算法

这个Pollard's Rho算法只将两个数放在内存中。我们不一下子随机生成k个数,反而一个一个生成检查相邻的2个数。反复执行这个步骤。

我们有一个Interesting的函数$f(x)=(x^2+a)\mod N$。这个a可以自己指定,也可以rand(),这个不是重点。

我们从$x_1=2$开始,让$x_2=f(x_1)$....$x_{i+1}=f(x_i)$。如果相邻的两个数的差与N的gcd大于1,返回他们的gcd。作为N的因子。但是这个可能会出现环(详见模域平方)。

有一个Floyed发明的神奇的算法。小明和小红在操场上跑圈,小明的速度恰好是小红的2倍,最后小明跑2圈时,小红恰好跑了1圈,小明从小红的后面追上小红,但是实际上小明比小红快整整一圈。

我们用floyed判环,我们如果失败了我们只需要重新选择一个a即可。

这就是P R算法。

upd 9.6

然后洛谷板子题真心坑人各种卡常,瞬间变身卡同学

然后这是代码(抄了最优解Rank1好多。。。)

外加Miller-Rabin素数测试,有空更新一下

(这份代码细节不太透彻,留个坑,NOIP后重新透彻一遍)

#include <bits/stdc++.h>

using namespace std;

long long const_A = 5;

unsigned long long myrand()
{
    static unsigned long long seed = 3589425289U;
    return seed * 74658973 + 348633;
}

void read(long long &x)
{
    x = 0;
    static char ch;
    ch = getchar();
    while (!isdigit(ch))
        ch = getchar();
    while (isdigit(ch))
    {
        x = x * 10 + ch - 48;
        ch = getchar();
    }
}

inline long long qpow(long long x, long long y, long long p)
{
    long long ans = 1;
    while (y > 0)
    {
        if (y & 1)
            ans = ((__int128)ans * (__int128) x) % p;
        x = ((__int128)x * (__int128)x) % p;
        y >>= 1;
    }
    return ans;
}

bool Miller_Rabin(long long N, long long a)
{
    if (N == a)
    	return true;
    if (N == 46856248255981LL || N < 2 || N % a == 0)
        return false; 
    long long d = N - 1;
    while (!(d & 1))
        d >>= 1;
    long long t = qpow(a, d, N);
    while (d != N - 1 && t != 1 && t != N - 1)
    {
        t = ((__int128)t * (__int128)t) % N;
        d <<= 1;
    }
    return (t == N - 1) || ((d & 1) == 1);
}

inline bool prime(long long x)
{
    return (Miller_Rabin(x, 2) && Miller_Rabin(x, 3) && Miller_Rabin(x, 7) && Miller_Rabin(x, 61) && Miller_Rabin(x, 24251));
}

inline long long gcd(long long a, long long b)
{
    if (!a) return b;
    if (!b) return a;
    int t = __builtin_ctzll(a | b);
    a >>= __builtin_ctzll(a);
    do
    {
        b >>= __builtin_ctzll(b);
        if (a > b)
        {
            long long tt = a;
            a = b;
            b = tt;
        }
        b -= a;
    }
    while(b != 0);
    return a << t;
}

long long Pollards_Rho(long long n)
{
    if (n % 2 == 0)
    	return 2;
    if (n % 3 == 0)
    	return 3;
    long long x = 0, y = x, t = 1, q = 1, c = (rand() % (n - 1)) + 1;
    for (register int k = 2; ; k <<= 1, y = x, q = 1)
    {
    	for (register int i = 1; i <= k; ++i)
    	{
    		x = (((__int128)x * (__int128)x) % n + c) % n;
    		q = ((__int128)q * (__int128)abs(x - y)) % n;
    		if ((i & 127) == 0)
            {
                t = gcd(q, n);
                if (t > 1)
                    break;
            }
    	}
    	if (t > 1 || (t = gcd(q, n)) > 1)
    		break;
    }
    return t;
}

void find(long long x, long long &ans)
{
//	printf("find %lld %lld\n", x, ans);
    if (x <= ans || x == 1)
        return;
    if (prime(x))
    {
        ans = ans > x ? ans : x;
        return;
    }
    long long t = x;
    while (t == x)
        t = Pollards_Rho(x);
    while (x % t == 0)
        x /= t;
    find(t, ans);
    find(x, ans);
}

int main()
{//freopen("d.txt", "r", stdin);
//	freopen("w.txt", "w", stdout);
    long long T;
    read(T);
    while (T--)
    {
        long long x;
        read(x);
        long long ans = 0;
        find(x, ans);
        if (x == ans)
            printf("Prime\n");
        else
            printf("%lld\n", ans);	
    }
    return 0;
}
o
粉丝 0
博文 500
码字总数 0
作品 0
私信 提问
加载中
请先登录后再评论。

暂无文章

hbase2.1.9 centos7 完全分布式 搭建随记

hbase2.1.9 centos7 完全分布式 搭建随记 这里是当初在三个ECS节点上搭建hadoop+zookeeper+hbase+solr的主要步骤,文章内容未经过润色,请参考的同学搭配其他博客一同使用,并记得根据实际情...

osc_4tfw1dxv
6分钟前
5
0
zookeeper3.5.5 centos7 完全分布式 搭建随记

zookeeper3.5.5 centos7 完全分布式 搭建随记 这里是当初在三个ECS节点上搭建hadoop+zookeeper+hbase+solr的主要步骤,文章内容未经过润色,请参考的同学搭配其他博客一同使用,并记得根据实...

osc_6jhxf9ab
8分钟前
0
0
steam夏日促销悄然开始,用Python爬取排行榜上的游戏打折信息

前言 很多人学习python,不知道从何学起。 很多人学习python,掌握了基本语法过后,不知道在哪里寻找案例上手。 很多已经做案例的人,却不知道如何去学习更加高深的知识。 那么针对这三类人,...

osc_ur9mmbck
8分钟前
0
0
python 里 certifi 库的作用

python 里 certifi 库的作用 安装了certifi之后,和requests库一样也有一个cacert.pem,可以用编辑器打开cacert.pem,里面包含了很多可信任知名公司的证书/公钥 库的路径,我这里是python2.7...

osc_1x6ycmfm
10分钟前
6
0
干掉"ZooKeeper",阿里为什么不用ZK做服务发现?

  20大进阶架构专题每日送达   链接:yq.aliyun.com/articles/601745   2020年Java面试题库连载中   !   正文   站在未来的路口,回望历史的迷途,常常会很有意思,因为我们会不...

osc_q5m9dzk0
12分钟前
9
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部