[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Jan 24:
; a dozen lines of code
(define sort #f)
(define merge #f)
(let ()
  (define dosort
    (lambda (pred? ls n)
      (if (= n 1)
          (list (car ls))
          (let ((i (quotient n 2)))
            (domerge pred?
                     (dosort pred? ls i)
                     (dosort pred? (list-tail ls i) (- n i)))))))
  (define domerge
    (lambda (pred? l1 l2)
      (cond
        ((null? l1) l2)
        ((null? l2) l1)
        ((pred? (car l2) (car l1))
         (cons (car l2) (domerge pred? l1 (cdr l2))))
        (else (cons (car l1) (domerge pred? (cdr l1) l2))))))
  (set! sort
    (lambda (pred? l)
      (if (null? l) l (dosort pred? l (length l)))))
  (set! merge
    (lambda (pred? l1 l2)
      (domerge pred? l1 l2))))

(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)))))

; primes n -- list of primes not greater than n in ascending order
(define (primes n) ; assumes n is an integer greater than one
  (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
    (let loop ((i 0) (p 3) (ps (list 2))) ; sieve of eratosthenes
      (cond ((< n (* p p))
              (do ((i i (+ i 1)) (p p (+ p 2))
                   (ps ps (if (vector-ref bits i) (cons p ps) ps)))
                  ((= i len) (reverse ps))))
            ((vector-ref bits i)
                (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
                    ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
                (vector-set! bits j #f)))
            (else (loop (+ i 1) (+ p 2) ps))))))

; prime? n -- #f if provably composite, else #t if probably prime
(define prime? (let ((ps (primes 100))) (lambda (n) ; integer n
  (define (spsp? n a) ; #f if n is provably composite, else #t
    (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
        ((odd? d)
          (if (= (expm a d n) 1) #t
            (do ((r 0 (+ r 1)))
                ((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
                  (< r s)))))))
  (if (member n ps) #t
    (do ((ps ps (cdr ps)))
        ((or (null? ps) (not (spsp? n (car ps)))) (null? ps)))))))

; factors n -- list of prime factors of n in ascending order
(define (factors n) ; assumes n is an integer, may be negative
  (if (<= -1 n 1) (list n) (if (< n 0) (cons -1 (factors (- n)))
    (let fact ((n n) (c 1) (fs (list))) ; pollard rho method
      (define (f x) (modulo (+ (* x x) c) n))
      (if (even? n) (fact (/ n 2) c (cons 2 fs)) (if (= n 1) fs
        (if (prime? n) (sort < (cons n fs))
          (let loop ((t 2) (h 2) (d 1))
            (if (= d 1) (let ((t (f t)) (h (f (f h))))
                          (loop t h (gcd (- t h) n)))
              (if (= d n) (fact n (+ c 1) fs)
                (if (prime? d) (fact (/ n d) (+ c 1) (cons d fs))
                  (fact n (+ c 1) fs))))))))))))

; project euler three -- largest prime factor of 600851475143
(display (car (reverse (factors 600851475143)))) (newline)

; project euler ten -- sum of primes less than two million
(display (apply + (primes 2000000))) (newline)


Output:
1
2
6857
142913828922


Create a new paste based on this one


Comments: