ROOT
ROOT
文章目录
  1. Miller_Rabin 算法
    1. 素数定理：
    2. 试除法
    3. 费马定理
    4. 伪素数测试过程
    5. Miller_Rabin 素数测试方法
    6. Miller_Rabin 素数测试的出错率
  2. Pollard_Rho 分解质因数算法
    1. 试除法分解合数
    2. 进一步优化
    3. 伪随机函数
    4. GCD
    5. 判环
    6. 最终代码
    7. 递归分解并储存所有质因数
  3. 习题
    1. POJ 1811
      1. 题目大意
      2. 思路
      3. AC 代码
    2. POJ 2429
      1. 题目大意
      2. 思路
      3. AC 代码
  4. 参考资料

蒙特卡罗法之随机素数判定算法 Miller_Rabin 与 Pollard_Rho 质因数分解算法

2015.12.06 Netcan 学习笔记  热度 ℃

昨天刷到 POJ 2429，没刷出来 = =，难度太大，翻Discuss 翻到了一题 POJ 1811，了解了两个算法：Miller_Rabin 素数判定与 Pollard_Rho 分解质因数算法。

Miller_Rabin 算法

这是我从 《算法导论》 和 Google 搜刮的资料上面整理出来的笔记，若有错误，请指出，谢谢。

素数定理：

素数分布函数 描述了小于等于 n 的素数数量。

在 非常小的情况下有 ，例如 次方时，，而，误差不超过。

随机选取一个数字 ，通过素数定理，该数字 为素数的概率为 ，期望值近似为，因此为了找出一个长度与 相同的素数，需要在 n 附近随机取 个整数。

试除法

假定 具有下列素数分解的因子：

当且仅当 并且 时，是素数。

所以可以用 去试除，最坏情况下，复杂度达到 ，若 非常大显然效率非常低，于是出现了一种高效判定大素数算法。。。

费马定理

若 是素数，则对任意，有

费 马定 理

证：

先引出一个结论：对任意 的数 除以 的余数正好为 1 到 的一个排列，例如 5 是素数，取，余数为，为 1 到 4 的一个排列。

反证法：若结论不成立的话，也就是至少有两个数 ，使得 与的余数相同，即

，做 差得

也就是说 可以整除 ，而 是素数且，与条件矛盾。

那么有

即

伪素数测试过程

根据费马定理，称 是一个基为 的伪素数。

费马定理意味着若 是素数，则对 都满足 费马 定理，令人惊讶的是它的逆命题也几乎成立（注意是几乎成立。。）

于是对 ，看看 是否满足 费马 定理 ，若不满足则返回COMPOSITE 声明 是合数，否则返回 PRIME，猜测 是素数（实际是只是 或者是素数，或者是基于 2 的伪素数）

if(mod_pow(a, n-1, n) != 1) return COMPOSITE; // 总是正确的
else return PRIME; // 可能错误，如果基于 a 的伪素数

这个过程可能产生错误，但是它判定合数总是正确的。如果判定素数，只有基于 2 的伪素数时才可能出错。

出错的几率非常小，但是不能无视它，在小于 10000 的 值中，有 22 个值会产生错误。可惜我们也不能通过取其他基数（例如 ）来消除所有的错误，因为还有一种称作Carmichael 的合数 ，对任意，满足 费马 定理。（真抓狂 233）

下一步是对素数测试方法改进，使得测试过程中不会将 Carmichael 当成素数。

Miller_Rabin 素数测试方法

Miller_Rabin 测试方法对上述测试过程做了改进，克服其存在的问题。

测试方法引进 witness，当且仅当 为合数 的证据时，返回 TRUE，witness 取测试的扩展：

来介绍下 witness 的构造过程，令 ，其中 且是奇数，即 的二进制表示是奇数 的二进制后面接 个。

则 ，所以可以通过先计算，然后 连续平方 次来计算最终的 ，期间若有一次 或者 ， 即则返回合数。

现在再来证明一下为什么 Math output error 时，为合数。

首先引出定理：

如果 是一个奇素数 () 且，则方程

仅有两个解，即。

证：

方程等价于 (| 是整除的意思。。)：

这说明了 或者 ，但不同时成立，否则 也能整除它们的差。

若，则，也就是。对称地，得出，综上可得。

所以 不断平方 ，期间Math output error，则可说明 是合数。

下面是 witness 实现过程

bool witness(ll a, ll n, ll u, ll t) { // a 随机数，u, t 于外部计算
	ll ret = mod_pow(a, u, n);
	ll last = ret;
	for(int i=1; i<=t; ++i) {ret = mod_mult(ret, ret, n);
		if(ret == 1 && last != 1 && last != n-1) return true;
		last = ret;
	}
	if(ret!=1) return true; // a^(n-1)%n!=1
	return false;
}

举个例子，比如 Carmichael 数，首先，随机取，这里采用从后往前的顺序计算。

  1. 先计算，即
  2. 再计算 ，直接说明了 是合数 = =

若按程序从前往后算，，，直到算到 即能检验出来。

整个 序列如下：

最后的 Miller_Rabin 函数就出炉了：

bool Miller_Rabin(ll n, int s) { // 随机取 s 次 a
	if(n==2) return true;
	if(n < 2 || !(n&1)) return false;
	ll u=n-1;
	ll t=0;
	while(!(u&1)) {u>>=1; ++t;} // 2^t*u = n-1
	for(int i=0; i<s; ++i) {ll a=rand()%(n-1)+1; // 若 n 为合数，证据为 a 的概率至少为 1/2
		if(witness(a, n, u, t)) return false; // 合数
	}
	return true;
}

Miller_Rabin 素数测试的出错率

如果 Miller_Rabin 返回 true(即PRIME)，则它还是有一种很小的可能性出错。然而概率非常低，仅取决于s 的大小和选取 a 值的运气。。而且每次都需要经过 witness 的严格判断，出错率应该极小，出错概率最多为 ，然而对于一个长度为1024 位的数字，大概需要 次，对于 64 位数，长度不超过 20，计算得 ，取 出错率可忽略不计。

Pollard_Rho 分解质因数算法

首先用 Miller_Rabin 判断一个数是否是质数，若不是质数则对其分解质因数。

试除法分解合数

试除法前面有介绍，效率较低，复杂度达。

试除法是一个一个试，可对其进行随机选取的方法提高效率。

简单起见，假设 是一个合数，并且只有两个真因子。

若从小于 的数随机取一个，取到真因子的概率为 ，若 比较大，效率也比较低。

进一步优化

假设从 1000 中随机去一个数，为 6 的概率为。

若换种方法，每次取两个数 ，问 是否为 6，概率是，概率提高了。。

试着将这个方法推广：随机选择 k 个数，并在其中两两比较。（对于 n，在 处有一半的命中率。）

所以直接推广，选出 个数并在之中两两比较，寻找到 （即 或）

可以问题没有这么简单，现在我们要储存 个数，还要进行 次比较，而且只有一半的概率找到真因子。

伪随机函数

显然，我们不能储存 个数，也不能进行 n 次比较。

如果我们只是在生成随机数时，比较生成时间上相邻的两个随机数，就可以同时解决上面的两个问题。

有伪随机函数 ，其中 N 是常数，a 是随机生成的常数。 并给出序列。（ 原文 没有解释为什么这个函数可以替换上述机制）

GCD

同时我们可以优化检测手段：我们不检测，而是检测

因为只有 整除，令，则 不 再内 ， 因为

共有 个。

判环

看似大功告成了，不过我们引入的伪随机函数 还有问题，生成的数的轨迹有可能有环（这也是算法名字中 Rho()的来历）。也就是对 的不断迭代可能陷入一个死循环。

为此我们需要在终点未知而空间有限的条件下判断是否陷入环中。

如果发现了有环，可以重新（随机）给出常数，再次运行算法。

最终代码

ll Pollard_rho(ll n, ll c) {// 伪随机函数 f(x)=x^2+c，c 为随机数
	ll i=1, k=2;
	ll x=rand()%n;
	ll y = x;
	while(1) {
		++i;
		x=(mod_mult(x, x, n)+c)%n; // 伪随机数
		ll d = gcd(y-x+n, n); // gcd 注意负数
		if(d!=1 && d!=n) return d;
		if(y==x) return n; // 遇环退出
		if(i==k) {y=x; k<<=1;}
	}
}

递归分解并储存所有质因数

这里再给出递归分解所有质因数的方法。

vector<ll> factors; // 储存质因数分解结果（无序）
void findfac(ll n) {if(Miller_Rabin(n, 20)) { // 素数
		factors.push_back(n); // 储存素因子
		return;
	}
	ll p = n;
	while(p>=n) p = Pollard_rho(p, rand()%(n-1)+1); // 找出当前合数的一个素因子
	findfac(p);
	findfac(n/p);
}

习题

POJ 1811

题目地址：http://poj.org/problem?id=1811

题目大意

给你一个 64 位大整数，判断是否为素数，若不是，则找出最小质因子。

思路

直接用 Miller_Rabin 进行判定，然后 findfac 求出所有质因子，最后输出最小的质因子。

AC 代码

/*************************************************************************
	> File Name: 1811.cpp
	  > Author: Netcan
	  > Blog: http://www.netcan.xyz
	  > Mail: 1469709759@qq.com
	  > Created Time: 2015-12-05 六 14:27:28 CST
 ************************************************************************/

#include <iostream>
#include <vector>
#include <string>
#include <queue>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <sstream>
#include <deque>
#include <functional>
#include <iterator>
#include <list>
#include <map>
#include <memory>
#include <stack>
#include <set>
#include <numeric>
#include <utility>
#include <cstring>
using namespace std;
typedef long long ll;

ll mod_mult(ll a, ll b, ll mod) { //return a*b%mod
	ll s = 0;
	while(b) {
		if(b&1) s=(s+a)%mod;
		a<<=1, b>>=1;
		if(a > mod) a %= mod;
	}
	return s;
}

ll mod_pow(ll x, ll n, ll mod) {  // a^b%mod
	ll res = 1;
	while(n) {
		if(n&1) res = mod_mult(res, x, mod);
		x=mod_mult(x, x, mod);
		n>>=1;
	}
	return res;
}

// 以 a 为基,n-1=u*2^t      a^(n-1)=1(mod n)  验证 n 是不是合数
// 一定是合数返回 true, 不一定返回 false
bool witness(ll a, ll n, ll u, ll t) { // a 随机数，u, t 于外部计算
	ll ret = mod_pow(a, u, n);
	ll last = ret;
	for(int i=1; i<=t; ++i) {
		ret = mod_mult(ret, ret, n);
		if(ret == 1 && last != 1 && last != n-1) return true;
		last = ret;
	}
	if(ret!=1) return true; // a^(n-1)!=1(mod n)
	return false;
}

// Miller_Rabin()算法素数判定
// 是素数返回 true.(可能是伪素数，但概率极小，P=(1/2)^s)
// 合数返回 false;
// s: 挑选 s 个随机数
bool Miller_Rabin(ll n, int s) {
	if(n==2) return true;
	if(n < 2 || !(n&1)) return false;
	ll u=n-1;
	ll t=0;
	while(!(u&1)) {u>>=1; ++t;} // 2^t*u = n-1
	for(int i=0; i<s; ++i) {
		ll a=rand()%(n-1)+1; // 若 n 为合数，证据为 a 的概率至少为 1/2
		if(witness(a, n, u, t)) return false; // 合数
	}
	return true;
}

//************************************************
//pollard_rho 算法进行质因数分解
//************************************************
vector<ll> factors; // 质因数分解结果（无序）

inline ll gcd(ll a, ll b) {
	return b==0?a:gcd(b, a%b);
}

ll Pollard_rho(ll n, ll c) { // 伪随机函数 f(x)=x^2+c，c 为随机数
	ll i=1, k=2;
	ll x=rand()%n;
	ll y = x;
	while(1) {
		++i;
		x=(mod_mult(x, x, n)+c)%n; // 伪随机数
		ll d = gcd(y-x+n, n); // gcd 注意负数
		if(d!=1 && d!=n) return d;
		if(y==x) return n; // 遇环退出
		if(i==k) {y=x; k<<=1;}
	}
}

// 对 n 进行递归素因子分解
void findfac(ll n) {
	if(Miller_Rabin(n, 20)) { // 素数
		factors.push_back(n); // 储存素因子
		return;
	}
	ll p = n;
	while(p>=n) p = Pollard_rho(p, rand()%(n-1)+1); // 找出当前合数的一个素因子
	findfac(p);
	findfac(n/p);
}

int T;
ll N;

void solve() {
	if(Miller_Rabin(N, 20)) {
		cout << "Prime" << endl;
	}
	else {
		factors.clear();
		findfac(N);
		cout << *min_element(factors.begin(), factors.end()) << endl;
	}
	// cout << "因子：\n";
	// for(vector<ll>::iterator i=factors.begin(); i!=factors.end(); ++i) {
		// printf("%lld", *i);
	// }
	// cout <<endl;
}

int main()
{
#ifdef Oj
	freopen("1811.in", "r", stdin);
#endif
	cin >> T;
	while(T--) {
		cin >> N;
		solve();
	}
	return 0;
}

POJ 2429

题目地址：http://poj.org/problem?id=2429

题目大意

给你两个数 ，求他们的最大公约数 和最小公倍数 比较简单；那么给出 和，求出 （升序），若有多解，求出 最小的那组解。

思路

题目给的 GCD, LCM 都是 64 位的，暴力是行不通的。

那么用 Pollard_Rho 分解 求出所有质因子，然后在合并相同的质因子，最终搜索最小 p+q，则

AC 代码

/*************************************************************************
	> File Name: 2429.cpp
	  > Author: Netcan
	  > Blog: http://www.netcan.xyz
	  > Mail: 1469709759@qq.com
	  > Created Time: 2015-12-07 一 21:11:13 CST
 ************************************************************************/

#include <iostream>
#include <vector>
#include <string>
#include <queue>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <sstream>
#include <deque>
#include <functional>
#include <iterator>
#include <list>
#include <map>
#include <memory>
#include <stack>
#include <set>
#include <numeric>
#include <utility>
#include <cstring>
using namespace std;
typedef long long ll;
ll GCD, LCM;

ll mod_mult(ll a, ll b, ll mod) { // return	a*b%mod
	ll s = 0;
	while(b) {
		if(b&1) s=(s+a)%mod;
		a<<=1, b>>=1;
		if(a > mod) a -= mod;
	}
	return s;
}

ll mod_pow(ll x, ll n, ll mod) { // return a^b%mod
	ll res = 1;
	while(n) {
		if(n&1) res = mod_mult(res, x, mod);
		x=mod_mult(x, x, mod);
		n>>=1;
	}
	return res;
}
// 以 a 为基,n-1=u*2^t      a^(n-1)=1(mod n)  验证 n 是不是合数
// 一定是合数返回 true, 不一定返回 false
bool witness(ll a, ll n, ll u, ll t) { // a 随机数，u, t 于外部计算
    ll ret = mod_pow(a, u, n);
    ll last = ret;
    for(int i=1; i<=t; ++i) {
        ret = mod_mult(ret, ret, n);
        if(ret == 1 && last != 1 && last != n-1) return true;
        last = ret;
    }
    if(ret!=1) return true; // a^(n-1)!=1(mod n)
    return false;
}

// Miller_Rabin()算法素数判定
// 是素数返回 true.(可能是伪素数，但概率极小，P=(1/2)^s)
// 合数返回 false;
// s: 挑选 s 个随机数
bool Miller_Rabin(ll n, int s) {
    if(n==2) return true;
    if(n < 2 || !(n&1)) return false;
    ll u=n-1;
    ll t=0;
    while(!(u&1)) {u>>=1; ++t;} // 2^t*u = n-1
    for(int i=0; i<s; ++i) {
        ll a=rand()%(n-1)+1; // 若 n 为合数，证据为 a 的概率至少为 1/2
        if(witness(a, n, u, t)) return false; // 合数
    }
    return true;
}

//************************************************
//pollard_rho 算法进行质因数分解
//************************************************
map<ll, int> factors; // 质因数分解结果（无序，每个质因子的 n 次方）

inline ll gcd(ll a, ll b) {
    return b==0?a:gcd(b, a%b);
}

ll Pollard_rho(ll n, ll c) { // 伪随机函数 f(x)=x^2+c，c 为随机数
    ll i=1, k=2;
    ll x=rand()%n;
    ll y = x;
    while(1) {
        ++i;
        x=(mod_mult(x, x, n)+c)%n; // 伪随机数
        ll d = gcd(y-x+n, n); // gcd 注意负数
        if(d!=1 && d!=n) return d;
        if(y==x) return n; // 遇环退出
        if(i==k) {y=x; k<<=1;}
    }
}

// 对 n 进行递归素因子分解
void findfac(ll n) {
    if(Miller_Rabin(n, 20)) { // 素数
		++factors[n]; // 储存素因子，用 map 后面方便处理
        return;
    }
    ll p = n;
    while(p>=n) p = Pollard_rho(p, rand()%(n-1)+1); // 找出当前合数的一个素因子
    findfac(p);
    findfac(n/p);
}

void dfs(const vector<ll> &nf,const ll d, int i, ll& sum,  ll x, ll &a) { // 搜索最小的 a+b，x 为当前搜索的 a，sum 为储存最小和，a 为最终答案
	if(i >= nf.size()) return;
	if(sum > d/x+x) { // 更新最小和
		a=x;
		sum = d/x+x;
	}
	dfs(nf, d, i+1, sum, x, a); // 不选 nf[i]并搜索
	x*=nf[i]; // 选择 nf[i]并搜索
	if(sum > d/x+x) { // 更新最小和
		a = x;
		sum = d/x+x;
	}
	dfs(nf, d, i+1, sum, x, a); // 这是选择 nf[i]的情况
}

void solve() {
	factors.clear();
	vector<ll> nfactors; // 最终因子
	ll d = LCM/GCD; // (a/gcd)(b/gcd)=lcm/gcd，然后分解 d，求出最终的(a/gcd), (b/gcd)
	if(d==1) {
		printf("%lld %lld\n", GCD, LCM);
		return;
	}
	findfac(d);

	for(map<ll, int>::iterator i = factors.begin(); i!=factors.end(); ++i)
		nfactors.push_back(mod_pow(i->first, i->second, (1ll<<63)-1)); // 对相同的质因子进行合并

	// for(vector<ll>::iterator i=nfactors.begin(); i!=nfactors.end(); ++i) {
		// printf("%lld", *i);
	// }
	// cout << endl;
	ll a = 1;
	ll sum = a+d/a;
	dfs(nfactors, d, 0, sum, a, a);

	// for(int i=0; i<(1<<nfactors.size()); ++i) {
		// ll aa = 1;
		// for(int j=0; j<nfactors.size(); ++j) // 神奇的搜索。。。
			// if(i & (1<<j)) aa*=nfactors[j];

		// if(d/aa + aa < min_sum) {
			// min_sum = d/aa + aa;
			// a = aa;
		// }
	// }
	a = min(a, d/a);

	printf("%lld %lld\n", a*GCD, d/a*GCD);
}

int main()
{
#ifdef Oj
	freopen("2429.in", "r", stdin);
#endif
	while(cin >> GCD >> LCM) {
		solve();
	}
	return 0;
}

参考资料

  1. 《算法导论》素数测试部分
  2. 数论部分第一节：素数与素性测试
  3. 笔记：A Quick Tutorial on Pollard’s Rho Algorithm
支持一下
扫一扫，支持Netcan
  • 微信扫一扫
  • 支付宝扫一扫