path: root/lib/scheme/base/80-number.csc

(export
  *
  +
  -
  /
  <
  <=
  =
  >
  >=
  abs
  ceiling
  denominator
  even?
  exact-integer-sqrt
  expt
  floor
  floor-quotient
  floor-remainder
  floor/
  gcd
  integer?
  lcm
  max
  min
  modulo
  negative?
  number?
  numerator
  odd?
  positive?
  quotient
  rational?
  remainder
  round
  square
  truncate
  truncate-quotient
  truncate-remainder
  truncate/
  zero?)
(import (only (csc builtins)
          call-builtin))
(begin


  ;; Integers


  (define-record-type <boxed-int>
    (make-boxed-int positive? digits)
    boxed-int?
    (positive? boxed-int-positive?)
    (digits boxed-int-digits))


  (define (small-int? obj)
    (call-builtin eq 1 (call-builtin typeof obj)))


  (define (integer? obj)
    (or (small-int? obj)
        (boxed-int? obj)))


  (define (int=? n1 n2)
    (cond
      ((and (small-int? n1) (small-int? n2))
        (call-builtin eq n1 n2))
      ((and (boxed-int? n1) (boxed-int? n2))
        (let ((n1-digits (boxed-int-digits n1))
              (n2-digits (boxed-int-digits n2)))
          (and (boolean=? (boxed-int-positive? n1) (boxed-int-positive? n2))
               (call-builtin eq (vector-length n1-digits) (vector-length n2-digits))
               (let loop ((i 0))
                 (if (call-builtin lt i (vector-length n1-digits))
                   (and (call-builtin eq (vector-ref n1-digits i) (vector-ref n2-digits i))
                        (loop (+ 1 i)))
                   #t)))))
      (else #f)))


  (define (int-positive? n)
    (or (and (small-int? n)
             (call-builtin lt 0 n))
        (and (boxed-int? n)
             (boxed-int-positive? n))))


  (define (int-negative? n)
    (or (and (small-int? n)
             (call-builtin lt n 0))
        (and (boxed-int? n)
             (not (boxed-int-positive? n)))))


  (define (int<? n1 n2)
    (cond
      ((and (small-int? n1) (small-int? n2))
        (call-builtin lt n1 n2))
      ((and (int-negative? n1) (int-negative? n2))
        (int<? (- n2) (- n1)))
      ((int-negative? n1) #t)
      ((int-negative? n2) #f)
      ((small-int? n1) #t)
      ((small-int? n2) #f)
      (else
        (let* ((n1-digits (boxed-int-digits n1))
               (n2-digits (boxed-int-digits n2))
               (n1-digits-len (vector-length n1-digits))
               (n2-digits-len (vector-length n2-digits)))
          (or (call-builtin lt n1-digits-len n2-digits-len)
              (and (int=? n1-digits-len n2-digits-len)
                   (let loop ((i (call-builtin sub n1-digits-len 1)))
                     (cond
                       ((call-builtin lt i 0) #f) ; n1 is equal to n2
                       ((call-builtin lt (vector-ref n1-digits i) (vector-ref n2-digits i)) #t)
                       ((int=? (vector-ref n1-digits i) (vector-ref n2-digits i))
                         (loop (- i 1)))
                       (else #f))))))))) ; n1 > n2


  (define (int>? n1 n2)
    (int<? n2 n1))


  (define (int<=? n1 n2)
    (or (int=? n1 n2)
        (int<? n1 n2)))


  (define (int>=? n1 n2)
    (or (int=? n1 n2)
        (int>? n1 n2)))


  (define (zero? obj)
    (call-builtin eq 0 obj))


  (define (odd? n)
    (cond
      ((small-int? n)
        (int=? 1 (call-builtin mod n 2)))
      (else
        (odd? (vector-ref (boxed-int-digits n) 0)))))


  (define (even? n)
    (cond
      ((small-int? n)
        (int=? 0 (call-builtin mod n 2)))
      (else
        (even? (vector-ref (boxed-int-digits n) 0)))))


  (define small-int-max #x3FFFFFFFFFFFFFFF)  ; 2^62 - 1


  (define small-int-min (- #x4000000000000000))  ; -(2^62)


  (define (digits n)
    (if (small-int? n)
      (vector n)
      (boxed-int-digits n)))


  ; adds two small positive ints.
  (define (add2 x y)
    (define z (call-builtin add x y))
    (if (negative? z)
      (values 1 (call-builtin add
                  1
                  (call-builtin add
                    z
                    small-int-max)))
      (values 0 z)))


  (define (big-int+ x y)
    (cond
      ((and (int-negative? x) (int-negative? y))
        (- (int+ (- x) (- y))))
      ((int-negative? x)
        (int- y (- x)))
      ((int-negative? y)
        (int- x (- y)))
      (else
        (let* ((x-digits (digits x))
               (y-digits (digits y))
               (x-ndigits (vector-length x-digits))
               (y-ndigits (vector-length y-digits))
               (out-digits (make-vector (max x-ndigits
                                             y-ndigits))))
          (let loop ((i 0)
                     (carry 0))
            (cond
              ((and (int<? i x-ndigits)
                    (int<? i y-ndigits))
                (let*-values (((c1 d1) (add2
                                         (vector-ref x-digits i)
                                         (vector-ref y-digits i)))
                              ((c2 d2) (add2
                                         d1
                                         carry)))
                  (vector-set! out-digits i d2)
                  (loop (call-builtin add 1 i) (call-builtin add c1 c2))))
              ((int<? i x-ndigits)
                (let-values (((c d) (add2
                                      carry
                                      (vector-ref x-digits i))))
                  (vector-set! out-digits i d)
                  (loop (call-builtin add 1 i) c)))
              ((int<? i y-ndigits)
                (let-values (((c d) (add2
                                      carry
                                      (vector-ref y-digits i))))
                  (vector-set! out-digits i d)
                  (loop (call-builtin add 1 i) c)))
              ((int-positive? carry)
                (set! out-digits (vector-append out-digits (vector carry))))))
          (make-boxed-int #t out-digits)))))


  (define (int+ x y)
    (cond
      ((and (small-int? x)
            (small-int? y))
        (let ((z (call-builtin add x y))
              (x-positive (int-positive? x)))
          (if (and (boolean=? x-positive (int-positive? y))
                   (not (boolean=? x-positive (int-positive? z))))  ; overflow
            (big-int+ x y)
            z)))
      (else (big-int+ x y))))


  (define (remove-leading-zeros n)
    (define digits (boxed-int-digits n))
    (define leading-zeros
      (let loop ((i (int- (vector-length digits) 1))
                 (n 0))
        (if (and (int>=? i 0)
                 (zero? (vector-ref digits i)))
          (loop (int- i 1) (int+ 1 n))
          n)))
    (make-boxed-int (boxed-int-positive? n) (vector-copy digits 0 (int- (vector-length digits) leading-zeros))))


  (define (normalize n)
    (set! n (remove-leading-zeros n))
    (define digits (boxed-int-digits n))
    (define len (vector-length digits))
    (cond
      ((call-builtin eq 0 len)
        0)
      ((call-builtin eq 1 len)
        (vector-ref digits 0))
      (else n)))


  (define (big-int- n1 n2)
    (cond
      ((and (int-negative? n1) (int-negative? n2))
        (- (big-int- (- n1) (- n2))))
      ((and (int-positive? n1) (int-negative? n2))
        (big-int+ n1 (- n2)))
      ((and (int-negative? n1) (int-positive? n2))
        (- (big-int+ (- n1) n2)))
      ((int<? n1 n2)
        (- (big-int- n2 n1)))
      (else
        (let* ((n1-digits (digits n1))
               (n2-digits (digits n2))
               (n1-ndigits (vector-length n1-digits))
               (n2-ndigits (vector-length n2-digits))
               (out-digits (make-vector n1-ndigits)))  ; n.b.: n1 is bigger
          (let loop ((i 0)
                     (carry 0))
            (cond 
              ((and (int<? i n1-ndigits)
                    (int<? i n2-ndigits))
                (let ((x (call-builtin add
                           carry
                           (call-builtin sub
                             (vector-ref n1-digits i)
                             (vector-ref n2-digits i)))))
                  (if (int-negative? x)
                    (begin
                      (vector-set! out-digits i
                        (call-builtin add
                          1
                          (call-builtin add x small-int-max)))
                      (loop (+ 1 i) -1))
                    (begin
                      (vector-set! out-digits i x)
                      (loop (+ 1 i) 0)))))
              ((int<? i n1-ndigits)
                (let ((x (call-builtin add
                           carry
                           (vector-ref n1-digits i))))
                  (if (int-negative? x)
                    (begin
                      (vector-set! out-digits i
                        (call-builtin add
                          1
                          (call-builtin add x small-int-max)))
                      (loop (+ 1 i) -1))
                    (begin
                      (vector-set! out-digits i x)
                      (loop (+ 1 i) 0)))))))

          (normalize (make-boxed-int #t out-digits))))))


  (define (int- n1 n2)
    (cond
      ((and (small-int? n1) (small-int? n2))
        (let ((x (call-builtin sub n1 n2))
              (n1-positive (int-positive? n1)))
          (if (or (boolean=? n1-positive (int-positive? n2))
                  (boolean=? n1-positive (int-positive? x)))
            x
            (big-int- n1 n2))))
      (else
        (big-int- n1 n2))))


  (define half-word-mask #x80000000)


  ; multiplies two small positive ints.
  (define (mul2 x y)
    ; Split each argument into half-words.
    (define x0 (call-builtin mod x half-word-mask))
    (define x1 (call-builtin div x half-word-mask))
    (define y0 (call-builtin mod y half-word-mask))
    (define y1 (call-builtin div y half-word-mask))

    ; Do the grade-school multiplication algorithm.
    (define z0 (call-builtin mul x0 y0))
    (define-values (z1-hi z1-lo)
      (add2 (call-builtin mul x1 y0)
            (call-builtin mul x0 y1)))
    (define z2 (call-builtin mul x1 y1))

    ; Shift z1 left by a half-word.
    (define z1*-lo (call-builtin mul
                     half-word-mask
                     (call-builtin mod z1-lo half-word-mask)))
    (define z1*-hi (call-builtin add
                     (call-builtin mul z1-hi half-word-mask)
                     (call-builtin div z1-lo half-word-mask)))

    ; Add z0 + z1* + z2*word-size
    (define-values (carry result-lo) (add2 z0 z1*-lo))
    (values
      (call-builtin add
        carry
        (call-builtin add
          z1*-hi
          z2))
      result-lo))


  (define (split-in-half n m)
    (define n-digits (digits n))
    (if (int<? m (vector-length n-digits))
      (values
        (normalize (make-boxed-int #t (vector-copy n-digits m)))
        (normalize (make-boxed-int #t (vector-copy n-digits 0 m))))
      (values
        0
        n)))


  (define (lshift n m)
    (let* ((old-digits (digits n))
           (new-digits (make-vector (int+ m (vector-length old-digits)))))
      (vector-fill! new-digits 0 0 m)
      (vector-copy! new-digits m old-digits)
      (normalize (make-boxed-int #t new-digits))))


  ; This is Karatsuba's algorithm.
  (define (big-int* x y)
    (define m (call-builtin div
                (max (vector-length (digits x)) (vector-length (digits y)))
                2))
    (define-values (x1 x0) (split-in-half x m))
    (define-values (y1 y0) (split-in-half y m))
    (define z0 (* x0 y0))
    (define z2 (* x1 y1))
    (define z1 (- (* (+ x1 x0)
                     (+ y1 y0))
                  z2
                  z0))
    (define ans
    (+ (lshift z2 (call-builtin mul 2 m))
       (lshift z1 m)
       z0))
    ans)


  (define (slow-int* x y)
    (cond
      ((and (int-negative? x) (int-negative? y))
        (int* (- x) (- y)))
      ((int-negative? x)
        (- (int* (- x) y)))
      ((int-negative? y)
        (- (int* x (- y))))
      ((and (small-int? x) (small-int? y))
        (let-values (((z1 z0) (mul2 x y)))
          (make-boxed-int #t (vector z0 z1))))
      (else
        (big-int* x y))))


  (define (int* n1 n2)
    (cond
      ((and (small-int? n1) (small-int? n2))
        (let ((x (call-builtin mul n1 n2)))
          (if (and (not (zero? n1))
                   (not (int=? (call-builtin div x n1)
                               n2))) ; overflow
            (slow-int* n1 n2)
            x)))
      (else
        (slow-int* n1 n2))))


  (define (left-index v i)
    (vector-ref v (- (vector-length v) 1 i)))


  (define (find-beta d m)
    (let loop ((lo 0)
               (hi small-int-max))
      (define guess (+ lo (call-builtin div (- hi lo) 2)))
      (define check (- d (* m guess)))
      (cond
        ((int>? lo hi)
          (error "binary search is hard"))
        ((negative? check)
          ; guess was too big
          (loop lo (- guess 1)))
        ((int<? check m)
          ; got it!
          guess)
        ((int<=? lo hi)
          ; guess was too small
          (loop (+ guess 1) hi))
        (else (error "binary search is broken" d m)))))


  (define (big-int/ n m)
    (define n-digits (digits n))
    (define m-digits (digits m))
    (define k (vector-length n-digits))
    (define l (vector-length m-digits))
    (if (int<? k l)
      (values 0 n)
      (let loop ((i (- l 1))
                 (q 0)
                 (r (normalize
                      (make-boxed-int #t (vector-copy n-digits (- k (- l 1)))))))  ; last l-1 digits of n
        (if (int<? i k)
          (let* ((d (+ (lshift r 1)
                       (left-index n-digits i)))
                 (beta (find-beta d m)))
            (loop (+ 1 i)
                  (+ (lshift q 1)
                     beta)
                  (- d (* m beta))))
          (values q r)))))


  (define (truncate/ n1 n2)
    (cond
      ((and (small-int? n1) (small-int? n2))
        (if (= -1 n2)
          ; To avoid dividing INT_MIN by -1,
          ; we just convert every division by -1 into a negation.
          (values (- n1) 0)
          (values
            (call-builtin div n1 n2)
            (call-builtin mod n1 n2))))
      ((and (int-negative? n1) (int-negative? n2))
        (let-values (((q r) (truncate/ (- n1) (- n2))))
          (values q (- r))))
      ((int-negative? n1)
        (let-values (((q r) (truncate/ (- n1) n2)))
          (values (- q) (- r))))  ; n.b.: remainder takes its sign from n1.
      ((int-negative? n2)
        (let-values (((q r) (truncate/ n1 (- n2))))
          (values (- q) r)))
      (else
        (big-int/ n1 n2))))


  (define (floor/ x y)
    (define-values (q r) (truncate/ x y))
    (cond
      ((or (zero? r)
           (boolean=? (positive? x) (positive? y)))
        (values q r))
      ((negative? x)
        (values (- q 1) (+ r y)))
      (else ; negative y
        (values (- q 1) (- r y)))))


  (define (floor-quotient n1 n2)
    (define-values (q r) (floor/ n1 n2))
    q)


  (define (floor-remainder n1 n2)
    (define-values (q r) (floor/ n1 n2))
    r)


  (define (truncate-quotient n1 n2)
    (define-values (q r) (truncate/ n1 n2))
    q)


  (define (truncate-remainder n1 n2)
    (define-values (q r) (truncate/ n1 n2))
    r)


  (define quotient truncate-quotient)


  (define remainder truncate-remainder)


  (define modulo floor-remainder)


  ; Credit to Euclid for this one.
  (define (gcd2 a b)
    (if (zero? b)
      (abs a)
      (gcd2 b (remainder a b))))


  (define gcd
    (case-lambda
      (() 0)
      ((n . ns)
        (let loop ((n n)
                   (ns ns))
          (if (null? ns)
            n
            (loop (gcd2 n (car ns))
                  (cdr ns)))))))


  (define lcm
    (case-lambda
      (() 1)
      ((n . ns)
        (define ns* (cons n ns))
        (quotient (abs (apply * ns*))
                  (expt (apply gcd ns*) (- (length ns*) 1))))))


  ;; Rationals


  (define-record-type <quotient>
    (make-quotient numerator denominator)
    quotient?
    (numerator quotient-numerator)
    (denominator quotient-denominator))


  (define (rational? obj)
    (or (integer? obj)
        (quotient? obj)))


  (define number? rational?)  ; only rational numbers for now.


  (define (numerator q)
    (cond
      ((integer? q) q)
      (else (quotient-numerator q))))


  (define (denominator q)
    (cond
      ((integer? q) 1)
      (else (quotient-denominator q))))


  (define (= z1 z2 . zs)
    (cond
      ((integer? z1)
        (let loop ((zs (cons z2 zs)))
          (or (null? zs)
              (and (int=? z1 (car zs))
                   (loop (cdr zs))))))
      (else
        (let ((n (numerator z1))
              (d (denominator z1)))
          (let loop ((zs (cons z2 zs)))
            (or (null? zs)
                (and (int=? n (numerator (car zs)))
                     (int=? d (denominator (car zs)))
                     (loop (cdr zs)))))))))


  (define (<2 x1 x2)
    (if (and (integer? x1) (integer? x2))
      (int<? x1 x2)
      (let* ((d1 (denominator x1))
             (d2 (denominator x2))
             (> (gcd d1 d2)))
        (int<? (* (numerator x1)
                  (quotient d2 >))
               (* (numerator x2)
                  (quotient d1 >))))))


  (define <
    (case-lambda
      ((x1 x2)
        (<2 x1 x2))
      ((x1 x2 . xs)
        (and (<2 x1 x2)
             (apply < x2 xs)))))


  (define (> x1 x2 . xs)
    (let loop ((x1 x1)
               (xs (cons x2 xs)))
      (or (null? xs)
          (and (<2 (car xs) x1)
               (loop (car xs)
                     (cdr xs))))))


  (define (<= x1 x2 . xs)
    (let loop ((x1 x1)
               (xs (cons x2 xs)))
      (or (null? xs)
          (and (or (= x1 (car xs))
                   (< x1 (car xs)))
               (loop (car xs) (cdr xs))))))


  (define (>= x1 x2 . xs)
    (let loop ((x1 x1)
               (xs (cons x2 xs)))
      (or (null? xs)
          (and (or (= x1 (car xs))
                   (> x1 (car xs)))
               (loop (car xs) (cdr xs))))))


  (define (positive? x)
    (> x 0))


  (define (negative? x)
    (< x 0))


  (define (max x1 . xs)
    (let loop ((xs xs)
               (m x1))
      (if (null? xs)
        m
        (loop (cdr xs)
              (if (> (car xs) m)
                (car xs)
                m)))))


  (define (min x1 . xs)
    (let loop ((xs xs)
               (m x1))
      (if (null? xs)
        m
        (loop (cdr xs)
              (if (< (car xs) m)
                (car xs)
                m)))))


  (define +
    (case-lambda
      ((x1 x2)
        (cond
          ((and (integer? x1) (integer? x2))
            (int+ x1 x2))
          (else
            (let* ((d1 (denominator x1))
                   (d2 (denominator x2))
                   (> (gcd d1 d2))
                   (s1 (quotient d2 >))
                   (s2 (quotient d1 >)))
              (/ (+ (* s1 (numerator x1))
                    (* s2 (numerator x2)))
                 (* d1 s1))))))
      (xs
        (let loop ((xs xs)
                   (s 0))
          (if (null? xs)
            s
            (loop (cdr xs)
                  (+ s (car xs))))))))


  (define *
    (case-lambda
      ((x1 x2)
        (cond
          ((and (integer? x1) (integer? x2))
            (int* x1 x2))
          (else
            (/ (* (numerator x1) (numerator x2))
               (* (denominator x1) (denominator x2))))))
      (xs
        (let loop ((xs xs)
                   (p 1))
          (if (null? xs)
            p
            (loop (cdr xs)
                  (* p (car xs))))))))


  (define -
    (case-lambda
      ((z)
        (cond
          ((int=? z small-int-min)
            #x4000000000000000)
          ((small-int? z)
            (call-builtin sub 0 z))
          ((boxed-int? z)
            (make-boxed-int (not (positive? z)) (boxed-int-digits z)))
          (else
            (/ (- (numerator z))
               (denominator z)))))
      ((z1 z2)
        (cond
          ((and (integer? z1) (integer? z2))
            (int- z1 z2))
          (else (+ z1 (- z2)))))
      ((z1 . zs)
        (let loop ((zs zs)
                   (d z1))
          (if (null? zs)
            d
            (loop (cdr zs)
                  (- d (car zs))))))))


  (define (normalize-quotient q)
    (when (negative? q)
      (set! q (make-quotient (- (numerator q))
                             (- (denominator q)))))
    (let* ((n (numerator q))
           (d (denominator q))
           (> (gcd n d)))
      (make-quotient (quotient n >)
                     (quotient d >))))


  (define /
    (case-lambda
      ((z)
        (/ 1 z))
      ((z1 z2)
        (cond
          ((and (integer? z1) (integer? z2))
            (when (negative? z2)
              (set! z1 (- z1))
              (set! z2 (- z2)))
            (let ((g (gcd z1 z2)))
              (if (= g z2)
                (quotient z1 g)
                (make-quotient (quotient z1 g)
                               (quotient z2 g)))))
          (else
            (* z1
               (/ (denominator z2)
                  (numerator z2))))))
      ((z1 . zs)
        (let loop ((zs zs)
                   (q z1))
          (if (null? zs)
            q
            (loop (cdr zs)
                  (/ q (car zs))))))))


  (define (abs x)
    (if (negative? x)
      (- x)
      x))


  (define (floor x)
    (if (integer? x)
      x
      (floor-quotient (numerator x)
                      (denominator x))))


  (define (ceiling x)
    (if (integer? x)
      x
      (+ 1 (floor x))))


  (define (truncate x)
    (if (integer? x)
      x
      (truncate-quotient (numerator x)
                         (denominator x))))


  (define (round x)
    (define d (denominator x))
    (define-values (q r) (floor-quotient (numerator x)
                                         d))
    (cond
      ((< (* 2 r) d)
        q)
      ((= (* 2 r) d) ; round to even
        (if (even? q)
          q
          (+ 1 q)))
      (else
        (+ 1 q))))


  (define (square z)
    (* z z))


  (define (exact-integer-sqrt k)
    (if (zero? k)
      (values 0 0)
      (let loop ((x (quotient k 2)))  ; initial estimate
        (define y (quotient (+ x (quotient k x))
                            2))
        (if (>= y x)
          (values x (- k (square x)))
          (loop y)))))


  (define (expt x y)
    (unless (integer? y)
      (error "only integer exponents are supported for now"))
    (cond
      ((zero? y) 1)
      ((even? y)
        (expt (square x) (quotient y 2)))
      (else
        (* x (expt x (- y 1)))))))