From 6e5cf18539ee98bf2e5d05f7e6513e79a9aca1f7 Mon Sep 17 00:00:00 2001
From: Rose Hogenson <rhogenson@posteo.net>
Date: Sun, 7 Aug 2022 15:52:39 -0700
Subject: Add numbers.

---
 lib/scheme/base/80-number.csc | 821 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 821 insertions(+)
 create mode 100644 lib/scheme/base/80-number.csc

(limited to 'lib/scheme/base')

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)))))))
-- 
cgit v1.3.1