质因数分解的rho以及miller-rabin

一、前言

质因数分解,是一个在算法竞赛里老生常谈的经典问题。我们在解决许多问题的时候需要用到质因数分解来辅助运算,而且质因数分解牵扯到许许多多经典高效的算法,例如miller-rabin判断素数算法,rho启发式搜索质因数分解算法等。在此文里,我要介绍的就是miller-rabin算法以及rho启发式搜索分解算法。

二、算术基本定理

首先,我们得知道,任意一个大于1的自然数都可以分解为有限个质数的乘积。这里因子均为质数,且为正整数。我们把这样的分解成为N的标准分解式。关于算数基本定理的应用有许多,例如可以证明素数无限,定义god,lcm等,在此不一一赘述。有了算数基本定理,我们可以发现计算数论函数,约数和等都是十分方便的,这大大的方便了我们的解题。接下来我们介绍如何来分解质因数。

三、质因数分解的算法

所谓质因数,是某自然数的因素而且这个因素还得是质数。

我们不难想到基本的枚举,即暴力枚举1~n所有数,判断是否能够被n整除且该数是否为质数。大概代码如下:

For i:=1 to n do

If n mod i=0 then

If check(i) then  //check为检验i是否为质数的子函数,返回值为boolean

Writeln(i);

Function check(n:longint):boolean;

Var i:longint;

begin

For i:=2 to n do

If n mod i=0 then exit(false); //如果一旦有比它小的数除他为0,那么该数不为质数

Exit(true); //否则显然为质数

End;

这个是最朴素的方法,虽然代码简洁思路简单,但是时间复杂度实在是太高了。

于是我们有了优化1,我们发现判断一个数是否为素数只需枚举到sqrt(i),为何呢?因为一个数若不是素数,那么它必定是两个数的乘积,这两个数显然有一个小于等于sqrt(i),否则sqr(sqrt(i))=i,那么显然为质数。

Function check(n:longint):boolean;

Var i:longint;

Begin

For i:=2 to sqrt(n) do  If n mod i=0 then exit(false);

Exit(true);

End;

但是我们发现时间复杂度还是比较高O(sqrt(n)*n)

于是我们有了优化2。我们不暴力判断n的每个因素,而是不断去试除。

何谓试除?设x是n的质因数,我们在分解的过程中每找到一个x就把它记录下来,并且让n:=n div x;那么我们在遇见x的倍数时,显然不可能会再有x的倍数了。那么我们就能通过“类似筛法”的思想,在O(sqrt(n))时间内,对n进行质因数分解,这也是现在比较常用的方法。

I:=2;

While n<>1 do

Begin

While n mod i=0 do

Begin

Writeln(i);

N:=n div i;

End;

Inc(i);

End;

我们在不断的div i中,我们把i的所有存在于n的因素里的倍数全部给丢掉了。

(即,n=12时,n=2*2*3=3*4 那么我们在i=2的时候,把4给筛掉了。)

在这里,我们发现每个i都会是质数,即i的变化总是2→3→5→7.......

于是我们有了优化3。如果我们能预先存储不大于sqrt(n)的所有质数,我们就可以快速调用了,不要小看这个优化,在n的质因数比较大的时候,O(n)的线性时间是我们耗费不起的。先来说下筛法(摘自百度百科)

用筛法求素数的基本思想是:把从1开始的、某一范围内的正整数从小到大顺序排列, 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

1不是素数,去掉。剩下的数中2最小,是素数,去掉2的倍数,余下的数是:

3 5 7 9 11 13 15 17 19 21 23 25 27 29

剩下的数中3最小,是素数,去掉3的倍数,如此下去直到所有的数都被筛完,求出的素数为:

2 3 5 7 11 13 17 19 23 29

procedure get(n:longint);  //Get函数为利用筛法求素数

var maxfactor,newp:longint;

begin

maxfactor:=trunc(sqrt(n));

fillchar(hash,sizeof(hash),true);

for i:=2 to maxfactor do

if hash[i] then

begin

newp:=i shl 1;  //等价于newp:=i*2;

while newp<n do

begin

hash[newp]:=false;

inc(newp,i);

end;

end;

hash[1]:=false;

for i:=2 to n do

if hash[i] then

begin

inc(sushut);

sushu[sushut]:=i;

end;

end;

procedure solve(n:longint);   inline;

var i:longint;

begin

i:=1;

while sushu[i]<=n do

begin

tot:=0;

if sushu[i]=0 then break;

while n mod sushu[i]=0 do

begin

n:=n div sushu[i];

inc(tot);

if tot=1 then

begin

inc(ans);

sum[ans].number:=sushu[i];

end;

sum[ans].k:=tot;

end;

if sushu[i]=2 then fuck:=tot;

inc(i);

end;

end;

可是这样的复杂度在n为超大整数类型且素因子十分大的时候也无能为力了,不仅会MLE,更会TLE。那么,我们来说一下这次的重头戏,也就是rho分解大法。

我们必须知道Fermat分解:

首先,对于任意的一个偶数,我们都可以提取出一个2的质因子,如果结果仍为偶数,则可继续该操作,直至将其化为一个奇数和2的多少次幂的乘积,那么我们可以假定这个奇数可以被表示成2*N+1,如果这个数是合数,那么一定可以写成N=c*d的形式,不难发现,式中的c和d都是奇数,不妨设c>d,我们可以令a=(c+d)/2,b=(c-d)/2,那么的可以得到N=a*a-b*b,而这正是Fermat整数分解的基础;由不等式的关系,我们又可以得到a>=sqrt(c*d)=sqrt(N),那么,我们就可以枚举大于N的完全平方数a*a,计算a*a-N的值,判断计算的结果是否为一个完全平方数,如果是,那么a,b都是N的因子,我们就可以将算法递归的进行下去,知道求出N的所有质因子。容易看出,Fermat分解大数的效率其实并不高,但是比起试除法要好了很多;而且每次的计算都是计算出N的一个因子,更加降低了其效率。这就让我们想着去尝试新的算法,那就是Pollard rho算法。

Pollard rho算法的原理就是通过某种方法得到两个整数a和b,而待分解的大整数为n,计算p=gcd(a-b,n),直到p不为1,或者a,b出现循环为止。然后再判断p是否为n,如果p=n成立,那么返回n是一个质数,否则返回p是n的一个因子,那么我们又可以递归的计算Pollard(p)和Pollard(n/p),这样,我们就可以求出n的所有质因子。

具体操作中,我们通常使用函数x2=x1*x1+c来计算逐步迭代计算a和b的值,实践中,通常取c为1,即b=a*a+1,在下一次计算中,将b的值赋给a,再次使用上式来计算新的b的值,当a,b出现循环时,即可退出进行判断。

在实际计算中,a和b的值最终肯定一出现一个循环,而将这些值用光滑的曲线连接起来的话,可以近似的看成是一个ρ型的。

对于Pollard rho,它可以在O(sqrt(p))的时间复杂度内找到n的一个小因子p,可见效率还是可以的,但是对于一个因子很少、因子值很大的大整数n来说,Pollard rho算法的效率仍然不是很好,那么,我们还得寻找更加的方法了。

有了上述预备知识还不够,我们还得知道miller-robin素性检验。

根据费马小定理,

若n是一个奇素数,a是任何整数(1≤ a≤n-1) ,则 a^(n-1)≡1(mod n)。

miller-robin的理论基础:

如果n是一个奇素数, 将n-1表示成2^s*r的形式(r是奇 数),a 是和n互素的任何整数, 那么ar≡1(mod n) 或者对某个j(0≤j ≤s -1, j∈Z) 等式 a2jr ≡-1(mod n)成立。 这个理论是通过一个事实经由Fermat定理推导而来: n是一个奇素数,则方程x2 ≡ 1 mod n只有±1两个解。

注意,miller-rabin算法只是一种素性判断的概率算法,关于miller-rabin,pascal贴吧里应该有相关的帖子介绍,而miller-rabin的判断错误的概率是相当低的,(0.25)^t,t为素性测试的轮数。Miller-rabin算法在Int64范围内无法解决的强伪素数特殊判断一下即可。

接下来我们上代码:(我没找到Pas的,于是自己写了一份很丑的Pas代码,大神勿喷)

const count=10;

pri:array [0..10] of longint=(2,3,5,7,11,13,17,19,23,29,31);

var n,total:int64;

i:longint;

ans:array [0..10000] of int64;

function multi(a,b,m:int64):int64;

//考虑到64位二进制数能表示的范围,只需取前9个素数为基

var ans:int64;

begin

ans:=0;

a:=a mod m;

while b<>0 do

begin

if (b and 1)=1 then

begin

ans:=(ans+a) mod m;

dec(b);

end;

b:=b>>1;

a:=(a+a) mod m;

end;

exit(ans);

end;

function gcd(x,y:int64):int64;

begin

if x mod y=0 then

exit(y)

else exit(gcd(y,x mod y));

end;

function quick_mod(a,b,m:int64):int64;

var ans:int64;

begin

ans:=1; a:=a mod m;

while b<>0 do

begin

if (b and 1)=1 then

begin

ans:=multi(ans,a,m);

dec(b);

end;

b:=b>>1;

a:=multi(a,a,m);

end;

exit(ans);

end;

function prime(n:int64):boolean;

var m,k,a,x,y:int64; i,j:longint;

begin

if n=2 then exit(true);

if (n<2) or ((n and 1)=0) then exit(false);

m:=n-1; k:=0;

while (m and 1)=0 do

begin

inc(k);

m:=m>>1;

end;

randomize;

for i:=0 to count do

begin

a:=random(n) mod (n-1)+1;

x:=quick_mod(a,m,n);

y:=0;

for j:=0 to k-1 do

begin

y:=multi(x,x,n);

if (y=1) and (x<>1) and (x<>n-1) then exit(false);

x:=y;

end;

if y<>1 then exit(false);

end;

exit(true);

end;

function pollard_rho(n,c:int64):int64;

var i,k,x,y,d:int64;

begin

i:=1; k:=2;

randomize;

x:=random(n) mod (n-1)+1;

y:=x;

while true do

begin

inc(i);

x:=(multi(x,x,n)+c) mod n;

d:=gcd(y-x,n);

if (d>1) and (d<n) then exit(d);

if (x=y) then exit(n);

if (i=k) then

begin

y:=x;

k:=k<<1;

end;

end;

end;

procedure find(n,c:int64);

var p,k:int64;

begin

if n=1 then exit;

if prime(n) then

begin

inc(total);

ans[total]:=n;

exit;

end

else

begin

p:=n;

k:=c;

randomize;

while p>=n do

begin

p:=pollard_rho(p,c);

dec(c);

end;

find(p,k);

find(n div p,k);

end;

end;

begin

read(n);

find(n,120);

for i:=1 to total do

writeln(ans[i]);

end.

至此,我们在Int64范围内解决了质因数分解的若干问题。

时间: 2024-10-03 22:54:06

质因数分解的rho以及miller-rabin的相关文章

Pollard rho算法+Miller Rabin算法 BZOJ 3668 Rabin-Miller算法

BZOJ 3667: Rabin-Miller算法 Time Limit: 60 Sec  Memory Limit: 512 MBSubmit: 1044  Solved: 322[Submit][Status][Discuss] Description Input 第一行:CAS,代表数据组数(不大于350),以下CAS行,每行一个数字,保证在64位长整形范围内,并且没有负数.你需要对于每个数字:第一,检验是否是质数,是质数就输出Prime 第二,如果不是质数,输出它最大的质因子是哪个. O

HDU 3864 D_num Miller Rabin 质数判断+Pollard Rho大整数分解

链接:http://acm.hdu.edu.cn/showproblem.php?pid=3864 题意:给出一个数N(1<=N<10^18),如果N只有四个约数,就输出除1外的三个约数. 思路:大数的质因数分解只能用随机算法Miller Rabin和Pollard_rho,在测试多的情况下正确率是由保证的. 代码: #include <iostream> #include <cstdio> #include <cstring> #include <c

POJ1811_Prime Test【Miller Rabin素数测试】【Pollar Rho整数分解】

Prime Test Time Limit: 6000MS Memory Limit: 65536K Total Submissions: 29193 Accepted: 7392 Case Time Limit: 4000MS Description Given a big integer number, you are required to find out whether it's a prime number. Input The first line contains the num

POJ2429_GCD &amp;amp; LCM Inverse【Miller Rabin素数測试】【Pollar Rho整数分解】

GCD & LCM Inverse Time Limit: 2000MS Memory Limit: 65536K Total Submissions: 9756Accepted: 1819 Description Given two positive integers a and b, we can easily calculate the greatest common divisor (GCD) and the least common multiple (LCM) of a and b.

HDU1164_Eddy&amp;#39;s research I【Miller Rabin素数测试】【Pollar Rho整数分解】

Eddy's research I Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others) Total Submission(s): 6664    Accepted Submission(s): 3997 Problem Description Eddy's interest is very extensive, recently he is interested in prime

POJ2429_GCD &amp; LCM Inverse【Miller Rabin素数测试】【Pollar Rho整数分解】

GCD & LCM Inverse Time Limit: 2000MS Memory Limit: 65536K Total Submissions: 9756Accepted: 1819 Description Given two positive integers a and b, we can easily calculate the greatest common divisor (GCD) and the least common multiple (LCM) of a and b.

POJ1811_Prime Test【Miller Rabin素数測试】【Pollar Rho整数分解】

Prime Test Time Limit: 6000MS Memory Limit: 65536K Total Submissions: 29193 Accepted: 7392 Case Time Limit: 4000MS Description Given a big integer number, you are required to find out whether it's a prime number. Input The first line contains the num

HDU1164_Eddy&#39;s research I【Miller Rabin素数测试】【Pollar Rho整数分解】

Eddy's research I Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others) Total Submission(s): 6664    Accepted Submission(s): 3997 Problem Description Eddy's interest is very extensive, recently he is interested in prime

51_1037最长循环节 (miller rabin算法 pollard rho算法 原根)

1037 最长的循环节 V2 基准时间限制:1 秒 空间限制:131072 KB 分值: 320 难度:7级算法题 收藏 关注 正整数k的倒数1/k,写为10进制的小数如果为无限循环小数,则存在一个循环节,求<=n的数中,倒数循环节长度最长的那个数. 1/6= 0.1(6) 循环节长度为1 1/7= 0.(142857) 循环节长度为6 1/9= 0.(1)  循环节长度为1 Input 输入n(10 <= n <= 10^18) Output 输出<=n的数中倒数循环节长度最长的