diff options
| author | Rose Hogenson <rhogenson@posteo.net> | 2022-08-07 15:52:39 -0700 |
|---|---|---|
| committer | Rose Hogenson <rhogenson@posteo.net> | 2022-08-07 15:52:39 -0700 |
| commit | 6e5cf18539ee98bf2e5d05f7e6513e79a9aca1f7 (patch) | |
| tree | b8434c910ad28c8b9647bc871512fc9f386e7c80 /lib/scheme/base/80-number.csc | |
| parent | 24be96c3fc269d84df73c62c3fc26f87165b5032 (diff) | |
| download | chromatopelma-6e5cf18539ee98bf2e5d05f7e6513e79a9aca1f7.tar.zst | |
Add numbers.
Diffstat (limited to 'lib/scheme/base/80-number.csc')
| -rw-r--r-- | lib/scheme/base/80-number.csc | 821 |
1 files changed, 821 insertions, 0 deletions
diff --git a/lib/scheme/base/80-number.csc b/lib/scheme/base/80-number.csc new file mode 100644 index 0000000..20f8db8 --- /dev/null +++ b/lib/scheme/base/80-number.csc @@ -0,0 +1,821 @@ +(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 int<? 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 int<? 0 n)) + (and (boxed-int? n) + (boxed-int-positive? n)))) + + + (define (int-negative? n) + (or (and (small-int? n) + (call-builtin int<? 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 int<? n1 n2)) + ((and (int-negative? n1) (int-positive? n2)) #t) + ((and (int-positive? n1) (int-negative? n2)) #f) + ((and (int-negative? n1) (int-negative? n2)) (int<? (- n2) (- n1))) + ((and (small-int? n1) (boxed-int? n2)) #t) + ((and (boxed-int? n1) (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 int<? 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 int<? i 0) #f) ; n1 is equal to n2 + ((call-builtin int<? (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 digit-mask #x2000000000000000) ; 2^61. Digits are 61-bit unsigned numbers. + + + (define (split-carry-bit n) + (values + (call-builtin div n digit-mask) + (call-builtin mod n digit-mask))) + + + (define (denormalized-int n) + (define-values (sgn n*) (if (int-positive? n) + (values #t n) + (values #f (- n)))) + (define-values (n2 n1) (split-carry-bit n*)) + (if (zero? n2) + (make-boxed-int sgn (vector n1)) + (make-boxed-int sgn (vector n1 n2)))) + + + (define (big-int+ n1 n2) + (cond + ((and (boxed-int-positive? n1) (boxed-int-negative? n2)) + (big-int- n1 (- n2))) + ((and (boxed-int-negative? n1) (boxed-int-positive? n2)) + (big-int- n2 (- n1))) + ((and (boxed-int-negative? n1) (boxed-int-negative? n2)) + (- (big-int+ (- n1) (- n2)))) + (else + (let* ((n1-digits (boxed-int-digits n1)) + (n2-digits (boxed-int-digits n2)) + (n1-ndigits (vector-length n1-digits)) + (n2-ndigits (vector-length n2-digits)) + (out-digits (make-vector (max n1-ndigits + n2-ndigits)))) + (let loop ((i 0) + (carry 0)) + (cond + ((and (int<? i n1-ndigits) + (int<? i n2-ndigits)) + (let-values (((c d) (split-carry-bit + (call-builtin add + carry + (call-builtin add + (vector-ref n1-digits i) + (vector-ref n2-digits i)))))) + (vector-set! out-digits i d) + (loop (call-builtin add 1 i) c))) + ((int<? i n1-ndigits) + (let-values (((c d) (split-carry-bit + (call-builtin add + carry + (vector-ref n1-digits i))))) + (vector-set! out-digits i d) + (loop (call-builtin add 1 i) c))) + ((int<? i n2-ndigits) + (let-values (((c d) (split-carry-bit + (call-builtin add + carry + (vector-ref n2-digits i))))))) + ((int-positive? carry) + (set! out-digits (vector-append out-digits (vector carry)))))) + (make-boxed-int #t out-digits))))) + + + (define (int+ n1 n2) + (cond + ((and (small-int? n1) + (small-int? n2)) + (let ((x (call-builtin add n1 n2)) + (n1-positive (int-positive? n1))) + (if (and (boolean=? n1-positive (int-positive? n2)) + (not (boolean=? n1-positive (int-positive? x)))) ; overflow + (big-int+ (denormalized-int n1) (denormalized-int n2)) + x))) + (else + (when (small-int? n1) + (set! n1 (denormalized-int n1))) + (when (small-int? n2) + (set! n2 (denormalized-int n2))) + (big-int+ n1 n2)))) + + + (define (remove-leading-zeros 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) + (define digits (boxed-int-digits n)) + (if (and (int<? (vector-length digits) 3) + (int<? (vector-ref digits 1) (if (int-positive? n) + 2 + 3))) ; int min requires special handling + (if (int-negative? n) + (- (+ (vector-ref digits 0) + (* digit-mask (vector-ref digits 1)))) + (+ (vector-ref digits 0) + (* digit-mask (vector-ref digits 1)))) + n)) ; I have no doubt that there is a bug in this procedure. + + + (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 (boxed-int-digits n1)) + (n2-digits (boxed-int-digits n2)) + (n1-ndigits (vector-length n1-digits)) + (n2-ndigits (vector-length n2-digits)) + (out-digits (make-vector (max n1-ndigits + n2-ndigits)))) + (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 x digit-mask)) + (loop (+ 1 i) -1)) + (begin + (vector-set! out-digits i x) + (loop (+ 1 i) 0))))))) + (normalize (remove-leading-zeros (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- (denormalized-int n1) (denormalized-int n2))))) + (else + (when (small-int? n1) + (set! n1 (denormalized-int n1))) + (when (small-int? n2) + (set! n2 (denormalized-int n2))) + (big-int- n1 n2)))) + + + (define (split-in-half n m) + (define digits (boxed-int-digits n)) + (if (int<? m (vector-length digits)) + (values + (normalize (make-boxed-int #t (vector-copy m))) + (normalize (make-boxed-int #t (vector-copy 0 m)))) + (values + 0 + n))) + + + (define (lshift n m) + (when (small-int? n) + (set! n (denormalized-int n))) + (let* ((old-digits (boxed-int-digits m)) + (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) + (cond + ((and (boxed-int-positive? x) (boxed-int-positive? y)) + (define m (call-builtin div + (max (vector-length (boxed-int-digits x)) (vector-length (boxed-int-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)) + (+ (lshift z2 (call-builtin mul 2 m)) + (lshift z1 m) + z0)) + ((and (boxed-int-negative? x) (boxed-int-negative? y)) + (big-int* (- x) (- y))) + ((boxed-int-negative? x) + (- (big-int* (- x) y))) + (else ; y is negative + (- (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 + (big-int* (denormalized-int n1) (denormalized-int n2)) + x))) + (else + (when (small-int? n1) + (set! n1 (denormalized-int n1))) + (when (small-int? n2) + (set! n2 (denormalized-int n2))) + (big-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 digit-mask)) + (define guess (call-builtin div (- hi lo) 2)) + (define check (- d (* m guess))) + (cond + ((negative? check) + ; guess was too big + (loop lo guess)) + ((int<? check m) + ; got it! + guess) + (else + ; guess was too small + (loop (+ guess 1) hi))))) + + + (define (big-int/ n m) + (define n-digits (boxed-int-digits n)) + (define m-digits (boxed-int-digits m)) + (define k (vector-length n-digits)) + (define l (vector-length m-digits)) + (if (int<? k l) + (values 0 (normalize n)) + (let loop ((i (- l 1)) + (q 0) + (r (normalize + (make-boxed-int #t (vector-copy (boxed-int-digits n) (- k l 1)))))) ; last l-1 digits of n + (if (int<? i k) + (let* ((d (+ (* digit-mask r) + (left-index n-digits i))) + (beta (find-beta d m))) + (loop (+ 1 i) + (+ (* digit-mask q) + 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 + (when (small-int? n1) + (set! n1 (denormalized-int n1))) + (when (small-int? n2) + (set! n2 (denormalized-int n2))) + ; big and positive + (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)) + (g (gcd d1 d2))) + (int<? (* (numerator x1) + (quotient d2 g)) + (* (numerator x2) + (quotient d1 g)))))) + + + (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 m1)) + (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)) + (g (gcd d1 d2)) + (s1 (quotient d2 g)) + (s2 (quotient d1 g))) + (/ (+ (* 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 + ((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? d) + (set! q (make-quotient (- (numerator q)) + (- (denominator q))))) + (let* ((n (numerator q)) + (d (denominator q)) + (g (gcd n d))) + (make-quotient (quotient n g) + (quotient d g)))) + + + (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))) + (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 (/ k 4))) ; initial estimate + (define y (/ (+ x (/ k x)) + 2)) ; 1/2 (x + k/x) + (if (> 1 (abs (- x y))) + (let ((s (truncate y))) + (values s (- k (* s s)))) + (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))))))) |
