给你…
Input: n > 3, an odd integer to be tested for primality;
Input: k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime
Write n − 1 as (2^s)·d with d odd by factoring powers of 2 from n − 1
WitnessLoop: repeat k times:
pick a random integer a in the range [2, n − 2]
x ← a^d mod n
if x = 1 or x = n − 1 then do next WitnessLoop
repeat s − 1 times:
x ← x^2 mod n
if x = 1 then return composite
if x = n − 1 then do next WitnessLoop
return composite
return probably prime
我从维基百科关于米勒-拉宾素数测试的文章中得到了这个。但是我一直没能理解......我不希望理解它背后的数学,而只是在一个程序中实现它。这个算法对我来说,有点令人困惑。一个更好,更简单的伪代码或在vb.net中实现它,将是有帮助的。
EDIT目前编写的代码:
Function Miller_Rabin(ByVal n As Integer) As Boolean
If n <= 3 Then : Return True
ElseIf n Mod 2 = 0 Then : Return False
Else
Dim k, s, a, d, x As Integer
k = 3
d = n - 1
While d Mod 2 = 0
d = d / 2
s += 1
End While
For c = 1 To k
a = Random(2, n - 1)
x = a ^ d Mod n
If x = 1 Or x = n - 1 Then GoTo skip
For r = 1 To s - 1
x = x ^ 2 Mod n
If x = 1 Then
Return False
Exit Function
Else
If x = n - 1 Then
GoTo skip
Else
Return False
Exit Function
End If
End If
Next
skip: Next
Return True
End If
End Function
Function Random(ByVal x As Integer, ByVal n As Integer) As Integer
Dim a As Integer = Now.Millisecond * Now.Second
skip:
a = (a ^ 2 + 1) Mod (n + 1)
If a < x Then
GoTo skip
Else
Return a
End If
End Function
下面是简单的伪代码:
function isStrongPseudoprime(n, a)
d := n - 1; s := 0
while d % 2 == 0
d := d / 2
s := s + 1
t := powerMod(a, d, n)
if t == 1 return ProbablyPrime
while s > 0
if t == n - 1
return ProbablyPrime
t := (t * t) % n
s := s - 1
return Composite
function isPrime(n)
for i from 1 to k
a := randInt(2, n-1)
if isStrongPseudoprime(n, a) == Composite
return Composite
return ProbablyPrime
function powerMod(b, e, m)
x := 1
while e > 0
if e % 2 == 1
x := (b * x) % m
b := (b * b) % m
e := e // 2 # integer division
return x
isStrongPseudoprime
函数检验了a是否证明n是复合的;注意,如果isStrongPseudoprime
返回Composite
,则数字肯定是复合的,但与之相反的是ProbablyPrime
,因为该数字有可能仍然是复合的。isPrime
功能测试k见证;通过设置k的值,您可以确定错误的可能性为4^k中的1次机会。大多数人使用的值k介于10和25之间。powerMod
函数通过平方执行求幂,如果您的语言不提供此功能,则提供此函数。
如果你想了解更多关于这个测试背后的数学知识,我在我的博客上谦虚地推荐这篇文章,它还包括五种语言的实现,尽管它们都不是VBA。
编辑:虽然他没有这么说,但最初的海报实际上想要做的是找到小于200万的质数的和,从而解决欧拉项目10。循环从2到n是一种非常低效的方法来求和小于n的质数;相反,推荐的方法是使用筛子。伪代码:
function sumPrimes(n)
sum := 0
sieve := makeArray(2..n, True)
for p from 2 to n step 1
if sieve[p]
sum := sum + p
for i from p * p to n step p
sieve[i] := False
return sum
这里使用的算法是埃拉托色尼筛法,由一位希腊数学家在两千多年前发明。同样,在我博客上的文章中有解释和代码。
关键思想与概念(p在这里代表素数):
- 费马小定理(a^(p-1) = 1 (mod p))
- 如果p是素数且x^2 = 1 (mod p),则x = +1或-1 (mod p)。
我们可以证明如下:
x^2 = 1 ( mod p )
x^2 - 1 = 0 ( mod p )
(x-1)(x+1) = 0 ( mod p )
如果p不能整除(x-1)和(x+1)而是整除它们的乘积,那么它不可能是素数,这是一个矛盾。因此,p要么除以(x-1)要么除以(x+1),所以x = +1或-1 (mod p)。
假设p - 1 = 2^d * s,其中s为奇数且d>= 0。如果p是素数,那么as = 1(对p取模)在这种情况下,as的重复平方总是得到1,所以(a^(p-1))%p等于1;或者a^(s*(2^r)) = -1 (mod p)对于某个r使得0 <= r
算法:
- 设p为我们要检验其是否为原数的给定数
- 首先我们把p-1写成(2^d)*s。(其中s为奇数且d>= 0)。
- 现在我们在[1,n-1]范围内选择一些a,然后检查as = 1 (mod p)或a^(s*(2^r)) = -1 (mod p)。
- 如果两者都失败,那么p肯定是复合的。否则p可能是素数。我们可以选择另一个a,重复相同的测试。
- 我们可以在一些固定的迭代次数后停止,并声称p肯定是合数,或者可能是素数。
小代码:Miller-Rabin素数检验,迭代表示检验的准确性
bool Miller(long long p,int iteration)
{
if(p<2)
return false;
if(p!=2 && p%2==0){
return false;
long long s=p-1;
while(s%2==0)
{
s/=2;
}
for(int i=0;i<iteration;i++)
{
long long a=rand()%(p-1)+1;
long long temp=s;
long long mod=modulo(a,temp,p);
while(temp!=p-1 && mod!=1 && mod!=p-1)
{
mod=mulmod(mod,mod,p);
temp *= 2;
}
if(mod!=p-1 && temp%2==0)
{
return false;
}
}
return true;
}
关于性能的几点:
可以证明,对于任何合数p,当在上述检验中选择"a"时,小于p的数中至少有(3/4)见证p是合数。这意味着如果我们做一次迭代,一个合数作为素数返回的概率是(1/4)。对于k次迭代,测试失败的概率是(1/4)k或4(-k)。与费马的测试相比,这个测试相对较慢,但它不会分解任何特定的合数,对于大多数应用程序来说,18-20次迭代是一个相当好的选择。
PS:这个函数计算(a*b)%c,考虑到a*b可能溢出,我在上面的MILLER RABBIN TEST中使用过。
long long mulmod(long long a,long long b,long long c)
{
long long x = 0,y=a%c;
while(b > 0)
{
if(b%2 == 1)
{
x = (x+y)%c;
}
y = (y*2)%c;
b /= 2;
}
return x%c;
}
VB实现使用十六进制转换函数来处理模求幂之前的大数。注释中提供的示例:
' USAGE:
' Example: strResult = mpModExp("3c", "03", "face")
' computes (0x3c)^3 mod 0xface = 0x5b56
' or, in decimal, 60^3 mod 64206 = 23382
' Parameters may be hex strings of any length subject to limitations
' of VB and your computer. May take a long time!