aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorRose Hogenson <rhogenson@posteo.net>2022-08-07 15:52:39 -0700
committerRose Hogenson <rhogenson@posteo.net>2022-08-07 15:52:39 -0700
commit6e5cf18539ee98bf2e5d05f7e6513e79a9aca1f7 (patch)
treeb8434c910ad28c8b9647bc871512fc9f386e7c80
parent24be96c3fc269d84df73c62c3fc26f87165b5032 (diff)
downloadchromatopelma-6e5cf18539ee98bf2e5d05f7e6513e79a9aca1f7.tar.zst
Add numbers.
-rw-r--r--lib/scheme/base/80-number.csc821
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)))))))