sicp-ex-1.28



<< Previous exercise (1.27) | sicp-solutions | Next exercise (1.29) >>


This exercise is based on sicp-ex-1.27.

  
 ;; ex 1.28 
  
 (define (square x) (* x x)) 
  
 (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 (miller-rabin-test n) 
   (define (try-it a) 
     (define (check-it x) 
       (and (not (= x 0)) (= x 1))) 
     (check-it (miller-rabin-expmod a (- n 1) n))) 
   (try-it (+ 1 (random (- n 1))))) 
  
 (define (fast-prime? n times) 
   (cond ((= times 0) true) 
         ((miller-rabin-test n) (fast-prime? n (- times 1))) 
         (else false))) 
  
 (define (prime? n)  
    ; Perform the test how many times?  
    ; Use 100 as an arbitrary value.  
    (fast-prime? n 100))  
  
 (define (report-prime n expected)  
   (define (report-result n result expected)  
     (newline)  
     (display n)  
     (display ": ")  
     (display result)  
     (display ": ")  
     (display (if (eq? result expected) "OK" "FOOLED")))  
   (report-result n (prime? n) expected))  
    
 (report-prime 2 true)  
 (report-prime 7 true)  
 (report-prime 13 true)  
 (report-prime 15 false) 
 (report-prime 37 true)  
 (report-prime 39 false) 
   
 (report-prime 561 false)  ; Carmichael number  
 (report-prime 1105 false) ; Carmichael number  
 (report-prime 1729 false) ; Carmichael number  
 (report-prime 2465 false) ; Carmichael number  
 (report-prime 2821 false) ; Carmichael number  
 (report-prime 6601 false) ; Carmichael number 
  

Another solution that avoids nested functions. Similar to http://wiki.drewhess.com/wiki/SICP_exercise_1.28. However, this version does not intermingle the non-trivial-sqrt property check with returning the correct result in expmod. Instead, the non-trivial-sqrt property check is a dedicated function which returns true or false, as it should be. Only works for n>0. Uses the 'let' special form not introduced at this point of the book.

  
 (define (miller-rabin n) 
   (miller-rabin-test (- n 1) n)) 
  
 (define (miller-rabin-test a n) 
   (cond ((= a 0) true) 
         ; expmod is congruent to 1 modulo n 
         ((= (expmod a (- n 1) n) 1) (miller-rabin-test (- a 1) n)) 
         (else false))) 
  
 (define (expmod base exp m) 
   (cond ((= exp 0) 1) 
         ((even? exp) 
          (let ((x (expmod base (/ exp 2) m))) 
            (if (non-trivial-sqrt? x m) 0 (remainder (square x) m)))) 
         (else 
          (remainder (* base (expmod base (- exp 1) m)) 
                     m)))) 
  
 (define (non-trivial-sqrt? n m) 
   (cond ((= n 1) false) 
         ((= n (- m 1)) false) 
         ; book reads: whose square is equal to 1 modulo n 
         ; however, what was meant is square is congruent 1 modulo n 
         (else (= (remainder (square n) m) 1)))) 
  

Another solution in GNU Guile:

 ;;;; Under Creative Commons Attribution-ShareAlike 4.0 
 ;;;; International. See 
 ;;;; <https://creativecommons.org/licenses/by-sa/4.0/>. 
  
 (define-module (net ricketyspace sicp one twentyeight) 
   #:use-module (srfi srfi-1) 
   #:export (miller-rabin-test prime? run-tests)) 
  
 (define (sqmod x m) 
   "Return x^2 if `x^2 mod m` is not equal to `1 mod m` 
 and x != m - 1 and x != 1; 0 otherwise." 
   (let ((square (* x x))) 
     (cond ((and  (= (remainder square m) 1) ; 1 mod m = 1 
                  (not (= x (1- m))) 
                  (not (= x 1))) 
            0) 
           (else square)))) 
  
 (define (expmod base exp m) 
   (cond ((= exp 0) 1) 
         ((even? exp) 
          (remainder (sqmod (expmod base (/ exp 2) m) m) 
                     m)) 
         (else 
          (remainder (* base (expmod base (- exp 1) m)) 
                     m)))) 
  
 (define (miller-rabin-test n) 
   (define (pass? a) 
     (= (expmod a (1- n) n) 1)) 
   (fold (lambda (a p) (and (pass? a) p)) #t (iota (1- n) 1))) 
  
 (define (prime? n) 
   (if (miller-rabin-test n) #t #f)) 
  
  
 ;;; Tests 
  
 (define (carmichael-numbers-pass?) 
   "Return #t if the sample carmichael numbers are detected as non-prime." 
   (let ((numbers '(561 1105 1729 2465 2821 6601))) 
     (cons "carmichael-numbers-pass?" 
           (fold (lambda (n p) (and (not (prime? n)) p)) #t numbers)))) 
  
 (define (prime-numbers-pass?) 
   "Return #t if the sample prime numbers are detected as prime" 
   (let ((numbers '(311 641 829 599 809 127 419 13 431 883))) 
     (cons "prime-numbers-pass?" 
           (fold (lambda (n p) (and (prime? n) p)) #t numbers)))) 
  
 (define (even-numbers-pass?) 
   "Return #t if the sample even numbers are detected as non-prime" 
   (let ((numbers '(302 640 828 594 804 128 414 12 436 888))) 
     (cons "prime-numbers-pass?" 
           (fold (lambda (n p) (and (not (prime? n)) p)) #t numbers)))) 
  
 (define (run-tests) 
   (map (lambda (test) (test)) (list carmichael-numbers-pass? 
                                     prime-numbers-pass? 
                                     even-numbers-pass?))) 
  
 ;;; Guile REPL 
 ;;; 
 ;;; scheme@(guile-user)> ,use (net ricketyspace sicp one twentyeight) 
 ;;; scheme@(guile-user)> (run-tests) 
 ;;; $18 = (("carmichael-numbers-pass?" . #t)  
 ;;;        ("prime-numbers-pass?" . #t)  
 ;;;        ("prime-numbers-pass?" . #t))