[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Jul 2:
; big numbers

; Represented as a list with the car containing the signed logarithm of the
; number to the internal base followed in the cdr by the unsigned digits to
; the internal base, with the least significant digit first. Thus, zero is
; (0), one is (1 1), negative one is (-1 1), and the internal base is (2 0 1).

; (integer->big int), (big->integer big), (big? x), (big-abs b),
; (big-negate b), (big-even? b), (big-odd? b), (big-positive? b),
; (big-negative? b), (big-zero? b), (big-sign b), (big-compare b1 b2),
; (big-eq? b1 b2), (big-ne? b1 b2), (big-lt? b1 b2), (big-gt? b1 b2),
; (big-le? b1 b2), (big-ge? b1 b2), (big-plus b1 b2), (big-minus b1 b2),
; (big-times b1 b2), (big-divide b1 b2) => (values quotient remainder),
; (big-quotient b1 b2), (big-remainder b1 b2), (big-modulo b1 b2),
; (big-gcd m n), (big-expt b e), (big-expm b e m), (big-sqrt n),
; (big-rand 'seed n) must be 0 < n < 2^32, (big-rand [first] past),
; (string->big str [radix]), (big->string b [radix])

(define (big-abs big)
  (if (positive? (car big)) big (cons (- (car big)) (cdr big))))
(define (big-negate big) (cons (* (car big) -1) (cdr big)))

(define (big-positive? big) (positive? (car big)))
(define (big-negative? big) (negative? (car big)))
(define (big-zero? big) (zero? (car big)))
(define (big-sign x) (if (big-negative? x) -1 (if (big-positive? x) 1 0)))

(define (big-even? big) (or (big-zero? big) (even? (cadr big))))
(define (big-odd? big) (not (or (big-zero? big) (even? (cadr big)))))

(define (big-compare big1 big2)
  ; big1 < big2 => -1 ; big1 = big2 => 0 ; big1 > big2 => 1
  (define (sign x) (if (negative? x) -1 (if (positive? x) 1 0)))
  (cond ((< (car big1) (car big2)) -1)
        ((< (car big2) (car big1)) 1)
        (else (let loop ((b1 (reverse (cdr big1))) (b2 (reverse (cdr big2))))
                (cond ((null? b1) 0)
                      ((< (car b1) (car b2)) (* (sign (car big1)) -1))
                      ((< (car b2) (car b1)) (* (sign (car big1)) 1))
                      (else (loop (cdr b1) (cdr b2))))))))

(define (big-eq? big1 big2) (zero? (big-compare big1 big2)))
(define (big-ne? big1 big2) (not (zero? (big-compare big1 big2))))
(define (big-lt? big1 big2) (negative? (big-compare big1 big2)))
(define (big-gt? big1 big2) (positive? (big-compare big1 big2)))
(define (big-le? big1 big2) (not (positive? (big-compare big1 big2))))
(define (big-ge? big1 big2) (not (negative? (big-compare big1 big2))))

(define big? #f)
(define big-plus #f)
(define big-minus #f)
(define big-times #f)
(define big-divide #f) ; => (values q r)
(define big-modulo #f)

(let ()
  (define (take n xs)
    (let loop ((n n) (xs xs) (ys '()))
      (if (or (zero? n) (null? xs)) (reverse ys)
        (loop (- n 1) (cdr xs) (cons (car xs) ys)))))
  (define (drop n xs)
    (let loop ((n n) (xs xs))
      (if (or (zero? n) (null? xs)) xs
        (loop (- n 1) (cdr xs)))))
  (define (all? pred? xs)
    (cond ((null? xs) #t)
          ((pred? (car xs)) (all? pred? (cdr xs)))
          (else #f)))
  (define base (expt 2 14))
  (define (sign x) (if (negative? x) -1 (if (positive? x) 1 0)))
  (define (lt? xs ys)
    (let ((xlen (length xs)) (ylen (length ys)))
      (if (< xlen ylen) #t (if (< ylen xlen) #f
        (let loop ((xs (reverse xs)) (ys (reverse ys)))
          (cond ((null? xs) #f) ; equal
                ((< (car xs) (car ys)) #t)
                ((< (car ys) (car xs)) #f)
                (else (loop (cdr xs) (cdr ys)))))))))
  (define (reduce xs)
    (if (null? (cdr xs)) xs
      (if (positive? (car xs)) xs
        (reduce (cdr xs)))))
  (define (add b1 b2)
    (let loop ((b1 b1) (b2 b2) (c 0) (bs '()))
      (cond ((null? b1)
              (if (zero? c) (reverse bs) (reverse (cons 1 bs))))
            ((null? b2)
              (let* ((sum (+ (car b1) c))
                     (new-c (if (<= base sum) 1 0)))
                (loop (cdr b1) b2 new-c
                      (cons (modulo sum base) bs))))
            (else (let* ((sum (+ (car b1) (car b2) c))
                          (new-c (if (<= base sum) 1 0)))
                    (loop (cdr b1) (cdr b2) new-c
                          (cons (modulo sum base) bs)))))))
  (define (sub b1 b2)
    (let loop ((b1 b1) (b2 b2) (c 0) (bs '()))
      (cond ((null? b1) (reverse (reduce bs)))
            ((null? b2)
              (let* ((diff (- (car b1) c))
                     (new-c (if (< diff 0) 1 0)))
                (loop (cdr b1) b2 new-c
                      (cons (modulo diff base) bs))))
            (else (let* ((diff (- (car b1) (car b2) c))
                         (new-c (if (< diff 0) 1 0)))
                    (loop (cdr b1) (cdr b2) new-c
                          (cons (modulo diff base) bs)))))))
  (define (times big1 big2)
    (let loop ((b1 big1) (b2 big2) (zs '())
               (c 0) (ps '()) (bs '()))
      (cond ((null? b1) bs)
            ((null? b2) (let ((zs (cons 0 zs)))
              (loop (cdr b1) big2 zs 0 zs
                (add (reverse (if (zero? c) ps (cons c ps))) bs))))
            (else (let* ((x (+ (* (car b1) (car b2)) c))
                         (c (quotient x base))
                         (p (modulo x base)))
                    (loop b1 (cdr b2) zs c (cons p ps) bs))))))
  (define (mul1 ns d)
    (let loop ((ns ns) (c 0) (ps '()))
      (if (null? ns) (reverse (if (zero? c) ps (cons c ps)))
        (let* ((x (+ (* (car ns) d) c))
               (c (quotient x base))
               (p (modulo x base)))
          (loop (cdr ns) c (cons p ps))))))
  (define (div1 ns d)
    (let loop ((rev-ns (reverse ns)) (qs '()) (r 0))
      (if (null? rev-ns) (values (reverse (reduce (reverse qs))) (list r))
        (let* ((x (+ (* r base) (car rev-ns)))
               (q (quotient x d)) (r (modulo x d)))
          (loop (cdr rev-ns) (cons q qs) r)))))
  (define (nextq x d0 d1 ns)
    (let loop ((q (quotient x d0)))
      (if (< (* q d1) (+ (* (- x (* q d0)) base) (caddr ns))) q
        (loop (- q 1)))))
  (define (nextn j n rev-ns ds*q)
    (let ((zs (append (reverse (sub (reverse (take (+ n 1) rev-ns)) ds*q))
                      (drop (+ n 1) rev-ns))))
      (if (< (length zs) (+ n j)) (cons 0 zs) zs)))
  (define (div ns ds)
    (if (lt? ns ds) (values '(0) ns)
      (let* ((n (length ds)) (m (- (length ns) n)))
        (if (= n 1) (div1 ns (car ds))
          (let* ((d (quotient base (+ (car (reverse ds)) 1)))
                 (rev-ns (reverse (mul1 ns d)))
                 (rev-ns (if (= (length ns) (length rev-ns))
                           (cons 0 rev-ns) rev-ns))
                 (ds (mul1 ds d))
                 (d0 (car (reverse ds)))
                 (d1 (cadr (reverse ds))))
            (let loop ((j m) (rev-ns rev-ns) (qs '()))
              (if (negative? j)
                    (lambda () (div1 (reverse rev-ns) d))
                    (lambda (q r) (values (reverse (reduce (reverse qs))) q)))
                  (let* ((x (+ (* (car rev-ns) base) (cadr rev-ns)))
                         (q (nextq x d0 d1 rev-ns))
                         (ds*q (mul1 ds q)))
                    (if (lt? (reverse (take (+ n 1) rev-ns)) ds*q)
                        (let* ((q (- q 1)) (ds*q (sub ds*q ds)))
                          (loop (- j 1) (nextn j n rev-ns ds*q) (cons q qs)))
                        (loop (- j 1) (nextn j n rev-ns ds*q) (cons q qs)))))))))))
  (set! big? (lambda (xs)
    (and (list? xs)
         (all? integer? xs)
         (all? (lambda (x) (< -1 x base)) xs)
         (= (abs (car xs)) (length (cdr xs))))))
  (set! big-plus (lambda (big1 big2)
    (if (zero? (car big1)) big2
      (if (zero? (car big2)) big1
        (let* ((b1 (cdr big1)) (b2 (cdr big2))
               (lt? (big-lt? (big-abs big1) (big-abs big2)))
               (s1 (if (positive? (car big1)) 1 -1))
               (s2 (if (positive? (car big2)) 1 -1))
               (x (apply (if (= s1 s2) add sub)
                         (if lt? (list b2 b1) (list b1 b2))))
               (len (length x)))
          (if (equal? x (list 0)) x
            (cons (* len (if (or (and lt? (= s2 1))
                             (and (not lt?) (= s1 1)))
                           1 -1))
  (set! big-minus (lambda (big1 big2)
    (big-plus big1 (big-negate big2))))
  (set! big-times (lambda (big1 big2)
    (if (or (big-zero? big1) (big-zero? big2)) (list 0)
      (let* ((b1 (cdr big1)) (b2 (cdr big2))
             (x (times b1 b2)) (len (length x)))
        (cons (* len (sign (* (car big1) (car big2)))) x)))))
  (set! big-divide (lambda (ns ds)
    (if (big-zero? ds) (error 'big-divide "divide by zero")
      (let ((sn (sign (car ns))) (sd (sign (car ds))))
          (lambda () (div (cdr (big-abs ns)) (cdr (big-abs ds))))
          (lambda (qs rs)
            (if (equal? qs '(0))
                (values '(0) (cons (* sn (length rs)) rs))
                (values (cons (* (if (= sn sd) 1 -1) (length qs)) qs)
                        (cons (* sn (length rs)) rs)))))))))
  (set! big-modulo (lambda (ns ds)
    (if (big-zero? ds) (error 'big-modulo "divide by zero")
      (let ((sn (sign (car ns))) (sd (sign (car ds))))
          (lambda () (div (cdr (big-abs ns)) (cdr (big-abs ds))))
          (lambda (qs rs)
            (if (equal? rs '(0)) rs
              (let ((r (cons (length rs) rs)))
                (if (= sn sd) r (big-plus r base)))))))))))

(define (big-quotient ns ds)
  (call-with-values (lambda () (big-divide ns ds)) (lambda (qs rs) qs)))
(define (big-remainder ns ds)
  (call-with-values (lambda () (big-divide ns ds)) (lambda (qs rs) rs)))

(define (string->big str . args)
  (define (char->digit c)
    (cond ((char-numeric? c) (- (char->integer c) 48))
          ((char-upper-case? c) (- (char->integer c) 55))
          ((char-lower-case? c) (- (char->integer c) 87))
          (else (error 'char->digit "illegal character"))))
  (if (string=? str "") '(0)
    (let* ((radix (list 1 (if (null? args) 10 (car args))))
           (s (if (char=? (string-ref str 0) #\-) -1 1))
           (ds (string->list str))
           (ds (if (positive? s) ds (cdr ds)))
           (ds (map char->digit ds)))
      (let loop ((ds (cdr ds)) (big (list 1 (car ds))))
        (if (null? ds)
            (if (equal? big '(1 0)) '(0) (cons (* s (car big)) (cdr big)))
            (loop (cdr ds)
                  (big-plus (big-times big radix)
                            (list 1 (car ds)))))))))

(define (big->string big . args)
  (define (digit->char d)
    (cond ((< d 10) (integer->char (+ d 48)))
          ((< d 36) (integer->char (+ d 55)))
          (else (error 'digit->char "illegal digit"))))
  (if (big-zero? big) "0"
    (if (big-negative? big) (string-append "-" (big->string (big-negate big)))
      (let ((radix (list 1 (if (null? args) 10 (car args)))))
        (let loop ((big big) (ds '()))
          (if (big-zero? big) (list->string ds)
              (lambda () (big-divide big radix))
              (lambda (q r) (loop q (cons (digit->char (cadr r)) ds))))))))))

(define (big-gcd m n)
  (let loop ((m (big-abs m)) (n (big-abs n)))
    (cond ((big-lt? m n) (loop m (big-minus n m)))
          ((big-lt? n m) (loop (big-minus m n) n))
          (else m))))

(define (big-expt b e)
  (if (big-zero? e) '(1 1)
    (let loop ((s b) (i e) (a '(1 1)))
      (let ((a (if (big-odd? i) (big-times a s) a))
            (i (big-quotient i '(1 2))))
        (if (big-zero? i) a (loop (big-times s s) i a))))))

(define (big-expm b e m)
  (define (m* x y) (big-modulo (big-times x y) m))
  (cond ((big-zero? e) '(1 1))
        ((big-even? e) (big-expm (m* b b) (big-quotient e '(1 2)) m))
        (else (m* b (big-expm (m* b b) (big-quotient (big-minus e '(1 1)) '(1 2)) m)))))

(define (big-sqrt n)
  (if (not (big-positive? n))
      (error 'big-sqrt "must be positive")
      (let loop ((x (cons (quotient (+ (car n) 1) 2) (cdr n))))
        (let ((y (big-quotient (big-plus x (big-quotient n x)) '(1 2))))
          (if (big-lt? y x) (loop y) x)))))

(define big-rand ; Park/Miller algorithm
  (let ((a (string->big "16807"))
        (m (string->big "2147483647"))
        (s (string->big "1043618065")))
    (lambda args
      (if (eq? (car args) 'seed)
          (begin (set! s (cadr args)) s)
          (let ((first (if (pair? (cdr args)) (car args) '(0)))
                (past (if (pair? (cdr args)) (cadr args) (car args))))
            (set! s (big-modulo (big-times a s) m))
            (big-plus first (big-quotient (big-times s (big-minus past first)) m)))))))

(define integer->big
  (let ((base (string->number (big->string '(2 0 1)))))
    (lambda (n)
      (if (zero? n) '(0)
        (if (negative? n) (big-negate (integer->big (- n)))
          (let loop ((n n) (bs '()))
            (if (zero? n) (cons (length bs) (reverse bs))
              (loop (quotient n base) (cons (modulo n base) bs)))))))))

(define big->integer
  (let ((base (string->number (big->string '(2 0 1)))))
    (lambda (b)
      (if (big-zero? b) 0
        (if (big-negative? b) (* -1 (big->integer (big-negate b)))
          (let loop ((bs (reverse (cdr b))) (n 0))
            (if (null? bs) n
              (loop (cdr bs) (+ (car bs) (* n base))))))))))

; examples

(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)
        ((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 (big-factorial n)
  (let loop ((n (- n 1)) (k '(1 1)) (f '(1 1)))
    (display (big->string k)) (display " ")
    (display (big->string f)) (newline)
    (when (positive? n)
      (let ((k (big-plus k '(1 1))))
        (loop (- n 1) k (big-times f k))))))

(big-factorial 50) (newline)

(define (big-factors n) ; trial division
  (if (big-even? n)
      (cons '(1 2) (factors (big-quotient n '(1 2))))
      (let loop ((n n) (f '(1 3)) (fs '()))
        (cond ((big-lt? n (big-times f f))
                (reverse (cons n fs)))
              ((big-zero? (big-modulo n f))
                (loop (big-quotient n f) f (cons f fs)))
              (else (loop n (big-plus f '(1 2)) fs))))))

(display (map big->integer (big-factors (integer->big 11111111111111111111))))

(define (big-prime? n)
  (define (witness? a n)
    (let loop ((r '(0)) (s (big-minus n '(1 1))))
      (if (big-even? s)
          (loop (big-plus r '(1 1)) (big-quotient s '(1 2)))
          (if (big-eq? (big-expm a s n) '(1 1)) #t
            (let loop ((j '(0)) (s s))
              (cond ((big-eq? j r) #f)
                    ((big-eq? (big-expm a s n) (big-minus n '(1 1))) #t)
                    (else (loop (big-plus j '(1 1)) (big-times s '(1 2))))))))))
  (cond ((big-lt? n '(1 2)) #f)
        ((big-even? n) (big-eq? n '(1 2)))
        (else (let loop ((k 20))
                (cond ((zero? k) #t)
                      ((not (witness? (big-rand '(1 1) n) n)) #f)
                      (else (loop (- k 1))))))))

(display (big-prime? (integer->big (- (expt 2 89) 1)))) (newline)

(define (big-rho n)
  (define (factor n c)
    (define (f x) (big-modulo (big-plus (big-times x x) c) n))
    (let loop ((x '(1 2)) (y '(1 2)) (d '(1 1)))
      (cond ((big-eq? d '(1 1))
              (let ((x (f x)) (y (f (f y))))
                (loop x y (big-gcd (big-minus x y) n))))
            ((big-eq? d n) (factor n (big-plus c '(1 1))))
            (else d))))
  (sort big-lt?
    (let fact ((n n) (fs '()))
      (cond ((big-eq? n '(1 1)) fs)
            ((big-even? n) (fact (big-quotient n '(1 2)) (cons '(1 2) fs)))
            ((big-prime? n) (cons n fs))
            (else (let ((f (factor n '(1 1))))
              (append fs (fact f '()) (fact (big-quotient n f) '()))))))))

(display (map big->integer (big-rho (integer->big 13290059)))) (newline)

1 1
2 2
3 6
4 24
5 120
6 720
7 5040
8 40320
9 362880
10 3628800
11 39916800
12 479001600
13 6227020800
14 87178291200
15 1307674368000
16 20922789888000
17 355687428096000
18 6402373705728000
19 121645100408832000
20 2432902008176640000
21 51090942171709440000
22 1124000727777607680000
23 25852016738884976640000
24 620448401733239439360000
25 15511210043330985984000000
26 403291461126605635584000000
27 10888869450418352160768000000
28 304888344611713860501504000000
29 8841761993739701954543616000000
30 265252859812191058636308480000000
31 8222838654177922817725562880000000
32 263130836933693530167218012160000000
33 8683317618811886495518194401280000000
34 295232799039604140847618609643520000000
35 10333147966386144929666651337523200000000
36 371993326789901217467999448150835200000000
37 13763753091226345046315979581580902400000000
38 523022617466601111760007224100074291200000000
39 20397882081197443358640281739902897356800000000
40 815915283247897734345611269596115894272000000000
41 33452526613163807108170062053440751665152000000000
42 1405006117752879898543142606244511569936384000000000
43 60415263063373835637355132068513997507264512000000000
44 2658271574788448768043625811014615890319638528000000000
45 119622220865480194561963161495657715064383733760000000000
46 5502622159812088949850305428800254892961651752960000000000
47 258623241511168180642964355153611979969197632389120000000000
48 12413915592536072670862289047373375038521486354677760000000000
49 608281864034267560872252163321295376887552831379210240000000000
50 30414093201713378043612608166064768844377641568960512000000000000

(11 41 101 271 3541 9091 27961)
(3119 4261)

Create a new paste based on this one