Permalink
Fetching contributors…
Cannot retrieve contributors at this time
11511 lines (10133 sloc) 410 KB
;;;============================================================================
;;; File: "_num.scm"
;;; Copyright (c) 1994-2016 by Marc Feeley, All Rights Reserved.
;;; Copyright (c) 2004-2016 by Brad Lucier, All Rights Reserved.
;;;============================================================================
(macro-case-target
((C)
(c-declare "#include \"mem.h\"")
(##define-macro (use-fast-bignum-algorithms) #t))
(else
(##define-macro (use-fast-bignum-algorithms) #f)))
;;;============================================================================
;;; Implementation of exceptions.
(implement-library-type-range-exception)
(define-prim (##raise-range-exception arg-num proc . args)
(##extract-procedure-and-arguments
proc
args
arg-num
#f
#f
(lambda (procedure arguments arg-num dummy1 dummy2)
(macro-raise
(macro-make-range-exception procedure arguments arg-num)))))
(implement-library-type-divide-by-zero-exception)
(define-prim (##raise-divide-by-zero-exception proc . args)
(##extract-procedure-and-arguments
proc
args
#f
#f
#f
(lambda (procedure arguments dummy1 dummy2 dummy3)
(macro-raise
(macro-make-divide-by-zero-exception procedure arguments)))))
(implement-library-type-fixnum-overflow-exception)
(define-prim (##raise-fixnum-overflow-exception proc . args)
(##extract-procedure-and-arguments
proc
args
#f
#f
#f
(lambda (procedure arguments dummy1 dummy2 dummy3)
(macro-raise
(macro-make-fixnum-overflow-exception procedure arguments)))))
;;;----------------------------------------------------------------------------
;;; Define type checking procedures.
(define-fail-check-type exact-signed-int8 'exact-signed-int8)
(define-fail-check-type exact-signed-int8-list 'exact-signed-int8-list)
(define-fail-check-type exact-unsigned-int8 'exact-unsigned-int8)
(define-fail-check-type exact-unsigned-int8-list 'exact-unsigned-int8-list)
(define-fail-check-type exact-signed-int16 'exact-signed-int16)
(define-fail-check-type exact-signed-int16-list 'exact-signed-int16-list)
(define-fail-check-type exact-unsigned-int16 'exact-unsigned-int16)
(define-fail-check-type exact-unsigned-int16-list 'exact-unsigned-int16-list)
(define-fail-check-type exact-signed-int32 'exact-signed-int32)
(define-fail-check-type exact-signed-int32-list 'exact-signed-int32-list)
(define-fail-check-type exact-unsigned-int32 'exact-unsigned-int32)
(define-fail-check-type exact-unsigned-int32-list 'exact-unsigned-int32-list)
(define-fail-check-type exact-signed-int64 'exact-signed-int64)
(define-fail-check-type exact-signed-int64-list 'exact-signed-int64-list)
(define-fail-check-type exact-unsigned-int64 'exact-unsigned-int64)
(define-fail-check-type exact-unsigned-int64-list 'exact-unsigned-int64-list)
(define-fail-check-type inexact-real 'inexact-real)
(define-fail-check-type inexact-real-list 'inexact-real-list)
(define-fail-check-type number 'number)
(define-fail-check-type real 'real)
(define-fail-check-type finite-real 'finite-real)
(define-fail-check-type rational 'rational)
(define-fail-check-type integer 'integer)
(define-fail-check-type exact-integer 'exact-integer)
(define-fail-check-type fixnum 'fixnum)
(define-fail-check-type flonum 'flonum)
;;;----------------------------------------------------------------------------
;;; Numerical type predicates.
(define-prim (##number? x)
(##complex? x))
(define-prim (##complex? x)
(macro-number-dispatch x #f
#t ;; x = fixnum
#t ;; x = bignum
#t ;; x = ratnum
#t ;; x = flonum
#t)) ;; x = cpxnum
(define-prim (number? x)
(macro-force-vars (x)
(##number? x)))
(define-prim (complex? x)
(macro-force-vars (x)
(##complex? x)))
(define-prim (##real? x)
(macro-number-dispatch x #f
#t ;; x = fixnum
#t ;; x = bignum
#t ;; x = ratnum
#t ;; x = flonum
(macro-cpxnum-real? x))) ;; x = cpxnum
(define-prim (real? x)
(macro-force-vars (x)
(##real? x)))
(define-prim (##rational? x)
(macro-number-dispatch x #f
#t ;; x = fixnum
#t ;; x = bignum
#t ;; x = ratnum
(macro-flonum-rational? x) ;; x = flonum
(macro-cpxnum-rational? x))) ;; x = cpxnum
(define-prim (rational? x)
(macro-force-vars (x)
(##rational? x)))
(define-prim (##integer? x)
(macro-number-dispatch x #f
#t ;; x = fixnum
#t ;; x = bignum
#f ;; x = ratnum
(macro-flonum-int? x) ;; x = flonum
(macro-cpxnum-int? x))) ;; x = cpxnum
(define-prim (integer? x)
(macro-force-vars (x)
(##integer? x)))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; Exactness predicates.
(define-prim (##exact? x)
(define (type-error) #f)
(macro-number-dispatch x (type-error)
#t ;; x = fixnum
#t ;; x = bignum
#t ;; x = ratnum
#f ;; x = flonum
(and (##not (##flonum? (macro-cpxnum-real x))) ;; x = cpxnum
(##not (##flonum? (macro-cpxnum-imag x))))))
(define-prim (exact? x)
(macro-force-vars (x)
(let ()
(define (type-error)
(##fail-check-number 1 exact? x))
(macro-number-dispatch x (type-error)
#t ;; x = fixnum
#t ;; x = bignum
#t ;; x = ratnum
#f ;; x = flonum
(and (##not (##flonum? (macro-cpxnum-real x))) ;; x = cpxnum
(##not (##flonum? (macro-cpxnum-imag x))))))))
(define-prim (##inexact? x)
(define (type-error) #f)
(macro-number-dispatch x (type-error)
#f ;; x = fixnum
#f ;; x = bignum
#f ;; x = ratnum
#t ;; x = flonum
(or (##flonum? (macro-cpxnum-real x)) ;; x = cpxnum
(##flonum? (macro-cpxnum-imag x)))))
(define-prim (inexact? x)
(macro-force-vars (x)
(let ()
(define (type-error)
(##fail-check-number 1 inexact? x))
(macro-number-dispatch x (type-error)
#f ;; x = fixnum
#f ;; x = bignum
#f ;; x = ratnum
#t ;; x = flonum
(or (##flonum? (macro-cpxnum-real x)) ;; x = cpxnum
(##flonum? (macro-cpxnum-imag x)))))))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; Comparison predicates.
(define-prim (##= x y)
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(macro-number-dispatch x (type-error-on-x)
(macro-number-dispatch y (type-error-on-y) ;; x = fixnum
(##fx= x y)
#f
#f
(if (##fixnum->flonum-exact? x)
(##fl= (##fixnum->flonum x) y)
(and (##flfinite? y)
(##ratnum.= (##exact-int->ratnum x) (##flonum->ratnum y))))
(##cpxnum.= (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = bignum
#f
(or (##eq? x y)
(##exact-int.= x y))
#f
(and (##flfinite? y)
(##ratnum.= (##exact-int->ratnum x) (##flonum->ratnum y)))
(##cpxnum.= (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = ratnum
#f
#f
(or (##eq? x y)
(##ratnum.= x y))
(and (##flfinite? y)
(##ratnum.= x (##flonum->ratnum y)))
(##cpxnum.= (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = flonum
(if (##fixnum->flonum-exact? y)
(##fl= x (##fixnum->flonum y))
(and (##flfinite? x)
(##ratnum.= (##flonum->ratnum x) (##exact-int->ratnum y))))
(and (##flfinite? x)
(##ratnum.= (##flonum->ratnum x) (##exact-int->ratnum y)))
(and (##flfinite? x)
(##ratnum.= (##flonum->ratnum x) y))
(##fl= x y)
(##cpxnum.= (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = cpxnum
(##cpxnum.= x (##noncpxnum->cpxnum y))
(##cpxnum.= x (##noncpxnum->cpxnum y))
(##cpxnum.= x (##noncpxnum->cpxnum y))
(##cpxnum.= x (##noncpxnum->cpxnum y))
(##cpxnum.= x y))))
(define-prim-nary-bool (= x y)
#t
(if (##number? x) #t '(1))
(##= x y)
macro-force-vars
macro-no-check
(##pair? ##fail-check-number))
(define-prim (##< x y #!optional (nan-result #f))
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(macro-number-dispatch x (type-error-on-x)
(macro-number-dispatch y (type-error-on-y) ;; x = fixnum
(##fx< x y)
(##not (##bignum.negative? y))
(##ratnum.< (##exact-int->ratnum x) y)
(cond ((##flfinite? y)
(if (##fixnum->flonum-exact? x)
(##fl< (##fixnum->flonum x) y)
(##ratnum.< (##exact-int->ratnum x) (##flonum->ratnum y))))
((##flnan? y)
nan-result)
(else
(##flpositive? y)))
(if (macro-cpxnum-real? y)
(##< x (macro-cpxnum-real y) nan-result)
(type-error-on-y)))
(macro-number-dispatch y (type-error-on-y) ;; x = bignum
(##bignum.negative? x)
(##exact-int.< x y)
(##ratnum.< (##exact-int->ratnum x) y)
(cond ((##flfinite? y)
(##ratnum.< (##exact-int->ratnum x) (##flonum->ratnum y)))
((##flnan? y)
nan-result)
(else
(##flpositive? y)))
(if (macro-cpxnum-real? y)
(##< x (macro-cpxnum-real y) nan-result)
(type-error-on-y)))
(macro-number-dispatch y (type-error-on-y) ;; x = ratnum
(##ratnum.< x (##exact-int->ratnum y))
(##ratnum.< x (##exact-int->ratnum y))
(##ratnum.< x y)
(cond ((##flfinite? y)
(##ratnum.< x (##flonum->ratnum y)))
((##flnan? y)
nan-result)
(else
(##flpositive? y)))
(if (macro-cpxnum-real? y)
(##< x (macro-cpxnum-real y) nan-result)
(type-error-on-y)))
(macro-number-dispatch y (type-error-on-y) ;; x = flonum
(cond ((##flfinite? x)
(if (##fixnum->flonum-exact? y)
(##fl< x (##fixnum->flonum y))
(##ratnum.< (##flonum->ratnum x) (##exact-int->ratnum y))))
((##flnan? x)
nan-result)
(else
(##flnegative? x)))
(cond ((##flfinite? x)
(##ratnum.< (##flonum->ratnum x) (##exact-int->ratnum y)))
((##flnan? x)
nan-result)
(else
(##flnegative? x)))
(cond ((##flfinite? x)
(##ratnum.< (##flonum->ratnum x) y))
((##flnan? x)
nan-result)
(else
(##flnegative? x)))
(if (or (##flnan? x) (##flnan? y))
nan-result
(##fl< x y))
(if (macro-cpxnum-real? y)
(##< x (macro-cpxnum-real y) nan-result)
(type-error-on-y)))
(if (macro-cpxnum-real? x) ;; x = cpxnum
(macro-number-dispatch y (type-error-on-y)
(##< (macro-cpxnum-real x) y nan-result)
(##< (macro-cpxnum-real x) y nan-result)
(##< (macro-cpxnum-real x) y nan-result)
(##< (macro-cpxnum-real x) y nan-result)
(if (macro-cpxnum-real? y)
(##< (macro-cpxnum-real x) (macro-cpxnum-real y) nan-result)
(type-error-on-y)))
(type-error-on-x))))
(define-prim-nary-bool (< x y)
#t
(if (##real? x) #t '(1))
(##< x y #f)
macro-force-vars
macro-no-check
(##pair? ##fail-check-real))
(define-prim-nary-bool (> x y)
#t
(if (##real? x) #t '(1))
(##< y x #f)
macro-force-vars
macro-no-check
(##pair? ##fail-check-real))
(define-prim-nary-bool (<= x y)
#t
(if (##real? x) #t '(1))
(##not (##< y x #t))
macro-force-vars
macro-no-check
(##pair? ##fail-check-real))
(define-prim-nary-bool (>= x y)
#t
(if (##real? x) #t '(1))
(##not (##< x y #t))
macro-force-vars
macro-no-check
(##pair? ##fail-check-real))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; Numerical property predicates.
(define-prim (##zero? x)
(define (type-error)
(##fail-check-number 1 zero? x))
(macro-number-dispatch x (type-error)
(##fxzero? x)
#f
#f
(##flzero? x)
(and (let ((imag (macro-cpxnum-imag x)))
(and (##flonum? imag) (##flzero? imag)))
(let ((real (macro-cpxnum-real x)))
(if (##fixnum? real)
(##fxzero? real)
(and (##flonum? real) (##flzero? real)))))))
(define-prim (zero? x)
(macro-force-vars (x)
(##zero? x)))
(define-prim (##positive? x)
(define (type-error)
(##fail-check-real 1 positive? x))
(macro-number-dispatch x (type-error)
(##fxpositive? x)
(##not (##bignum.negative? x))
(##positive? (macro-ratnum-numerator x))
(##flpositive? x)
(if (macro-cpxnum-real? x)
(##positive? (macro-cpxnum-real x))
(type-error))))
(define-prim (positive? x)
(macro-force-vars (x)
(##positive? x)))
(define-prim (##negative? x)
(define (type-error)
(##fail-check-real 1 negative? x))
(macro-number-dispatch x (type-error)
(##fxnegative? x)
(##bignum.negative? x)
(##negative? (macro-ratnum-numerator x))
(##flnegative? x)
(if (macro-cpxnum-real? x)
(##negative? (macro-cpxnum-real x))
(type-error))))
(define-prim (negative? x)
(macro-force-vars (x)
(##negative? x)))
(define-prim (##odd? x)
(define (type-error)
(##fail-check-integer 1 odd? x))
(macro-number-dispatch x (type-error)
(##fxodd? x)
(macro-bignum-odd? x)
(type-error)
(if (macro-flonum-int? x)
(##odd? (##flonum->exact-int x))
(type-error))
(if (macro-cpxnum-int? x)
(##odd? (##inexact->exact (macro-cpxnum-real x)))
(type-error))))
(define-prim (odd? x)
(macro-force-vars (x)
(##odd? x)))
(define-prim (##even? x)
(define (type-error)
(##fail-check-integer 1 even? x))
(macro-number-dispatch x (type-error)
(##not (##fxodd? x))
(##not (macro-bignum-odd? x))
(type-error)
(if (macro-flonum-int? x)
(##even? (##flonum->exact-int x))
(type-error))
(if (macro-cpxnum-int? x)
(##even? (##inexact->exact (macro-cpxnum-real x)))
(type-error))))
(define-prim (even? x)
(macro-force-vars (x)
(##even? x)))
(define-prim (##finite? x)
(define (type-error)
(##fail-check-real 1 finite? x))
(macro-number-dispatch x (type-error)
#t
#t
#t
(##flfinite? x)
(if (macro-cpxnum-real? x)
(let ((real (macro-cpxnum-real x)))
(or (##not (##flonum? real))
(##flfinite? real)))
(type-error))))
(define-prim (finite? x)
(macro-force-vars (x)
(##finite? x)))
(define-prim (##infinite? x)
(define (type-error)
(##fail-check-real 1 infinite? x))
(macro-number-dispatch x (type-error)
#f
#f
#f
(##flinfinite? x)
(if (macro-cpxnum-real? x)
(let ((real (macro-cpxnum-real x)))
(and (##flonum? real)
(##flinfinite? real)))
(type-error))))
(define-prim (infinite? x)
(macro-force-vars (x)
(##infinite? x)))
(define-prim (##nan? x)
(define (type-error)
(##fail-check-real 1 nan? x))
(macro-number-dispatch x (type-error)
#f
#f
#f
(##flnan? x)
(if (macro-cpxnum-real? x)
(let ((real (macro-cpxnum-real x)))
(and (##flonum? real)
(##flnan? real)))
(type-error))))
(define-prim (nan? x)
(macro-force-vars (x)
(##nan? x)))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; Max and min.
(define-prim (##max x y)
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(macro-number-dispatch x (type-error-on-x)
(macro-number-dispatch y (type-error-on-y) ;; x = fixnum
(##fxmax x y)
(if (##< x y) y x)
(if (##< x y) y x)
(##flmax (##fixnum->flonum x) y)
(if (macro-cpxnum-real? y)
(##max x (macro-cpxnum-real y))
(type-error-on-y)))
(macro-number-dispatch y (type-error-on-y) ;; x = bignum
(if (##< x y) y x)
(if (##< x y) y x)
(if (##< x y) y x)
(##flmax (##exact-int->flonum x) y)
(if (macro-cpxnum-real? y)
(##max x (macro-cpxnum-real y))
(type-error-on-y)))
(macro-number-dispatch y (type-error-on-y) ;; x = ratnum
(if (##< x y) y x)
(if (##< x y) y x)
(if (##< x y) y x)
(##flmax (##ratnum->flonum x) y)
(if (macro-cpxnum-real? y)
(##max x (macro-cpxnum-real y))
(type-error-on-y)))
(macro-number-dispatch y (type-error-on-y) ;; x = flonum
(##flmax x (##fixnum->flonum y))
(##flmax x (##exact-int->flonum y))
(##flmax x (##ratnum->flonum y))
(##flmax x y)
(if (macro-cpxnum-real? y)
(##max x (macro-cpxnum-real y))
(type-error-on-y)))
(if (macro-cpxnum-real? x) ;; x = cpxnum
(macro-number-dispatch y (type-error-on-y)
(##max (macro-cpxnum-real x) y)
(##max (macro-cpxnum-real x) y)
(##max (macro-cpxnum-real x) y)
(##max (macro-cpxnum-real x) y)
(if (macro-cpxnum-real? y)
(##max (macro-cpxnum-real x) (macro-cpxnum-real y))
(type-error-on-y)))
(type-error-on-x))))
(define-prim-nary (max x y)
()
(if (##real? x) x '(1))
(##max x y)
macro-force-vars
macro-no-check
(##pair? ##fail-check-real))
(define-prim (##min x y)
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(macro-number-dispatch x (type-error-on-x)
(macro-number-dispatch y (type-error-on-y) ;; x = fixnum
(##fxmin x y)
(if (##< x y) x y)
(if (##< x y) x y)
(##flmin (##fixnum->flonum x) y)
(if (macro-cpxnum-real? y)
(##min x (macro-cpxnum-real y))
(type-error-on-y)))
(macro-number-dispatch y (type-error-on-y) ;; x = bignum
(if (##< x y) x y)
(if (##< x y) x y)
(if (##< x y) x y)
(##flmin (##exact-int->flonum x) y)
(if (macro-cpxnum-real? y)
(##min x (macro-cpxnum-real y))
(type-error-on-y)))
(macro-number-dispatch y (type-error-on-y) ;; x = ratnum
(if (##< x y) x y)
(if (##< x y) x y)
(if (##< x y) x y)
(##flmin (##ratnum->flonum x) y)
(if (macro-cpxnum-real? y)
(##min x (macro-cpxnum-real y))
(type-error-on-y)))
(macro-number-dispatch y (type-error-on-y) ;; x = flonum
(##flmin x (##fixnum->flonum y))
(##flmin x (##exact-int->flonum y))
(##flmin x (##ratnum->flonum y))
(##flmin x y)
(if (macro-cpxnum-real? y)
(##min x (macro-cpxnum-real y))
(type-error-on-y)))
(if (macro-cpxnum-real? x) ;; x = cpxnum
(macro-number-dispatch y (type-error-on-y)
(##min (macro-cpxnum-real x) y)
(##min (macro-cpxnum-real x) y)
(##min (macro-cpxnum-real x) y)
(##min (macro-cpxnum-real x) y)
(if (macro-cpxnum-real? y)
(##min (macro-cpxnum-real x) (macro-cpxnum-real y))
(type-error-on-y)))
(type-error-on-x))))
(define-prim-nary (min x y)
()
(if (##real? x) x '(1))
(##min x y)
macro-force-vars
macro-no-check
(##pair? ##fail-check-real))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; +, *, -, /
(define-prim (##+ x y)
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(macro-number-dispatch x (type-error-on-x)
(macro-number-dispatch y (type-error-on-y) ;; x = fixnum
(or (##fx+? x y)
(##bignum.+ (##fixnum->bignum x) (##fixnum->bignum y)))
(if (##fxzero? x)
y
(##bignum.+ (##fixnum->bignum x) y))
(if (##fxzero? x)
y
(##ratnum.+ (##exact-int->ratnum x) y))
(if (and (macro-special-case-exact-zero?) (##fxzero? x))
y
(##fl+ (##fixnum->flonum x) y))
(if (and (macro-special-case-exact-zero?) (##fxzero? x))
y
(##cpxnum.+ (##noncpxnum->cpxnum x) y)))
(macro-number-dispatch y (type-error-on-y) ;; x = bignum
(if (##fxzero? y)
x
(##bignum.+ x (##fixnum->bignum y)))
(##bignum.+ x y)
(##ratnum.+ (##exact-int->ratnum x) y)
(##fl+ (##exact-int->flonum x) y)
(##cpxnum.+ (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = ratnum
(if (##fxzero? y)
x
(##ratnum.+ x (##exact-int->ratnum y)))
(##ratnum.+ x (##exact-int->ratnum y))
(##ratnum.+ x y)
(##fl+ (##ratnum->flonum x) y)
(##cpxnum.+ (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = flonum
(if (and (macro-special-case-exact-zero?) (##fxzero? y))
x
(##fl+ x (##fixnum->flonum y)))
(##fl+ x (##exact-int->flonum y))
(##fl+ x (##ratnum->flonum y))
(##fl+ x y)
(##cpxnum.+ (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = cpxnum
(if (and (macro-special-case-exact-zero?) (##fxzero? y))
x
(##cpxnum.+ x (##noncpxnum->cpxnum y)))
(##cpxnum.+ x (##noncpxnum->cpxnum y))
(##cpxnum.+ x (##noncpxnum->cpxnum y))
(##cpxnum.+ x (##noncpxnum->cpxnum y))
(##cpxnum.+ x y))))
(define-prim-nary (+ x y)
0
(if (##number? x) x '(1))
(##+ x y)
macro-force-vars
macro-no-check
(##pair? ##fail-check-number))
(define-prim (##* x y)
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(macro-number-dispatch x (type-error-on-x)
(macro-number-dispatch y (type-error-on-y) ;; x = fixnum
(cond ((##fx= y 0)
0)
((if (##fx= y -1)
(##fx-? x)
(##fx*? x y))
=> (lambda (result) result))
(else
(##bignum.* (##fixnum->bignum x) (##fixnum->bignum y))))
(cond ((##fxzero? x)
0)
((##fx= x 1)
y)
((##fx= x -1)
(##negate y))
(else
(##bignum.* (##fixnum->bignum x) y)))
(cond ((##fxzero? x)
0)
((##fx= x 1)
y)
((##fx= x -1)
(##negate y))
(else
(##ratnum.* (##exact-int->ratnum x) y)))
(cond ((and (macro-special-case-exact-zero?)
(##fxzero? x))
0)
((##fx= x 1)
y)
(else
(##fl* (##fixnum->flonum x) y)))
(cond ((and (macro-special-case-exact-zero?)
(##fxzero? x))
0)
((##fx= x 1)
y)
(else
(##cpxnum.* (##noncpxnum->cpxnum x) y))))
(macro-number-dispatch y (type-error-on-y) ;; x = bignum
(cond ((##fx= y 0)
0)
((##fx= y 1)
x)
((##fx= y -1)
(##negate x))
(else
(##bignum.* x (##fixnum->bignum y))))
(##bignum.* x y)
(##ratnum.* (##exact-int->ratnum x) y)
(##fl* (##exact-int->flonum x) y)
(##cpxnum.* (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = ratnum
(cond ((##fxzero? y)
0)
((##fx= y 1)
x)
((##fx= y -1)
(##negate x))
(else
(##ratnum.* x (##exact-int->ratnum y))))
(##ratnum.* x (##exact-int->ratnum y))
(##ratnum.* x y)
(##fl* (##ratnum->flonum x) y)
(##cpxnum.* (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = flonum
(cond ((and (macro-special-case-exact-zero?) (##fxzero? y))
0)
((##fx= y 1)
x)
(else
(##fl* x (##fixnum->flonum y))))
(##fl* x (##exact-int->flonum y))
(##fl* x (##ratnum->flonum y))
(##fl* x y)
(##cpxnum.* (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = cpxnum
(cond ((and (macro-special-case-exact-zero?) (##fxzero? y))
0)
((##fx= y 1)
x)
(else
(##cpxnum.* x (##noncpxnum->cpxnum y))))
(##cpxnum.* x (##noncpxnum->cpxnum y))
(##cpxnum.* x (##noncpxnum->cpxnum y))
(##cpxnum.* x (##noncpxnum->cpxnum y))
(##cpxnum.* x y))))
(define-prim-nary (* x y)
1
(if (##number? x) x '(1))
(##* x y)
macro-force-vars
macro-no-check
(##pair? ##fail-check-number))
(define-prim (##square x)
(define (type-error)
(##fail-check-number 1 square x))
(macro-number-dispatch x (type-error)
(or (##fxsquare? x)
(let ((x (##fixnum->bignum x)))
(##bignum.* x x)))
(##bignum.* x x)
(##ratnum.* x x)
(##fl* x x)
(let ((real (macro-cpxnum-real x))
(imag (macro-cpxnum-imag x)))
(if (or (##eqv? real 0)
(##exact? x))
(##make-rectangular (##* (##- real imag) (##+ real imag))
(##* 2 (##* real imag)))
(##csquare (##exact->inexact x))))))
(define-prim (square x)
(macro-force-vars (x)
(##square x)))
(define-prim (##negate x)
(##define-macro (type-error) `'(1))
(macro-number-dispatch x (type-error)
(or (##fx-? x)
(##bignum.- (##fixnum->bignum 0) (##fixnum->bignum ##min-fixnum)))
(##bignum.- (##fixnum->bignum 0) x)
(macro-ratnum-make (##negate (macro-ratnum-numerator x))
(macro-ratnum-denominator x))
(##fl- x)
(##make-rectangular (##negate (macro-cpxnum-real x))
(##negate (macro-cpxnum-imag x)))))
(define-prim (##- x y)
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(macro-number-dispatch x (type-error-on-x)
(macro-number-dispatch y (type-error-on-y) ;; x = fixnum
(or (##fx-? x y)
(##bignum.- (##fixnum->bignum x) (##fixnum->bignum y)))
(##bignum.- (##fixnum->bignum x) y)
(if (##fxzero? x)
(##negate y)
(##ratnum.- (##exact-int->ratnum x) y))
(if (and (macro-special-case-exact-zero?) (##fxzero? x))
(##fl- y)
(##fl- (##fixnum->flonum x) y))
(##cpxnum.- (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = bignum
(if (##fxzero? y)
x
(##bignum.- x (##fixnum->bignum y)))
(##bignum.- x y)
(##ratnum.- (##exact-int->ratnum x) y)
(##fl- (##exact-int->flonum x) y)
(##cpxnum.- (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = ratnum
(if (##fxzero? y)
x
(##ratnum.- x (##exact-int->ratnum y)))
(##ratnum.- x (##exact-int->ratnum y))
(##ratnum.- x y)
(##fl- (##ratnum->flonum x) y)
(##cpxnum.- (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = flonum
(if (and (macro-special-case-exact-zero?) (##fxzero? y))
x
(##fl- x (##fixnum->flonum y)))
(##fl- x (##exact-int->flonum y))
(##fl- x (##ratnum->flonum y))
(##fl- x y)
(##cpxnum.- (##noncpxnum->cpxnum x) y))
(macro-number-dispatch y (type-error-on-y) ;; x = cpxnum
(if (and (macro-special-case-exact-zero?) (##fxzero? y))
x
(##cpxnum.- x (##noncpxnum->cpxnum y)))
(##cpxnum.- x (##noncpxnum->cpxnum y))
(##cpxnum.- x (##noncpxnum->cpxnum y))
(##cpxnum.- x (##noncpxnum->cpxnum y))
(##cpxnum.- x y))))
(define-prim-nary (- x y)
()
(##negate x)
(##- x y)
macro-force-vars
macro-no-check
(##pair? ##fail-check-number))
(define-prim (##inverse x)
(##define-macro (type-error) `'(1))
(define (divide-by-zero-error) #f)
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
(divide-by-zero-error)
(if (##fxnegative? x)
(if (##fx= x -1)
x
(macro-ratnum-make -1 (##negate x)))
(if (##fx= x 1)
x
(macro-ratnum-make 1 x))))
(if (##bignum.negative? x)
(macro-ratnum-make -1 (##negate x))
(macro-ratnum-make 1 x))
(let ((num (macro-ratnum-numerator x))
(den (macro-ratnum-denominator x)))
(cond ((##fx= num 1)
den)
((##fx= num -1)
(##negate den))
(else
(if (##negative? num)
(macro-ratnum-make (##negate den) (##negate num))
(macro-ratnum-make den num)))))
(##fl/ (macro-inexact-+1) x)
(##cpxnum./ (##noncpxnum->cpxnum 1) x)))
(define-prim (##/ x y)
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(define (divide-by-zero-error) #f)
(macro-number-dispatch y (type-error-on-y)
(macro-number-dispatch x (type-error-on-x) ;; y = fixnum
(cond ((##fxzero? y)
(divide-by-zero-error))
((##fx= y 1)
x)
((##fx= y -1)
(##negate x))
((##fxzero? x)
0)
((##fx= x 1)
(##inverse y))
(else
(##ratnum./ (##exact-int->ratnum x) (##exact-int->ratnum y))))
(cond ((##fxzero? y)
(divide-by-zero-error))
((##fx= y 1)
x)
((##fx= y -1)
(##negate x))
(else
(##ratnum./ (##exact-int->ratnum x) (##exact-int->ratnum y))))
(cond ((##fxzero? y)
(divide-by-zero-error))
((##fx= y 1)
x)
((##fx= y -1)
(##negate x))
(else
(##ratnum./ x (##exact-int->ratnum y))))
(cond ((##fxzero? y)
(divide-by-zero-error))
((##fx= y 1)
x)
(else
(##fl/ x (##fixnum->flonum y))))
(cond ((##fxzero? y)
(divide-by-zero-error))
((##fx= y 1)
x)
(else
(##cpxnum./ x (##noncpxnum->cpxnum y)))))
(macro-number-dispatch x (type-error-on-x) ;; y = bignum
(cond ((##fxzero? x)
0)
((##fx= x 1)
(##inverse y))
(else
(##ratnum./ (##exact-int->ratnum x) (##exact-int->ratnum y))))
(##ratnum./ (##exact-int->ratnum x) (##exact-int->ratnum y))
(##ratnum./ x (##exact-int->ratnum y))
(##fl/ x (##exact-int->flonum y))
(##cpxnum./ x (##noncpxnum->cpxnum y)))
(macro-number-dispatch x (type-error-on-x) ;; y = ratnum
(cond ((##fxzero? x)
0)
((##fx= x 1)
(##inverse y))
(else
(##ratnum./ (##exact-int->ratnum x) y)))
(##ratnum./ (##exact-int->ratnum x) y)
(##ratnum./ x y)
(##fl/ x (##ratnum->flonum y))
(##cpxnum./ x (##noncpxnum->cpxnum y)))
(macro-number-dispatch x (type-error-on-x) ;; y = flonum, no error possible
(##fl/ (##fixnum->flonum x) y)
(##fl/ (##exact-int->flonum x) y)
(##fl/ (##ratnum->flonum x) y)
(##fl/ x y)
(##cpxnum./ x (##noncpxnum->cpxnum y)))
(macro-number-dispatch x (type-error-on-x) ;; y = cpxnum
(##cpxnum./ (##noncpxnum->cpxnum x) y)
(##cpxnum./ (##noncpxnum->cpxnum x) y)
(##cpxnum./ (##noncpxnum->cpxnum x) y)
(##cpxnum./ (##noncpxnum->cpxnum x) y)
(##cpxnum./ x y))))
(define-prim-nary (/ x y)
()
(##inverse x)
(##/ x y)
macro-force-vars
macro-no-check
(##pair? ##fail-check-number)
(##not ##raise-divide-by-zero-exception))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; abs
(define-prim (##exact-int.negative? x)
(if (##fixnum? x)
(##fxnegative? x)
(##bignum.negative? x)))
(define-prim (##abs x)
(define (type-error)
(##fail-check-real 1 abs x))
(macro-number-dispatch x (type-error)
(if (##fxnegative? x) (##negate x) x)
(if (##bignum.negative? x) (##negate x) x)
(if (##exact-int.negative? (macro-ratnum-numerator x))
(macro-ratnum-make (##negate (macro-ratnum-numerator x))
(macro-ratnum-denominator x))
x)
(##flabs x)
(if (macro-cpxnum-real? x)
(##make-rectangular (##abs (macro-cpxnum-real x))
(##abs (macro-cpxnum-imag x)))
(type-error))))
(define-prim (abs x)
(macro-force-vars (x)
(##abs x)))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; quotient, remainder, modulo
(define-prim (##quotient x y)
(define (type-error-on-x)
(##fail-check-integer 1 quotient x y))
(define (type-error-on-y)
(##fail-check-integer 2 quotient x y))
(define (divide-by-zero-error)
(##raise-divide-by-zero-exception quotient x y))
(define (exact-quotient x y)
(##car (##exact-int.div x y)))
(define (inexact-quotient x y)
(let ((exact-y (##inexact->exact y)))
(if (##eqv? exact-y 0)
(divide-by-zero-error)
(##exact->inexact
(##quotient (##inexact->exact x) exact-y)))))
(macro-number-dispatch y (type-error-on-y)
(macro-number-dispatch x (type-error-on-x) ;; y = fixnum
(cond ((##fx= y 0)
(divide-by-zero-error))
((##fx= y -1) ;; (quotient ##min-fixnum -1) is a bignum
(##negate x))
(else
(##fxquotient x y)))
(cond ((##fx= y 0)
(divide-by-zero-error))
(else
(exact-quotient x y)))
(type-error-on-x)
(if (macro-flonum-int? x)
(inexact-quotient x y)
(type-error-on-x))
(if (macro-cpxnum-int? x)
(inexact-quotient x y)
(type-error-on-x)))
(macro-number-dispatch x (type-error-on-x) ;; y = bignum
(exact-quotient x y)
(exact-quotient x y)
(type-error-on-x)
(if (macro-flonum-int? x)
(inexact-quotient x y)
(type-error-on-x))
(if (macro-cpxnum-int? x)
(inexact-quotient x y)
(type-error-on-x)))
(type-error-on-y) ;; y = ratnum
(macro-number-dispatch x (type-error-on-x) ;; y = flonum
(if (macro-flonum-int? y)
(inexact-quotient x y)
(type-error-on-y))
(if (macro-flonum-int? y)
(inexact-quotient x y)
(type-error-on-y))
(type-error-on-x)
(if (macro-flonum-int? x)
(if (macro-flonum-int? y)
(inexact-quotient x y)
(type-error-on-y))
(type-error-on-x))
(if (macro-cpxnum-int? x)
(if (macro-flonum-int? y)
(inexact-quotient x y)
(type-error-on-y))
(type-error-on-x)))
(if (macro-cpxnum-int? y) ;; y = cpxnum
(macro-number-dispatch x (type-error-on-x)
(inexact-quotient x y)
(inexact-quotient x y)
(type-error-on-x)
(if (macro-flonum-int? x)
(inexact-quotient x y)
(type-error-on-x))
(if (macro-cpxnum-int? x)
(inexact-quotient x y)
(type-error-on-x)))
(type-error-on-y))))
(define-prim (quotient x y)
(macro-force-vars (x y)
(##quotient x y)))
(define-prim (##remainder x y)
(define (type-error-on-x)
(##fail-check-integer 1 remainder x y))
(define (type-error-on-y)
(##fail-check-integer 2 remainder x y))
(define (divide-by-zero-error)
(##raise-divide-by-zero-exception remainder x y))
(define (exact-remainder x y)
(##cdr (##exact-int.div x y
#f ;; need-quotient?
#t ;; keep-dividend?
)))
(define (inexact-remainder x y)
(let ((exact-y (##inexact->exact y)))
(if (##eqv? exact-y 0)
(divide-by-zero-error)
(##exact->inexact
(##remainder (##inexact->exact x) exact-y)))))
(macro-number-dispatch y (type-error-on-y)
(macro-number-dispatch x (type-error-on-x) ;; y = fixnum
(cond ((##fx= y 0)
(divide-by-zero-error))
(else
(##fxremainder x y)))
(cond ((##fx= y 0)
(divide-by-zero-error))
(else
(exact-remainder x y)))
(type-error-on-x)
(if (macro-flonum-int? x)
(inexact-remainder x y)
(type-error-on-x))
(if (macro-cpxnum-int? x)
(inexact-remainder x y)
(type-error-on-x)))
(macro-number-dispatch x (type-error-on-x) ;; y = bignum
(exact-remainder x y)
(exact-remainder x y)
(type-error-on-x)
(if (macro-flonum-int? x)
(inexact-remainder x y)
(type-error-on-x))
(if (macro-cpxnum-int? x)
(inexact-remainder x y)
(type-error-on-x)))
(type-error-on-y) ;; y = ratnum
(macro-number-dispatch x (type-error-on-x) ;; y = flonum
(if (macro-flonum-int? y)
(inexact-remainder x y)
(type-error-on-y))
(if (macro-flonum-int? y)
(inexact-remainder x y)
(type-error-on-y))
(type-error-on-x)
(if (macro-flonum-int? x)
(if (macro-flonum-int? y)
(inexact-remainder x y)
(type-error-on-y))
(type-error-on-x))
(if (macro-cpxnum-int? x)
(if (macro-flonum-int? y)
(inexact-remainder x y)
(type-error-on-y))
(type-error-on-x)))
(if (macro-cpxnum-int? y) ;; y = cpxnum
(macro-number-dispatch x (type-error-on-x)
(inexact-remainder x y)
(inexact-remainder x y)
(type-error-on-x)
(if (macro-flonum-int? x)
(inexact-remainder x y)
(type-error-on-x))
(if (macro-cpxnum-int? x)
(inexact-remainder x y)
(type-error-on-x)))
(type-error-on-y))))
(define-prim (remainder x y)
(macro-force-vars (x y)
(##remainder x y)))
(define-prim (##modulo x y)
(define (type-error-on-x)
(##fail-check-integer 1 modulo x y))
(define (type-error-on-y)
(##fail-check-integer 2 modulo x y))
(define (divide-by-zero-error)
(##raise-divide-by-zero-exception modulo x y))
(define (exact-modulo x y)
(let ((r (##cdr (##exact-int.div x
y
#f ;; need-quotient?
#t ;; keep-dividend?
))))
(if (##eqv? r 0)
0
(if (##eq? (##negative? x) (##negative? y))
r
(##+ r y)))))
(define (inexact-modulo x y)
(let ((exact-y (##inexact->exact y)))
(if (##eqv? exact-y 0)
(divide-by-zero-error)
(##exact->inexact
(##modulo (##inexact->exact x) exact-y)))))
(macro-number-dispatch y (type-error-on-y)
(macro-number-dispatch x (type-error-on-x) ;; y = fixnum
(cond ((##fx= y 0)
(divide-by-zero-error))
(else
(##fxmodulo x y)))
(cond ((##fx= y 0)
(divide-by-zero-error))
(else
(exact-modulo x y)))
(type-error-on-x)
(if (macro-flonum-int? x)
(inexact-modulo x y)
(type-error-on-x))
(if (macro-cpxnum-int? x)
(inexact-modulo x y)
(type-error-on-x)))
(macro-number-dispatch x (type-error-on-x) ;; y = bignum
(exact-modulo x y)
(exact-modulo x y)
(type-error-on-x)
(if (macro-flonum-int? x)
(inexact-modulo x y)
(type-error-on-x))
(if (macro-cpxnum-int? x)
(inexact-modulo x y)
(type-error-on-x)))
(type-error-on-y) ;; y = ratnum
(macro-number-dispatch x (type-error-on-x) ;; y = flonum
(if (macro-flonum-int? y)
(inexact-modulo x y)
(type-error-on-y))
(if (macro-flonum-int? y)
(inexact-modulo x y)
(type-error-on-y))
(type-error-on-x)
(if (macro-flonum-int? x)
(if (macro-flonum-int? y)
(inexact-modulo x y)
(type-error-on-y))
(type-error-on-x))
(if (macro-cpxnum-int? x)
(if (macro-flonum-int? y)
(inexact-modulo x y)
(type-error-on-y))
(type-error-on-x)))
(if (macro-cpxnum-int? y) ;; y = cpxnum
(macro-number-dispatch x (type-error-on-x)
(inexact-modulo x y)
(inexact-modulo x y)
(type-error-on-x)
(if (macro-flonum-int? x)
(inexact-modulo x y)
(type-error-on-x))
(if (macro-cpxnum-int? x)
(inexact-modulo x y)
(type-error-on-x)))
(type-error-on-y))))
(define-prim (modulo x y)
(macro-force-vars (x y)
(##modulo x y)))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; gcd, lcm
(define-prim (##gcd x y)
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(define (##fast-gcd u v)
;; See the paper "Fast Reduction and Composition of Binary
;; Quadratic Forms" by Arnold Schoenhage. His algorithm and proof
;; are derived from, and basically the same for, his Controlled
;; Euclidean Descent algorithm for gcd, which he has never
;; published. This algorithm has complexity log N times a
;; constant times the complexity of a multiplication of the same
;; size. We don't use it until we get to about 6800 bits. Note
;; that this is the same place that we start using FFT
;; multiplication and fast division with Newton's method for
;; finding inverses.
;; Niels Mo"ller has written two papers about an improved version
;; of this algorithm.
;; assumes u and v are nonnegative exact ints
(define (make-gcd-matrix A_11 A_12
A_21 A_22)
(##vector A_11 A_12
A_21 A_22))
(define (gcd-matrix_11 A)
(##vector-ref A 0))
(define (gcd-matrix_12 A)
(##vector-ref A 1))
(define (gcd-matrix_21 A)
(##vector-ref A 2))
(define (gcd-matrix_22 A)
(##vector-ref A 3))
(define (make-gcd-vector v_1 v_2)
(##vector v_1 v_2))
(define (gcd-vector_1 v)
(##vector-ref v 0))
(define (gcd-vector_2 v)
(##vector-ref v 1))
(define gcd-matrix-identity '#(1 0
0 1))
(define (gcd-matrix-multiply A B)
(cond ((##eq? A gcd-matrix-identity)
B)
((##eq? B gcd-matrix-identity)
A)
(else
(let ((A_11 (gcd-matrix_11 A)) (A_12 (gcd-matrix_12 A))
(A_21 (gcd-matrix_21 A)) (A_22 (gcd-matrix_22 A))
(B_11 (gcd-matrix_11 B)) (B_12 (gcd-matrix_12 B))
(B_21 (gcd-matrix_21 B)) (B_22 (gcd-matrix_22 B)))
(make-gcd-matrix (##+ (##* A_11 B_11)
(##* A_12 B_21))
(##+ (##* A_11 B_12)
(##* A_12 B_22))
(##+ (##* A_21 B_11)
(##* A_22 B_21))
(##+ (##* A_21 B_12)
(##* A_22 B_22)))))))
(define (gcd-matrix-multiply-strassen A B)
;; from http://mathworld.wolfram.com/StrassenFormulas.html
(cond ((##eq? A gcd-matrix-identity)
B)
((##eq? B gcd-matrix-identity)
A)
(else
(let ((A_11 (gcd-matrix_11 A)) (A_12 (gcd-matrix_12 A))
(A_21 (gcd-matrix_21 A)) (A_22 (gcd-matrix_22 A))
(B_11 (gcd-matrix_11 B)) (B_12 (gcd-matrix_12 B))
(B_21 (gcd-matrix_21 B)) (B_22 (gcd-matrix_22 B)))
(let ((Q_1 (##* (##+ A_11 A_22) (##+ B_11 B_22)))
(Q_2 (##* (##+ A_21 A_22) B_11))
(Q_3 (##* A_11 (##- B_12 B_22)))
(Q_4 (##* A_22 (##- B_21 B_11)))
(Q_5 (##* (##+ A_11 A_12) B_22))
(Q_6 (##* (##- A_21 A_11) (##+ B_11 B_12)))
(Q_7 (##* (##- A_12 A_22) (##+ B_21 B_22))))
(make-gcd-matrix (##+ (##+ Q_1 Q_4) (##- Q_7 Q_5))
(##+ Q_3 Q_5)
(##+ Q_2 Q_4)
(##+ (##+ Q_1 Q_3) (##- Q_6 Q_2))))))))
(define (gcd-matrix-solve A y)
(let ((y_1 (gcd-vector_1 y))
(y_2 (gcd-vector_2 y)))
(make-gcd-vector (##- (##* y_1 (gcd-matrix_22 A))
(##* y_2 (gcd-matrix_12 A)))
(##- (##* y_2 (gcd-matrix_11 A))
(##* y_1 (gcd-matrix_21 A))))))
(define (x>=2^n x n)
(##fx< n (##integer-length x)))
(define (determined-minimal? u v s)
;; assumes 2^s <= u , v; s>= 0 fixnum
;; returns #t if we can determine that |u-v|<2^s
;; at least one of u and v is a bignum
(let ((u (if (##fixnum? u) (##fixnum->bignum u) u))
(v (if (##fixnum? v) (##fixnum->bignum v) v)))
(let ((u-length (##bignum.mdigit-length u)))
(and (##fx= u-length (##bignum.mdigit-length v))
(let loop ((i (##fx- u-length 1)))
(let ((v-digit (##bignum.mdigit-ref v i))
(u-digit (##bignum.mdigit-ref u i)))
(if (and (##fxzero? u-digit)
(##fxzero? v-digit))
(loop (##fx- i 1))
(and (##fx= (##fxquotient s ##bignum.mdigit-width)
i)
(##fx< (##fxmax (##fx- u-digit v-digit)
(##fx- v-digit u-digit))
(##fxarithmetic-shift-left
1
(##fxremainder s ##bignum.mdigit-width)))))))))))
(define (gcd-small-step cont M u v s)
;; u, v >= 2^s
;; M is the matrix product of the partial sums of
;; the continued fraction representation of a/b so far
;; returns updated M, u, v, and a truth value
;; u, v >= 2^s and
;; if last return value is #t, we know that
;; (- (max u v) (min u v)) < 2^s, i.e, u, v are minimal above 2^s
(define (gcd-matrix-multiply-low M q)
(let ((M_11 (gcd-matrix_11 M))
(M_12 (gcd-matrix_12 M))
(M_21 (gcd-matrix_21 M))
(M_22 (gcd-matrix_22 M)))
(make-gcd-matrix (##+ M_11 (##* q M_12)) M_12
(##+ M_21 (##* q M_22)) M_22)))
(define (gcd-matrix-multiply-high M q)
(let ((M_11 (gcd-matrix_11 M))
(M_12 (gcd-matrix_12 M))
(M_21 (gcd-matrix_21 M))
(M_22 (gcd-matrix_22 M)))
(make-gcd-matrix M_11 (##+ (##* q M_11) M_12)
M_21 (##+ (##* q M_21) M_22))))
(if (or (##bignum? u)
(##bignum? v))
;; if u and v are nearly equal bignums, the two ##<
;; following this condition could take O(N) time to compute.
;; When this happens, however, it will be likely that
;; determined-minimal? will return true.
(cond ((determined-minimal? u v s)
(cont M
u
v
#t))
((##< u v)
(let* ((qr (##exact-int.div v u))
(q (##car qr))
(r (##cdr qr)))
(cond ((x>=2^n r s)
(cont (gcd-matrix-multiply-low M q)
u
r
#f))
((##eqv? q 1)
(cont M
u
v
#t))
(else
(cont (gcd-matrix-multiply-low M (##- q 1))
u
(##+ r u)
#t)))))
((##< v u)
(let* ((qr (##exact-int.div u v))
(q (##car qr))
(r (##cdr qr)))
(cond ((x>=2^n r s)
(cont (gcd-matrix-multiply-high M q)
r
v
#f))
((##eqv? q 1)
(cont M
u
v
#t))
(else
(cont (gcd-matrix-multiply-high M (##- q 1))
(##+ r v)
v
#t)))))
(else
(cont M
u
v
#t)))
;; here u and v are fixnums, so 2^s, which is <= u and v, is
;; also a fixnum
(let ((two^s (##fxarithmetic-shift-left 1 s)))
(if (##fx< u v)
(if (##fx< (##fx- v u) two^s)
(cont M
u
v
#t)
(let ((r (##fxremainder v u))
(q (##fxquotient v u)))
(if (##fx>= r two^s)
(cont (gcd-matrix-multiply-low M q)
u
r
#f)
;; the case when q is one and the remainder is < two^s
;; is covered in the first test
(cont (gcd-matrix-multiply-low M (##fx- q 1))
u
(##fx+ r u)
#t))))
;; here u >= v, but the case u = v is covered by the first test
(if (##fx< (##fx- u v) two^s)
(cont M
u
v
#t)
(let ((r (##fxremainder u v))
(q (##fxquotient u v)))
(if (##fx>= r two^s)
(cont (gcd-matrix-multiply-high M q)
r
v
#f)
;; the case when q is one and the remainder is < two^s
;; is covered in the first test
(cont (gcd-matrix-multiply-high M (##fx- q 1))
(##fx+ r v)
v
#t))))))))
(define (gcd-middle-step cont a b h m-prime cont-needs-M?)
((lambda (cont)
(if (and (x>=2^n a h)
(x>=2^n b h))
(MR cont a b h cont-needs-M?)
(cont gcd-matrix-identity a b)))
(lambda (M x y)
(let loop ((M M)
(x x)
(y y))
(if (or (x>=2^n x h)
(x>=2^n y h))
((lambda (cont) (gcd-small-step cont M x y m-prime))
(lambda (M x y minimal?)
(if minimal?
(cont M x y)
(loop M x y))))
((lambda (cont) (MR cont x y m-prime cont-needs-M?))
(lambda (M-prime alpha beta)
(cont (if cont-needs-M?
(if (##fx> (##fx- h m-prime) 1024)
;; here we trade off 1 multiplication
;; for 21 additions
(gcd-matrix-multiply-strassen M M-prime)
(gcd-matrix-multiply M M-prime))
gcd-matrix-identity)
alpha
beta))))))))
(define (MR cont a b m cont-needs-M?)
((lambda (cont)
(if (and (x>=2^n a (##fx+ m 2))
(x>=2^n b (##fx+ m 2)))
(let ((n (##fx- (##fxmax (##integer-length a)
(##integer-length b))
m)))
((lambda (cont)
(if (##fx<= m n)
(cont m 0)
(cont n (##fx- (##fx+ m 1) n))))
(lambda (m-prime p)
(let ((h (##fx+ m-prime (##fxquotient n 2))))
(if (##fx< 0 p)
(let ((a (##arithmetic-shift a (##fx- p)))
(b (##arithmetic-shift b (##fx- p)))
(a_0 (##extract-bit-field p 0 a))
(b_0 (##extract-bit-field p 0 b)))
((lambda (cont)
(gcd-middle-step cont a b h m-prime #t))
(lambda (M alpha beta)
(let ((M-inverse-v_0 (gcd-matrix-solve M (make-gcd-vector a_0 b_0))))
(cont (if cont-needs-M? M gcd-matrix-identity)
(##+ (##arithmetic-shift alpha p)
(gcd-vector_1 M-inverse-v_0))
(##+ (##arithmetic-shift beta p)
(gcd-vector_2 M-inverse-v_0)))))))
(gcd-middle-step cont a b h m-prime cont-needs-M?))))))
(cont gcd-matrix-identity
a
b)))
(lambda (M alpha beta)
(let loop ((M M)
(alpha alpha)
(beta beta)
(minimal? #f))
(if minimal?
(cont M alpha beta)
(gcd-small-step loop M alpha beta m))))))
((lambda (cont)
(if (and (use-fast-bignum-algorithms)
(##bignum? u)
(##bignum? v)
(x>=2^n u ##bignum.fast-gcd-size)
(x>=2^n v ##bignum.fast-gcd-size))
(MR cont u v ##bignum.fast-gcd-size #f)
(cont 0 u v)))
(lambda (M a b)
(general-base a b))))
(define (general-base a b)
(if (##eqv? b 0)
a
(let ((rem (cdr (##exact-int.div a b ;; calculate (remainder a b)
#f ;; need-quotient?
#f ;; keep-dividend?
))))
(if (##fixnum? b)
(fixnum-base b rem)
(general-base b rem)))))
(define (fixnum-base a b)
(##declare (not interrupts-enabled))
(if (##eqv? b 0)
a
(let ((a b)
(b (##fxremainder a b)))
(if (##eqv? b 0)
a
(fixnum-base b (##fxremainder a b))))))
(define (exact-gcd x y)
;; always returns an exact result, even with inexact arguments.
(let ((x (cond ((##inexact? x)
(##inexact->exact (##flabs x)))
((##negative? x)
(##negate x))
((##bignum? x)
(##bignum.copy x))
(else ;; nonnegative fixnum
x)))
(y (cond ((##inexact? y)
(##inexact->exact (##flabs y)))
((##negative? y)
(##negate y))
((##bignum? y)
(##bignum.copy y))
(else ;; nonnegative fixnum
y))))
;; now x and y are newly allocated, so we can overwrite them if
;; necessary in general-base
(cond ((##eqv? x 0)
y)
((##eqv? y 0)
x)
((and (##fixnum? x) (##fixnum? y))
(fixnum-base x y))
(else
(##fast-gcd x y)))))
(cond ((##not (##integer? x))
(type-error-on-x))
((##not (##integer? y))
(type-error-on-y))
((##eq? x y)
(##abs x))
(else
(if (and (##exact? x) (##exact? y))
(exact-gcd x y)
(##exact->inexact (exact-gcd x y))))))
(define-prim-nary (gcd x y)
0
(if (##integer? x) (##abs x) '(1))
(##gcd x y)
macro-force-vars
macro-no-check
(##pair? ##fail-check-integer))
(define-prim (##lcm x y)
(##define-macro (type-error-on-x) `'(1))
(##define-macro (type-error-on-y) `'(2))
(define (exact-lcm x y)
(if (or (##eqv? x 0) (##eqv? y 0))
0
(##abs (##* (##quotient x (##gcd x y))
y))))
(define (inexact-lcm x y)
(##exact->inexact
(exact-lcm (##inexact->exact x)
(##inexact->exact y))))
(cond ((##not (##integer? x))
(type-error-on-x))
((##not (##integer? y))
(type-error-on-y))
(else
(if (and (##exact? x) (##exact? y))
(exact-lcm x y)
(inexact-lcm x y)))))
(define-prim-nary (lcm x y)
1
(if (##integer? x) (##abs x) '(1))
(##lcm x y)
macro-force-vars
macro-no-check
(##pair? ##fail-check-integer))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; numerator, denominator
(define-prim (##numerator x)
(define (type-error)
(##fail-check-rational 1 numerator x))
(macro-number-dispatch x (type-error)
x
x
(macro-ratnum-numerator x)
(cond ((##flzero? x)
x)
((macro-flonum-rational? x)
(##exact->inexact (##numerator (##flonum->exact x))))
(else
(type-error)))
(if (macro-cpxnum-rational? x)
(##numerator (macro-cpxnum-real x))
(type-error))))
(define-prim (numerator x)
(macro-force-vars (x)
(##numerator x)))
(define-prim (##denominator x)
(define (type-error)
(##fail-check-rational 1 denominator x))
(macro-number-dispatch x (type-error)
1
1
(macro-ratnum-denominator x)
(if (macro-flonum-rational? x)
(##exact->inexact (##denominator (##flonum->exact x)))
(type-error))
(if (macro-cpxnum-rational? x)
(##denominator (macro-cpxnum-real x))
(type-error))))
(define-prim (denominator x)
(macro-force-vars (x)
(##denominator x)))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; floor, ceiling, truncate, round
(define-prim (##floor x)
(define (type-error)
(##fail-check-finite-real 1 floor x))
(macro-number-dispatch x (type-error)
x
x
(let ((num (macro-ratnum-numerator x))
(den (macro-ratnum-denominator x)))
(if (##negative? num)
(##quotient (##- num (##- den 1)) den)
(##quotient num den)))
(if (##flfinite? x)
(##flfloor x)
(type-error))
(if (macro-cpxnum-real? x)
(##floor (macro-cpxnum-real x))
(type-error))))
(define-prim (floor x)
(macro-force-vars (x)
(##floor x)))
(define-prim (##ceiling x)
(define (type-error)
(##fail-check-finite-real 1 ceiling x))
(macro-number-dispatch x (type-error)
x
x
(let ((num (macro-ratnum-numerator x))
(den (macro-ratnum-denominator x)))
(if (##negative? num)
(##quotient num den)
(##quotient (##+ num (##- den 1)) den)))
(if (##flfinite? x)
(##flceiling x)
(type-error))
(if (macro-cpxnum-real? x)
(##ceiling (macro-cpxnum-real x))
(type-error))))
(define-prim (ceiling x)
(macro-force-vars (x)
(##ceiling x)))
(define-prim (##truncate x)
(define (type-error)
(##fail-check-finite-real 1 truncate x))
(macro-number-dispatch x (type-error)
x
x
(##quotient (macro-ratnum-numerator x)
(macro-ratnum-denominator x))
(if (##flfinite? x)
(##fltruncate x)
(type-error))
(if (macro-cpxnum-real? x)
(##truncate (macro-cpxnum-real x))
(type-error))))
(define-prim (truncate x)
(macro-force-vars (x)
(##truncate x)))
(define-prim (##round x)
(define (type-error)
(##fail-check-finite-real 1 round x))
(macro-number-dispatch x (type-error)
x
x
(##ratnum.round x)
(if (##flfinite? x)
(##flround x)
(type-error))
(if (macro-cpxnum-real? x)
(##round (macro-cpxnum-real x))
(type-error))))
(define-prim (round x)
(macro-force-vars (x)
(##round x)))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; rationalize
(define-prim (##rationalize x y)
(define (simplest-rational1 x y)
(if (##< y x)
(simplest-rational2 y x)
(simplest-rational2 x y)))
(define (simplest-rational2 x y)
(cond ((##not (##< x y))
x)
((##positive? x)
(simplest-rational3 x y))
((##negative? y)
(##negate (simplest-rational3 (##negate y) (##negate x))))
(else
0)))
(define (simplest-rational3 x y)
(let ((fx (##floor x))
(fy (##floor y)))
(cond ((##not (##< fx x))
fx)
((##= fx fy)
(##+ fx
(##inverse
(simplest-rational3
(##inverse (##- y fy))
(##inverse (##- x fx))))))
(else
(##+ fx 1)))))
(cond ((##not (##rational? x))
(##fail-check-finite-real 1 rationalize x y))
((and (##flonum? y)
(##fl= y (macro-inexact-+inf)))
(macro-inexact-+0))
((##not (##rational? y))
(##fail-check-real 2 rationalize x y))
((##negative? y)
(##raise-range-exception 2 rationalize x y))
((and (##exact? x) (##exact? y))
(simplest-rational1 (##- x y) (##+ x y)))
(else
(let ((exact-x (##inexact->exact x))
(exact-y (##inexact->exact y)))
(##exact->inexact
(simplest-rational1 (##- exact-x exact-y)
(##+ exact-x exact-y)))))))
(define-prim (rationalize x y)
(macro-force-vars (x y)
(##rationalize x y)))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; trigonometry and complex numbers
#|
The next functions are from
Functions from
Branch Cuts for Complex Elementary Functions
or
Much Ado About Nothing's Sign Bit
by W. Kahan
Full reference:
Kahan, W: Branch cuts for complex elementary functions; or, Much ado about nothing’s sign bit. In Iserles, A., and Powell, M. (eds.), The state of the art in numerical analysis. Clarendon Press (1987) pp 165-211.
Note that Kahan's paper contains two treatments of branch cuts---Section 4, which deals with arithmetic with signed zeros (like IEEE arithmetic) and Section 5, which deals with arithmetic with only unsigned zeros. The codes in the paper are only for IEEE-style arithmetic.
Gambit Scheme is in a funny position, as it allows mixed-exactness complex numbers. We'll consider inexact real zeros (+0., -0.) as signed (of course), but we'll interpret exact zero (0) as unsigned.
The branch cuts of all the functions considered here lie on the exact real axis or the exact imaginary axis.
All of the inverse functions are defined in terms of log and sqrt, and the side of the continuity at the branch cuts is determined by the sides of continuity of those two functions.
I believe that this is the same as the continuity rules that the CLHS gives for atan, asin, acos, etc., along branch cuts.
Thanks to Raymond Toy for email discussions and for the code for cmucl, which gets this stuff right in the Common Lisp context.
See
http://140.177.205.23/InverseHyperbolicFunctions.html
for a discussion of branch cuts.
|#
(define-prim (##cabs z)
;; As far as I can tell, this is just magic. It works, and I'm not
;; going to touch it.
#|
Code to compute the constants using my computable reals package.
(load "exact-reals")
(define r2-exact
(computable-sqrt (exact->computable 2)))
(define r2p1-exact
(computable-+ r2-exact (exact->computable 1)))
(define r2p1
(computable->inexact r2p1-exact))
(define t2p1-exact
(computable-- r2p1-exact (exact->computable (inexact->exact r2p1))))
(define r2
(computable->inexact r2-exact))
(define t2p1
(computable->inexact t2p1-exact))
(for-each pretty-print
`((define r2 ,r2) (define r2p1 ,r2p1) (define t2p1 ,t2p1)))
|#
(define r2 1.4142135623730951)
(define r2p1 2.414213562373095)
(define t2p1 1.2537167179050217e-16)
(let ((x (##flabs (macro-cpxnum-real z)))
(y (##flabs (macro-cpxnum-imag z))))
(define (continue x y)
(let* ((x (if (##flinfinite? y) y x))
(t (##fl- x y)))
(if (and (##not (##fl= x +inf.0))
(##not (##fl= t x)))
(if (##fl> t y)
(let* ((s (##fl/ x y))
(s (##fl+ s (##flsqrt (##fl+ 1.0 (##fl* s s))))))
(##fl+ x (##fl/ y s)))
(let* ((s (##fl/ t y))
(t (##fl* (##fl+ 2.0 s) s))
(s (##fl+ r2p1
(##fl+ s
(##fl+ t2p1
(##fl/ t
(##fl+ r2 (##flsqrt (##fl+ 2.0 t)))))))))
(##fl+ x (##fl/ y s))))
x)))
(if (##fl< x y)
(continue y x)
(continue x y))))
(define-prim (##carg z)
(##angle z))
(define-prim (##csquare xi+ieta)
(let ((xi (macro-cpxnum-real xi+ieta))
(eta (macro-cpxnum-imag xi+ieta)))
(let ((x (##fl* (##fl- xi eta) (##fl+ xi eta)))
(y (##fl* 2.0 xi eta)))
(cond ((##flnan? x)
(cond ((##flinfinite? y)
(macro-cpxnum-make (##flcopysign (macro-inexact-+0) xi) y))
((##flinfinite? eta)
(macro-cpxnum-make (macro-inexact--inf) y))
((##flinfinite? xi)
(macro-cpxnum-make (macro-inexact-+inf) y))
(else
(macro-cpxnum-make x y))))
((and (##flnan? y)
(##flinfinite? x))
(macro-cpxnum-make x (##flcopysign (macro-inexact-+0) y)))
(else
(macro-cpxnum-make x y))))))
(define-prim (##cssqs x+iy)
(let ((x (macro-cpxnum-real x+iy))
(y (macro-cpxnum-imag x+iy)))
(cond ((or (##flinfinite? x)
(##flinfinite? y))
(##cons (macro-inexact-+inf) 0))
((and (##flzero? x)
(##flzero? y))
(##cons 0. 0))
(else
;; from now on, neither x nor y are infinite, and one is non-zero
(let* ((x^2 (##flsquare x))
(y^2 (##flsquare y))
(rho (##fl+ x^2 y^2)))
(if (or (##flinfinite? rho) ;; if rho is NaN, this is false
(and (or (##fl< x^2 (macro-inexact-lambda)) ;; poor man's way to see whether underflow flag was set
(##fl< y^2 (macro-inexact-lambda)))
(##fl< rho (##fl/ (macro-inexact-lambda) (macro-inexact-epsilon))))) ;; if rho is NaN, this is false
;; rho is not NaN, so x and y are not NaN, and x and y are not infinite. Whew.
(let ((k (##flilogb (##flmax (##flabs x) (##flabs y)))))
(##cons (##fl+ (##flsquare (##flscalbn x (##fx- k)))
(##flsquare (##flscalbn y (##fx- k))))
k))
(##cons rho 0)))))))
(define-prim (##csqrt x+iy)
(let* ((x (macro-cpxnum-real x+iy))
(y (macro-cpxnum-imag x+iy))
(rho+ik (##cssqs x+iy))
(rho (##car rho+ik))
(k (##cdr rho+ik))
(rho (if (##flnan? x)
rho
(##fl+ (##flscalbn (##flabs x) (##fx- k))
(##flsqrt rho))))
(rho (if (##fxodd? k)
(##flscalbn (##flsqrt rho) (##fxquotient (##fx- k 1) 2))
(##flscalbn (##flsqrt (##fl* 2.0 rho)) (##fx- (##fxquotient k 2) 1))))
(xi rho)
(eta y))
(if (##not (##fl= rho 0.0))
(let ((eta (if (##not (##flinfinite? (##flabs eta)))
(##fl/ (##fl/ eta rho) 2.0)
eta)))
(if (##flnegative? x)
(macro-cpxnum-make (##flabs eta) (##flcopysign rho y))
(macro-cpxnum-make xi eta)))
(macro-cpxnum-make xi eta))))
(define-prim (##cacos z)
(##- (macro-inexact-+pi/2) (##casin z)))
(define-prim (##cacosh z)
(let ((sqrt-z-1 (##sqrt (##- z 1)))
(sqrt-z+1 (##sqrt (##+ z 1))))
;; if z is real and > 1, then the imaginary part of the next expression can be
;; inexact 0, but that's OK because this routine is not called in this case.
(##make-rectangular (##asinh (##real-part (##* (##conjugate sqrt-z-1) sqrt-z+1)))
(##* 2 (##atan2 (##imag-part sqrt-z-1) (##real-part sqrt-z+1))))))
(define-prim (##casin z)
;; if (##real-part z) is exact zero, then there is a correlation of errors in sqrt-1-z and sqrt-1+z that
;; allows the next substitution
(let ((x (##real-part z)))
(if (##eqv? x 0)
(##make-rectangular 0 (##asinh (##imag-part z)))
(let ((sqrt-1-z (##sqrt (##- 1 z)))
(sqrt-1+z (##sqrt (##+ 1 z))))
(##make-rectangular (##atan2 x (##real-part (##* sqrt-1-z sqrt-1+z)))
(##asinh (##imag-part (##* (##conjugate sqrt-1-z) sqrt-1+z))))))))
(define-prim (##casinh z)
(##* -i (##casin (##* +i z))))
(define-prim (##catanh x+iy)
(define (x/x^2+y^2 x y)
(if (##fl< (##flabs y) (##flabs x))
(##fl/ 1. (##fl+ x (##fl* (##fl/ y x) y)))
(let ((x/y (##fl/ x y)))
(##fl/ x/y (##fl+ (##fl* x x/y) y)))))
(define (##->exact-sign x)
;; returns an exact number with the same sign as x, returns 1 if x is exact zero
(if (##flonum? x)
(##inexact->exact (##flcopysign 1. x))
(if (##negative? x) -1 1)))
(let* ((pi/2 (##* 2 (##atan 1)))
(theta (##fl/ (##flsqrt (macro-inexact-omega)) 4.))
(rho (##fl/ theta))
(beta (##->exact-sign (##real-part x+iy))) ;; beta is exact
(x+iy (##* beta (##conjugate x+iy)))
(x (##real-part x+iy))
(y (##imag-part x+iy))
(inexact-x (##exact->inexact x))
(inexact-y (##exact->inexact y))
(abs-y (##flabs inexact-y))
(zeta (cond ((or (##fl< theta inexact-x)
(##fl< theta abs-y))
(macro-cpxnum-make (##exact->inexact (x/x^2+y^2 inexact-x inexact-y))
(##flcopysign pi/2 inexact-y)))
((##fl= inexact-x 1.)
(macro-cpxnum-make (##fllog (##fl/ (##flsqrt (##flsqrt (##fl+ 4. (##flsquare abs-y))))
(##flsqrt (##fl+ abs-y rho))))
(##fl/ (##flcopysign (##fl+ pi/2 (##flatan (##fl/ (##fl+ abs-y rho) 2.0)))
inexact-y)
2.)))
(else
(macro-cpxnum-make (if (##eqv? x 0)
;; if rho and abs-y were exact in the next expression (no matter their values)
;; then the argument to fllog1p would be exact 0, so the result would be exact 0.
0
(##fl/ (##fllog1p (##fl/ (##fl* 4. inexact-x) ;; was (##* 4 x) originally
(##fl+ (##flsquare (##fl- 1. inexact-x))
(##flsquare (##fl+ abs-y rho)))))
4.))
(##fl/ (##carg (macro-cpxnum-make (##fl- (##fl* (##fl- 1. inexact-x)
(##fl+ 1. inexact-x))
(##flsquare (##fl+ abs-y rho)))
(##fl* 2. inexact-y)))
2.0))))))
(##* beta (##conjugate zeta))))
(define-prim (##ctanh xi+ieta)
;; we assume that neither xi nor eta can be exact 0
(let* ((xi (macro-cpxnum-real xi+ieta))
(eta (macro-cpxnum-imag xi+ieta)))
(if (##< (##fl/ (##flasinh (macro-inexact-omega)) 4.)
(##abs xi))
(macro-cpxnum-make (##flcopysign 1. (##exact->inexact xi)) ;; xi cannot be exact 0
(##flcopysign 0. (##exact->inexact eta))) ;; eta cannot be exact 0
(let* ((t (##tan eta)) ;; sin(eta)/cos(eta), can't be exact 0, so can't be exact
(beta (##fl+ 1. (##flsquare t))) ;; 1/cos^2(eta), can't be exact
(s (##sinh xi)) ;; sinh(xi), can't be exact zero, so can't be exact
(rho (##flsqrt (##fl+ 1. (##flsquare s))))) ;; cosh(xi), can't be exact
(if (##infinite? t) ;; if sin(eta)/cos(eta) = infinity (how, I don't know)
(macro-cpxnum-make (##fl/ rho s)
(##fl/ t))
(let ((one+beta*s^2 (##fl+ 1. (##fl* beta (##flsquare s)))))
(macro-cpxnum-make (##fl/ (##fl* beta (##fl* rho s))
one+beta*s^2)
(##fl/ t
one+beta*s^2))))))))
(define-prim (##ctan zeta)
(##* -i (##ctanh (##* +i zeta))))
;;; End of Kahan's functions
(define-prim (##conjugate x)
(define (type-error)
(##fail-check-number 1 conjugate x))
(macro-number-dispatch x (type-error)
x x x x (macro-cpxnum-make (macro-cpxnum-real x)
(##negate (macro-cpxnum-imag x)))))
(define-prim (conjugate x)
(macro-force-vars (x)
(##conjugate x)))
(define-prim (##exp x)
(define (type-error)
(##fail-check-number 1 exp x))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
1
(##flexp (##fixnum->flonum x)))
(##flexp (##exact-int->flonum x))
(##flexp (##ratnum->flonum x))
(##flexp x)
(##make-polar (##exp (macro-cpxnum-real x))
(macro-cpxnum-imag x))))
(define-prim (exp x)
(macro-force-vars (x)
(##exp x)))
(define-prim (##flonum-full-precision? x)
(let ((y (##flabs x)))
(and (##fl< y (macro-inexact-+inf))
(##fl<= (macro-flonum-min-normal) y))))
(define-prim (##log x)
(define (type-error)
(##fail-check-number 1 log x))
(define (range-error)
(##raise-range-exception 1 log x))
(define (negative-log x)
(##make-rectangular (##log (##negate x)) (macro-inexact-+pi)))
(define (exact-log x)
;; x is positive, x is not 1.
;; There are three places where just converting to a flonum and
;; taking the flonum logarithm doesn't work well.
;; 1. Overflow in the conversion
;; 2. Underflow in the conversion (or even loss of precision
;; because of a denormalized conversion result)
;; 3. When the number is close to 1.
(let ((float-x (##exact->inexact x)))
(cond ((##= x float-x)
(##fllog float-x)) ;; first, we trust the builtin flonum log
((##not (##flonum-full-precision? float-x))
;; direct conversion to flonum could incur massive relative
;; rounding errors, or would just lead to an infinite result
;; so we tolerate more than one rounding error in the calculation
(let* ((wn (##integer-length (##numerator x)))
(wd (##integer-length (##denominator x)))
(p (##fx- wn wd))
(float-p (##fixnum->flonum p))
(partial-result (##fllog
(##exact->inexact
(##* x (##expt 2 (##fx- p)))))))
(##fl+ (##fl* float-p
(macro-inexact-log-2))
partial-result)))
((or (##fl< (macro-inexact-exp-+1/2) float-x)
(##fl< float-x (macro-inexact-exp--1/2)))
;; here the absolute value of the logarithm is at least 0.5,
;; so there is less rounding error in the final result.
(##fllog float-x))
(else
;; use ln1p for arguments near one.
(##fllog1p (##exact->inexact (##- x 1)))))))
(define (complex-log-magnitude x)
(define (log-mag a b)
;; both are finite, 0 <= a <= b, b is nonzero
(let* ((c (##/ a b))
(approx-mag (##* b (##sqrt (##+ 1 (##* c c))))))
(if (or (##exact? approx-mag)
(and (##flonum-full-precision? approx-mag)
(or (##fl< (macro-inexact-exp-+1/2) approx-mag)
(##fl< approx-mag (macro-inexact-exp--1/2)))))
;; log composed with magnitude will compute a relatively accurate answer
(##log approx-mag)
(let ((a (##inexact->exact a))
(b (##inexact->exact b)))
(##* 1/2 (exact-log (##+ (##* a a) (##* b b))))))))
(let ((abs-r (##abs (##real-part x)))
(abs-i (##abs (##imag-part x))))
;; abs-i is not exact 0
(cond ((or (and (##flonum? abs-r)
(##fl= abs-r (macro-inexact-+inf)))
(and (##flonum? abs-i)
(##fl= abs-i (macro-inexact-+inf))))
(macro-inexact-+inf))
;; neither abs-r or abs-i is infinite
((and (##flonum? abs-r)
(##flnan? abs-r))
abs-r)
;; abs-r is not a NaN
((and (##flonum? abs-i)
(##flnan? abs-i))
abs-i)
;; abs-i is not a NaN
((##eqv? abs-r 0)
(##log abs-i))
;; abs-r is not exact 0
((and (##zero? abs-r)
(##zero? abs-i))
(macro-inexact--inf))
;; abs-i and abs-r are not both zero
(else
(if (##< abs-r abs-i)
(log-mag abs-r abs-i)
(log-mag abs-i abs-r))))))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
(range-error)
(if (##fxnegative? x)
(negative-log x)
(if (##eqv? x 1)
0
(exact-log x))))
(if (##bignum.negative? x)
(negative-log x)
(exact-log x))
(if (##negative? (macro-ratnum-numerator x))
(negative-log x)
(exact-log x))
(if (or (##flnan? x)
(##not (##flnegative?
(##flcopysign (macro-inexact-+1) x))))
(##fllog x)
(negative-log x))
(##make-rectangular (complex-log-magnitude x) (##angle x))))
(define-prim (log x)
(macro-force-vars (x)
(##log x)))
(define-prim (##sin x)
(define (type-error)
(##fail-check-number 1 sin x))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
0
(##flsin (##fixnum->flonum x)))
(##flsin (##exact-int->flonum x))
(##flsin (##ratnum->flonum x))
(##flsin x)
(let ((real (macro-cpxnum-real x))
(imag (macro-cpxnum-imag x)))
(if (and (##flonum? real)
(##flonum? imag))
;; fast path for flonums case
(macro-cpxnum-make (##fl* (##flsin real) (##flcosh imag))
(##fl* (##flcos real) (##flsinh imag)))
(##make-rectangular (##* (##sin real) (##cosh imag))
(##* (##cos real) (##sinh imag)))))))
(define-prim (sin x)
(macro-force-vars (x)
(##sin x)))
(define-prim (##cos x)
(define (type-error)
(##fail-check-number 1 cos x))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
1
(##flcos (##fixnum->flonum x)))
(##flcos (##exact-int->flonum x))
(##flcos (##ratnum->flonum x))
(##flcos x)
(let ((real (macro-cpxnum-real x))
(imag (macro-cpxnum-imag x)))
(if (and (##flonum? real)
(##flonum? imag))
;; fast path for flonums case
(macro-cpxnum-make (##fl* (##flcos real) (##flcosh imag))
(##fl- (##fl* (##flsin real) (##flsinh imag))))
(##make-rectangular (##* (##cos real) (##cosh imag))
(##negate (##* (##sin real) (##sinh imag))))))))
(define-prim (cos x)
(macro-force-vars (x)
(##cos x)))
(define-prim (##tan x)
(define (type-error)
(##fail-check-number 1 tan x))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
0
(##fltan (##fixnum->flonum x)))
(##fltan (##exact-int->flonum x))
(##fltan (##ratnum->flonum x))
(##fltan x)
;; complex ##tanh is the basic one here.
(##* -i (##tanh (##* +i x)))))
(define-prim (tan x)
(macro-force-vars (x)
(##tan x)))
(define-prim (##asin x)
(define (type-error)
(##fail-check-number 1 asin x))
(define (real-case x)
(if (or (##< 1 x)
(##< x -1))
(##casin (macro-cpxnum-make x 0))
(##flasin (##exact->inexact x))))
(macro-number-dispatch x (type-error)
(if (##eqv? x 0)
0
(real-case x))
(real-case x)
(real-case x)
(real-case x)
(##casin x)))
(define-prim (asin x)
(macro-force-vars (x)
(##asin x)))
(define-prim (##acos x)
(define (type-error)
(##fail-check-number 1 acos x))
(define (real-case x)
(if (or (##< 1 x)
(##< x -1))
(##cacos (macro-cpxnum-make x 0))
(##flacos (##exact->inexact x))))
(macro-number-dispatch x (type-error)
(if (##eqv? x 1)
0
(real-case x))
(real-case x)
(real-case x)
(real-case x)
(##cacos x)))
(define-prim (acos x)
(macro-force-vars (x)
(##acos x)))
(define-prim (##atan x)
(define (type-error)
(##fail-check-number 1 atan x))
(define (range-error)
(##raise-range-exception 1 atan x))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
0
(##flatan (##fixnum->flonum x)))
(##flatan (##exact-int->flonum x))
(##flatan (##ratnum->flonum x))
(##flatan x)
(let ((real (macro-cpxnum-real x))
(imag (macro-cpxnum-imag x)))
(if (and (##eqv? real 0)
(or (##eqv? imag 1)
(##eqv? imag -1)))
(range-error)
(##* -i (##atanh (##* +i x)))))))
(define-prim (##atan2 y x)
(cond ((or (and (##flonum? x) (##flnan? x))
(and (##flonum? y) (##flnan? y)))
+nan.0)
((##eqv? 0 y)
(if (##exact? x)
(if (##negative? x)
(macro-inexact-+pi)
0)
(if (##negative? (##flcopysign (macro-inexact-+1) x))
(macro-inexact-+pi)
0.)))
((and (##not (##finite? x))
(##not (##finite? y)))
(if (##positive? x)
(##flcopysign (macro-inexact-+pi/4) y)
(##flcopysign (macro-inexact-+3pi/4) y)))
(else
(let ((inexact-x (##exact->inexact x))
(inexact-y (##exact->inexact y)))
(if (and (or (##flonum? x)
(##flonum-full-precision? inexact-x)
(##= x inexact-x))
(or (##flonum? y)
(##flonum-full-precision? inexact-y)
(##= y inexact-y)))
(##flatan inexact-y inexact-x)
;; at least one of x or y is nonzero
;; and at least one of them is not a flonum
(let* ((exact-x (##inexact->exact x))
(exact-y (##inexact->exact y))
(max-arg (##max (##abs exact-x)
(##abs exact-y)))
(normalizer (##expt 2 (##- (##integer-length (##denominator max-arg))
(##integer-length (##numerator max-arg))))))
;; now the largest argument will be about 1.
(##flatan (##exact->inexact (##* normalizer exact-y))
(##exact->inexact (##* normalizer exact-x)))))))))
(define-prim (atan x #!optional (y (macro-absent-obj)))
(macro-force-vars (x)
(if (##eq? y (macro-absent-obj))
(##atan x)
(macro-force-vars (y)
(cond ((##not (##real? x))
(##fail-check-real 1 atan x y))
((##not (##real? y))
(##fail-check-real 2 atan x y))
(else
(##atan2 x y)))))))
;;; Hyperbolic functions
(define-prim (##sinh x)
(define (type-error)
(##fail-check-number 1 sinh x))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
0
(##flsinh (##fixnum->flonum x)))
(##flsinh (##exact-int->flonum x))
(##flsinh (##ratnum->flonum x))
(##flsinh x)
(let ((real (macro-cpxnum-real x))
(imag (macro-cpxnum-imag x)))
(if (and (##flonum? real) (##flonum? imag))
;; fast path for flonum case
(macro-cpxnum-make (##fl* (##flsinh real) (##flcos imag))
(##fl* (##flcosh real) (##flsin imag)))
(macro-cpxnum-make (##* (##sinh real) (##cos imag))
(##* (##cosh real) (##sin imag)))))))
(define-prim (sinh x)
(macro-force-vars (x)
(##sinh x)))
(define-prim (##cosh x)
(define (type-error)
(##fail-check-number 1 cosh x))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
1
(##flcosh (##fixnum->flonum x)))
(##flcosh (##exact-int->flonum x))
(##flcosh (##ratnum->flonum x))
(##flcosh x)
(let ((real (macro-cpxnum-real x))
(imag (macro-cpxnum-imag x)))
(if (and (##flonum? real) (##flonum? imag))
;; fast path for flonum case
(macro-cpxnum-make (##fl* (##flcosh real) (##flcos imag))
(##fl* (##flsinh real) (##flsin imag)))
(macro-cpxnum-make (##* (##cosh real) (##cos imag))
(##* (##sinh real) (##sin imag)))))))
(define-prim (cosh x)
(macro-force-vars (x)
(##cosh x)))
(define-prim (##tanh x)
(define (type-error)
(##fail-check-number 1 tanh x))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
0
(##fltanh (##fixnum->flonum x)))
(##fltanh (##exact-int->flonum x))
(##fltanh (##ratnum->flonum x))
(##fltanh x)
(let ((real (macro-cpxnum-real x))
(imag (macro-cpxnum-imag x)))
(if (##eqv? real 0)
;; the argument of the next ##tan is real
;; (##* +i (##tan (##* -i x)))
(macro-cpxnum-make 0 (##tan imag))
(##ctanh x)))))
(define-prim (tanh x)
(macro-force-vars (x)
(##tanh x)))
;;; Inverse hyperbolic functions
(define-prim (##asinh x)
(define (type-error)
(##fail-check-number 1 asinh x))
(macro-number-dispatch x (type-error)
(if (##fxzero? x)
0
(##flasinh (##fixnum->flonum x)))
(##flasinh (##exact-int->flonum x))
(##flasinh (##ratnum->flonum x))
(##flasinh x)
(##casinh x)))
(define-prim (asinh x)
(macro-force-vars (x)
(##asinh x)))
(define-prim (##acosh x)
(define (type-error)
(##fail-check-number 1 acosh x))
(define (real-case x)
(if (##< x 1)
(##cacosh (macro-cpxnum-make x 0))
(##flacosh (##exact->inexact x))))
(macro-number-dispatch x (type-error)
(if (##fx= x 1)
0
(real-case x))
(real-case x)
(real-case x)
(real-case x)
(##cacosh x)))
(define-prim (acosh x)
(macro-force-vars (x)
(##acosh x)))
(define-prim (##atanh x)
(define (type-error)
(##fail-check-number 1 atanh x))
(define (range-error)
(##raise-range-exception 1 atanh x))
(define (real-case x)
(cond ((##< 1 x)
(##negate (real-case (##negate x))))
((##< x -1)
(##make-rectangular (##fl/ (##fllog1p (##exact->inexact (##/ (##* 4 x)
(##square (##- x 1)))))
4.)
(macro-inexact-+pi/2)))
(else
(##flatanh (##exact->inexact x)))))
(macro-number-dispatch x (type-error)
(case x
((0)
0)
((-1 1)
(range-error))
(else
(real-case x)))
(real-case x)
(real-case x)
(real-case x)
(##catanh x)))
(define-prim (atanh x)
(macro-force-vars (x)
(##atanh x)))
(define-prim (##sqrt x)
(define (type-error)
(##fail-check-number 1 sqrt x))
(define (exact-int-sqrt x)
(if (##negative? x)
(##make-rectangular 0 (exact-int-sqrt (##negate x)))
(let ((y (##exact-int.sqrt x)))
(cond ((##eqv? (##cdr y) 0)
(##car y))
((##not (##exact-int.< (macro-flonum-+m-max-plus-1) x))
;; 0 <= x <= (macro-flonum-+m-max-plus-1), can be
;; converted to flonum exactly so avoids double
;; rounding in next expression. This has a relatively
;; fast path for small integers.
(##flsqrt (##exact-int->flonum x)))
((##not (##< (##car y) (macro-flonum-+m-max-plus-1)))
;; ##exact-int->flonum uses second argument correctly
(##exact-int->flonum (##car y) #t))
(else
;; The integer part of y does not have enough bits accuracy
;; to round it correctly to a flonum, so to
;; make sure (##car y) is big enough in the next call we
;; multiply by (expt 2 (macro-flonum-m-bits-plus-1*2)),
;; which is somewhat extravagant;
;; (expt 2 (+ 1 (macro-flonum-m-bits-plus-1))) should
;; work fine.
(##fl* (macro-flonum-inverse-+m-max-plus-1-inexact)
(exact-int-sqrt
(##arithmetic-shift
x
(macro-flonum-m-bits-plus-1*2)))))))))
(define (ratnum-sqrt x)
(if (##negative? x)
(##make-rectangular 0 (ratnum-sqrt (##negate x)))
(let ((p (macro-ratnum-numerator x))
(q (macro-ratnum-denominator x)))
(let ((sqrt-p (##exact-int.sqrt p))
(sqrt-q (##exact-int.sqrt q)))
(if (and (##zero? (##cdr sqrt-p))
(##zero? (##cdr sqrt-q)))
;; both (abs p) and q are perfect squares and
;; their square roots do not have any common factors
(macro-ratnum-make (##car sqrt-p)
(##car sqrt-q))
(let ((wp (##integer-length p))
(wq (##integer-length q)))
;; for IEEE 754 double precision, we need at least
;; 53 or 54 (I can't seem to work it out) of the
;; leading bits of (sqrt (/ p q)). Here we get
;; about 64 leading bits. We just shift p (either
;; right or left) until it is about 128 bits longer
;; than q (shift must be even), then take the
;; integer square root of the result.
(let* ((shift
(##fxarithmetic-shift-left
(##fxarithmetic-shift-right
(##fx- 128 (##fx- wp wq))
1)
1))
(leading-bits
(##car
(##exact-int.sqrt
(##quotient
(##arithmetic-shift p shift)
q))))
(pre-rounded-result
(if (##fxnegative? shift)
(##arithmetic-shift
leading-bits
(##fx-
(##fxarithmetic-shift-right
shift
1)))
(##ratnum.normalize
leading-bits
(##arithmetic-shift
1
(##fxarithmetic-shift-right
shift
1))))))
(if (##ratnum? pre-rounded-result)
(##ratnum->flonum pre-rounded-result #t)
(##exact-int->flonum pre-rounded-result #t)))))))))
(macro-number-dispatch x (type-error)
(exact-int-sqrt x)
(exact-int-sqrt x)
(ratnum-sqrt x)
(if (##flnegative? x)
(##make-rectangular 0 (##flsqrt (##fl- x)))
(##flsqrt x))
(let ((real (##real-part x))
(imag (##imag-part x)))
(cond ((and (##exact? real)
(##exact? imag)
(let ((discriminant (##sqrt (##+ (##* real real)
(##* imag imag)))))
(and (##exact? discriminant)
(let ((result-real (##sqrt (##/ (##+ real discriminant) 2))))
(and (##exact? result-real)
(##make-rectangular result-real (##/ imag (##* 2 result-real))))))))
=>
values)
(else
(##csqrt (##exact->inexact x)))))))
(define-prim (sqrt x)
(macro-force-vars (x)
(##sqrt x)))
(define-prim (##expt x y)
(define (exact-int-expt x y)
(define (positive-int-expt x y)
;; x is an exact number and y is a positive exact integer
(define (square x)
(##* x x))
(define (expt-aux x y)
;; x is an exact integer (not 0 or 1) and y is a nonzero exact integer
(if (##eqv? y 1)
x
(let ((temp (square (expt-aux x (##arithmetic-shift y -1)))))
(if (##even? y)
temp
(##* x temp)))))
(cond ((or (##eqv? x 0)
(##eqv? x 1))
x)
((eqv? x -1)
(if (##odd? y)
-1
1))
((##ratnum? x)
(macro-ratnum-make
(exact-int-expt (macro-ratnum-numerator x) y)
(exact-int-expt (macro-ratnum-denominator x) y)))
(else
(expt-aux x y))))
(define (invert z)
;; z is exact
(let ((result (##inverse z)))
(if (##not result)
(##raise-range-exception 1 expt x y)
result)))
(if (##negative? y)
(invert (positive-int-expt x (##negate y)))
(positive-int-expt x y)))
(define (complex-expt x y)
(##exp (##* (##log x) y)))
(define (ratnum-expt x y)
;; x is exact-int or ratnum
(cond ((##eqv? x 0)
(if (##negative? y)
(##raise-range-exception 1 expt x y)
0))
((##eqv? x 1)
1)
((##negative? x)
;; We'll do some nice multiples of angles of pi carefully
(case (macro-ratnum-denominator y)
((2)
(##* (##expt (##negate x) y)
(case (##modulo (macro-ratnum-numerator y) 4)
((1)
(macro-cpxnum-+i))
(else ;; (3)
(macro-cpxnum--i)))))
((3)
(##* (##expt (##negate x) y)
(case (##modulo (macro-ratnum-numerator y) 6)
((1)
(macro-cpxnum-+1/2+sqrt3/2i))
((2)
(macro-cpxnum--1/2+sqrt3/2i))
((4)
(macro-cpxnum--1/2-sqrt3/2i))
(else ;; (5)
(macro-cpxnum-+1/2-sqrt3/2i)))))
((6)
(##* (##expt (##negate x) y)
(case (##modulo (macro-ratnum-numerator y) 12)
((1)
(macro-cpxnum-+sqrt3/2+1/2i))
((5)
(macro-cpxnum--sqrt3/2+1/2i))
((7)
(macro-cpxnum--sqrt3/2-1/2i))
(else ;; (11)
(macro-cpxnum-+sqrt3/2-1/2i)))))
;; otherwise, we punt
(else
(complex-expt x y))))
((or (##fixnum? x)
(##bignum? x))
(let* ((y-den (macro-ratnum-denominator y))
(temp (##exact-int.nth-root x y-den)))
(if (##= x (exact-int-expt temp y-den))
(exact-int-expt temp (macro-ratnum-numerator y))
(##flexpt (##exact-int->flonum x)
(##ratnum->flonum y)))))
(else
;; x is a ratnum
(let ((x-num (macro-ratnum-numerator x))
(x-den (macro-ratnum-denominator x))
(y-num (macro-ratnum-numerator y))
(y-den (macro-ratnum-denominator y)))
(let ((temp-num (##exact-int.nth-root x-num y-den)))
(if (##= (exact-int-expt temp-num y-den) x-num)
(let ((temp-den (##exact-int.nth-root x-den y-den)))
(if (##= (exact-int-expt temp-den y-den) x-den)
(exact-int-expt (macro-ratnum-make temp-num temp-den)
y-num)
(##flexpt (##ratnum->flonum x)
(##ratnum->flonum y))))
(##flexpt (##ratnum->flonum x)
(##ratnum->flonum y))))))))
(macro-number-dispatch y (##fail-check-number 2 expt x y)
(macro-number-dispatch x (##fail-check-number 1 expt x y) ;; y a fixnum
(if (##fx= y 0)
1
(exact-int-expt x y))
(if (##fx= y 0)
1
(exact-int-expt x y))
(if (##fx= y 0)
1
(exact-int-expt x y))
(cond ((##fx= y 0)
1)
((##flnan? x)
x)
((##flnegative? x)
;; we do this because (##fixnum->flonum y) is always
;; even for large enough y on 64-bit machines
(let ((abs-result
(##flexpt (##fl- x) (##fixnum->flonum y))))
(if (##fxodd? y)
(##fl- abs-result)
abs-result)))
(else
(##flexpt x (##fixnum->flonum y))))
(cond ((##fx= y 0)
1)
((##fx= y 1)
x)
((##exact? x)
(exact-int-expt x y))
(else
(complex-expt x y))))
(macro-number-dispatch x (##fail-check-number 1 expt x y) ;; y a bignum
(exact-int-expt x y)
(exact-int-expt x y)
(exact-int-expt x y)
(cond ((##flnan? x)
x)
((##flnegative? x)
;; we do this because (##exact-int->flonum y) is always
;; even for large enough y
(let ((abs-result
(##flexpt (##fl- x) (##exact-int->flonum y))))
(if (##odd? y)
(##fl- abs-result)
abs-result)))
(else
(##flexpt x (##exact-int->flonum y))))
(if (##exact? x)
(exact-int-expt x y)
(complex-expt x y)))
(macro-number-dispatch x (##fail-check-number 1 expt x y) ;; y a ratnum
(ratnum-expt x y)
(ratnum-expt x y)
(ratnum-expt x y)
(cond ((##flnan? x)
x)
((##flnegative? x)
(if (##eqv? 2 (macro-ratnum-denominator y))
(let ((magnitude (##flexpt (##fl- x) (##ratnum->flonum y))))
(if (##eqv? 1 (##modulo (macro-ratnum-numerator y) 4))
;; multiple of i
(macro-cpxnum-make 0 magnitude)
;; multiple of -i
(macro-cpxnum-make 0 (##fl- magnitude))))
(complex-expt x y)))
(else
(##flexpt x (##ratnum->flonum y))))
(or (and (##eqv? 2 (macro-ratnum-denominator y))
(or (and (##eqv? 1 (macro-ratnum-numerator y))
(##sqrt x))
(and (##exact? x)
(let ((sqrt-x (##sqrt x)))
(and (##exact? sqrt-x)
(##* sqrt-x (##expt x (##quotient (##- (macro-ratnum-numerator y) 1) 2))))))))
(complex-expt x y)))
(macro-number-dispatch x (##fail-check-number 1 expt x y) ;; y a flonum
(cond ((##flnan? y)
y)
((##eqv? x 0)
(if (##flnegative? y)
(##raise-range-exception 1 expt x y)
0.))
((or (##fxpositive? x)
(macro-flonum-int? y))
(##flexpt (##fixnum->flonum x) y))
(else
(complex-expt x y)))
(cond ((##flnan? y)
y)
((or (##positive? x)
(macro-flonum-int? y))
(##flexpt (##exact-int->flonum x) y))
(else
(complex-expt x y)))
(cond ((##flnan? y)
y)
((or (##positive? x)
(macro-flonum-int? y))
(##flexpt (##ratnum->flonum x) y))
(else
(complex-expt x y)))
(cond ((##flnan? x)
x)
((##flnan? y)
y)
((or (##flpositive? x)
(macro-flonum-int? y))
(##flexpt x y))
(else
(complex-expt x y)))
(cond ((##flnan? y)
y)
(else
(complex-expt x y))))
(macro-number-dispatch x (##fail-check-number 1 expt x y) ;; y a cpxnum
(if (##eqv? x 0)
(let ((real (##real-part y)))
(if (##positive? real)
0
;; If we call (complex-expt 0 y),
;; we'll try to take (##log 0) in complex-expt,
;; so we raise the exception here.
(##raise-range-exception 1 expt x y)))
(complex-expt x y))
(complex-expt x y)
(complex-expt x y)
(complex-expt x y)
(complex-expt x y))))
(define-prim (expt x y)
(macro-force-vars (x y)
(##expt x y)))
(define-prim (##make-rectangular x y)
(cond ((##not (##real? x))
(##fail-check-real 1 make-rectangular x y))
((##not (##real? y))
(##fail-check-real 2 make-rectangular x y))
(else
(let ((real (##real-part x))
(imag (##real-part y)))
(if (##eqv? imag 0)
real
(macro-cpxnum-make real imag))))))
(define-prim (make-rectangular x y)
(macro-force-vars (x y)
(##make-rectangular x y)))
(define-prim (##make-polar x y)
(cond ((##not (##real? x))
(##fail-check-real 1 make-polar x y))
((##not (##real? y))
(##fail-check-real 2 make-polar x y))
(else
(let ((real-x (##real-part x))
(real-y (##real-part y)))
(##make-rectangular (##* real-x (##cos real-y))
(##* real-x (##sin real-y)))))))
(define-prim (make-polar x y)
(macro-force-vars (x y)
(##make-polar x y)))
(define-prim (##real-part x)
(define (type-error)
(##fail-check-number 1 real-part x))
(macro-number-dispatch x (type-error)
x x x x (macro-cpxnum-real x)))
(define-prim (real-part x)
(macro-force-vars (x)
(##real-part x)))
(define-prim (##imag-part x)
(define (type-error)
(##fail-check-number 1 imag-part x))
(macro-number-dispatch x (type-error)
0 0 0 0 (macro-cpxnum-imag x)))
(define-prim (imag-part x)
(macro-force-vars (x)
(##imag-part x)))
(define-prim (##magnitude x)
(define (type-error)
(##fail-check-number 1 magnitude x))
(macro-number-dispatch x (type-error)
(if (##fxnegative? x) (##negate x) x)
(if (##bignum.negative? x) (##negate x) x)
(if (##exact-int.negative? (macro-ratnum-numerator x))
(macro-ratnum-make (##negate (macro-ratnum-numerator x))
(macro-ratnum-denominator x))
x)
(##flabs x)
(let ((real (macro-cpxnum-real x))
(imag (macro-cpxnum-imag x)))
(cond ((and (##flonum? real) (##flonum? imag))
(##cabs x))
;; at least one of real or imag is exact
((and (##exact? real) (##exact? imag))
(##sqrt (##+ (##square real) (##square imag))))
;; one is exact, other is inexact
((and (##finite? real) (##finite? imag))
(##exact->inexact (##sqrt (##+ (##square (##inexact->exact real))
(##square (##inexact->exact imag))))))
;; one is exact, other is not finite inexact
(else
(##cabs (macro-cpxnum-make (##exact->inexact real)
(##exact->inexact imag))))))))
(define-prim (magnitude x)
(macro-force-vars (x)
(##magnitude x)))
(define-prim (##angle x)
(define (type-error)
(##fail-check-number 1 angle x))
(macro-number-dispatch x (type-error)
(if (##fxnegative? x)
(macro-inexact-+pi)
0)
(if (##bignum.negative? x)
(macro-inexact-+pi)
0)
(if (##negative? (macro-ratnum-numerator x))
(macro-inexact-+pi)
0)
(if (##flnegative? (##flcopysign (macro-inexact-+1) x))
(macro-inexact-+pi)
(macro-inexact-+0))
(##atan2 (macro-cpxnum-imag x) (macro-cpxnum-real x))))
(define-prim (angle x)
(macro-force-vars (x)
(##angle x)))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; exact->inexact, inexact->exact
(define-prim (##exact->inexact x)
(define (type-error)
(##fail-check-number 1 exact->inexact x))
(macro-number-dispatch x (type-error)
(##fixnum->flonum x)
(##exact-int->flonum x)
(##ratnum->flonum x)
x
(##make-rectangular (##exact->inexact (macro-cpxnum-real x))
(##exact->inexact (macro-cpxnum-imag x)))))
(define-prim (exact->inexact x)
(macro-force-vars (x)
(##exact->inexact x)))
(define-prim (##inexact->exact x)
(define (type-error)
(##fail-check-number 1 inexact->exact x))
(define (range-error)
(##raise-range-exception 1 inexact->exact x))
(macro-number-dispatch x (type-error)
x
x
x
(if (macro-flonum-rational? x)
(##flonum->exact x)
(range-error))
(let ((real (macro-cpxnum-real x))
(imag (macro-cpxnum-imag x)))
(if (and (macro-noncpxnum-rational? real)
(macro-noncpxnum-rational? imag))
(##make-rectangular (##inexact->exact real)
(##inexact->exact imag))
(range-error)))))
(define-prim (inexact->exact x)
(macro-force-vars (x)
(##inexact->exact x)))
;;; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;;; number->string, string->number
(define-prim (##exact-int->string x rad force-sign?)
(##define-macro (macro-make-block-size)
(let* ((max-rad 16)
(t (make-vector (+ max-rad 1) 0)))
(define max-fixnum 536870911) ;; OK to be conservative
(define (block-size-for rad)
(let loop ((i 0) (rad^i 1))
(let ((new-rad^i (* rad^i rad)))
(if (<= new-rad^i max-fixnum)
(loop (+ i 1) new-rad^i)
i))))
(let loop ((i max-rad))
(if (< 1 i)
(begin
(vector-set! t i (block-size-for i))
(loop (- i 1)))))
`',t))
(define block-size (macro-make-block-size))
(##define-macro (macro-make-rad^block-size)
(let* ((max-rad 16)
(t (make-vector (+ max-rad 1) 0)))
(define max-fixnum 536870911) ;; OK to be conservative
(define (rad^block-size-for rad)
(let loop ((i 0) (rad^i 1))
(let ((new-rad^i (* rad^i rad)))
(if (<= new-rad^i max-fixnum)
(loop (+ i 1) new-rad^i)
rad^i))))
(let loop ((i max-rad))
(if (< 1 i)
(begin
(vector-set! t i (rad^block-size-for i))
(loop (- i 1)))))
`',t))
(define rad^block-size (macro-make-rad^block-size))
(define (make-string-from-last-fixnum rad x len pos)
(let loop ((x x) (len len) (pos pos))
(if (##fx= x 0)
(##make-string len)
(let* ((new-pos
(##fx+ pos 1))
(s
(loop (##fxquotient x rad)
(##fx+ len 1)
new-pos)))
(##string-set!
s
(##fx- (##string-length s) new-pos)
(##string-ref ##digit-to-char-table
(##fx- (##fxremainder x rad))))
s))))
(define (convert-non-last-fixnum s rad x pos)
(let loop ((x x)
(size (##vector-ref block-size rad))
(i (##fx- (##string-length s) pos)))
(if (##fx< 0 size)
(let ((new-i (##fx- i 1)))
(##string-set!
s
new-i
(##string-ref ##digit-to-char-table
(##fxremainder x rad)))
(loop (##fxquotient x rad)
(##fx- size 1)
new-i)))))
(define (make-string-from-fixnums rad lst len pos)
(let loop ((lst lst) (pos pos))
(let ((new-lst (##cdr lst)))
(if (##null? new-lst)
(make-string-from-last-fixnum
rad
(##fx- (##car lst))
(##fx+ len pos)
pos)
(let* ((size
(##vector-ref block-size rad))
(new-pos
(##fx+ pos size))
(s
(loop new-lst new-pos)))
(convert-non-last-fixnum s rad (##car lst) pos)
s)))))
(define (uinteger->fixnums level sqs x lst)
(cond ((and (##null? lst) (##eqv? x 0))
lst)
((##fx= level 0)
(##cons x lst))
(else
(let* ((qr (##exact-int.div x (##car sqs)))
(new-level (##fx- level 1))
(new-sqs (##cdr sqs))
(q (##car qr))
(r (##cdr qr)))
(uinteger->fixnums
new-level
new-sqs
r
(uinteger->fixnums new-level new-sqs q lst))))))
(define (uinteger->string x rad len)
(make-string-from-fixnums
rad
(let ((rad^size
(##vector-ref rad^block-size rad))
(x-length
(##integer-length x)))
(let loop ((level 0)
(sqs '())
(rad^size^2^level rad^size))
(let ((new-level
(##fx+ level 1))
(new-sqs
(##cons rad^size^2^level sqs)))
(if (##fx< x-length
(##fx-
(##fx* (##integer-length rad^size^2^level) 2)
1))
(uinteger->fixnums new-level new-sqs x '())
(let ((new-rad^size^2^level
(##exact-int.square rad^size^2^level)))
(if (##< x new-rad^size^2^level)
(uinteger->fixnums new-level new-sqs x '())
(loop new-level
new-sqs
new-rad^size^2^level)))))))
len
0))
(if (##fixnum? x)
(cond ((##fxnegative? x)
(let ((s (make-string-from-last-fixnum rad x 1 0)))
(##string-set! s 0 #\-)
s))
((##fxzero? x)
(if force-sign?
(##string #\+ #\0)
(##string #\0)))
(else
(if force-sign?
(let ((s (make-string-from-last-fixnum rad (##fx- x) 1 0)))
(##string-set! s 0 #\+)
s)
(make-string-from-last-fixnum rad (##fx- x) 0 0))))
(cond ((##bignum.negative? x)
(let ((s (uinteger->string (##negate x) rad 1)))
(##string-set! s 0 #\-)
s))
(else
(if force-sign?
(let ((s (uinteger->string x rad 1)))
(##string-set! s 0 #\+)
s)
(uinteger->string x rad 0))))))
(define ##digit-to-char-table "0123456789abcdefghijklmnopqrstuvwxyz")
(define-prim (##ratnum->string x rad force-sign?)
(##string-append
(##exact-int->string (macro-ratnum-numerator x) rad force-sign?)
"/"
(##exact-int->string (macro-ratnum-denominator x) rad #f)))
(##define-macro (macro-r6rs-fp-syntax) #t)
(##define-macro (macro-chez-fp-syntax) #f)
(##define-macro (macro-make-10^constants)
(define n 326)
(let ((v (make-vector n)))
(let loop ((i 0) (x 1))
(if (< i n)
(begin
(vector-set! v i x)
(loop (+ i 1) (* x 10)))))
`',v))
(define ##10^-constants
(if (use-fast-bignum-algorithms)
(macro-make-10^constants)
#f))
(define-prim (##flonum-printout v sign-prefix)
;; This algorithm is derived from the paper "Printing Floating-Point
;; Numbers Quickly and Accurately" by Robert G. Burger and R. Kent Dybvig,
;; SIGPLAN'96 Conference on Programming Language Design an Implementation.
;; v is a flonum
;; f is an exact integer (fixnum or bignum)
;; e is an exact integer (fixnum only)
(define (10^ n) ;; 0 <= n < 326
(if (use-fast-bignum-algorithms)
(##vector-ref ##10^-constants n)
(##expt 10 n)))
(define (base-10-log x)
(##define-macro (1/log10) `',(/ (log 10)))
(##fl* (##fllog x) (1/log10)))
(##define-macro (epsilon)
1e-10)
(define (scale r s m+ m- round? v)
;; r is an exact integer (fixnum or bignum)
;; s is an exact integer (fixnum or bignum)
;; m+ is an exact integer (fixnum or bignum)
;; m- is an exact integer (fixnum or bignum)
;; round? is a boolean
;; v is a flonum
(let ((est
(##flonum->fixnum
(##flceiling (##fl- (base-10-log v) (epsilon))))))
(if (##fxnegative? est)
(let ((factor (10^ (##fx- est))))
(fixup (##* r factor)
s
(##* m+ factor)
(##* m- factor)
est
round?))
(let ((factor (10^ est)))
(fixup r
(##* s factor)
m+
m-
est
round?)))))
(define (fixup r s m+ m- k round?)
(if (if round?
(##not (##< (##+ r m+) s))
(##< s (##+ r m+)))
(##cons (##fx+ k 1)
(generate r
s
m+
m-
round?
0))
(##cons k
(generate (##* r 10)
s
(##* m+ 10)
(##* m- 10)
round?
0))))
(define (generate r s m+ m- round? n)
(let* ((dr (##exact-int.div r s))
(d (##car dr))
(r (##cdr dr))
(tc (if round?
(##not (##< (##+ r m+) s))
(##< s (##+ r m+)))))
(if (if round? (##not (##< m- r)) (##< r m-))
(let* ((last-digit
(if tc
(let ((r*2 (##arithmetic-shift r 1)))
(if (or (and (##fxeven? d)
(##= r*2 s)) ;; tie, round d to even
(##< r*2 s))
d
(##fx+ d 1)))
d))
(str
(##make-string (##fx+ n 1))))
(##string-set!
str
n
(##string-ref ##digit-to-char-table last-digit))
str)
(if tc
(let ((str
(##make-string (##fx+ n 1))))
(##string-set!
str
n
(##string-ref ##digit-to-char-table (##fx+ d 1)))
str)
(let ((str
(generate (##* r 10)
s
(##* m+ 10)
(##* m- 10)
round?
(##fx+ n 1))))
(##string-set!
str
n
(##string-ref ##digit-to-char-table d))
str)))))
(define (flonum->exponent-and-digits v)
(let* ((x (##flonum->exact-exponential-format v))
(f (##vector-ref x 0))
(e (##vector-ref x 1))
(round? (##not (##odd? f))))
(if (##fxnegative? e)
(if (and (##not (##fx= e (macro-flonum-e-min)))
(##= f (macro-flonum-+m-min)))
(scale (##arithmetic-shift f 2)
(##arithmetic-shift 1 (##fx- 2 e))
2
1
round?
v)
(scale (##arithmetic-shift f 1)
(##arithmetic-shift 1 (##fx- 1 e))
1
1
round?
v))
(let ((2^e (##arithmetic-shift 1 e)))
(if (##= f (macro-flonum-+m-min))
(scale (##arithmetic-shift f (##fx+ e 2))
4
(##arithmetic-shift 1 (##fx+ e 1))
2^e
round?
v)
(scale (##arithmetic-shift f (##fx+ e 1))
2
2^e
2^e
round?
v))))))
(let* ((x (flonum->exponent-and-digits v))
(e (##car x))
(d (##cdr x)) ;; d = digits
(n (##string-length d))) ;; n = number of digits
(cond ((and (##not (##fx< e 0)) ;; 0<=e<=10
(##not (##fx< 10 e)))
(cond ((##fx= e 0) ;; e=0
;; Format 1: .DDD (0.DDD in chez-fp-syntax)
(##string-append sign-prefix
(if (macro-chez-fp-syntax) "0." ".")
d))
((##fx< e n) ;; e<n
;; Format 2: D.DDD up to DDD.D
(##string-append sign-prefix
(##substring d 0 e)
"."
(##substring d e n)))
((##fx= e n) ;; e=n
;; Format 3: DDD. (DDD.0 in chez-fp-syntax)
(##string-append sign-prefix
d
(if (macro-chez-fp-syntax) ".0" ".")))
(else ;; e>n
;; Format 4: DDD000000. (DDD000000.0 in chez-fp-syntax)
(##string-append sign-prefix
d
(##make-string (##fx- e n) #\0)
(if (macro-chez-fp-syntax) ".0" ".")))))
((and (##not (##fx< e -2)) ;; -2<=e<=-1
(##not (##fx< -1 e)))
;; Format 5: .0DDD or .00DDD (0.0DDD or 0.00DDD in chez-fp-syntax)
(##string-append sign-prefix
(if (macro-chez-fp-syntax) "0." ".")
(##make-string (##fx- e) #\0)
d))
(else
;; Format 6: D.DDDeEEE
;;
;; This is the most general format. We insert a period after
;; the first digit (unless there is only one digit) and add
;; an exponent.
(##string-append sign-prefix
(##substring d 0 1)
(if (##fx= n 1) "" ".")
(##substring d 1 n)
"e"
(##number->string (##fx- e 1) 10))))))
(define-prim (##flonum->string x rad force-sign?)
(define (non-neg-num->str x rad sign-prefix)
(if (##flzero? x)
(##string-append sign-prefix (if (macro-chez-fp-syntax) "0.0" "0."))
(##flonum-printout x sign-prefix)))
(cond ((##flnan? x)
(##string-copy (if (or (macro-r6rs-fp-syntax)
(macro-chez-fp-syntax))
"+nan.0"
"+nan.")))
((##flnegative? (##flcopysign (macro-inexact-+1) x))
(let ((abs-x (##flcopysign x (macro-inexact-+1))))
(cond ((##fl= abs-x (macro-inexact-+inf))
(##string-copy (if (or (macro-r6rs-fp-syntax)
(macro-chez-fp-syntax))
"-inf.0"
"-inf.")))
(else
(non-neg-num->str abs-x rad "-")))))
(else
(cond ((##fl= x (macro-inexact-+inf))
(##string-copy (if (or (macro-r6rs-fp-syntax)
(macro-chez-fp-syntax))
"+inf.0"
"+inf.")))
(force-sign?
(non-neg-num->str x rad "+"))
(else
(non-neg-num->str x rad ""))))))
(define-prim (##cpxnum->string x rad force-sign?)
(let* ((real
(macro-cpxnum-real x))
(real-str
(if (##eqv? real 0) "" (##number->string real rad force-sign?))))
(let ((imag (macro-cpxnum-imag x)))
(cond ((##eqv? imag 1)
(##string-append real-str "+i"))
((##eqv? imag -1)
(##string-append real-str "-i"))
(else
(##string-append real-str
(##number->string imag rad #t)
"i"))))))
(define-prim (##number->string x #!optional (rad 10) (force-sign? #f))
(macro-number-dispatch x '()
(##exact-int->string x rad force-sign?)
(##exact-int->string x rad force-sign?)
(##ratnum->string x rad force-sign?)
(##flonum->string x rad force-sign?)
(##cpxnum->string x rad force-sign?)))
(define-prim (number->string n #!optional (r (macro-absent-obj)))
(macro-force-vars (n r)
(let ((rad (if (##eq? r (macro-absent-obj)) 10 r)))
(if (macro-exact-int? rad)
(if (or (##eqv? rad 2)
(##eqv? rad 8)
(##eqv? rad 10)
(##eqv? rad 16))
(let ((result (##number->string n rad #f)))
(if (##null? result)
(##fail-check-number 1 number->string n r)
result))
(##raise-range-exception 2 number->string n r))
(##fail-check-exact-integer 2 number->string n r)))))
(##define-macro (macro-make-char-to-digit-table)
(let ((t (make-vector 128 99)))
(vector-set! t (char->integer #\#) 0) ;; #\# counts as 0
(let loop1 ((i 9))
(if (not (< i 0))
(begin
(vector-set! t (+ (char->integer #\0) i) i)
(loop1 (- i 1)))))
(let loop2 ((i 25))
(if (not (< i 0))
(begin
(vector-set! t (+ (char->integer #\A) i) (+ i 10))
(vector-set! t (+ (char->integer #\a) i) (+ i 10))
(loop2 (- i 1)))))
`',(list->u8vector (vector->list t))))
(define ##char-to-digit-table (macro-make-char-to-digit-table))
(define-prim (##string->number str #!optional (rad 10) (check-only? #f))
;; The number grammar parsed by this procedure is:
;;
;; <num R E> : <prefix R E> <complex R E>
;;
;; <complex R E> : <real R E>
;; | <real R E> @ <real R E>
;; | <real R E> <sign> <ureal R> i
;; | <real R E> <sign-inf-nan R E> i
;; | <real R E> <sign> i
;; | <sign> <ureal R> i
;; | <sign-inf-nan R E> i
;; | <sign> i
;;
;; <real R E> : <ureal R>
;; | <sign> <ureal R>
;; | <sign-inf-nan R E>
;;
;; <sign-inf-nan R i> : +inf.0
;; | -inf.0
;; | +nan.0
;; <sign-inf-nan R empty> : <sign-inf-nan R i>
;;
;; <ureal R> : <uinteger R>
;; | <uinteger R> / <uinteger R>
;; | <decimal R>
;;
;; <decimal 10> : <uinteger 10> <suffix>
;; | . <digit 10>+ #* <suffix>
;; | <digit 10>+ . <digit 10>* #* <suffix>
;; | <digit 10>+ #+ . #* <suffix>
;;
;; <uinteger R> : <digit R>+ #*
;;
;; <prefix R E> : <radix R E> <exactness E>
;; | <exactness E> <radix R E>
;;
;; <suffix> : <empty>
;; | <exponent marker> <digit 10>+
;; | <exponent marker> <sign> <digit 10>+
;;
;; <exponent marker> : e | s | f | d | l
;; <sign> : + | -
;; <exactness empty> : <empty>
;; <exactness i> : #i
;; <exactness e> : #e
;; <radix 2> : #b
;; <radix 8> : #o
;; <radix 10> : <empty> | #d
;; <radix 16> : #x
;; <digit 2> : 0 | 1
;; <digit 8> : 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7
;; <digit 10> : 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
;; <digit 16> : <digit 10> | a | b | c | d | e | f
(##define-macro (macro-make-exact-10^n-table)
(define max-exact-power-of-10 22) ;; (floor (inexact->exact (/ (log (expt 2 (macro-flonum-m-bits-plus-1))) (log 5))))
(let ((t (make-vector (+ max-exact-power-of-10 1))))
(let loop ((i max-exact-power-of-10))
(if (not (< i 0))
(begin
(vector-set! t i (exact->inexact (expt 10 i)))
(loop (- i 1)))))
`',(list->f64vector (vector->list t))))
(define exact-10^n-table (macro-make-exact-10^n-table))
(##define-macro (macro-make-block-size)
(let* ((max-rad 16)
(t (make-vector (+ max-rad 1) 0)))
(define max-fixnum 536870911) ;; OK to be conservative
(define (block-size-for rad)
(let loop ((i 0) (rad^i 1))
(let ((new-rad^i (* rad^i rad)))
(if (<= new-rad^i max-fixnum)
(loop (+ i 1) new-rad^i)
i))))
(let loop ((i max-rad))
(if (< 1 i)
(begin
(vector-set! t i (block-size-for i))
(loop (- i 1)))))
`',t))
(define block-size (macro-make-block-size))
(##define-macro (macro-make-rad^block-size)
(let* ((max-rad 16)
(t (make-vector (+ max-rad 1) 0)))
(define max-fixnum 536870911) ;; OK to be conservative
(define (rad^block-size-for rad)
(let loop ((i 0) (rad^i 1))
(let ((new-rad^i (* rad^i rad)))
(if (<= new-rad^i max-fixnum)
(loop (+ i 1) new-rad^i)
rad^i))))
(let loop ((i max-rad))
(if (< 1 i)
(begin
(vector-set! t i (rad^block-size-for i))
(loop (- i 1)))))
`',t))
(define rad^block-size (macro-make-rad^block-size))
(define (substring->uinteger-fixnum str rad i j)
;; Simple case: result is known to fit in a fixnum.
(let loop ((i i) (n 0))
(if (##fx< i j)
(let* ((c (##string-ref str i))
(ic (##char->integer c)))
(loop (##fx+ i 1)
(##fx+ (##fx* n rad)