4

这里是...

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

我从关于Miller-Rabin primality test的维基百科文章中得到了这个。但我一直无法理解它......我不想理解它背后的数学，而只是在程序中实现它。在我看来，这个算法有点令人困惑。一个更好、更简单的伪代码或在 vb.net 中的实现会很有帮助。

编辑到目前为止编写的代码：

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
4

3 回答 3

5

根据要求，这是简单的伪代码：

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 次机会。大多数人使用介于 10 和 25 之间的kpowerMod值。该函数通过平方来执行幂运算，并且在您的语言没有为您提供它的情况下提供。

如果你想了解更多关于这个测试背后的数学知识，我谦虚地在我的博客上推荐这篇文章，其中还包括五种语言的实现，尽管它们都不是 VBA。

编辑：虽然他没有这么说，但原始发帖人实际上想要做的是找到小于 200 万的素数之和，从而解决 Project Euler 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

这里使用的算法是两千多年前由一位希腊数学家发明的埃拉托色尼筛法。同样，解释和代码在我博客的文章中。

于 2013-06-13T03:40:09.507 回答
4

关键思想和概念（p在这里代表素数）：

  1. 费马小定理。( a^(p-1) = 1 ( mod p ))
  2. 如果 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 ( mod p )，从 as 重复平方总是会产生 1 ，所以 (a^(p-1))%p 将为 1；或 a^(s*(2^r)) = -1 ( mod p ) 对于某些 r 使得 0 <= r < d，因为从它重复平方将始终产生 1，最后 a^(p-1) = 1 ( 模 p )。如果这些都不成立，那么对于任何素数 a，a^(p-1) 都不会为 1（否则将与事实 #2 相矛盾）。

算法 ：

  1. 让 p 是我们必须测试素数的给定数字。
  2. 首先我们将 p-1 重写为 (2^d)*s。（其中 s 为奇数且 d >= 0）。
  3. 现在我们在 [1,n-1] 范围内选择一些 a，然后检查是 as = 1 ( mod p ) 还是 a^(s*(2^r)) = -1 ( mod p )。
  4. 如果它们都失败了，那么 p 肯定是复合的。否则 p 可能是素数。我们可以选择另一个 a 并重复相同的测试。
  5. 我们可以在一些固定次数的迭代后停止并声称 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”时，至少有 (3/4) 小于 p 的数字将证明 p 是合数。这意味着如果我们进行 1 次迭代，复合数作为素数返回的概率是 (1/4)。通过 k 次迭代，测试失败的概率为 (1/4)k 或 4(-k)。与 Fermat 的测试相比，该测试相对较慢，但它不会分解任何特定的合数，并且 18-20 次迭代对于大多数应用程序来说是一个相当不错的选择。

PS：此函数计算 (a*b)%c 时考虑到 a*b 可能会溢出，我在 MILLER RABBIN 测试中使用过。

   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;
    }
于 2013-06-12T12:09:54.233 回答
2

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!
于 2013-06-12T16:04:40.187 回答