; 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).

; (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),
; (big-quotient b1 b2), (big-remainder b1 b2), (big-modulo b1 b2),
; (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)
                  (call-with-values
                    (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))
                  x)))))))
  (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))))
        (call-with-values
          (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))))
        (call-with-values
          (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)
            (call-with-values
              (lambda () (big-divide big radix))
              (lambda (q r) (loop q (cons (digit->char (cadr r)) ds))))))))))

;;;;; testing

(define-syntax assert
  (syntax-rules ()
    ((assert expr result)
      (if (not (equal? expr result))
          (for-each display `(
            #\newline "failed assertion:" #\newline
            expr #\newline "expected: " ,result
            #\newline "returned: " ,expr #\newline))))))

(define next ; park-miller generator, suitable only for testing
  (let ((a (string->big "16807"))
        (m (string->big "2147483647"))
        (seed (string->big "1043618065")))
    (lambda (n)
      (let* ((minus-n (big-negate n))
             (span (big-plus (big-times n '(1 2)) '(1 1))))
        (set! seed (big-modulo (big-times a seed) m))
        (big-plus minus-n (big-quotient (big-times seed span) m))))))

(define (test-big) ; no news is good news
  (let ((zero '(0)) (one '(1 1)) (minus-one '(-1 1))
        (base '(2 0 1)) (minus-base '(-2 0 1))
        (b1 (string->big "737")) (b2 (string->big "13290059"))
        (ten25 (string->big "10000000000000000000000000"))
        (ten24 (string->big "1000000000000000000000000")))
    (assert (string->big "0") zero)
    (assert (big->string (string->big "0")) "0")
    (assert (big->string (string->big "1")) "1")
    (assert (big->string (string->big "2")) "2")
    (assert (big->string (string->big "-1")) "-1")
    (assert (big->string (string->big "42")) "42")
    (assert (big->string (string->big "-42")) "-42")
    (assert (big->string (string->big "737")) "737")
    (assert (big->string (string->big "13290059")) "13290059")
    (assert (big-abs one) one)
    (assert (big-abs minus-one) one)
    (assert (big-negate one) minus-one)
    (assert (big-negate minus-one) one)
    (assert (big-negate base) minus-base)
    (assert (big->string (big-negate b1)) "-737")
    (assert (big->string (big-negate b2)) "-13290059")
    (assert (big-zero? zero) #t)
    (assert (big-zero? b2) #f)
    (assert (big-zero? minus-one) #f)
    (assert (big-positive? zero) #f)
    (assert (big-positive? b2) #t)
    (assert (big-positive? minus-one) #f)
    (assert (big-negative? zero) #f)
    (assert (big-negative? b2) #f)
    (assert (big-negative? minus-one) #t)
    (assert (big-even? zero) #t)
    (assert (big-odd? zero) #f)
    (assert (big-even? b2) #f)
    (assert (big-odd? b2) #t)
    (assert (big-even? minus-base) #t)
    (assert (big-odd? minus-base) #f)
    (assert (big-lt? b1 b2) #t)
    (assert (big-lt? b1 (big-negate b2)) #f)
    (assert (big-le? b1 b1) #t)
    (assert (big-le? b2 b1) #f)
    (assert (big-lt? zero one) #t)
    (assert (big-le? zero one) #t)
    (assert (big-lt? b1 b1) #f)
    (assert (big-eq? b1 b1) #t)
    (assert (big-eq? b1 b2) #f)
    (assert (big-ne? b1 b2) #t)
    (assert (big-ne? b1 b1) #f)
    (assert (big-gt? b1 b2) #f)
    (assert (big-gt? b1 (big-negate b2)) #t)
    (assert (big-ge? b1 b1) #t)
    (assert (big-ge? b2 b1) #t)
    (assert (big-gt? zero one) #f)
    (assert (big-ge? zero one) #f)
    (assert (big-gt? b1 b1) #f)
    (do ((i 0 (+ i 1))
         (n (next ten25) (next ten25))
         (d (next ten24) (next ten24)))
        ((= i 1000))
      (call-with-values
        (lambda () (big-divide n d))
        (lambda (q r)
          (if (big-negative? n)
              (begin (assert (big-le? r '(0)) #t)
                     (assert (big-lt? (big-negate (big-abs d)) r) #t))
              (begin (assert (big-ge? r '(0)) #t)
                     (assert (big-lt? r (big-abs d)) #t)))
          (assert (big-minus (big-times q d) (big-negate r)) n))))))

(test-big)