首页 > 代码库 > 素数检测杂谈

素数检测杂谈

最开始我们学到的是这种朴素的写法:

 

(define (prime? n)
  (cond
   [(or (= n 2) (= n 3)) #t]
   [(even? n) #f]
   [else (prime-test-loop n)]))

(define (prime-test-loop n)
  (let ((top (ceiling (sqrt n))))
    (let iter ((start 3))
      (cond [(> start top) #t]
            [(= (mod n start) 0) #f]
            [else (iter (+ start 2))]))))

 

虽然很笨拙,但其实这也是经过“优化”的,首先,除了 2 以外的偶数就不用判断了。其次,试除迭代从 3 起步,每次 +2,同样避开了偶数。最后,循环的结束条件是 (sqrt n). 

后来又看到一种生成素数序列的方法,大体思想是维护一个已知的素数表,每一个新的数字 n 都用已知素数表里的数去试除。如果都不能整除,说明 n 是素数,将 n 添加到已知的素数表里,进入下一轮迭代。原文是用 C 实现的,事先搞一个大数组,每一轮迭代都遍历数组,将已知素数的整数倍所对应的数组下标都划掉,最后留在数组中的就是素数。下面是我写的 Scheme 版:

 

(define (prime-list n)
  (let iter ((start 3) (primes ‘(2)))
    (cond [(> start n) primes]
          [(andmap (lambda (p) (not (= (mod start p) 0)))
                   primes)
           (iter (+ start 2) (cons start primes))]
          [else (iter (+ start 2) primes)])))

 

后来知道有一种基于概率的检测方法叫费马检查,SICP 上就有一个 Scheme 实现:

(define (square n) (* n n))

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (mod (square (expmod base (/ exp 2) m))
              m))
        (else
         (mod (* base (expmod base (- exp 1) m))
              m))))

(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) #t)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (else #f)))

这个算法很不可靠,因为它所依据的理论是个伪命题,费马小定理断言:如果 n 是一个素数,那么小于 n 的任意正整数 a 的 n 次方再对 n 取模,结果依然为 a。

它的逆命题是,a 是任意一个小于 n 的正整数,  如果 a ^ n % n == a,那么 n 就是一个素数。很遗憾,这个逆命题不成立,因为有一类数叫做 Carmichael 数,它符合  a ^ n % n == a 这个特性,但却不是素数。直接用这段程序检查一个数是不是素数有很大的机会出错。这是已知 Carmichael 数的一部分:

561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361, 101101, 115921, 126217, 162401, 172081, 188461, 252601, 278545, 294409, 314821, 334153, 340561, 399001, 410041, 449065, 488881, 512461

可以看到,最小的数是561。如果你准备写一个程序,生成一个素数表,然而它迭代了500多次以后就开始胡言乱语,那玩笑就有点开大了。SICP 提到一个费马检测的变形,叫做 Miller-Rabin 检测,它能够避免 Carmichael 数的愚弄。

这个检测依据的断言是:如果 n 是一个素数,a 是任意一个小于 n 的正整数,那么 a ^ (n-1) % n == 1 % n。

要用 米勒-拉宾 检测算法测试一个数字是不是素数,我们应当选择一个小于 n 的随机正整数 a, 然后利用 expmod 函数计算 a ^ (n-1) % n . 但是当执行到 expmod 里的 square 那一步时,我们检测是不是找到了一个 (1 % n) 的非平凡的平方根,换言之,看一看是不是有这样一个数字,它既不是 1, 也不是 n-1,但是它的平方等于 1 % n。很容易证明,如果这样一个 1 的非平凡的平方根存在,那么 n 肯定不是素数。同样可以证明,如果 n 是一个并非素数的奇数,在小于 n 的正整数中,至少有一半在用这种方法计算 a ^ (n-1) 时能够找到这种非凡平方根。

中文版的 SICP 看不太明白,我自己翻译英文题目,还是觉得看不太明白。但大概意思明白了,如果一个小于 n 的正整数 a,它符合 a ^ (n-1) & n == 1,但是 a 既不是 1, 也不是 n-1,那么这个 n 肯定不是素数。对一个伪素数来说,只做一次测试,正确与错误的概率对半分,基本和算命差不多。如果做两次测试,那就很有可能发现它的真面目。重复的次数越多,结果越可信。如果重复多次依然通过了测试,那个这个数 n 大概真的是个素数了。这就是它能避免 Carmicheal 数愚弄的原因。

 

SICP 留了一个作业 1.28 ,要求实现 Miller-Rabin 检测,下面就是实现:

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (check-nontrivial-sqrt (expmod base (/ exp 2) m) m)) ;; look here
        (else
         (mod (* base (expmod base (- exp 1) m)) m))))

(define (check-nontrivial-sqrt n m)
  (let ((x (mod (square n) m)))
    (if (and (not (= n 1))
             (not (= n (- m 1)))
             (= x 1))
        0
        x)))

(define (miller-rabin-test n)
  (define (try-it a)
    (= (expmod a (- n 1) n) 1))
  (try-it (+ 1 (random (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) #t)
        ((miller-rabin-test n) (fast-prime? n (- times 1)))
        (else #f)))

;; 搞了个包装函数,偶数就不必判断了吧 (define (prime? n) (cond ((= n 2) #t) ((even? n) #f) (else (fast-prime? n 20))))

好吧,写到是写出来了,但是测试发现,它依然胡言乱语。。。

等有空再调试一下。

素数检测杂谈