SICP 练习 1.28 - Miller-Rabin - "at least half the numbers will reveal a nontrivial square root of 1 modulo n"

SICP Exercise 1.28 - Miller-Rabin - "at least half the numbers will reveal a nontrivial square root of 1 modulo n"

SICP 练习 1.28

https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-11.html#%_thm_1.28

One variant of the Fermat test that cannot be fooled is called the Miller-Rabin test (Miller 1976; Rabin 1980). This starts from an alternate form of Fermat's Little Theorem, which states that if n is a prime number and a is any positive integer less than n, then a raised to the (n - 1)st power is congruent to 1 modulo n. To test the primality of a number n by the Miller-Rabin test, we pick a random number a < n and raise a to the (n - 1)st power modulo n using the expmod procedure. However, whenever we perform the squaring step in expmod, we check to see if we have discovered a ``nontrivial square root of 1 modulo n,'' that is, a number not equal to 1 or n - 1 whose square is equal to 1 modulo n. It is possible to prove that if such a nontrivial square root of 1 exists, then n is not prime. It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing a^(n-1) in this way will reveal a nontrivial square root of 1 modulo n. (This is why the Miller-Rabin test cannot be fooled.) Modify the expmod procedure to signal if it discovers a nontrivial square root of 1, and use this to implement the Miller-Rabin test with a procedure analogous to fermat-test. Check your procedure by testing various known primes and non-primes. Hint: One convenient way to make expmod signal is to have it return 0.

我写了自己的解决方案,其结果与这里提供的解决方案一致:

http://community.schemewiki.org/?sicp-ex-1.28

15 是一个非质数的奇数,因此对于从 114 的至少一半数字 a,我希望计算 expmod(a, 14, 15) 将显示 1 模 n 的非平凡平方根,由 expmod 返回 0.

表示

然而,这些是我得到的结果:

(expmod 1 14 15)
> 1
(expmod 2 14 15)
> 4
(expmod 3 14 15)
> 9
(expmod 4 14 15)
> 0
(expmod 5 14 15)
> 10
(expmod 6 14 15)
> 6
(expmod 7 14 15)
> 4
(expmod 8 14 15)
> 4
(expmod 9 14 15)
> 6
(expmod 10 14 15)
> 10
(expmod 11 14 15)
> 0
(expmod 12 14 15)
> 9
(expmod 13 14 15)
> 4
(expmod 14 14 15)
> 1

可以看出,这些结果中只有 2 个 0,与预期的至少 7 个相差甚远。

我是不是误解了这个说法?我是个彻头彻尾的白痴吗?代码有错吗? SICP错了吗?非常感谢。

编辑 1:要求我提供我正在使用的确切代码。就是这样,虽然我基本上只是复制我链接到的解决方案,并将 remainder 别名为 mod 因为我的解释器就是这样称呼它的。

 (define (square x) (* x x))

 (define remainder mod)

 (define (miller-rabin-expmod base exp m) 
   (define (squaremod-with-check x) 
     (define (check-nontrivial-sqrt1 x square) 
       (if (and (= square 1) 
                (not (= x 1)) 
                (not (= x (- m 1)))) 
           0 
           square)) 
     (check-nontrivial-sqrt1 x (remainder (square x) m))) 
   (cond ((= exp 0) 1) 
         ((even? exp) (squaremod-with-check 
                       (miller-rabin-expmod base (/ exp 2) m))) 
         (else 
          (remainder (* base (miller-rabin-expmod base (- exp 1) m)) 
                     m))))

(define expmod miller-rabin-expmod)

(print (expmod 1 14 15))
(print (expmod 2 14 15))
(print (expmod 3 14 15))
(print (expmod 4 14 15))
(print (expmod 5 14 15))
(print (expmod 6 14 15))
(print (expmod 7 14 15))
(print (expmod 8 14 15))
(print (expmod 9 14 15))
(print (expmod 10 14 15))
(print (expmod 11 14 15))
(print (expmod 12 14 15))
(print (expmod 13 14 15))
(print (expmod 14 14 15))

编辑 2:我现在还手动计算了 expmod(a, 14, 15) 的步骤(它总是通过 exp = 14exp = 7exp = 6exp = 3 递归, exp = 2, exp = 1, exp = 0), 对于 a 从 1 到 14 的所有值,我确定只有 a = 4a = 11 遇到 1 的非平凡平方根。所以我倾向于认为 SICP 在这方面要么是错误的,要么没有清楚地表达自己。

One variant of the Fermat test that cannot be fooled is called the Miller-Rabin test (Miller 1976; Rabin 1980). This starts from an alternate form of Fermat's Little Theorem, which states that if n is a prime number and a is any positive integer less than n, then a raised to the (n - 1)st power is congruent to 1 modulo n.

  • n是素数,a < n ==> a^(n-1) = 1 (mod n)

To test the primality of a number n by the Miller-Rabin test, we pick a random number a < n and raise a to the (n - 1)st power modulo n using the expmod procedure.

所以随机选择一个并检查是否 a^(n-1) = 1 (mod n)。如果它 不是 那么你知道 n 不是素数。

However, whenever we perform the squaring step in expmod, we check to see if we have discovered a ``nontrivial square root of 1 modulo n,'' that is, a number not equal to 1 or n - 1 whose square is equal to 1 modulo n.

这是关于在 expmod 函数中添加的额外检查。这可能是你忽略的事情。

It is possible to prove that if such a nontrivial square root of 1 exists, then n is not prime.

让我们详细了解一下。 1 的非平凡平方根将是一个数字 x,使得 x^2 = 1 (mod n)。并且 x 不是 1 或 -1。

为什么其中之一表明 n 不是质数?

我们知道 x^2 - 1 = (x - 1)(x + 1) 和 (working modulo n) x-1 和 x+1 都不为零,但它们的乘积为零。这意味着我们有一个复合 modulus,您可以通过计算这两个值的 GCD 来分解它。

It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing a^(n-1) in this way will reveal a nontrivial square root of 1 modulo n. (This is why the Miller-Rabin test cannot be fooled.)

这又是在谈论将内部测试添加到 expmod 函数的平方分支。

Modify the expmod procedure to signal if it discovers a nontrivial square root of 1, and use this to implement the Miller-Rabin test with a procedure analogous to fermat-test. Check your procedure by testing various known primes and non-primes. Hint: One convenient way to make expmod signal is to have it return 0.

希望对您有所帮助!询问您是否需要更多指导。

我找到了一篇涵盖此测试的论文和特定结果的证明,即 2 和 n-2 之间的一半以上的值将导致 1 的非平凡平方根。(定理 4.1)

我制作这段代码是为了仔细检查它

(define (print x) (display x) (newline))

(define (assert p) (unless p (error 'assert-failed)))

(define (power-of-two-split m)
  ;; write m = 2^e k
  (let loop ((e 0) (k m))
    (if (even? k)
        (loop (+ e 1) (quotient k 2))
        (cons e k))))

(define (exp/mod a k n)
  ;; a^k (mod n)
  (cond ((= k 0) 1)
        ((= k 1) (modulo a n))
        (else (modulo (* a (exp/mod a (- k 1) n)) n))))

(define (miller-rabin a n)
  (assert (odd? n))
  (assert (= 3 (modulo n 4))) ;; only handles e=1 case, need to use power-of-two-split for full test
  (let ((k (quotient (- n 1) 2)))
    (exp/mod a k n)))

(define (test n)
  (for ((i (in-range 2 (- n 2))))
    (let ((m (miller-rabin i n)))
      (print `(,i -> ,m squared ,(exp/mod m 2 n))))))

(test 15)

它打印出以下结果

(2 -> 8 squared 4)
(3 -> 12 squared 9)
(4 -> 4 squared 1)
(5 -> 5 squared 10)
(6 -> 6 squared 6)
(7 -> 13 squared 4)
(8 -> 2 squared 4)
(9 -> 9 squared 6)
(10 -> 10 squared 10)
(11 -> 11 squared 1)
(12 -> 3 squared 9)

所以对照米勒-拉宾证人的正式定义,他们实际上是 all wintesses:

Definition 2.3. For odd n > 1, write n − 1 = 2^ek with k odd and pick a ∈ {1, . . . , n − 1}. We say a is a Miller–Rabin witness for n if all of the congruences in are false: * a^k = 1 mod n and a^{(2^i)k} = −1 mod n for all i ∈ {0, . . . , e − 1}.

你可以看到 'm' 列中的 none 个值是 1,平方列中的 none 个值是 14。所以他们都是见证人,所以>50% 是。

执行 "check-nontrivial-sqrt1" 的代码与 n=3 (mod 4) 的这种特殊情况无关,因为在这种情况下 e=1。


更新:

我刚刚意识到为什么我们有很多证人但我们并不总能从他们身上找到一个平方根的原因:

the idea behind Miller–Rabin witnesses is to find an unexpected square root of 1 mod n. This is not always what we actually find, since the premise a n−1 ≡ 1 mod n for prime n might not actually be true for composite n.

这里有一个 table 的 a^(n-1) (mod n) for n=15

(2 -> 4)
(3 -> 9)
(4 -> 1)
(5 -> 10)
(6 -> 6)
(7 -> 4)
(8 -> 4)
(9 -> 6)
(10 -> 10)
(11 -> 1)
(12 -> 9)

如您所见,实际上只有两次同余 a^(n-1) = 1 (mod n) 成立。

SICP 是错误的,因为它使用了错误的 Miller-Rabin 证人定义(参见 Keith Conrad,The Miller–Rabin Test)。特别是,以下行是错误的:

Wrong claim. It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a < n, computing a^{n - 1} in this way will reveal a nontrivial square root of 1 modulo n.

你可以验证这是假的,例如当n = 9时。

根据上述参考资料中的定理 2.9,正确的说法应该是这样的:

Correct claim. Let n > 1 be an odd number that is not prime. Then we can write n - 1 as (2^e)k such that e ≥ 1 and k is odd. (For example, if n = 21, we can write 21 - 1 = 20 = (2^2)·5, so that e = 2 ≥ 1 and k = 5 is odd.) It is possible to prove that for at least half the numbers a < n, we will have a^k ≢ 1 and a^k ≢ n - 1 and a^{2k} ≢ n - 1 and ... and a^{2^{e-1}k} ≢ n - 1 modulo n.

因此,对于 n = 21,我们可以证明对于至少一半的数字 a < 21,我们将有 a^5 ≢ 1 和 a^5 ≢ 20 以及 a^10 ≢ 20。我们得到以下 table(所有值均以 21 为模):

+----+-----+-------+
| a  | a^5 |  a^10 | MILLER–RABIN WITNESS? 
+----+-----+-------+
|  1 |   1 |     1 | NO, a^5 ≡ 1
|  2 |  11 |    16 | YES
|  3 |  12 |    18 | YES
|  4 |  16 |     4 | YES
|  5 |  17 |    16 | YES
|  6 |   6 |    15 | YES
|  7 |   7 |     7 | YES
|  8 |   1 |     1 | NO, a^5 ≡ 1
|  9 |  18 |     9 | YES
| 10 |  19 |     4 | YES
| 11 |   2 |     4 | YES
| 12 |   3 |     9 | YES
| 13 |  13 |     1 | YES
| 14 |  14 |     7 | YES
| 15 |  15 |    15 | YES
| 16 |   4 |    16 | YES
| 17 |   5 |     4 | YES
| 18 |   9 |    18 | YES
| 19 |  10 |    16 | YES
| 20 |  20 |     1 | NO, a^5 ≡ 20
+----+-----+-------+

果然,超过一半的 a < 21(事实上,超过 75%)满足所有三个同余 a^5 ≢ 1、a^5 ≢ 20 和 a^10 ≢ 20。 (我们称他们为 Miller–Rabin 证人;因为他们见证了 n 不是素数的事实。一般来说,许多素数测试依赖于一些 属性 对所有素数都成立——如果你证明这样的 属性 某些数字失败,则该数字不能是素数。目击者越多,素数测试效果越好。)

编辑。 举个例子来说明素数,这里是 n = 13 的 table。自然不会有任何 Miller-Rabin 见证 13 作为它是质数;没有非素数可证。由于 n = 13,我们有 n - 1 = 12 = (2^2)·3,因此 e = 2 ≥ 1 且 k = 3 为奇数。因此,正如基思·康拉德 (Keith Conrad) 的说明性论文第 1 页所论证的那样,所有 a < 13 都将满足三个同余项中的至少一个:a^3 ≡ 1、a^3 ≡ 12、a^6 ≡ 12。果然:

+----+-----+-------+
| a  | a^3 |   a^6 | MILLER–RABIN WITNESS? 
+----+-----+-------+
|  1 |   1 |     1 | NO, a^3 ≡ 1
|  2 |   8 |    12 | NO, a^6 ≡ 12
|  3 |   1 |     1 | NO, a^3 ≡ 1
|  4 |  12 |     1 | NO, a^3 ≡ 12
|  5 |   8 |    12 | NO, a^6 ≡ 12
|  6 |   8 |    12 | NO, a^6 ≡ 12
|  7 |   5 |    12 | NO, a^6 ≡ 12
|  8 |   5 |    12 | NO, a^6 ≡ 12
|  9 |   1 |     1 | NO, a^3 ≡ 1
| 10 |  12 |     1 | NO, a^3 ≡ 12
| 11 |   5 |    12 | NO, a^6 ≡ 12
| 12 |  12 |     1 | NO, a^3 ≡ 12
+----+-----+-------+

根据@Memes 的回答,我已经为它添加了方案代码:

(define (display-all . vs)
  (for-each display vs))

(define (find-e-k n)
  (define (find-e-k-iter possible-k possible-e)
    (if (= (remainder possible-k 2) 0)
        (find-e-k-iter (/ possible-k 2) (+ possible-e 1))
        (values possible-e possible-k)))
  (find-e-k-iter (- n 1) 0))

; first-witness-case-test: (a ^ k) mod n # 1
(define (first-witness-case-test a k n)
  (not (= (expmod a k n) 1)))

; second-witness-case-test: all a ^ ((2 ^ i) * k) (with i = {0..e-1}) mod n # (n - 1)
(define (second-witness-case-test a e k n)
  (define (second-witness-case-test-iter a i k n)
    (cond ((= i -1) true)
          (else (let ()
                 (define witness (not (= (expmod a (* (fast-expt 2 i) k) n) (- n 1))))
                 (if witness
                 (second-witness-case-test-iter a (- i 1) k n)
                 false)))))
  (second-witness-case-test-iter a (- e 1) k n))

(define (miller-rabin-test n)
  (define (try-it a e k)
    (if (and (first-witness-case-test a k n) (second-witness-case-test a e k n))
        (display-all "is not prime, with a = " a "\n")
        (if (< a (- n 1))
            (try-it (+ a 1) e k)
            (display "is prime\n"))))
  (cond ((< n 2) (display "not prime"))
        ((= (remainder n 2) 0) (display "not prime\n"))
        (else (let ()
               (define-values (e k) (find-e-k n))
               (try-it 1 e k)))))