[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/UudIdev3    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Dec 27:
; carmichael numbers

(define (expm b e m)
  (define (m* x y) (modulo (* x y) m))
  (cond ((zero? e) 1)
        ((even? e) (expm (m* b b) (/ e 2) m))
        (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))

(define (primes n)
  (let* ((max-index (quotient (- n 3) 2))
         (v (make-vector (+ 1 max-index) #t)))
    (let loop ((i 0) (ps '(2)))
      (let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
        (cond ((>= (* p p) n)
               (let loop ((j i) (ps ps))
                  (cond ((> j max-index) (reverse ps))
                        ((vector-ref v j)
                          (loop (+ j 1) (cons (+ j j 3) ps)))
                        (else (loop (+ j 1) ps)))))
              ((vector-ref v i)
                (let loop ((j startj))
                  (if (<= j max-index)
                      (begin (vector-set! v j #f)
                             (loop (+ j p)))))10
                      (loop (+ 1 i) (cons p ps)))
              (else (loop (+ 1 i) ps)))))))

(define (fermat-pseudo-prime? a n)
  (and (= (gcd a n) 1)
       (= (expm a (- n 1) n) 1)))

(display (fermat-pseudo-prime? 2 561)) (newline)
(display (fermat-pseudo-prime? 3 561)) (newline)

(define (strong-pseudo-prime? a n)
  (let loop ((r 0) (s (- n 1)))
    (if (even? s) (loop (+ r 1) (/ s 2))
      (if (= (expm a s n) 1) #t
        (let loop ((j 0) (s s))
          (cond ((= j r) #f)
                ((= (expm a s n) (- n 1)) #t)
                (else (loop (+ j 1) (* s 2)))))))))

(display (strong-pseudo-prime? 2 561)) (newline)
(display (strong-pseudo-prime? 3 561)) (newline)

(define pseudo-prime?
  (let ((as (primes 200)))
    (lambda (method? n)
      (if (member n as) #t
        (let loop ((as as))
          (cond ((null? as) #t)
                ((not (method? (car as) n)) #f)
                (else (loop (cdr as)))))))))

(define (fermat-prime? n)
  (pseudo-prime? fermat-pseudo-prime? n))

(define (strong-prime? n)
  (pseudo-prime? strong-pseudo-prime? n))

(display (fermat-prime? 561)) (newline)
(display (strong-prime? 561)) (newline)

(define (carmichael? n)
  (let ((ps (primes n)))
    (if (= (car (reverse ps)) n) #f
      (let loop ((ps (primes n)))
        (cond ((null? ps) #t)
              ((not (= (gcd (car ps) n) 1)) (loop (cdr ps)))
              ((not (fermat-pseudo-prime? (car ps) n)) #f)
              (else (loop (cdr ps))))))))

(display
  (let loop ((c 4000) (cs '()))
    (cond ((= c 2) cs)
          ((carmichael? c)
            (loop (- c 1) (cons c cs)))
          (else (loop (- c 1) cs)))))

(newline)

(define (factors n)
  (if (even? n) (cons 2 (factors (/ n 2)))
    (let loop ((n n) (f 3) (fs '()))
      (cond ((< n (* f f)) (reverse (cons n fs)))
            ((zero? (modulo n f)) (loop (/ n f) f (cons f fs)))
            (else (loop n (+ f 2) fs))))))

(define (carmichael? n)
  (let ((fs (factors n)))
    (and (odd? n) (pair? (cdr fs)) (apply < fs)
         (let loop ((fs fs))
           (cond ((null? fs) #t)
                 ((not (zero? (modulo (- n 1) (- (car fs) 1)))) #f)
                 (else (loop (cdr fs))))))))

(display
  (let loop ((c 10000) (cs '()))
    (cond ((= c 2) cs)
          ((carmichael? c)
            (loop (- c 1) (cons c cs)))
          (else (loop (- c 1) cs)))))


Output:
1
2
3
4
5
6
7
8
#t
#f
#f
#f
#f
#f
(561 1105 1729 2465 2821)
(561 1105 1729 2465 2821 6601 8911)


Create a new paste based on this one


Comments: