;;;; numbers.scm
;
; Copyright (c) 2008-2011 The CHICKEN Team
; Copyright (c) 2000-2007, Felix L. Winkelmann
; All rights reserved.
;
; Redistribution and use in source and binary forms, with or without
; modification, are permitted provided that the following conditions
; are met:
; 1. Redistributions of source code must retain the above copyright
;    notice, this list of conditions and the following disclaimer.
; 2. Redistributions in binary form must reproduce the above copyright
;    notice, this list of conditions and the following disclaimer in the
;    documentation and/or other materials provided with the distribution.
; 3. The name of the authors may not be used to endorse or promote products
;    derived from this software without specific prior written permission.
;
; THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
; IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
; IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT,
; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
; NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
; DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
; THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
; THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


(declare 
  (disable-interrupts)
  (no-bound-checks)
  (no-procedure-checks))


(module numbers 
    (+ - * / = > < >= <=
       add1 sub1 signum number->string string->number 
       bitwise-and bitwise-ior bitwise-xor bitwise-not arithmetic-shift
       eqv? equal? ; From scheme. Structural & bytevector comparisons Just Work
       exp log sin cos tan atan acos asin expt sqrt conj
       quotient modulo remainder quotient&modulo quotient&remainder
       numerator denominator
       abs max min gcd lcm
       positive? negative? odd? even? zero?
       exact? inexact?
       rationalize
       random randomize
       floor ceiling truncate round
       inexact->exact exact->inexact
       number? complex? real? rational? integer?
       make-rectangular make-polar real-part imag-part magnitude angle
       bignum? ratnum? cflonum? rectnum? compnum? cintnum? cplxnum?
       ;;; XXX These are deprecated. You can do this with the module system
       numbers:bitwise-and numbers:bitwise-ior numbers:bitwise-xor 
       numbers:bitwise-not
       numbers:+ numbers:- numbers:> numbers:< numbers:=
       numbers:>= numbers:<=)

  (import (except scheme
		  + - * / = > < >= <=
		  number->string string->number 
		  exp log sin cos tan atan acos asin expt sqrt
		  quotient modulo remainder
		  numerator denominator
		  abs max min gcd lcm
		  positive? negative? odd? even? zero?
		  exact? inexact?
		  rationalize
		  floor ceiling truncate round
		  inexact->exact exact->inexact
		  number? complex? real? rational? integer?
		  make-rectangular make-polar real-part imag-part magnitude angle)
	  (except chicken add1 sub1 random randomize conj signum
		  bitwise-and bitwise-ior bitwise-xor bitwise-not arithmetic-shift)
	  foreign regex)
  (require-library regex)

(foreign-declare "#include \"numbers-c.h\"")
(foreign-declare "#include \"numbers-c.c\"")

(define-foreign-variable FIX integer)
(define-foreign-variable FLO integer)
(define-foreign-variable BIG integer)
(define-foreign-variable RAT integer)
(define-foreign-variable COMP integer)
(define-foreign-variable NONE integer)

(define-foreign-variable bignum_and_op integer)
(define-foreign-variable bignum_ior_op integer)
(define-foreign-variable bignum_xor_op integer)

;;; Error handling

(define (bad-number loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a number" x))
(define (bad-real loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a real number" x))
(define (bad-ratnum loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a rational number" x))
(define (bad-integer loc x) (##sys#signal-hook #:type-error loc "bad argument type - not an integer" x))
(define (bad-complex/o loc x) (##sys#signal-hook #:type-error loc "bad argument type - complex number has no ordering" x))
(define (bad-base loc x) (##sys#signal-hook #:type-error loc "bad argument type - not a valid base" x))
(define (bad-inexact loc x) (##sys#signal-hook #:type-error loc "bad argument type - inexact number has no exact representation" x))
(define (div/0 loc x y) (##sys#signal-hook #:arithmetic-error loc "division by zero" x y))

(define-inline (%init-tags tagvec) (##core#inline "init_tags" tagvec))
(define-inline (%check-number x) (##core#inline "check_number" x))

(define-inline (assert-number x loc)
  (when (eq? NONE (%check-number x))
    (bad-number loc x) ) )

(define-inline (fix-div/0 x y loc)
  (if (eq? y 0)
      (div/0 loc x y)
      y) )

(define-inline (flo-div/0 x y loc)
  (if (##core#inline "C_flonum_equalp" y 0.0)
      (div/0 loc x y)
      y) )

;;; Constants

(define-constant PI 3.141592653589793)

;;; Primitives

(define-inline (fp/ x y) (##core#inline_allocate ("C_a_i_flonum_quotient" 4) x y))
(define-inline (fp+ x y) (##core#inline_allocate ("C_a_i_flonum_plus" 4) x y))
(define-inline (fp- x y) (##core#inline_allocate ("C_a_i_flonum_difference" 4) x y))
(define-inline (fp* x y) (##core#inline_allocate ("C_a_i_flonum_times" 4) x y))
(define-inline (fp= x y) (##core#inline "C_flonum_equalp" x y))
(define-inline (fp> x y) (##core#inline "C_flonum_greaterp" x y))
(define-inline (fp< x y) (##core#inline "C_flonum_lessp" x y))

(define-inline (%flonum? x) (##core#inline "flonump" x))
(define-inline (%flo-integer? x) (##core#inline "C_i_integerp" x))

(define-inline (complex-real c) (##sys#slot c 1))
(define-inline (complex-imag c) (##sys#slot c 2))
(define-inline (%make-complex r i) (##sys#make-structure 'compnum r i))

(define-inline (rat-numerator c) (##sys#slot c 1))
(define-inline (rat-denominator c) (##sys#slot c 2))
(define-inline (%make-rat r i) (##sys#make-structure 'ratnum r i))

(define-inline (%fix->flo n) (##core#inline_allocate ("fix_to_flo" 4) n))
(define-inline (%big->flo n) (##core#inline_allocate ("big_to_flo" 4) n))

(define %fix+fix (##core#primitive "fix_plus_fix"))
(define %fix+big (##core#primitive "fix_plus_big"))
(define %big+big (##core#primitive "big_plus_big"))

(define %big-neg (##core#primitive "big_neg"))
;; Can't use fxneg because that breaks in the edge case of negating
;; the most negative fixnum.  Yes, 2's complement is fun!
(define %fix-neg (##core#primitive "fix_neg"))

(define %fix-big (##core#primitive "fix_minus_big"))
(define %big-fix (##core#primitive "big_minus_fix"))
(define %big-big (##core#primitive "big_minus_big"))

(define %fix*fix (##core#primitive "fix_times_fix"))
(define %fix*big (##core#primitive "fix_times_big"))
(define %big*big (##core#primitive "big_times_big"))

(define %big-quotient-fix (##core#primitive "big_quotient_fix"))
(define %big-quotient-big (##core#primitive "big_quotient_big"))

(define %big-remainder-fix (##core#primitive "big_remainder_fix"))
(define %big-remainder-big (##core#primitive "big_remainder_big"))

(define %big-divrem-fix (##core#primitive "big_divrem_fix"))
(define %big-divrem-big (##core#primitive "big_divrem_big"))

;; This one should really be part of Chicken, hence the name
(define fxgcd (##core#primitive "C_fixnum_gcd"))
(define biggcd (##core#primitive "big_gcd"))
(define (fpgcd x y) (##core#inline_allocate ("C_a_i_flonum_gcd" 4) x y))

(define-inline (%big-comp-big x y) (##core#inline "big_comp_big" x y))

(define %big-abs (##core#primitive "big_abs"))

(define-inline (%big-odd? x) (##core#inline "big_oddp" x))
(define-inline (%big-negative? x) (##core#inline "big_negp" x))

(define %%expt-0 (##core#primitive "C_expt"))

(define (%expt-0 a b)
  (if (and (negative? a) (not (##sys#integer? b)))
      (%* (%%expt-0 (- a) b)
          (exp (make-complex 0.0 (%* PI b))))
      (%%expt-0 a b)))

(define %quotient-0 (##core#primitive "C_quotient"))

(define %flo->integer (##core#primitive "flo_to_int"))

(define %int-bitwise-int (##core#primitive "int_bitwise_int"))
(define %int-not (##core#primitive "int_not"))
(define %int-shift-fix (##core#primitive "int_shift_fix"))

(define string->number-0 (##core#primitive "C_string_to_number"))
(define number->string-0 (##core#primitive "C_number_to_string"))

(define %big->string (##core#primitive "big_to_string"))
(define %string->big (##core#primitive "string_to_big"))

(define-inline (%subchar s i) (##core#inline "C_subchar" s i))

(define-inline (%fix-randomize n) (##core#inline "fix_randomize" n))
(define-inline (%big-randomize n) (##core#inline "big_randomize" n))

(define-inline (%fix-random n) (##core#inline "fix_random" n))
(define %big-random (##core#primitive "big_random"))


;;; Support macros

(define-syntax switchq
  (syntax-rules (else)
    ((_ "aux" _) (##core#undefined))
    ((_ "aux" _ (else body ...))
     (begin body ...))
    ((_ "aux" tmp (val body ...) more ...)
     (if (eq? tmp val)
	 (begin body ...)
	 (switchq "aux" tmp more ...)))
    ((_ exp body ...)
     (let ((tmp exp))
       (switchq "aux" tmp body ...)))))


;;; Setup

(%init-tags
 (vector 'bignum			; BIG_TAG
	 'ratnum			; RAT_TAG
	 'compnum))			; COMP_TAG

(##sys#gc #f)				; move tag-vector into 2nd generation


;;; Basic arithmetic:

(define (+ . args)
  (if (null? args) 
      0
      (let ((x (##sys#slot args 0))
	    (rest (##sys#slot args 1)))
	(cond ((null? rest)
	       (assert-number x '+) 
	       x)
	      (else
	       (let loop ((args rest) (x x))
		 (if (null? args)
		     x
		     (loop (##sys#slot args 1) (%+ x (##sys#slot args 0))) ) ) ) ) ) ) )

(define (%+ x y)
  (switchq (%check-number x)
    [FIX 
     (switchq (%check-number y)
       [FIX (%fix+fix x y)]
       [FLO (fp+ (%fix->flo x) y)]
       [BIG (%fix+big x y)]
       ;; a/b + c/d = (a*d + b*c)/(b*d)  [with b = 1]
       [RAT (let ((d (rat-denominator y)))
              (%/ (%+ (%* x d) (rat-numerator y)) d))]
       [COMP (%comp+comp (%make-complex x 0) y)]
       [else (bad-number '+ y)] ) ]
    [FLO 
     (switchq (%check-number y)
       [FIX (fp+ x (%fix->flo y))]
       [FLO (fp+ x y)]
       [BIG (fp+ x (%big->flo y))]
       ;; a/b + c/d = (a*d + b*c)/(b*d)  [with b = 1]
       [RAT (let ((d (rat-denominator y)))
              (%/ (%+ (%* x d) (rat-numerator y)) d))]
       [COMP (%comp+comp (%make-complex x 0) y)]
       [else (bad-number '+ y)] ) ]
    [BIG
     (switchq (%check-number y)
       [FIX (%fix+big y x)]
       [FLO (fp+ (%big->flo x) y)]
       [BIG (%big+big x y)]
       ;; a/b + c/d = (a*d + b*c)/(b*d)  [with b = 1]
       [RAT (let ((d (rat-denominator y)))
              (%/ (%+ (%* x d) (rat-numerator y)) d))]
       [COMP (%comp+comp (%make-complex x 0) y)]
       [else (bad-number '+ y)] ) ]
    [RAT
     (switchq (%check-number y)
       [RAT (rat+/- x y %+)]
       [COMP (%comp+comp (%make-complex x 0) y)]
       [NONE (bad-number '+ y)]
       ;; a/b + c/d = (a*d + b*c)/(b*d)  [with d = 1]
       [else (let ((b (rat-denominator x)))
              (%/ (%+ (rat-numerator x) (%* b y)) b))] ) ]
    [COMP
     (switchq (%check-number y)
       [COMP (%comp+comp x y)]
       [NONE (bad-number '+ y)]
       [else (%comp+comp x (%make-complex y 0))] ) ]
    [else (bad-number '+ x)] ) )

(define numbers:+ %+)

(define (%comp+comp x y)
  (let ([r (%+ (complex-real x) (complex-real y))]
	[i (%+ (complex-imag x) (complex-imag y))] )
    (make-complex r i) ) )

(define (- arg1 . args)
  (if (null? args) 
      (switchq (%check-number arg1)
	[FIX (%fix-neg arg1)]
	[FLO (fpneg arg1)]
	[BIG (%big-neg arg1)]
	[RAT (%make-rat (- (rat-numerator arg1)) (rat-denominator arg1))]
	[COMP (%make-complex (%- 0 (complex-real arg1)) (complex-imag arg1))]
	[else (bad-number '- arg1)] )
      (let loop ([args (##sys#slot args 1)] [x (%- arg1 (##sys#slot args 0))])
	(if (null? args)
	    x
	    (loop (##sys#slot args 1) (%- x (##sys#slot args 0))) ) ) ) )

(define (%- x y)
  (switchq (%check-number x)
    [FIX 
     (switchq (%check-number y)
       [FIX (let ((n (%fix-neg y)))
              (if (= (%check-number n) BIG) ;; fix-neg(most negative) => bignum
                  (%fix+big x n)
                  (%fix+fix x n)))]
       [FLO (fp- (%fix->flo x) y)]
       [BIG (%fix-big x y)]
       ;; a/b - c/d = (a*d - b*c)/(b*d)  [with b = 1]
       [RAT (let ((d (rat-denominator y)))
              (%/ (%- (%* x d) (rat-numerator y)) d))]
       [COMP (%comp-comp (%make-complex x 0) y)]
       [else (bad-number '- y)] ) ]
    [FLO 
     (switchq (%check-number y)
       [FIX (fp- x (%fix->flo y))]
       [FLO (fp- x y)]
       [BIG (fp- x (%big->flo y))]
       ;; a/b - c/d = (a*d - b*c)/(b*d)  [with b = 1]
       [RAT (let ((d (rat-denominator y)))
              (%/ (%- (%* x d) (rat-numerator y)) d))]
       [COMP (%comp-comp (%make-complex x 0) y)]
       [else (bad-number '- y)] ) ]
    [BIG
     (switchq (%check-number y)
       [FIX (%big-fix x y)]
       [FLO (fp- (%big->flo x) y)]
       [BIG (%big-big x y)]		
       ;; a/b - c/d = (a*d - b*c)/(b*d)  [with b = 1]
       [RAT (let ((d (rat-denominator y)))
              (%/ (%- (%* x d) (rat-numerator y)) d))]
       [COMP (%comp-comp (%make-complex x 0) y)]
       [else (bad-number '- y)] ) ]
    [RAT
     (switchq (%check-number y)
       [RAT (rat+/- x y %-)]
       [COMP (%comp+comp (%make-complex x 0) y)]
       [NONE (bad-number '- y)]
       ;; a/b - c/d = (a*d - b*c)/(b*d)  [with d = 1]
       [else (let ((b (rat-denominator x)))
               (%/ (%- (rat-numerator x) (%* b y)) b))] ) ]
    [COMP
     (switchq (%check-number y)
       [COMP (%comp-comp x y)]
       [NONE (bad-number '- y)] 
       [else (%comp-comp x (%make-complex y 0))] ) ]
    [else (bad-number '- x)] ) )

(define numbers:- %-)

(define (%comp-comp x y)
  (let ([r (%- (complex-real x) (complex-real y))]
	[i (%- (complex-imag x) (complex-imag y))] )
    (make-complex r i) ) )

(define (* . args)
  (if (null? args) 
      1
      (let ((x (##sys#slot args 0))
	    (rest (##sys#slot args 1)))
	(cond ((null? rest)
	       (assert-number x '*) 
	       x)
	      (else
	       (let loop ((args rest) (x x))
		 (if (null? args)
		     x
		     (loop (##sys#slot args 1) (%* x (##sys#slot args 0))) ) ) ) ) ) ) )

(define-inline (%nonrat*rat x y)
  ;; a/b * c/d = a*c / b*d  [with b = 1]
  ;;  =  ((a / g) * c) / (d / g)
  ;; With   g = gcd(a, d)   and  a = x   [Knuth, 4.5.1]
  (let* ((d (rat-denominator y))
         (g (%gcd-0 x d)))
    (ratnum (%* (%quotient x g) (rat-numerator y))
            (%quotient d g))))

(define (%* x y)
  (switchq (%check-number x)
    [FIX
     (switchq (%check-number y)
       [FIX (%fix*fix x y)]
       [FLO (fp* (%fix->flo x) y)]
       [BIG (%fix*big x y)]
       [RAT (%nonrat*rat x y)]
       [COMP (%comp*comp (%make-complex x 0) y)]
       [else (bad-number '* y)] ) ]
    [FLO
     (switchq (%check-number y)
       [FIX (fp* x (%fix->flo y))]
       [FLO (fp* x y)]
       [BIG (fp* x (%big->flo y))]
       ;; TODO: This can be incorrect when the ratnum consists of bignums
       [RAT (fp* x (%exact->inexact y))]
       [COMP (%comp*comp (%make-complex x 0) y)]
       [else (bad-number '* y)] ) ]
    [BIG
     (switchq (%check-number y)
       [FIX (%fix*big y x)]
       [FLO (fp* (%big->flo x) y)]
       [BIG (%big*big x y)]
       [RAT (%nonrat*rat x y)]
       [COMP (%comp*comp (%make-complex x 0) y)]
       [else (bad-number '* y)] ) ]
    [RAT
     (switchq (%check-number y)
       ;; a/b * c/d = a*c / b*d  [generic]
       ;;   = ((a / g1) * (c / g2)) / ((b / g2) * (d / g1))
       ;; With   g1 = gcd(a, d)   and    g2 = gcd(b, c) [Knuth, 4.5.1]
       [RAT (let* ((a (rat-numerator x)) (b (rat-denominator x))
                   (c (rat-numerator y)) (d (rat-denominator y))
                   (g1 (%gcd-0 a d)) (g2 (%gcd-0 b c)))
              (ratnum (%* (%quotient a g1) (%quotient c g2))
                      (%* (%quotient b g2) (%quotient d g1))))]
       [COMP (%comp*comp (%make-complex x 0) y)]
       ;; TODO: This can be incorrect when the ratnum consists of bignums
       [FLO (fp* y (%exact->inexact x))]
       [NONE (bad-number '* y)]
       [else (%nonrat*rat y x)]) ]
    [COMP
     (switchq (%check-number y)
       [COMP (%comp*comp x y)]
       [NONE (bad-number '* y)] 
       [else (%comp*comp x (%make-complex y 0))] ) ]
    [else (bad-number '* x)] ) )

(define (%comp*comp x y)
  (let* ([a (complex-real x)]
	 [b (complex-imag x)]
	 [c (complex-real y)]
	 [d (complex-imag y)] 
	 [r (%- (%* a c) (%* b d))]
	 [i (%+ (%* a d) (%* b c))] )
    (make-complex r i) ) )

(define (/ arg1 . args)
  (if (null? args) 
      (%/ 1 arg1)
      (let loop ([args (##sys#slot args 1)] [x (%/ arg1 (##sys#slot args 0))])
	(if (null? args)
	    x
	    (loop (##sys#slot args 1) (%/ x (##sys#slot args 0))) ) ) ) )

(define-inline (%nonrat/rat x y)
  ;; a/b / c/d = a*d / b*c  [with b = 1]
  ;;   = ((a / g1) * d * sign(a)) / abs(c / g1)
  ;; With   g1 = gcd(a, c)   and   a = x  [Knuth, 4.5.1 ex. 4]
  (let* ((c (rat-numerator y))
         (g (%gcd-0 x c)))
    (%/ (%* (%quotient x g) (rat-denominator y))
        (%quotient c g))))

(define (%/ x y)
  (switchq (%check-number x)
    [FIX 
     (switchq (%check-number y)
       [FIX (let ((g (fxgcd x (fix-div/0 x y '/))))
              (ratnum (fx/ x g) (fx/ y g)))]
       [FLO (fp/ (%fix->flo x) (flo-div/0 x y '/))]
       [BIG (let ((g (%gcd-0 x y)))
              (ratnum (%quotient x g) (%quotient y g)))]
       [RAT (%nonrat/rat x y)]
       [COMP (%comp/comp (%make-complex x 0) y)]
       [else (bad-number '/ y)] ) ]
    [FLO 
     (switchq (%check-number y)
       [FIX (fp/ x (%fix->flo (fix-div/0 x y '/)))]
       [FLO (fp/ x (flo-div/0 x y '/))]
       [BIG (fp/ x (%big->flo y))]
       ;; TODO: This can be incorrect when the ratnum consists of bignums
       [RAT (fp/ x (%exact->inexact y))]
       [COMP (%comp/comp (%make-complex x 0) y)]
       [else (bad-number '/ y)] ) ]
    [BIG
     (switchq (%check-number y)
       [FIX (let ((g (%gcd-0 x (fix-div/0 x y '/))))
              (ratnum (%quotient x g) (%quotient y g)))]
       [FLO (fp/ (%big->flo x) (flo-div/0 x y '/))]
       [BIG (let ((g (%gcd-0 x y)))
              (ratnum (%quotient x g) (%quotient y g)))]
       [RAT (%nonrat/rat x y)]
       [COMP (%comp/comp (%make-complex x 0) y)]
       [else (bad-number '/ y)] ) ]
    [RAT
     (switchq (%check-number y)
       ;; a/b / c/d = a*d / b*c  [generic]
       ;;   = ((a / g1) * (d / g2) * sign(a)) / abs((b / g2) * (c / g1))
       ;; With   g1 = gcd(a, c)   and    g2 = gcd(b, d) [Knuth, 4.5.1 ex. 4]
       [RAT (let* ((a (rat-numerator x)) (b (rat-denominator x))
                   (c (rat-numerator y)) (d (rat-denominator y))
                   (g1 (%gcd-0 a c)) (g2 (%gcd-0 b d)))
              (%/ (%* (%quotient a g1) (%quotient d g2))
                  (%* (%quotient b g2) (%quotient c g1))))]
       [COMP (%comp-comp (%make-complex x 0) y)]
       ;; TODO: This can be incorrect when the ratnum consists of bignums
       [FLO (fp/ (%exact->inexact x) y)]
       [NONE (bad-number '/ y)]
       ;; a/b / c/d = a*d / b*c  [with d = 1]
       ;;   = ((a / g) * sign(a)) / abs(b * (c / g))
       ;; With   g = gcd(a, c)   and  c = y  [Knuth, 4.5.1 ex. 4]
       [else (let* ((a (rat-numerator x))
                    (g (%gcd-0 a y))) ;; TODO: Improve error message if /0
               (%/ (%quotient a g)
                   (%* (rat-denominator x) (%quotient y g))))] ) ]
    [COMP
     (switchq (%check-number y)
       [COMP (%comp/comp x y)]
       [NONE (bad-number '/ y)] 
       [else (%comp/comp x (%make-complex y 0))] ) ]
    [else (bad-number '/ x)] ) )

(define (%comp/comp p q)
  (let* ([a (complex-real p)]
	 [b (complex-imag p)]
	 [c (complex-real q)]
	 [d (complex-imag q)]
	 [r (%+ (%* c c) (%* d d))]
	 [x (%/ (%+ (%* a c) (%* b d)) r)]
	 [y (%/ (%- (%* b c) (%* a d)) r)] )
    (make-complex x y) ) )


;;; Comparisons:

(define (%= x y)
  (switchq (%check-number x)
    [FIX 
     (switchq (%check-number y)
       [FIX (fx= x y)]
       [FLO (fp= (%fix->flo x) y)]
       [BIG #f] ;; Needs bignum representation?  Can't be equal to a fixnum!
       [RAT #f] ;; Rats are never x/1, because those are normalised to just x
       [COMP #f] ;; Comps are only ever equal to other comps
       [else (bad-number '= y)] ) ]
    [FLO 
     (switchq (%check-number y)
       [FIX (fp= x (%fix->flo y))]
       [FLO (fp= x y)]
       [BIG (and (%flo-integer? x) (= (%flo->integer x) y))]
       [RAT (and (not (or (fp= x +inf) (fp= x -inf)))
                 (%= (%flo->rat '= x) y))] ; Compare as ratnums
       [COMP #f] ;; Comps are only ever equal to other comps
       [else (bad-number '= y)] ) ]
    [BIG
     (switchq (%check-number y)
       [FIX #f]  ;; Needs bignum representation?  Can't be equal to a fixnum!
       [FLO (and (%flo-integer? y) (= x (%flo->integer y)))]
       [BIG (fx= (%big-comp-big x y) 0)]
       [RAT #f] ;; Rats are never x/1, because those are normalised to just x
       [COMP #f] ;; Comps are only ever equal to other comps
       [else (bad-number '= y)] ) ]
    [RAT
     (switchq (%check-number y)
       [FIX #f] ;; Rats are never x/1, because those are normalised to just x
       [FLO (and (not (or (fp= y +inf) (fp= y -inf)))
                 (%= x (%flo->rat '= y)))] ; Compare as ratnums
       [BIG #f] ;; Rats are never x/1, because those are normalised to just x
       [RAT (and (%= (rat-numerator x) (rat-numerator y))
                 (%= (rat-denominator x) (rat-denominator y)))]
       [COMP #f] ;; Comps are only ever equal to other comps
       [else (bad-number '= y)] ) ]
    [COMP
     (switchq (%check-number y)
       [COMP (and (= (complex-real x) (complex-real y))
                  (= (complex-imag x) (complex-imag y)))]
       [NONE (bad-number '= y)]
       [else #f] ) ]
    [else (bad-number '= x)] ))

(define numbers:= %=)
(define (= x1 x2 . xs)
  (and (%= x1 x2)
       (let loop ([x x2] [xs xs])
	 (or (null? xs)
	     (let ([h (##sys#slot xs 0)])
	       (and (%= x h)
		    (loop h (##sys#slot xs 1)) ) ) ) ) ))

(define (> x1 x2 . xs)
  (and (%> x1 x2 '>)
       (let loop ([x x2] [xs xs])
	 (or (null? xs)
	     (let ([h (##sys#slot xs 0)])
	       (and (%> x h '>)
		    (loop h (##sys#slot xs 1)) ) ) ) ) ) )

(define (%> x y loc)
  (switchq (%check-number x)
    (FIX 
     (switchq (%check-number y)
       (FIX (fx> x y))
       (FLO (fp> (%fix->flo x) y))
       ;;   x neg?   y neg?   x > y?   reason
       ;;  ---------------------------------------------------------------
       ;;     no       no       no     abs(y) > abs(x), both positive
       ;;     no      yes      yes     (a > b)  true if  (a > 0) & (b < 0)
       ;;    yes       no       no     (a > b)  false if (a < 0) & (b > 0)
       ;;    yes      yes      yes     abs(y) > abs(x), both negative
       ;;
       ;; It follows that x is only bigger than y when y is negative
       (BIG (%big-negative? y))
       ;; a/b > c/d  when  a*d > b*c  [with b = 1]
       (RAT (%> (%* x (rat-denominator y))
                (rat-numerator y) loc))
       (COMP (bad-complex/o loc y))
       (else (bad-number loc y)) ) )
    (FLO
     (switchq (%check-number y)
       (FIX (fp> x (%fix->flo y)))
       (FLO (fp> x y))
       (BIG (or (fp= x +inf)
                (and (not (fp= x -inf))
                     (%> (%flo->rat loc x) y loc)))) ; Compare as ratnums
       ;; a/b > c/d  when  a*d > b*c  [with b = 1]
       (RAT (%> (%* x (rat-denominator y))
                (rat-numerator y) loc))
       (COMP (bad-complex/o loc y))
       (else (bad-number loc y)) ) )
    (BIG 
     (switchq (%check-number y)
       ;;   x neg?   y neg?   x > y?   reason
       ;;  ---------------------------------------------------------------
       ;;     no       no      yes     abs(x) > abs(y), both positive
       ;;     no      yes      yes     (a > b)  true if  (a > 0) & (b < 0)
       ;;    yes       no       no     (a > b)  false if (a < 0) & (b > 0)
       ;;    yes      yes       no     abs(x) > abs(y), both negative
       ;;
       ;; It follows that x is only bigger than y when x is not negative
       (FIX (not (%big-negative? x)))
       (FLO (or (fp= y -inf)
                (and (not (fp= y +inf))
                     (%> x (%flo->rat loc y) loc)))) ; Compare as ratnums
       (BIG (fx> (%big-comp-big x y) 0))
       ;; a/b > c/d  when  a*d > b*c  [with b = 1]
       (RAT (%> (%* x (rat-denominator y))
                (rat-numerator y) loc))
       (COMP (bad-complex/o loc y))
       (else (bad-number loc y)) ) )
    (RAT
     (switchq (%check-number y)
       ;; a/b > c/d  when  a*d > b*c  [generic]
       (RAT (%> (%* (rat-numerator x) (rat-denominator y))
                (%* (rat-denominator x) (rat-numerator y)) loc))
       (COMP (bad-complex/o loc y))
       (NONE (bad-number loc y))
       ;; a/b > c/d  when  a*d > b*c  [with d = 1]
       (else (%> (rat-numerator x)
                 (%* (rat-denominator x) y) loc)) ) )
    (COMP (bad-complex/o loc x))
    (else (bad-number loc x)) ) )

(define (numbers:> x y) (%> x y '>))
(define (numbers:<= x y) (not (%> x y '<=)))

(define (< x1 x2 . xs)
  (and (%< x1 x2 '<)
       (let loop ([x x2] [xs xs])
	 (or (null? xs)
	     (let ([h (##sys#slot xs 0)])
	       (and (%< x h '<)
		    (loop h (##sys#slot xs 1)) ) ) ) ) ) )

(define (%< x y loc)
  (switchq (%check-number x)
    (FIX 
     (switchq (%check-number y)
       (FIX (fx< x y))
       (FLO (fp< (%fix->flo x) y))
       ;;   x neg?   y neg?   x < y?   reason
       ;;  ---------------------------------------------------------------
       ;;     no       no      yes     abs(x) < abs(y), both positive
       ;;     no      yes       no     (a < b)  false if  (a > 0) & (b < 0)
       ;;    yes       no      yes     (a < b)  true if (a < 0) & (b > 0)
       ;;    yes      yes       no     abs(x) < abs(y), both negative
       ;;
       ;; It follows that x is only smaller than y when y is not negative
       (BIG (not (%big-negative? y)))
       ;; a/b < c/d  when  a*d < b*c  [with b = 1]
       (RAT (%< (%* x (rat-denominator y))
                (rat-numerator y) loc))
       (COMP (bad-complex/o loc y))
       (else (bad-number loc y)) ) )
    (FLO
     (switchq (%check-number y)
       (FIX (fp< x (%fix->flo y)))
       (FLO (fp< x y))
       (BIG (or (fp= x -inf)
                (and (not (fp= x +inf))
                     (%< (%flo->rat loc x) y loc)))) ; Compare as ratnums
       ;; a/b < c/d  when  a*d < b*c  [with b = 1]
       (RAT (%< (%* x (rat-denominator y))
                (rat-numerator y) loc))
       (COMP (bad-complex/o loc y))
       (else (bad-number loc y)) ) )
    (BIG 
     (switchq (%check-number y)
       ;;   x neg?   y neg?   x < y?   reason
       ;;  ---------------------------------------------------------------
       ;;     no       no       no     abs(y) < abs(x), both positive
       ;;     no      yes       no     (a < b)  false if  (a > 0) & (b < 0)
       ;;    yes       no      yes     (a < b)  true if (a < 0) & (b > 0)
       ;;    yes      yes      yes     abs(y) < abs(x), both negative
       ;;
       ;; It follows that x is only smaller than y when x is negative
       (FIX (%big-negative? x))
       (FLO (or (fp= y +inf)
                (and (not (fp= y -inf))
                     (%< x (%flo->rat loc y) loc)))) ; Compare as ratnums
       (BIG (fx< (%big-comp-big x y) 0))
       ;; a/b < c/d  when  a*d < b*c  [with b = 1]
       (RAT (%< (%* x (rat-denominator y))
                (rat-numerator y) loc))
       (COMP (bad-complex/o loc y))
       (else (bad-number loc y)) ) )
    (RAT
     (switchq (%check-number y)
       ;; a/b < c/d  when  a*d < b*c  [generic]
       (RAT (%< (%* (rat-numerator x) (rat-denominator y))
                (%* (rat-denominator x) (rat-numerator y)) loc))
       (COMP (bad-complex/o loc y))
       (NONE (bad-number loc y))
       ;; a/b < c/d  when  a*d < b*c  [with d = 1]
       (else (%< (rat-numerator x)
                 (%* (rat-denominator x) y) loc)) ) )
    (COMP (bad-complex/o loc x))
    (else (bad-number loc x)) ) )

(define (numbers:< x y) (%< x y '<))
(define (numbers:>= x y) (not (%< x y '>=)))

(define (>= x1 x2 . xs)
  (and (not (%< x1 x2 '>=))
       (let loop ([x x2] [xs xs])
	 (or (null? xs)
	     (let ([h (##sys#slot xs 0)])
	       (and (not (%< x h '>=))
		    (loop h (##sys#slot xs 1)) ) ) ) ) ) )

(define (<= x1 x2 . xs)
  (and (not (%> x1 x2 '<=))
       (let loop ([x x2] [xs xs])
	 (or (null? xs)
	     (let ([h (##sys#slot xs 0)])
	       (and (not (%> x h '<=))
		    (loop h (##sys#slot xs 1)) ) ) ) ) ) )


;;; Complex numbers

(define (make-complex r i)
  (if (or (eq? i 0) (and (%flonum? i) (fp= i 0.0)))
      r
      (%make-complex r i) ) )

(define (make-rectangular r i)
  (switchq (%check-number r)
    (COMP (bad-real 'make-rectangular r))
    (NONE (bad-number 'make-rectangular r)) )
  (switchq (%check-number i)
    (COMP (bad-real 'make-rectangular i))
    (NONE (bad-number 'make-rectangular i)) )
  (make-complex r i) )

(define (%make-polar r phi)
  (switchq (%check-number r)
    (COMP (bad-real 'make-polar r))
    (NONE (bad-number 'make-polar r)) )
  (switchq (%check-number phi)
    (COMP (bad-real 'make-polar phi))
    (NONE (bad-number 'make-polar phi)) )
  (make-complex (%* r (##core#inline_allocate ("C_a_i_cos" 4) phi))
		(%* r (##core#inline_allocate ("C_a_i_sin" 4) phi))))

(define make-polar %make-polar)

(define (real-part x)
  (switchq (%check-number x)
    (COMP (complex-real x))
    (NONE (bad-number 'real-part x))
    (else x) ) )

(define (imag-part x)
  (switchq (%check-number x)
    (COMP (complex-imag x))
    (NONE (bad-number 'imag-part x))
    (else 0) ) )

(define (%magnitude x)
  (switchq (%check-number x)
    (COMP (##core#inline_allocate 
	   ("C_a_i_sqrt" 4) 
	   (let ((r (complex-real x))
		 (i (complex-imag x)) )
	     (%+ (%* r r) (%* i i)) ) ) )
    (NONE (bad-number 'magnitude x))
    (else (%abs x)) ) )

(define magnitude %magnitude)

(define (%angle x)
  (switchq (%check-number x)
    (COMP (##core#inline_allocate ("C_a_i_atan2" 4) (complex-imag x) (complex-real x)))
    (NONE (bad-number 'angle x))
    (else (##core#inline_allocate ("C_a_i_atan2" 4) 0 x)) ) )

(define angle %angle)


;;; Rationals

(define (ratnum m n)
  (cond
   ((eq? n 1) m)
   ((eq? n -1) (- m))
   ((negative? n) (%make-rat (- m) (- n)))
   (else (%make-rat m n))))

;; Knuth, 4.5.1
(define (rat+/- x y op)
  (let ((a (rat-numerator x)) (b (rat-denominator x))
        (c (rat-numerator y)) (d (rat-denominator y)))
    (let ((g1 (%gcd-0 b d)))
      (cond
       ((eq? g1 1) (%make-rat (op (%* a d) (%* b c)) (%* b d)))
       ;; Save a quotient and multiplication if the gcd is equal
       ;; to one of the denominators since quotient of b or d and g1 = 1
       ;; TODO: Check properties of the gcd to see if g2 and t are needed
       ((%= g1 b) (let* ((t (op (%* a (%quotient d g1)) c))
                         (g2 (%gcd-0 t g1)))
                    (ratnum (%quotient t g2) (%quotient d g2))))
       ((%= g1 d) (let* ((b/g1 (%quotient b g1))
                         (t (op a (%* c b/g1))) ;; Is this worth it?
                         (g2 (%gcd-0 t g1)))
                    (ratnum (%quotient t g2)
                            (%* b/g1 (%quotient d g2)))))
       (else (let* ((b/g1 (%quotient b g1))
                    (t (op (%* a (%quotient d g1))
                           (%* c b/g1)))
                    (g2 (%gcd-0 t g1)))
               (%make-rat (%quotient t g2)
                          (%* b/g1 (%quotient d g2)))))))))

(define (numerator x)
  (switchq (%check-number x)
    (FIX x)
    (FLO (if (%flo-integer? x) x (bad-ratnum 'numerator x)) )
    (BIG x)
    (RAT (rat-numerator x))
    (COMP (bad-ratnum 'numerator x))
    (else (bad-number 'numerator x)) ) )

(define (denominator x)
  (switchq (%check-number x)
    (FIX 1)
    (FLO (if (%flo-integer? x) 1 (bad-ratnum 'denominator x)) )
    (BIG 1)
    (RAT (rat-denominator x))
    (COMP (bad-ratnum 'denominator x))
    (else (bad-number 'denominator x)) ) )


;;; Enhanced versions of other standard procedures

(define (%abs x)
  (switchq (%check-number x)
    (FIX (if (fx< x 0) (%fix-neg x) x))
    (FLO (##core#inline_allocate ("C_a_i_abs" 4) x))
    (BIG (%big-abs x))
    (RAT (%make-rat (%abs (rat-numerator x)) (rat-denominator x)))
    (COMP (##sys#signal-hook #:type-error 'abs "can not compute absolute value of complex number" x))
    (NONE (bad-number 'abs x)) ) )

(define abs %abs)

(define (number? x)
  (switchq (%check-number x)
    (NONE #f)
    (else #t) ) )

(set! ##sys#number? number?)
(define complex? number?)

(define (real? x)
  (switchq (%check-number x)
    (COMP #f)
    (NONE #f)
    (else #t) ) )

(define rational? real?)

(define (%integer? x)
  (switchq (%check-number x)
    (FIX #t)
    (FLO (%flo-integer? x))
    (BIG #t)
    (else #f) ) )

(set! ##sys#integer? %integer?)
(define integer? %integer?)

(define (%exact? x)
  (switchq (%check-number x)
    (FLO #f)
    (COMP (and (%exact? (complex-real x)) (%exact? (complex-imag x))))
    (NONE (bad-number 'exact? x))
    (else #t) ) )

(define exact? %exact?)
(define ##sys#exact? %exact?)

(define (%inexact? x)
  (switchq (%check-number x)
    (FLO #t)
    (COMP (and (%inexact? (complex-real x)) (%inexact? (complex-imag x))))
    (NONE (bad-number 'inexact? x))
    (else #f) ) )

(define inexact? %inexact?)
(define ##sys#inexact? %inexact?)

(define (positive? x) (%> x 0 'positive?))
(define (negative? x) (%< x 0 'negative?))

(define (%zero? x)
  (switchq (%check-number x)
    (FIX (eq? x 0))
    (FLO (fp= x 0.0))
    (NONE (bad-number 'zero? x))
    (else #f) ) )

(define zero? %zero?)

(define (odd? x)
  (switchq (%check-number x)
    (FIX (##core#inline "C_i_oddp" x))
    (FLO (##core#inline "C_i_oddp" x))
    (BIG (%big-odd? x))
    (else (bad-integer 'odd? x)) ) )

(define (even? x)
  (switchq (%check-number x)
    (FIX (##core#inline "C_i_evenp" x))
    (FLO (##core#inline "C_i_evenp" x))
    (BIG (not (%big-odd? x)))
    (else (bad-integer 'even? x)) ) )

(define (max x1 . xs)
  (let ((i (%flonum? x1)))
    (let loop ((m x1) (xs xs))
      (if (null? xs)
	  (if i (%exact->inexact m) m)
	  (let ((h (##sys#slot xs 0)))
	    (switchq (%check-number h)
	      (FLO (set! i #t))
	      (COMP (bad-complex/o 'max h)) )
	    (loop (if (%> h m 'max) h m) (##sys#slot xs 1)) ) ) ) ) )

(define (min x1 . xs)
  (let ((i (%flonum? x1)))
    (let loop ((m x1) (xs xs))
      (if (null? xs)
	  (if i (%exact->inexact m) m)
	  (let ((h (##sys#slot xs 0)))
	    (switchq (%check-number h)
	      (FLO (set! i #t))
	      (COMP (bad-complex/o 'min h)) )
	    (loop (if (%< h m 'min) h m) (##sys#slot xs 1)) ) ) ) ) )

(define (%quotient x y)
  (switchq (%check-number x)
    (FIX (switchq (%check-number y)
           ;; fix-quotient-big always returns 0 since abs(x) < abs(y)
           ;; But take care of MOST_NEGATIVE_FIXNUM (grrr!)
           (BIG (if (bignum? (- x)) -1 0))
           (RAT (bad-integer 'quotient y)) ; Perhaps convert to flonum?
           (NONE (bad-number 'quotient y))
           (else (%quotient-0 x y))))
    (BIG (switchq (%check-number y)
           (FIX (%big-quotient-fix x (fix-div/0 x y 'quotient)))
           (BIG (%big-quotient-big x y))
           (FLO (if (not (%flo-integer? y))
                    (%quotient-0 (%big->flo x) y) ; Could overflow
                    (%exact->inexact
                     (%quotient x (%flo->integer (flo-div/0 x y 'quotient))))))
           (NONE (bad-number 'quotient y))
           (else (bad-integer 'quotient y))))
    (FLO (switchq (%check-number y)
           (BIG (if (%flo-integer? x)
                    (%exact->inexact (%quotient (%flo->integer x) y))
                    (%quotient-0 x (%big->flo y)))) ; Will probably overflow
           (RAT (bad-integer 'quotient y))
           (NONE (bad-number 'quotient y))
           (else (%quotient-0 x y))))
    (NONE (bad-number 'quotient x))
    (else (bad-integer 'quotient x))))

(define quotient %quotient)

(define (%remainder x y)
  (switchq (%check-number x)
    [FIX (switchq (%check-number y)
           [FIX (fx- x (fx* (fx/ x y) y))]
           [FLO (let ((flx (%fix->flo x)))
                  (fp- flx (fp* (##sys#floor (fp/ flx (flo-div/0 x y 'remainder))) y)))]
           ;; If abs(x) < abs(y), then remainder is always just x
           ;; But again, take care of MOST_NEGATIVE_FIXNUM
           [BIG (if (bignum? (- x)) 0 x)]
           [else (bad-integer 'remainder y)])]
    [FLO (switchq (%check-number y)
           [FLO (fp- x (fp* (##sys#floor (fp/ x (flo-div/0 x y 'remainder))) y))]
           [FIX (let ((fly (%fix->flo (fix-div/0 x y 'remainder))))
                  (fp- x (fp* (##sys#floor (fp/ x fly)) fly)))]
           [BIG (if (%flo-integer? x)
                    (%exact->inexact (%remainder (%flo->integer x) y))
                    (%remainder x (%big->flo y)))] ; Could overflow
           [else (bad-integer 'remainder y)])]
    [BIG (switchq (%check-number y)
           [FIX (%big-remainder-fix x (fix-div/0 x y 'remainder))]
           [FLO (if (%flo-integer? y)
                    (%exact->inexact (%remainder x (%flo->integer y)))
                    (%remainder (%big->flo x) y))] ; Could overflow
           [BIG (%big-remainder-big x y)]
           [else (bad-integer 'remainder y)])]
    [else (bad-integer 'remainder x)]) )

(define remainder %remainder)

;; Modulo's sign follows y  (whereas remainder's sign follows x)
(define (modulo x y)
   (let ((r (%remainder x y)))
      (if (%> y 0 'modulo)
	  (if (%< r 0 'modulo)
	      (%+ r y)
	      r)
	  (if (%> r 0 'modulo)
	      (%+ r y)
	      r))))

(define (quotient&remainder x y)
  (switchq (%check-number x)
    [FIX (switchq (%check-number y)
           [FIX (values (fx/ x y) (remainder x y))]
           [FLO (values (quotient x y) (remainder x y))]
           ;; If abs(x) < abs(y), then remainder is always just x
           ;; But again, take care of MOST_NEGATIVE_FIXNUM
           [BIG (if (bignum? (- x)) (values -1 0) (values 0 x))]
           [else (bad-integer 'quotient&remainder y)])]
    [FLO (switchq (%check-number y)
           [FLO (values (quotient x y) (remainder x y))]
           [FIX (values (quotient x y) (remainder x y))]
           [BIG (if (%flo-integer? x)
                    (receive (div rem)
                      (quotient&remainder (%flo->integer x) y)
                      (values (%exact->inexact div) (%exact->inexact rem)))
                    (quotient&remainder x (%big->flo y)))] ; Probably overflows
           [else (bad-integer 'quotient&remainder y)])]
    [BIG (switchq (%check-number y)
           [FIX (%big-divrem-fix x (fix-div/0 x y 'remainder))]
           [FLO (if (%flo-integer? y)
                    (receive (div rem)
                      (quotient&remainder x (%flo->integer y))
                      (values (%exact->inexact div) (%exact->inexact rem)))
                    (quotient&remainder (%big->flo x) y))] ; Probably overflows
           [BIG (%big-divrem-big x y)]
           [else (bad-integer 'quotient&remainder y)])]
    [else (bad-integer 'quotient&remainder x)]))

;; Modulo's sign follows y  (whereas remainder's sign follows x)
(define (quotient&modulo x y)
   (receive (div rem)
     (quotient&remainder x y)
     (if (%> y 0 'modulo)
         (if (%< rem 0 'modulo)
             (values div (%+ rem y))
             (values div rem))
         (if (%> rem 0 'modulo)
             (values div (%+ rem y))
             (values div rem)))))

;; Try to multiply by two until we reach an integer
(define (%float-fraction-length x)
  (do ((x x (fp* x 2.0))
       (i 0 (+ i 1)))
      ((%flo-integer? x) i)
    ;; 3000 -> %float-precision?
    (if (> i 3000) (error "I'm bored." x))))

(define (%flo->rat loc x)
  (define (deliver y d)
    (let ((q (expt 2 (%float-fraction-length y))))
      (if (%exact? q) ; XXX Ever untrue? float-fraction-length returns natnums
          (let ((e (%/ (%/ (%flo->integer
                            (%* y (%exact->inexact q)))
                           q)
                       d)))
             (if (%exact? e)
                e
                (bad-inexact loc x)))
          (bad-inexact loc x))))
  (if (and (< x (%fix->flo 1))    ; watch out for denormalized numbers
           (> x (%fix->flo -1)))
      (deliver (%* x (expt (%fix->flo 2) flonum-precision))
               ;; Can be bignum (is on 32-bit), so must wait until after init.
               ;; We shouldn't need to calculate this every single time, tho..
               (expt 2 flonum-precision))
      (deliver x 1)))

(define (%inexact->exact x)
  (switchq (%check-number x)
    (FIX x)
    (FLO (cond
          ((or (= x +inf) (= x -inf)) (bad-inexact 'inexact->exact x))
          ((%flo-integer? x) (%flo->integer x))
          (else (%flo->rat 'inexact->exact x))))
    (BIG x)
    (RAT x)
    (COMP (make-complex (%inexact->exact (complex-real x))
                        (%inexact->exact (complex-imag x))))
    (NONE (bad-number 'inexact->exact x)) ) )

(define inexact->exact %inexact->exact)
(define ##sys#inexact->exact %inexact->exact)

(define (%exact->inexact x)
  (switchq (%check-number x)
    (FIX (%fix->flo x))
    (FLO x)
    (BIG (%big->flo x))
    (RAT (%/ (%exact->inexact (rat-numerator x))
             (%exact->inexact (rat-denominator x))))
    (COMP (make-complex (%exact->inexact (complex-real x)) (%exact->inexact (complex-imag x))))
    (NONE (bad-number 'exact->inexact x)) ) )

(define exact->inexact %exact->inexact)
(define ##sys#exact->inexact %exact->inexact)

(define (%gcd-0 x y)
  (switchq (%check-number x)
    [FIX (switchq (%check-number y)
           [FIX (fxgcd x y)]
           [FLO (if (%flo-integer? y)
                    (fpgcd (%fix->flo x) y)
                    (bad-integer 'gcd y))]
           [BIG (if (eq? x 0) y (fxgcd x (%remainder y x)))]
           [else (bad-integer 'gcd y)])]
    [FLO (switchq (%check-number y)
           [FIX (if (%flo-integer? x)
                    (fpgcd x (%fix->flo y))
                    (bad-integer 'gcd x))]
           [FLO (if (%flo-integer? x)
                    (if (%flo-integer? y)
                        (fpgcd x (%fix->flo y))
                        (bad-integer 'gcd x))
                    (bad-integer 'gcd x))]
           [BIG (if (%zero? x) y (fpgcd x (%remainder y x)))]
           [else (bad-integer 'gcd y)])]
    [BIG (switchq (%check-number y)
           [FIX (if (eq? y 0) x (fxgcd y (%remainder x y)))]
           [FLO (if (%zero? y) x (fpgcd y (%remainder x y)))]
           [BIG (biggcd x y)]
           [else (bad-integer 'gcd y)])]
    [else (bad-integer 'gcd x)]) )

(define (gcd . ns)
  (if (eq? ns '())
      0
      (let loop ([ns ns] [f #t])
	(let ([head (##sys#slot ns 0)]
	      [next (##sys#slot ns 1)] )
	  (if (null? next)
	      (%abs head)
	      (let ([n2 (##sys#slot next 0)])
		(loop (cons (%gcd-0 head n2) (##sys#slot next 1)) #f) ) ) ) ) ) )

(define (%lcm-0 x y)
  (%quotient (%* x y) (%gcd-0 x y)) )

(define (lcm . ns)
  (if (null? ns)
      1
      (let loop ([ns ns] [f #t])
	(let ([head (##sys#slot ns 0)]
	      [next (##sys#slot ns 1)] )
	  (if (null? next)
	      (%abs head)
	      (let ([n2 (##sys#slot next 0)])
		(loop (cons (%lcm-0 head (##sys#slot next 0)) (##sys#slot next 1)) #f) ) ) ) ) ) )

(define (%floor x)
  (switchq (%check-number x)
    (FIX x)
    (FLO (##sys#floor x))
    (BIG x)
    ;; (floor x) = greatest integer <= x
    (RAT (let* ((n (rat-numerator x))
                (q (quotient n (rat-denominator x))))
           (if (>= n 0)
               q
               (%- q 1))))
    (else (bad-real 'floor x))) )

(define floor %floor)

(define (ceiling x)
  (switchq (%check-number x)
    (FIX x)
    (FLO (##sys#ceiling x))
    (BIG x)
    ;; (ceiling x) = smallest integer >= x
    (RAT (let* ((n (rat-numerator x))
                (q (quotient n (rat-denominator x))))
           (if (>= n 0)
               (%+ q 1)
               q)))
    (else (bad-real 'ceiling x))) )

(define (truncate x)
  (switchq (%check-number x)
    (FIX x)
    (FLO (##sys#truncate x))
    (BIG x)
    ;; (rational-truncate x) = integer of largest magnitude <= (abs x)
    (RAT (%quotient (rat-numerator x) (rat-denominator x)))
    (else (bad-real 'truncate x))) )

(define (round x)
  (switchq (%check-number x)
    (FIX x)
    (FLO (##core#inline_allocate ("C_a_i_flonum_round_proper" 4) x))
    (BIG x)
    (RAT (let* ((x+1/2 (%+ x (%make-rat 1 2)))
                (r (%floor x+1/2)))
           (if (and (%= r x+1/2)
                    (odd? r))
               (%- r 1)
               r)))
    (else (bad-real 'round x)) ) )

(define (find-ratio-between x y)
  (define (sr x y)
    (let ((fx (%inexact->exact (%floor x))) 
	  (fy (%inexact->exact (%floor y))))
      (cond ((not (%< fx x 'rationalize)) (list fx 1))
	    ((%= fx fy) 
	     (let ((rat (sr (%/ 1 (%- y fy)) (%/ 1 (%- x fx)))))
	       (list (%+ (cadr rat) (%* fx (car rat))) (car rat))))
	    (else (list (%+ 1 fx) 1)))))
  (cond ((%< y x 'rationalize)
	 (find-ratio-between y x))
	((not (%< x y 'rationalize))
	 (list x 1))
	((%> x 0 'rationalize)
	 (sr x y))
	((%< y 0 'rationalize) 
	 (let ((rat (sr (%- 0 y) (%- 0 x))))
	   (list (%- 0 (car rat)) (cadr rat))))
	(else '(0 1))))

(define (find-ratio x e) (find-ratio-between (%- x e) (%+ x e)))

(define (rationalize x e) (apply %/ (find-ratio x e))) ; doesn't preserve exactness

(define (%exp n)
  (switchq (%check-number n)
    (NONE (bad-number 'exp n))
    (COMP (%* (##core#inline_allocate ("C_a_i_exp" 4) (complex-real n))
	      (let ((p (complex-imag n)))
		(make-complex
		 (##core#inline_allocate ("C_a_i_cos" 4) p)
		 (##core#inline_allocate ("C_a_i_sin" 4) p) ) ) ) )
    (else (##core#inline_allocate ("C_a_i_exp" 4) (%exact->inexact n)) ) ))

(define exp %exp)

(define (%log n)
  (switchq (%check-number n)
    (NONE (bad-number 'log n))
    (COMP (make-complex (##core#inline_allocate ("C_a_i_log" 4) (%magnitude n)) (%angle n)))
    (else (##core#inline_allocate ("C_a_i_log" 4) (%exact->inexact n)) ) ) )

(define log %log)

(define %i (%make-complex 0 1))
(define %-i (%make-complex 0 -1))
(define %i2 (%make-complex 0 2))

(define (%sin n)
  (switchq (%check-number n)
    (NONE (bad-number 'sin n))
    (COMP (let ((in (%* %i n)))
	    (%/ (%- (%exp in) (%exp (%- 0 in))) %i2)))
    (else (##core#inline_allocate ("C_a_i_sin" 4) (%exact->inexact n)) ) ))

(define sin %sin)

(define (%cos n)
  (switchq (%check-number n)
    (NONE (bad-number 'cos n))
    (COMP (let ((in (%* %i n)))
	    (%/ (%+ (%exp in) (%exp (%- 0 in))) 2) ) )
    (else (##core#inline_allocate ("C_a_i_cos" 4) (%exact->inexact n)) ) ) )

(define cos %cos)

(define (tan n)
  (switchq (%check-number n)
    (NONE (bad-number 'tan n))
    (COMP (%/ (%sin n) (%cos n)))
    (else (##core#inline_allocate ("C_a_i_tan" 4) (%exact->inexact n)) ) ))

(define (%asin n)
  (switchq (%check-number n)
    (NONE (bad-number 'asin n))
    (COMP (%* %-i (%log (%+ (%* %i n) (%sqrt (%- 1 (%* n n)))))))
    (else (##core#inline_allocate ("C_a_i_asin" 4) (%exact->inexact n)) ) ))

(define asin %asin)

(define acos
  (let ((asin1 (##core#inline_allocate ("C_a_i_asin" 4) 1)))
    (lambda (n)
      (switchq (%check-number n)
        (NONE (bad-number 'asin n))
	(COMP (%- asin1 (%asin n)))
	(else (##core#inline_allocate ("C_a_i_acos" 4) (%exact->inexact n)) ) ) ) ) )

(define (atan n #!optional b)
  (switchq (%check-number n)
    (NONE (bad-number 'atan n))
    (COMP (if b
	      (bad-real 'atan n)
	      (let ((in (%* %i n)))
		(%/ (%- (%log (%+ 1 in)) (%log (%- 1 in))) %i2) ) ) )
    (BIG (set! n (%big->flo n)))
    (RAT (set! n (%exact->inexact n))) )
  (if b
      (##core#inline_allocate ("C_a_i_atan2" 4) n b)
      (##core#inline_allocate ("C_a_i_atan" 4) n) ) )

(define (%sqrt n)
  (switchq (%check-number n)
    (NONE (bad-number 'sqrt n))
    (COMP (let ((p (%/ (%angle n) 2))
		(m (##core#inline_allocate ("C_a_i_sqrt" 4) (%magnitude n))) )
	    (make-complex (%* m (%cos p)) (%* m (%sin p)) ) ) )
    (else
     (if (negative? n)
       (make-complex
        0.0
        (##core#inline_allocate ("C_a_i_sqrt" 4) (%exact->inexact (- n))))
       (##core#inline_allocate ("C_a_i_sqrt" 4) (%exact->inexact n)) ) )))

(define sqrt %sqrt)

(define (square x) (%* x x))

(define (%integer-power base e)
  (if (negative? e)
      (ratnum 1 (%integer-power base (- e)))
      (let lp ((res 1) (e2 e))
        (cond
          ((%zero? e2) res)
          ((even? e2) ; recursion is faster than iteration here
           (%* res (square (lp 1 (arithmetic-shift e2 -1)))))
          (else
           (lp (%* res base) (%- e2 1)))))))

(define (expt a b)
  (define (slow-expt a b) (%exp (%* b (%log a))))
  (let ((ta (%check-number a))
	(tb (%check-number b)) )
    (cond ((eq? NONE ta) (bad-number 'expt a))
	  ((eq? NONE tb) (bad-number 'expt b))
	  ((eq? FLO ta)
	   (switchq tb
	     (FIX (%expt-0 a b))
	     (FLO (%expt-0 a b))
	     (BIG (%expt-0 a (%big->flo b)))
	     (RAT (%expt-0 a (%exact->inexact b)))
	     (else (slow-expt a b)) ) )
	  ((eq? FLO tb)
	   (switchq ta
	     (FIX (%expt-0 a b))
	     (FLO (%expt-0 a b))
	     (BIG (%expt-0 (%big->flo a) b))
	     (RAT (%expt-0 (%exact->inexact a) b))
	     (else (slow-expt a b)) ) )
	  ;; is there a better way
	  ((eq? RAT tb)
	   (let ((e (%exact->inexact b)))
	     (switchq ta
	       (FIX (%expt-0 a e))
	       (FLO (%expt-0 a e))
	       (BIG (%expt-0 (%big->flo a) e))
	       (RAT (%expt-0 (%exact->inexact a) e))
	       (else (slow-expt a b)))))
	  ((or (eq? COMP ta) (eq? COMP tb)) (slow-expt a b))
	  ;; this doesn't work that well, yet...
	  (else (%integer-power a b)) ) ) )

(define (conj n)
  (switchq (%check-number n)
    (NONE (bad-number 'conj n))
    (COMP (make-complex (complex-real n) (%- 0 (complex-imag n))))
    (else n) ) )

(define (add1 n) (%+ n 1))
(define (sub1 n) (%- n 1))

(define (signum n)
  (switchq (%check-number n)
    (FIX (cond ((eq? 0 n) 0)
	       ((fx< n 0) -1)
	       (else 1) ) )
    (FLO (cond ((fp= n 0.0) 0.0)
	       ((fp< n 0.0) -1.0)
	       (else 1.0) ) )
    (BIG (if (%big-negative? n) -1 1)) ; Can't be 0; it would be a fixnum then
    (RAT (signum (rat-numerator n)))
    (COMP (bad-complex/o 'signum n))
    (else (bad-number 'signum n)) ) )

(define (%->integer loc n)
  (switchq (%check-number n)
    (FIX n)
    (FLO (if (%integer? n)
             (%flo->integer n)
	     (bad-integer loc n)))
    (BIG n)
    (else (bad-integer loc n)) ) )

(define (numbers:bitwise-and . xs)
  (let loop ((x -1) (xs xs))
    (if (null? xs)
	x
	(let ((xi (##sys#slot xs 0)))
	  (loop
	   (%int-bitwise-int bignum_and_op x (%->integer 'bitwise-and xi))
	   (##sys#slot xs 1) ) ) ) ) )

(define (numbers:bitwise-ior . xs)
  (let loop ((x 0) (xs xs))
    (if (null? xs)
	x
	(let ((xi (##sys#slot xs 0)))
	  (loop
	   (%int-bitwise-int bignum_ior_op x (%->integer 'bitwise-ior xi))
	   (##sys#slot xs 1) ) ) ) ) )

(define (numbers:bitwise-xor . xs)
  (let loop ((x 0) (xs xs))
    (if (null? xs)
	x
	(let ((xi (##sys#slot xs 0)))
	  (loop
	   (%int-bitwise-int bignum_xor_op x (%->integer 'bitwise-xor xi))
	   (##sys#slot xs 1) ) ) ) ) )

(define (numbers:bitwise-not n)
  (%int-not (%->integer 'bitwise-not n)) )

(define bitwise-and numbers:bitwise-and)
(define bitwise-ior numbers:bitwise-ior)
(define bitwise-xor numbers:bitwise-xor)
(define bitwise-not numbers:bitwise-not)

(define (arithmetic-shift n m)
  (let ((n (%->integer 'arithmetic-shift n))
        (m (%->integer 'arithmetic-shift m)))
    (if (bignum? m)
        (##sys#signal-hook #:type-error 'arithmetic-shift
                           "can not shift by bignum amounts" n m)
        (%int-shift-fix n m))) )

(define %number->string
  (let ((string-append string-append))
    (lambda (n #!optional (base 10))
      (unless (memq base '(2 8 10 16)) (bad-base 'number->string base))
      (let numstr ((n n))
	(switchq (%check-number n)
          (FIX (number->string-0 n base))
	  (FLO (number->string-0 n base))
	  (BIG (%big->string n base))
	  (RAT (string-append (numstr (rat-numerator n))
                              "/"
                              (numstr (rat-denominator n))))
	  (COMP (let ((r (complex-real n))
		      (i (complex-imag n)) )
		  (string-append (numstr r) (if (%> i 0 'number->string) "+" "") (numstr i) "i") ) )
	  (else (bad-number 'number->string n)) ) ) ) ) )

(define number->string %number->string)
(define ##sys#number->string %number->string) ; for printer

(define %string->number
  (let ((copy string-copy)
	(string-match-positions string-match-positions)
	(rxp (regexp "([-+0-9A-Fa-f#./]+)@([-+0-9A-Fa-f#./]+)"))
	(rxr0 (regexp "([-+][-+0-9A-Fa-f#./]+)i"))
	(rxr (regexp "([-+0-9A-Fa-f#./]+)([-+][-+0-9A-Fa-f#./]*)i")) )
    (lambda (str #!optional (base 10))
      (##sys#check-string str 'string->number)
      (##sys#check-exact base 'string->number)
      (let ((e 0)
	    (str (copy str)) 
	    (len (##sys#size str)) )
	(call/cc
	 (lambda (return)
	   (define (real str start end)
	     (let ((rat #f))
	       (let loop ((i start))
		 (if (fx>= i end)
		     (if rat
			 (or (and-let* ((r1 (real str start rat))
					(r2 (real str (%+ rat 1) end))
					((not (zero? r2))))
			       (%/ r1 r2))
			     (return #f))
			 (%string->big
			  (##sys#make-c-string (##sys#substring str start end))
			  base))
		     (let ((c (%subchar str i)))
		       (case c
			 ((#\#)
			  (set! e #f)
			  (##core#inline "C_setsubchar" str i #\0) 
			  (loop (fx+ i 1)) )
			 ((#\.) (string->number-0 (##sys#substring str start end) base))
			 ((#\+ #\-)
			  (if (fx> i start)
			      (string->number-0 (##sys#substring str start end) base)
			      (loop (fx+ i 1)) ) )
			 ((#\e #\E) 
			  (if (eq? base 16)
			      (loop (fx+ i 1))
			      (string->number-0 (##sys#substring str start end) base) ) )
			 ((#\/)
			  (if rat #f (begin (set! rat i) (loop (fx+ i 1)))) )
			 (else (loop (fx+ i 1))) ) ) ) ) ) )
	   (define (fin n)
	     (and n
		  (cond ((eq? e 0) n)
			(e (%inexact->exact n))
			(else (%exact->inexact n)) ) ) )
	   (if (string=? "#" str)
	       0.0
	       (and (fx> len 0)
		    (let ((start
			   (let loop ((i 0))
			     (if (fx< i len)
				 (let ((c (%subchar str i)))
				   (if (eq? c #\#)
				       (let* ((i (fx+ i 1))
					      (c (%subchar str i)) )
					 (case c
					   ((#\e)
					    (set! e #t)
					    (loop (fx+ i 1)) )
					   ((#\i)
					    (set! e #f)
					    (loop (fx+ i 1)) )
					   ((#\x)
					    (set! base 16)
					    (loop (fx+ i 1)) )
					   ((#\d)
					    (set! base 10)
					    (loop (fx+ i 1)) )
					   ((#\o)
					    (set! base 8)
					    (loop (fx+ i 1)) )
					   ((#\b)
					    (set! base 2)
					    (loop (fx+ i 1)) )
					   (else (fx- i 1)) ) )
				       i) )
				 i) ) ) )
		      (let ((sub (##sys#substring str start len)))
			(cond ((string=? sub "+i") (fin (make-complex 0 1)))
			      ((string=? sub "-i") (fin (make-complex 0 -1)))
			      (else
			       (let ((m (string-match-positions rxp sub)))
				 (if (and m (= 3 (length m))
					  (pair? (cadr m))
					  (pair? (caddr m)))
				     (and-let* ((a (real sub (caadr m) (cadadr m)))
						(b (real sub (caaddr m) (cadadr (cdr m)))))
				       (fin (%make-polar a b) ) ) 
				     (let* ((m (string-match-positions rxr sub))
					    (lm (and m (length m))))
				       (cond ((and lm (= 3 lm)
						   (pair? (cadr m))
						   (not (caddr m)))
					      (and-let* ((a (real sub (caadr m) (cadadr m))))
						(fin (make-complex 0 a)) ) )
					     ((and lm (= 3 lm) 
						   (pair? (cadr m))
						   (pair? (caddr m)))
					      (let ((r1 (caadr m))
						    (r2 (cadadr m))
						    (i1 (caaddr m))
						    (i2 (cadadr (cdr m))))
						(and-let* ((rp (real sub r1 r2))
							   (ip (if (eq? i2 (fx+ i1 1))
								   (case (%subchar sub i1)
								     ((#\-) -1)
								     ((#\+) 1)
								     (else #f) )
								   (real sub i1 i2))))
						  (fin (make-complex rp ip)) ) ) )
					     (else
					      (let ((m (string-match-positions rxr0 sub)))
						(if (and m (pair? (cdr m)) (pair? (cadr m)))
						    (fin (make-complex 0 (real sub (caadr m) (cadadr m))))
						    (fin (or (real str start len)
							     (string->number-0 str base) )) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) )

(define (randomize #!optional (seed (##sys#fudge 2)))
  (switchq (%check-number seed)
    (FIX (%fix-randomize seed))
    (BIG (%big-randomize seed))
    (else (bad-integer 'randomize seed)) ) )

(define (random n)
  (switchq (%check-number n)
    (FIX (%fix-random n))
    (BIG (%big-random n))
    (else (bad-integer 'random n)) ) )

(define string->number %string->number)
(define ##sys#string->number %string->number) ; for reader


;;; Non-standard type procedures

(define (bignum? x) (eq? (%check-number x) BIG)) ; big number
(define (ratnum? x) (eq? (%check-number x) RAT)) ; rational number
(define (cplxnum? x) (eq? (%check-number x) COMP)) ; complex number

(define (rectnum? x)    ; "exact" complex number
  (and (eq? (%check-number x) COMP)
       (%integer? (complex-real x))
       (%integer? (complex-imag x))))

(define (compnum? x)    ; inexact complex number
  (and (eq? (%check-number x) COMP)
       (%inexact? (complex-real x))
       (%inexact? (complex-imag x))))

(define (cintnum? x)    ; integer number
  (switchq (%check-number x)
    (FIX #t)
    (BIG #t)
    (FLO (%flo-integer? x))
    (COMP (and (%integer? (complex-real x)) (%integer? (complex-imag x))))
    (else #f) ) )

(define (cflonum? x)    ; floatingpoint number
  (switchq (%check-number x)
    (FLO #t)
    (COMP (and (%flonum? (complex-real x)) (%flonum? (complex-imag x))))
    (else #f) ) )

;;; What we provide

(register-feature! #:full-numeric-tower)

)