;;;; numbers.scm
;
; Copyright (c) 2008-2014 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 
  (no-bound-checks)
  (no-procedure-checks))

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

       ;; Specialization hooks, to be removed when integrating into core
       @fixnum-2-plus @integer-2-plus @basic-2-plus @bignum-2-plus
       @fixnum-negate @integer-negate @basic-negate @bignum-negate
       @fixnum-abs @integer-abs @basic-abs @bignum-abs
       @fixnum-signum @flonum-signum @integer-signum @basic-signum @bignum-signum
       @fixnum-2-minus @integer-2-minus @basic-2-minus @bignum-2-minus
       @fixnum-2-times @integer-2-times @basic-2-times @bignum-2-times
       @fixnum-quotient @flonum-quotient @integer-quotient @basic-quotient @bignum-quotient
       @fixnum-remainder @flonum-remainder @integer-remainder @basic-remainder @bignum-remainder
       @fixnum-quotient&remainder @flonum-quotient&remainder
       @integer-quotient&remainder @basic-quotient&remainder
       @bignum-quotient&remainder
       @fixnum-2-gcd @flonum-2-gcd @integer-2-gcd @basic-2-gcd @bignum-2-gcd
       @basic-2-= @integer-2-= @bignum-2-=
       @basic-2-< @integer-2-< @bignum-2-<
       @basic-2-> @integer-2-> @bignum-2->
       @integer-even? @integer-odd? @bignum-even? @bignum-odd?
       @basic-nan? @flonum-nan?
       @basic-finite? @flonum-finite? @basic-infinite? @flonum-infinite?
       @basic-number->string @fixnum->string @flonum->string @integer->string

       @integer-2-bitwise-and @integer-2-bitwise-ior @integer-2-bitwise-xor
       @integer-arithmetic-shift @fixnum-length @integer-length @bignum-length

       ;; These must stay exported because they're called as hooks from C.
       ;; We might later use them in the types db, but it won't help much?
       @extended-2-= @extended-2-> @extended-2-<
       @extended-2-plus @extended-2-minus @extended-2-times
       @extended-abs @extended-signum @extended-negate
       @extended-number->string)

  (import (except scheme
		  + - * / = > < >= <= eqv?
		  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 finite?
		  bitwise-and bitwise-ior bitwise-xor bitwise-not arithmetic-shift)
	  foreign)

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

;;; THESE SERVE AS SPECIALIZATION ENTRIES.
;;; Once this is integrated into core, they can be killed and the
;;; C functions inlined directly.  Remember to fix the allocations!
(define @basic-2-plus (##core#primitive "C_2_basic_plus"))
(define (@fixnum-2-plus a b) (##core#inline_allocate ("C_a_u_i_2_fixnum_plus" 8) a b))
(define @integer-2-plus (##core#primitive "C_u_2_integer_plus"))
(define @bignum-2-plus (##core#primitive "C_u_2_bignum_plus"))

(define @basic-abs (##core#primitive "C_basic_abs"))
(define (@fixnum-abs x) (##core#inline_allocate ("C_a_u_i_fixnum_abs" 8) x))
(define @integer-abs (##core#primitive "C_u_integer_abs"))
(define @bignum-abs (##core#primitive "C_u_bignum_abs"))

(define @basic-signum (##core#primitive "C_basic_signum"))
(define (@fixnum-signum x) (##core#inline "C_u_i_fixnum_signum" x))
(define (@flonum-signum x) (##core#inline_allocate ("C_a_u_i_flonum_signum" 4) x))
(define (@integer-signum x) (##core#inline "C_u_i_integer_signum" x))
(define (@bignum-signum x) (##core#inline "C_u_i_bignum_signum" x))

(define @basic-negate (##core#primitive "C_basic_negate"))
(define (@fixnum-negate x) (##core#inline_allocate ("C_a_u_i_fixnum_negate" 8) x))
(define @integer-negate (##core#primitive "C_u_integer_negate"))
(define @bignum-negate (##core#primitive "C_u_bignum_negate"))

(define @basic-2-minus (##core#primitive "C_2_basic_minus"))
(define (@fixnum-2-minus a b) (##core#inline_allocate ("C_a_u_i_2_fixnum_minus" 8) a b))
(define @integer-2-minus (##core#primitive "C_u_2_integer_minus"))
(define @bignum-2-minus (##core#primitive "C_u_2_bignum_minus"))

(define @basic-2-times (##core#primitive "C_2_basic_times"))
(define (@fixnum-2-times a b) (##core#inline_allocate ("C_a_u_i_2_fixnum_times" 32) a b)) ; Over-estimated
(define @integer-2-times (##core#primitive "C_u_2_integer_times"))
(define @bignum-2-times (##core#primitive "C_u_2_bignum_times"))

(define @basic-quotient (##core#primitive "C_basic_quotient"))
(define (@fixnum-quotient loc a b) (##core#inline_allocate ("C_a_u_i_fixnum_quotient_checked_loc" 8) loc a b))
(define (@flonum-quotient loc a b) (##core#inline_allocate ("C_a_i_flonum_quotient_checked_loc" 4) loc a b))
(define @integer-quotient (##core#primitive "C_u_integer_quotient"))
(define @bignum-quotient (##core#primitive "C_u_bignum_quotient"))

(define @basic-remainder (##core#primitive "C_basic_remainder"))
(define (@fixnum-remainder loc a b) (##core#inline "C_u_i_fixnum_remainder_checked_loc" loc a b))
(define (@flonum-remainder loc a b) (##core#inline_allocate ("C_a_i_flonum_remainder_checked_loc" 4) loc a b))
(define @integer-remainder (##core#primitive "C_u_integer_remainder"))
(define @bignum-remainder (##core#primitive "C_u_bignum_remainder"))

(define @basic-quotient&remainder (##core#primitive "C_basic_divrem"))
(define (@fixnum-quotient&remainder loc a b)
  (values (##core#inline_allocate ("C_a_u_i_fixnum_quotient_checked_loc" 8) loc a b)
          (##core#inline "C_u_i_fixnum_remainder_checked_loc" loc a b)))
(define (@flonum-quotient&remainder loc a b)
  (values (##core#inline_allocate ("C_a_i_flonum_quotient_checked_loc" 4) loc a b)
          (##core#inline_allocate ("C_a_i_flonum_remainder_checked_loc" 4) loc a b)))
(define @integer-quotient&remainder (##core#primitive "C_u_integer_divrem"))
(define @bignum-quotient&remainder (##core#primitive "C_u_bignum_divrem"))

;; There is no specialization for flonums because of the integer check.
;; Perhaps that should change?
(define @basic-2-gcd (##core#primitive "C_2_basic_gcd"))
(define (@fixnum-2-gcd a b) (##core#inline "C_u_i_2_fixnum_gcd" a b))
(define (@flonum-2-gcd a b) (##core#inline_allocate ("C_a_u_i_2_flonum_gcd" 4) a b))
(define @integer-2-gcd (##core#primitive "C_u_2_integer_gcd"))
(define @bignum-2-gcd (##core#primitive "C_u_2_bignum_gcd"))

(define @basic-2-= (##core#primitive "C_2_basic_equalp"))
(define (@integer-2-= a b) (##core#inline "C_u_i_2_integer_equalp" a b))
(define (@bignum-2-= a b) (##core#inline "C_u_i_2_bignum_equalp" a b))

(define @basic-2-< (##core#primitive "C_2_basic_lessp"))
(define (@integer-2-< a b) (##core#inline "C_u_i_2_integer_lessp" a b))
(define (@bignum-2-< a b) (##core#inline "C_u_i_2_bignum_lessp" a b))

(define @basic-2-> (##core#primitive "C_2_basic_greaterp"))
(define (@integer-2-> a b) (##core#inline "C_u_i_2_integer_greaterp" a b))
(define (@bignum-2-> a b) (##core#inline "C_u_i_2_bignum_greaterp" a b))

(define (@integer-even? x) (##core#inline "C_u_i_integer_evenp" x))
(define (@integer-odd? x) (##core#inline "C_u_i_integer_oddp" x))
(define (@fixnum-even? x) (##core#inline "C_i_fixnum_evenp" x))
(define (@fixnum-odd? x) (##core#inline "C_i_fixnum_oddp" x))
(define (@bignum-even? x) (##core#inline "C_u_i_bignum_evenp" x))
(define (@bignum-odd? x) (##core#inline "C_u_i_bignum_oddp" x))

(define @integer-2-bitwise-and (##core#primitive "C_u_2_integer_bitwise_and"))
(define @integer-2-bitwise-ior (##core#primitive "C_u_2_integer_bitwise_ior"))
(define @integer-2-bitwise-xor (##core#primitive "C_u_2_integer_bitwise_xor"))
(define @integer-arithmetic-shift (##core#primitive "C_u_integer_shift"))

(define (@fixnum-length x) (##core#inline "C_u_i_fixnum_length" x))
(define (@bignum-length x) (##core#inline "C_u_i_bignum_length" x))
(define (@integer-length x) (##core#inline "C_u_i_integer_length" x))

(define (@flonum-nan? x) (##core#inline "C_u_i_flonum_nanp" x))
;; Seems silly, but in core's types.db it can improve things considerably
(define (@basic-nan? x) (##core#inline "C_i_nanp" x))
(define (@flonum-finite? x) (##core#inline "C_u_i_flonum_finitep" x))
;; Seems silly, but in core's types.db it can improve things considerably
(define (@basic-finite? x) (##core#inline "C_i_numbers_finitep" x))
(define (@flonum-infinite? x) (##core#inline "C_u_i_flonum_infinitep" x))
;; Seems silly, but in core's types.db it can improve things considerably
(define (@basic-infinite? x) (##core#inline "C_i_numbers_infinitep" x))

(define @basic-number->string (##core#primitive "C_basic_number_to_string"))
(define @fixnum->string (##core#primitive "C_u_fixnum_to_string"))
(define @flonum->string (##core#primitive "C_u_flonum_to_string"))
(define @integer->string (##core#primitive "C_u_integer_to_string"))

(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)

;;; 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-natural loc x) (##sys#signal-hook #:type-error loc "bad argument type - must be an nonnegative integer" x))
(define (bad-positive loc x) (##sys#signal-hook #:type-error loc "bad argument type - must be a positive (non-zero) 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 (bad-exact loc x) (##sys#signal-hook #:type-error loc "bad argument type - must be an exact number" x))
(define (log0 loc x) (##sys#signal-hook #:arithmetic-error loc "log of exact 0 is undefined" x))
(define (expt0 loc x y) (##sys#signal-hook #:arithmetic-error loc "exponent of exact 0 with complex argument is undefined" x y))
(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))

;;; 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_times" 4) x y))
(define-inline (fp= x y) (##core#inline "C_flonum_equalp" x y))

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

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

;;; 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)))
	(if (null? rest)
            (if (number? x) x (bad-number '+ x))
            (let loop ((args rest) (x x))
              (if (null? args)
                  x
                  (loop (##sys#slot args 1) (%+ x (##sys#slot args 0))) ) )  ) ) ) )

(define %+ (##core#primitive "C_2_basic_plus"))

(define (@extended-2-plus x y)
  (cond ((or (cplxnum? x) (cplxnum? y))
         ;; Just add real and imag parts together
         (let ((r (%+ (real-part x) (real-part y)))
               (i (%+ (imag-part x) (imag-part y))) )
           (make-complex r i) ))
        ((ratnum? x)
         (if (ratnum? y)
             (rat+/- '+ %+ x y)
             ;; a/b + c/d = (a*d + b*c)/(b*d)  [with d = 1]
             (let ((b (ratnum-denominator x)))
               (%/ (%+ (ratnum-numerator x) (%* b y)) b))))
        ((ratnum? y)
         ;; a/b + c/d = (a*d + b*c)/(b*d)  [with b = 1]
         (let ((d (ratnum-denominator y)))
           (%/ (%+ (%* x d) (ratnum-numerator y)) d)))
        (else (bad-number '+ y)) ) )

(define (- arg1 . args)
  (if (null? args)
      ((##core#primitive "C_basic_negate") 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 (@extended-negate x)
  (cond ((ratnum? x)
         (%make-ratnum ((##core#primitive "C_u_integer_negate")
                        (ratnum-numerator x))
                       (ratnum-denominator x)))
        ((cplxnum? x)
         (%make-complex ((##core#primitive "C_basic_negate") (compnum-real x))
                        ((##core#primitive "C_basic_negate") (compnum-imag x))))
        (else (bad-number '- x)) ) ) ; loc?

(define %- (##core#primitive "C_2_basic_minus"))

(define (@extended-2-minus x y)
  (cond ((or (cplxnum? x) (cplxnum? y))
         ;; Just subtract real and imag parts from eachother
         (let ((r (%- (real-part x) (real-part y)))
               (i (%- (imag-part x) (imag-part y))))
           (make-complex r i) ))
        ((ratnum? x)
         (if (ratnum? y)
             (rat+/- '- %- x y)
             ;; a/b - c/d = (a*d - b*c)/(b*d)  [with d = 1]
             (let ((b (ratnum-denominator x)))
               (%/ (%- (ratnum-numerator x) (%* b y)) b))))
        ((ratnum? y)
         ;; a/b - c/d = (a*d - b*c)/(b*d)  [with b = 1]
         (let ((d (ratnum-denominator y)))
           (%/ (%- (%* x d) (ratnum-numerator y)) d)))
        (else (bad-number '- y)) ) )

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

(define %* (##core#primitive "C_2_basic_times"))

(define (@extended-2-times x y)
  (define (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 (ratnum-denominator y))
           (g ((##core#primitive "C_2_basic_gcd") '* x d)))
      (ratnum (%* (%quotient '* x g) (ratnum-numerator y))
              (%quotient '* d g))))

  (cond ((or (cplxnum? x) (cplxnum? y))
         (let* ((a (real-part x)) (b (imag-part x))
                (c (real-part y)) (d (imag-part y))
                (r (%- (%* a c) (%* b d)))
                (i (%+ (%* a d) (%* b c))) )
           (make-complex r i) ) )
        ((or (##core#inline "C_i_flonump" x) (##core#inline "C_i_flonump" y))
         ;; This may be incorrect when one is a ratnum consisting of bignums
         (fp* (exact->inexact y) (exact->inexact x))) ; loc?
        ((ratnum? x)
         (if (ratnum? 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]
             (let* ((a (ratnum-numerator x)) (b (ratnum-denominator x))
                    (c (ratnum-numerator y)) (d (ratnum-denominator y))
                    (g1 ((##core#primitive "C_u_2_integer_gcd") a d))
                    (g2 ((##core#primitive "C_u_2_integer_gcd") b c)))
               (ratnum (%* (%quotient '* a g1) (%quotient '* c g2))
                       (%* (%quotient '* b g2) (%quotient '* d g1))))
             (nonrat*rat y x)))
        ((ratnum? y) (nonrat*rat x y))
        (else (bad-number '* x))))

(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 (%/ x y)
  (cond ((and (exact-integer? x) (exact-integer? y))
         (receive (q r) ((##core#primitive "C_u_integer_divrem") '/ x y)
           (if (eq? r 0)
               q
               (let ((g ((##core#primitive "C_u_2_integer_gcd") y r)))
                 (ratnum ((##core#primitive "C_u_integer_quotient") '/ x g)
                         ((##core#primitive "C_u_integer_quotient") '/ y g))))))
        ;; Compnum *must* be checked first
        ((or (cplxnum? x) (cplxnum? y))
         (let* ((a (real-part x)) (b (imag-part x))
                (c (real-part y)) (d (imag-part y))
                (r (%+ (%* c c) (%* d d)))
                (x (%/ (%+ (%* a c) (%* b d)) r))
                (y (%/ (%- (%* b c) (%* a d)) r)) )
           (make-complex x y) ))
        ((or (##core#inline "C_i_flonump" x) (##core#inline "C_i_flonump" y))
         ;; This may be incorrect when one is a ratnum consisting of bignums
         (fp/ (exact->inexact x)
              (if (eq? y 0) (div/0 '/ x y) (exact->inexact y))))
        ((ratnum? x)
         (if (ratnum? 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]
             (let* ((a (ratnum-numerator x)) (b (ratnum-denominator x))
                    (c (ratnum-numerator y)) (d (ratnum-denominator y))
                    (g1 ((##core#primitive "C_u_2_integer_gcd") a c))
                    (g2 ((##core#primitive "C_u_2_integer_gcd") b d)))
               ;; TODO: can't we make-rat here?
               (%/ (%* (%quotient '/ a g1) (%quotient '/ d g2))
                   (%* (%quotient '/ b g2) (%quotient '/ c g1))))
             ;; 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]
             (let* ((a (ratnum-numerator x))
                    (g ((##core#primitive "C_2_basic_gcd") '/ a y)))
               ;; TODO: can't we make-rat here?
               (%/ (%quotient '/ a g)
                   (%* (ratnum-denominator x) (%quotient '/ y g))))))
        ((ratnum? 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 (ratnum-numerator y))
                (g ((##core#primitive "C_2_basic_gcd") '/ x c)))
           ;; TODO: can't we make-rat here?
           (%/ (%* (%quotient '/ x g) (ratnum-denominator y))
               (%quotient '/ c g))))
        ((not (number? x)) (bad-number '/ x))
        (else (bad-number '/ y))) )


;;; Comparisons:

(define (@extended-2-= x y)
  (cond ((or (cplxnum? x) (cplxnum? y))
         (and (%= (real-part x) (real-part y))
              (%= (imag-part x) (imag-part y))))
        ((ratnum? x)
         (cond ((##core#inline "C_i_flonump" y)
                (and (##core#inline "C_u_i_flonum_finitep" x)
                     (%= x (%flo->rat '= y))))
               ((ratnum? y)
                (and (%= (ratnum-numerator x) (ratnum-numerator y))
                     (%= (ratnum-denominator x) (ratnum-denominator y))))
               ((number? y) #f)
               (else (bad-number '= y))))
        ((ratnum? y)
         (cond ((##core#inline "C_i_flonump" x)
                (and (##core#inline "C_u_i_flonum_finitep" x)
                     (%= (%flo->rat '= x) y)))
               ((number? x) #f)
               (else (bad-number '= x))))
        (else (bad-number '= x))))

(define %= (##core#primitive "C_2_basic_equalp"))

(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 (eqv? a b) (##core#inline "C_i_numbers_eqvp" a b))

(define (@extended-2-> loc x y)
  (cond ((cplxnum? x) (bad-complex/o loc x))
        ((cplxnum? y) (bad-complex/o loc y))
        ((##core#inline "C_i_flonump" x)
         (or (fp= x +inf.0)
             (and (##core#inline "C_u_i_flonum_finitep" x)
                  (%> loc (%flo->rat loc x) y))))
        ((##core#inline "C_i_flonump" y)
         (or (fp= y -inf.0)
             (and (##core#inline "C_u_i_flonum_finitep" x)
                  (%> loc x (%flo->rat loc y)))))
        ((and (number? x) (number? y))  ; Compare as two ratnums
         ;; a/b > c/d  when  a*d > b*c  [generic]
         (%> loc
             (%* (numerator x) (denominator y))
             (%* (denominator x) (numerator y))))
        (else (bad-number loc x))))

(define %> (##core#primitive "C_2_basic_greaterp"))

(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 (@extended-2-< loc x y)
  (cond ((cplxnum? x) (bad-complex/o loc x))
        ((cplxnum? y) (bad-complex/o loc y))
        ((##core#inline "C_i_flonump" x)
         (or (fp= x -inf.0)
             (and (##core#inline "C_u_i_flonum_finitep" x)
                  (%< loc (%flo->rat loc x) y))))
        ((##core#inline "C_i_flonump" y)
         (or (fp= y +inf.0)
             (and (##core#inline "C_u_i_flonum_finitep" x)
                  (%< loc x (%flo->rat loc y)))))
        ((and (number? x) (number? y))  ; Compare as two ratnums
         ;; a/b < c/d  when  a*d < b*c  [generic]
         (%< loc
             (%* (numerator x) (denominator y))
             (%* (denominator x) (numerator y))))
        (else (bad-number loc x))))

(define %< (##core#primitive "C_2_basic_lessp"))

(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)
  ;; Could use nan? here, but that would give an error message with bad loc :(
  (and (or (not (##core#inline "C_i_flonump" x1))
           (not (##core#inline "C_u_i_flonum_nanp" x1)))
       (or (not (##core#inline "C_i_flonump" x2))
           (not (##core#inline "C_u_i_flonum_nanp" x2)))
       (not (%< '>= x1 x2))
       (let loop ((x x2) (xs xs))
         (or (null? xs)
             (let ((h (##sys#slot xs 0)))
               (and (or (not (##core#inline "C_i_flonump" h))
                        (not (##core#inline "C_u_i_flonum_nanp" h)))
                    (not (%< '>= x h))
                    (loop h (##sys#slot xs 1)) ) ) ) ) ) )

(define (<= x1 x2 . xs)
  ;; Could use nan? here, but that would give an error message with bad loc :(
  (and (or (not (##core#inline "C_i_flonump" x1))
           (not (##core#inline "C_u_i_flonum_nanp" x1)))
       (or (not (##core#inline "C_i_flonump" x2))
           (not (##core#inline "C_u_i_flonum_nanp" x2)))
       (not (%> '<= x1 x2))
       (let loop ((x x2) (xs xs))
         (or (null? xs)
             (let ((h (##sys#slot xs 0)))
               (and (or (not (##core#inline "C_i_flonump" h))
                        (not (##core#inline "C_u_i_flonum_nanp" h)))
                    (not (%> '<= x h))
                    (loop h (##sys#slot xs 1)) ) ) ) ) ) )


;;; Complex numbers

(define (make-complex r i)
  (if (or (eq? i 0) (and (##core#inline "C_i_flonump" i) (fp= i 0.0)))
      r
      (%make-complex (if (inexact? i) (exact->inexact r) r)
                     (if (inexact? r) (exact->inexact i) i)) ) )

(define (make-rectangular r i)
  (unless (real? r) (bad-real 'make-rectangular r))
  (unless (real? i) (bad-real 'make-rectangular i))
  (make-complex r i) )

(define (make-polar r phi)
  (unless (real? r) (bad-real 'make-rectangular r))
  (unless (real? phi) (bad-real 'make-rectangular phi))
  (let ((fphi (exact->inexact phi)))
    (make-complex (%* r (##core#inline_allocate ("C_a_i_cos" 4) fphi))
                  (%* r (##core#inline_allocate ("C_a_i_sin" 4) fphi)))))

(define (real-part x)
  (cond ((cplxnum? x) (compnum-real x))
        ((number? x) x)
        (else (bad-number 'real-part x))))

(define (imag-part x)
  (cond ((cplxnum? x) (compnum-imag x))
        ((##core#inline "C_i_flonump" x) 0.0)
        ((number? x) 0)
        (else (bad-number 'imag-part x))))

(define (magnitude x)
  (cond ((cplxnum? x)
         (let ((r (compnum-real x))
               (i (compnum-imag x)) )
           (%sqrt 'magnitude (%+ (%* r r) (%* i i))) ))
        ;; Avoid bad error message (location is hardcoded in abs)
        ((number? x) ((##core#primitive "C_basic_abs") x))
        (else (bad-number 'magnitude x))))

(define (angle x)
  (if (number? x)
      (##core#inline_allocate
       ("C_a_i_atan2" 4)
       (exact->inexact (imag-part x))
       (exact->inexact (real-part x)))
      (else (bad-number 'angle x))))

;;; Rationals

(define (ratnum m n)
  (cond
   ((eq? n 1) m)
   ((eq? n -1) ((##core#primitive "C_u_integer_negate") m))
   ((negative? n)
    (%make-ratnum ((##core#primitive "C_u_integer_negate") m)
                  ((##core#primitive "C_u_integer_negate") n)))
   (else (%make-ratnum m n))))

;; Knuth, 4.5.1
;; TODO: Use integer-quotient here
(define (rat+/- loc op x y)
  (let ((a (ratnum-numerator x)) (b (ratnum-denominator x))
        (c (ratnum-numerator y)) (d (ratnum-denominator y)))
    (let ((g1 ((##core#primitive "C_u_2_integer_gcd") b d)))
      (cond
       ((eq? g1 1) (%make-ratnum (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 loc d g1)) c))
                         (g2 ((##core#primitive "C_u_2_integer_gcd") t g1)))
                    (ratnum (%quotient loc t g2) (%quotient loc d g2))))
       ((%= g1 d) (let* ((b/g1 (%quotient loc b g1))
                         (t (op a (%* c b/g1))) ;; Is this worth it?
                         (g2 ((##core#primitive "C_u_2_integer_gcd") t g1)))
                    (ratnum (%quotient loc t g2)
                            (%* b/g1 (%quotient loc d g2)))))
       (else (let* ((b/g1 (%quotient loc b g1))
                    (t (op (%* a (%quotient loc d g1))
                           (%* c b/g1)))
                    (g2 ((##core#primitive "C_u_2_integer_gcd") t g1)))
               (%make-ratnum (%quotient loc t g2)
                             (%* b/g1 (%quotient loc d g2)))))))))

(define (numerator x)
  (cond ((##core#inline "C_i_basic_numberp" x)
         (cond ((not (##core#inline "C_i_flonump" x)) x)
               ((not (finite? x)) (bad-inexact 'numerator x))
               ((##core#inline "C_u_i_fpintegerp" x) x)
               ;; XXX modf?
               (else (exact->inexact (numerator (%flo->rat 'numerator x))))))
        ((ratnum? x) (ratnum-numerator x))
        (else (bad-ratnum 'numerator x))))

(define (denominator x)
  (cond ((##core#inline "C_i_basic_numberp" x)
         (cond ((not (##core#inline "C_i_flonump" x)) 1)
               ((not (finite? x)) (bad-inexact 'denominator x))
               ((##core#inline "C_u_i_fpintegerp" x) 1.0)
               ;; XXX modf?
               (else (exact->inexact (denominator (%flo->rat 'denominator x))))))
        ((ratnum? x) (ratnum-denominator x))
        (else (bad-ratnum 'denominator x))))

;;; Enhanced versions of other standard procedures

(define (@extended-abs x)
  (cond ((ratnum? x)
         (%make-ratnum
          ((##core#primitive "C_u_integer_abs") (ratnum-numerator x))
          (ratnum-denominator x)))
        ((cplxnum? x)
         (##sys#signal-hook
          #:type-error 'abs
          "can not compute absolute value of complex number" x))
        (else (bad-number 'abs x))))

(define (abs x) ((##core#primitive "C_basic_abs") x))

;; TODO: Should this be one thing, with the basic_number calling into
;; the extended number checking procedure?  Makes no real sense...
(define (number? x)
  (or (##core#inline "C_i_basic_numberp" x) (extended-number? x)))

;; TODO: Extend C_i_integerp
(define (integer? x)
  (and (##core#inline "C_i_basic_numberp" x)
       (or (not (##core#inline "C_i_flonump" x))
           (##core#inline "C_u_i_fpintegerp" x))))

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

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

;; All numbers are real, except for compnums
(define (real? x)
  (or (##core#inline "C_i_basic_numberp" x)
      (##sys#structure? x 'ratnum) ) )

(define (rational? x) (and (real? x) (finite? x)))

(define (exact-integer? x)
  (or (##core#inline "C_fixnump" x) (##core#inline "C_i_bignump" x)) )


(define (exact? x)
  (cond ((##core#inline "C_fixnump" x))
        ((##core#inline "C_i_bignump" x))
        ((##core#inline "C_i_flonump" x) #f)
        ((ratnum? x))
        ((cplxnum? x)
         (and (exact? (compnum-real x)) (exact? (compnum-imag x))))
        (else (bad-number 'exact? x))))

(define ##sys#exact? exact?)

(define (inexact? x)
  (cond ((##core#inline "C_fixnump" x) #f)
        ((##core#inline "C_i_bignump" x) #f)
        ((##core#inline "C_i_flonump" x))
        ((ratnum? x) #f)
        ((cplxnum? x)
         (or (inexact? (compnum-real x)) (inexact? (compnum-imag x))))
        (else (bad-number 'inexact? x))))

(define ##sys#inexact? inexact?)

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

(define (zero? x)
  (or (eq? x 0) (and (##core#inline "C_i_flonump" x) (fp= x 0.0))
      (and (not (number? x)) (bad-number 'zero? x))))

(define (odd? x) (##core#inline "C_i_basic_oddp" x))
(define (even? x) (##core#inline "C_i_basic_evenp" x))

;; TODO: Move to C?
(define (max x1 . xs)
  (let loop ((i (##core#inline "C_i_flonump" x1)) (m x1) (xs xs))
    (if (null? xs)
        (if i (exact->inexact m) m)
        (let ((h (##sys#slot xs 0)))
          (loop (or i (##core#inline "C_i_flonump" h))
                (if (%> 'max h m) h m)
                (##sys#slot xs 1)) ) ) )  )

;; TODO: Move to C?
(define (min x1 . xs)
  (let loop ((i (##core#inline "C_i_flonump" x1)) (m x1) (xs xs))
    (if (null? xs)
        (if i (exact->inexact m) m)
        (let ((h (##sys#slot xs 0)))
          (loop (or i (##core#inline "C_i_flonump" h))
                (if (%< 'min h m) h m)
                (##sys#slot xs 1)) ) ) )  )

;; TODO: Inline and remove
(define %quotient (##core#primitive "C_basic_quotient"))

(define (quotient x y)
  ((##core#primitive "C_basic_quotient") 'quotient x y))
(define (truncate-quotient x y)
  ((##core#primitive "C_basic_quotient") 'truncate-quotient x y))

(define (remainder x y)
  ((##core#primitive "C_basic_remainder") 'remainder x y))
(define (truncate-remainder x y)
  ((##core#primitive "C_basic_remainder") 'truncate-remainder x y))

;; Modulo's sign follows y  (whereas remainder's sign follows x)
(define (%modulo loc x y)
  (let ((r ((##core#primitive "C_basic_remainder") loc x y)))
    (if (%> loc y 0)
        (if (%< loc r 0)
            (%+ r y)
            r)
        (if (%> loc r 0)
            (%+ r y)
            r))))

(define (modulo x y) (%modulo 'modulo x y))
(define (floor-remainder x y) (%modulo 'floor-remainder x y))

(define %quotient&remainder )

(define (quotient&remainder x y)
  ((##core#primitive "C_basic_divrem") 'quotient&remainder x y))
(define (truncate/ x y)
  ((##core#primitive "C_basic_divrem") 'truncate/ x y))

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

;; Same as above, but quotient gets adjusted along with the remainder
(define (floor/ x y)
  (receive (div rem)
      ((##core#primitive "C_basic_divrem") 'floor/ x y)
    (if (%> 'floor/ y 0)
        (if (%< 'floor/ rem 0)
            (values (%- div 1) (%+ rem y))
            (values div rem))
        (if (%> 'floor/ rem 0)
            (values (%- div 1) (%+ rem y))
            (values div rem)))))

;; XXX This is a bad bad bad definition; very inefficient But to
;; improve it we would need to provide another implementation of the
;; quotient procedure which floors instead of truncates.
(define (floor-quotient x y)
  (receive (div rem) (floor/ x y) div))


(define (%flo->rat loc x)
  ;; Try to multiply by two until we reach an integer
  (define (float-fraction-length x)
    (do ((x x (fp* x 2.0))
         (i 0 (fx+ i 1)))
        ((##core#inline "C_u_i_fpintegerp" x) i)
      ;; 3000 -> %float-precision?
      ;; XXX TODO: If this gets passed +nan.0, we hit this check
      (if (fx> i 3000) (error "I'm bored." 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 (%/ (%/ ((##core#primitive "C_u_flo_to_int")
                            loc (%* y (exact->inexact q)))
                           q)
                       d)))
            (if (exact? e)
                e
                (bad-inexact loc x)))
          (bad-inexact loc x))))

  (if (and (%< loc x 1.0)         ; watch out for denormalized numbers
           (%> loc x -1.0))
      (deliver (%* x (expt 2.0 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)
  (cond ((exact? x) x)
        ((##core#inline "C_i_flonump" x)
         (if (##core#inline "C_u_i_fpintegerp" x)
             ((##core#primitive "C_u_flo_to_int") 'inexact->exact x)
             (%flo->rat 'inexact->exact x)))
        ((cplxnum? x)
         (make-complex (inexact->exact (compnum-real x))
                       (inexact->exact (compnum-imag x))))
        (else (bad-number 'inexact->exact x))))

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

;; Exponent of the lowest allowed flonum; if we get any lower we get zero.
;; In other words, this is the first (smallest) flonum after 0.
;; Equal to (expt 2.0 (- flonum-minimum-exponent flonum-precision))
(define minimum-denorm-flonum-expt (fx- flonum-minimum-exponent flonum-precision))

(define (exact->inexact x)
  (cond ((##core#inline "C_fixnump" x)
         (##core#inline_allocate ("C_a_i_fix_to_flo" 4) x))
        ((##core#inline "C_i_flonump" x) x)
        ((##core#inline "C_i_bignump" x)
         (##core#inline_allocate ("C_a_u_i_big_to_flo" 4) x))
        ((ratnum? x)
         ;; This tries to keep the numbers within representable ranges
         ;; and tries to drop as few significant digits as possible by
         ;; bringing the two numbers to within the same powers of two.
         ;; See algorithms M & N in Knuth, 4.2.1
         ;;
         ;; TODO: Use (fp/ n d) if both are finite after conversion to
         ;; flonums
         (let* ((n1 (ratnum-numerator x))
                (an ((##core#primitive "C_u_integer_abs") n1))
                (d1 (ratnum-denominator x))
                ;; Approximate distance between the numbers in powers
                ;; of 2 ie, 2^e-1 < n/d < 2^e+1 (e is the *un*biased
                ;; value of e_w in M2)
                ;; XXX: What if b != 2 (ie, flonum-radix is not 2)?
                (e (fx- (integer-length an) (integer-length d1)))
                (rnd (lambda (n d e) ; Here, 1 <= n/d < 2  (normalized) [N5]
                       ;; Cannot shift above the available precision,
                       ;; and can't have an exponent that's below the
                       ;; minimum flonum exponent.
                       (let* ((s (min (fx- flonum-precision 1)
                                      (fx- e minimum-denorm-flonum-expt)))
                              (normalized (%/ (arithmetic-shift n s) d))
                              (r (round normalized))
                              (fraction (exact->inexact r))
                              (exp (fx- e s)))
                         (let ((res (fp* fraction (expt 2.0 exp))))
                           (if (negative? n1) (%- 0 res) res)))))
                (scale (lambda (n d)  ; Here, 1/2 <= n/d < 2   [N3]
                         (if (%< 'exact->inexact n d) ; n/d < 1?
                             ;; Scale left [N3]; only needed once (see note in M3)
                             (rnd (arithmetic-shift n 1) d (fx- e 1))
                             ;; Already normalized
                             (rnd n d e)))))
           ;; After this step, which shifts the smaller number to
           ;; align with the larger, "f" in algorithm N is represented
           ;; in the procedures above by n/d.
           (if (negative? e)
               (scale (arithmetic-shift an (%- 0 e)) d1)
               (scale an (arithmetic-shift d1 e)))))
        ((cplxnum? x)
         (make-complex (exact->inexact (compnum-real x))
                       (exact->inexact (compnum-imag x))))
        (else (bad-number 'exact->inexact x))))

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

(define (gcd . ns)
  (if (eq? ns '())
      0
      (let loop ((head (##sys#slot ns 0))
                 (next (##sys#slot ns 1)))
        (if (null? next)
            (if (integer? head)
                ((##core#primitive "C_basic_abs") head)
                (bad-integer 'gcd head))
            (let ((n2 (##sys#slot next 0)))
              (loop ((##core#primitive "C_2_basic_gcd") 'gcd head n2)
                    (##sys#slot next 1)) ) )  ) ) )

(define (lcm . ns)
  (if (null? ns)
      1
      (let loop ((head (##sys#slot ns 0))
                 (next (##sys#slot ns 1)))
        (if (null? next)
            (if (integer? head)
                ((##core#primitive "C_basic_abs") head)
                (bad-integer 'gcd head))
            (let ((n2 (##sys#slot next 0)))
              (loop (%quotient
                     'lcm (%* head n2)
                     ((##core#primitive "C_2_basic_gcd") 'lcm head n2))
                    (##sys#slot next 1)) ) )  ) ) )

(define (floor x)
  (cond ((exact-integer? x) x)
        ((##core#inline "C_i_flonump" x) (fpfloor x))
        ;; (floor x) = greatest integer <= x
        ((ratnum? x) (let* ((n (ratnum-numerator x))
                            (q (%quotient 'floor n (ratnum-denominator x))))
                       (if (>= n 0) q (%- q 1))))
        (else (bad-real 'floor x))))

(define (ceiling x)
  (cond ((exact-integer? x) x)
        ((##core#inline "C_i_flonump" x) (fpceiling x))
        ;; (ceiling x) = smallest integer >= x
        ((ratnum? x) (let* ((n (ratnum-numerator x))
                            (q (%quotient 'ceiling n (ratnum-denominator x))))
                       (if (>= n 0) (%+ q 1) q)))
        (else (bad-real 'ceiling x))) )

(define (truncate x)
  (cond ((exact-integer? x) x)
        ((##core#inline "C_i_flonump" x) (fptruncate x))
        ;; (rational-truncate x) = integer of largest magnitude <= (abs x)
        ((ratnum? x) (%quotient 'truncate (ratnum-numerator x)
                                (ratnum-denominator x)))
        (else (bad-real 'truncate x))))

(define (round x)
  (cond ((exact-integer? x) x)
        ((##core#inline "C_i_flonump" x)
         (##core#inline_allocate ("C_a_i_flonum_round_proper" 4) x))
        ((ratnum? x) (let* ((x+1/2 (%+ x (%make-ratnum 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 (%< 'rationalize fx x)) (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 ((%< 'rationalize y x)
	 (find-ratio-between y x))
	((not (%< 'rationalize x y))
	 (list x 1))
	((%> 'rationalize x 0)
	 (sr x y))
	((%< 'rationalize y 0) 
	 (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)
  (let ((result (apply %/ (find-ratio x e))))
    (if (or (inexact? x) (inexact? e))
        (exact->inexact result)
        result)))

(define (exp n)
  (cond ((not (number? n)) (bad-number 'exp n))
        ((cplxnum? n)
         (%* (##core#inline_allocate ("C_a_i_exp" 4) (compnum-real n))
             (let ((p (compnum-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 (%log x)
  (cond
   ((eq? x 0)                        ; Exact zero?  That's undefined
    (log0 'log x))
   ;; avoid calling inexact->exact on X here (to avoid overflow?)
   ((or (cplxnum? x) (%< 'log x 0.0)) ; General case
    (%+ (%log (magnitude x)) (* (make-complex 0 1) (angle x))))
   (else ; Real number case (< already ensured the argument type is a number)
    (##core#inline_allocate ("C_a_i_log" 4) (exact->inexact x)))))

(define (log a #!optional b)
  (if b (%/ (%log a) (%log b)) (%log a)))

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

(define (sin n)
  (cond ((not (number? n)) (bad-number 'sin n))
        ((cplxnum? n) (let ((in (%* %i n)))
                        (%/ (%- (exp in) (exp (%- 0 in))) %i2)))
        (else (##core#inline_allocate ("C_a_i_sin" 4) (exact->inexact n)))))

(define (cos n)
  (cond ((not (number? n)) (bad-number 'cos n))
        ((cplxnum? n) (let ((in (%* %i n)))
                        (%/ (%+ (exp in) (exp (%- 0 in))) 2) ))
        (else (##core#inline_allocate ("C_a_i_cos" 4) (exact->inexact n)))))

(define (tan n)
  (cond ((not (number? n)) (bad-number 'tan n))
        ((cplxnum? n) (%/ (sin n) (cos n)))
        (else (##core#inline_allocate ("C_a_i_tan" 4) (exact->inexact n)))))

;; General case: sin^{-1}(z) = -i\ln(iz + \sqrt{1-z^2})
(define (asin n)
  (cond ((not (number? n)) (bad-number 'asin n))
        ((and (##core#inline "C_i_flonump" n) (fp>= n -1.0) (fp<= n 1.0))
         (##core#inline_allocate ("C_a_i_asin" 4) n))
        ((and (##core#inline "C_fixnump" n) (fx>= n -1) (fx<= n 1))
         (##core#inline_allocate ("C_a_i_asin" 4)
                                 (##core#inline_allocate
                                  ("C_a_i_fix_to_flo" 4) n)))
        ;; General definition can return compnums
        (else (%* %-i (%log (%+ (%* %i n) (%sqrt 'asin (%- 1 (%* n n)))))))))

;; General case:
;; cos^{-1}(z) = 1/2\pi + i\ln(iz + \sqrt{1-z^2}) = 1/2\pi - sin^{-1}(z) = sin(1) - sin(z)
(define acos
  (let ((asin1 (##core#inline_allocate ("C_a_i_asin" 4) 1)))
    (lambda (n)
      (cond ((not (number? n)) (bad-number 'acos n))
            ((and (##core#inline "C_i_flonump" n) (fp>= n -1.0) (fp<= n 1.0))
             (##core#inline_allocate ("C_a_i_acos" 4) n))
            ((and (##core#inline "C_fixnump" n) (fx>= n -1) (fx<= n 1))
             (##core#inline_allocate ("C_a_i_acos" 4)
                                     (##core#inline_allocate
                                      ("C_a_i_fix_to_flo" 4) n)))
            ;; General definition can return compnums
            (else (%- asin1 (asin n)))))))

(define (atan n #!optional b)
  (cond ((not (number? n)) (bad-number 'atan n))
        ((cplxnum? n) (if b
                          (bad-real 'atan n)
                          (let ((in (%* %i n)))
                            (%/ (%- (%log (%+ 1 in)) (%log (%- 1 in))) %i2))))
        (else (if b
                  (##core#inline_allocate
                   ("C_a_i_atan2" 4) (exact->inexact n) (exact->inexact b))
                  (##core#inline_allocate
                   ("C_a_i_atan" 4) (exact->inexact n))))))

;; TODO: This could be moved into C, and do the same trick we do for
;; gcd: allocate one number and keep scaling it down.  After that's
;; done, we may perhaps create C_basic_sqrt() which may call this one.
(define (%exact-integer-sqrt loc k)
  (if (or (eq? 0 k) (eq? 1 k))
      (values k 0)
      ;; Hacker's Delight, figure 11-1 (Newton's method - see also SICP 1.1.7)
      ;; Initial guess is about the smallest x^2 which is bigger than sqrt(k)
      ;; We could subtract one from ilen(k) to avoid overshooting by 1, but
      ;; that's one more operation.
      (let* ((len (fx+ (quotient (integer-length k) 2) 1))
             (g0 (arithmetic-shift 1 len)))
        (let lp ((g0 g0)
                 (g1 (arithmetic-shift
                      (%+ g0 (arithmetic-shift k (fxneg len))) -1)))
          (if (%< loc g1 g0)
              (lp g1 (arithmetic-shift (%+ g1 (quotient k g1)) -1))
              (values g0 (%- k (%* g0 g0))))))))

(define (exact-integer-sqrt x)
  (if (or (not (exact-integer? x)) (##core#inline "C_u_i_2_integer_lessp" x 0))
      (bad-natural 'exact-integer-sqrt x)
      (%exact-integer-sqrt 'exact-integer-sqrt x)))

;; This procedure is so large because it tries very hard to compute
;; exact results if at all possible.
(define (%sqrt loc n)
  (cond ((cplxnum? n)     ; Must be checked before we call "negative?"
         (let ((p (%/ (angle n) 2))
               (m (##core#inline_allocate ("C_a_i_sqrt" 4) (magnitude n))) )
           (make-complex (%* m (cos p)) (%* m (sin p)) ) ))
        ((negative? n)
         (make-complex
          0.0
          (##core#inline_allocate
           ("C_a_i_sqrt" 4)
           (exact->inexact ((##core#primitive "C_basic_negate") n)))))
        ((exact-integer? n)
         (receive (s^2 r)
             (%exact-integer-sqrt loc n)
           (if (eq? 0 r)
               s^2
               (##core#inline_allocate ("C_a_i_sqrt" 4) (exact->inexact n)))))
        ((ratnum? n)                    ; Try to compute exact sqrt
         (receive (ns^2 nr)             ; We already know n is positive!
             (%exact-integer-sqrt loc (ratnum-numerator n))
           (if (eq? nr 0)
               (receive (ds^2 dr)
                   (%exact-integer-sqrt loc (ratnum-denominator n))
                 (if (eq? dr 0)
                     (%/ ns^2 ds^2)
                     (%sqrt loc (exact->inexact n))))
               (%sqrt loc (exact->inexact n)))))
        (else (##core#inline_allocate ("C_a_i_sqrt" 4) (exact->inexact n)))))

(define (sqrt x) (%sqrt 'sqrt x))

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

;; Generalized Newton's algorithm for positive integers, with a little help
;; from Wikipedia ;)  https://en.wikipedia.org/wiki/Nth_root_algorithm
;; TODO: This can also be moved to C, if exact-integer-sqrt is successful?
(define (%exact-integer-nth-root loc k n)
  (cond
   ((or (eq? 0 k) (eq? 1 k) (eq? 1 n)) ; Maybe call exact-integer-sqrt on n=2?
    (values k 0))
   ((%< loc n 1)
    (bad-positive loc n))
   (else
    (let ((len (integer-length k)))
      (if (%< loc len n)      ; Idea from Gambit: 2^{len-1} <= k < 2^{len}
          (values 1 (%- k 1)) ; Since we know x >= 2, we know x^{n} can't exist
          ;; Set initial guess to (at least) 2^ceil(ceil(log2(k))/n)
          (let* ((shift-amount (inexact->exact (ceiling (/ (fx+ len 1) n))))
                 (g0 (arithmetic-shift 1 shift-amount))
                 (n-1 (%- n 1)))
            (let lp ((g0 g0)
                     (g1 (%quotient loc (%+ (%* n-1 g0) (%quotient loc k (%integer-power g0 n-1))) n)))
              (if (%< loc g1 g0)
                  (lp g1 (%quotient loc (%+ (%* n-1 g1) (%quotient loc k (%integer-power g1 n-1))) n))
                  (values g0 (%- k (%integer-power g0 n)))))))))))

(define (exact-integer-nth-root k n)
  (cond ((or (not (exact-integer? k))
             (##core#inline "C_u_i_2_integer_lessp" k 0))
         (bad-natural 'exact-integer-nth-root k))
        ((not (exact-integer? n))
         (bad-integer 'exact-integer-nth-root n))
        (else (%exact-integer-nth-root 'exact-integer-nth-root k n))))

;; TODO: Why does this have to be basic_negate and not C_u_integer_negate?
(define (%integer-power base e)
  (if (negative? e)
      (%/ 1 (%integer-power base ((##core#primitive "C_basic_negate") e)))
      (let lp ((res 1) (e2 e))
        (cond
         ((eq? e2 0) res)
         ((even? e2) ; recursion is faster than iteration here
          (%* res (square (lp 1 (arithmetic-shift (inexact->exact e2) -1)))))
         (else
          (lp (%* res base) (%- e2 1)))))))

(define (expt a b)
  (define (log-expt a b)
    (exp (%* b (%log a))))
  (define (slow-expt a b)
    (if (eq? 0 a)
        (expt0 'expt a b)
        (exp (%* b (%log a)))))
  (cond ((not (number? a)) (bad-number 'expt a))
        ((not (number? b)) (bad-number 'expt b))
        ((and (ratnum? a) (not (inexact? b)))
         ;; (n*d)^b = n^b * d^b = n^b * x^{-b}  | x = 1/b
         ;; Hopefully faster than integer-power
         (%* (expt (ratnum-numerator a) b)
             (expt (ratnum-denominator a)
                   ((##core#primitive "C_basic_negate") b))))
        ((ratnum? b)
         ;; x^{a/b} = (x^{1/b})^a
         (cond
          ((exact-integer? a)
           (if (negative? a)
               (log-expt (exact->inexact a) (exact->inexact b))
               (receive (ds^n r)
                   (%exact-integer-nth-root 'expt a (ratnum-denominator b))
                 (if (eq? r 0)
                     (expt ds^n (ratnum-numerator b))
                     ((##core#primitive "C_expt")
                      (exact->inexact a) (exact->inexact b))))))
          ((##core#inline "C_i_flonump" a)
           (log-expt a (exact->inexact b)))
          (else (slow-expt a b))))
        ((or (cplxnum? b) (and (cplxnum? a) (not (integer? b))))
         (slow-expt a b))
        ((and (##core#inline "C_i_flonump" b)
              (not (##core#inline "C_u_i_fpintegerp" b)))
         (if (negative? a)
             (log-expt (exact->inexact a) (exact->inexact b))
             ((##core#primitive "C_expt") (exact->inexact a) (exact->inexact b))))
        ((##core#inline "C_i_flonump" a)
         ((##core#primitive "C_expt") (exact->inexact a) (exact->inexact b)))
        ;; this doesn't work that well, yet...
        ;; (XXX: What does this mean? why not? I do know this is ugly... :P)
        (else (if (or (inexact? a) (inexact? b))
                  (exact->inexact (%integer-power a b))
                  (%integer-power a b)))) )

(define (conj x)                        ; Nonstandard, not part of CHICKEN
  (cond ((##core#inline "C_i_basic_numberp" x) x)
        ((ratnum? x) x)
        ((cplxnum? x)
         (make-complex (compnum-real x)
                       ((##core#primitive "C_basic_negate") (compnum-imag x))))
        (else (bad-number 'conj x)) ) )

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

(define (@extended-signum x)
  (cond ((ratnum? x) (##core#inline "C_u_i_integer_signum" (ratnum-numerator x)))
        ((cplxnum? x) (make-polar 1 (angle x)))
        (else (bad-number 'signum x))))

(define (signum x) ((##core#primitive "C_basic_signum") x))

;; From SRFI-60 (oddly not in Chicken core) TODO: Maybe get rid of
;; ->integer, move this to C and error out if we get passed anything
;; non-integerlike?
(define (integer-length x)
  (unless (exact-integer? x) (bad-exact 'integer-length x))
  (##core#inline "C_u_i_integer_length" x))

(define (bitwise-and . xs)
  (let loop ((x -1) (xs xs))
    (if (null? xs)
	x
	(let ((xi (##sys#slot xs 0)))
          (unless (exact-integer? xi) (bad-exact 'bitwise-and xi))
	  (loop
	   ((##core#primitive "C_u_2_integer_bitwise_and") x xi)
	   (##sys#slot xs 1) ) ) ) ) )

(define (bitwise-ior . xs)
  (let loop ((x 0) (xs xs))
    (if (null? xs)
	x
	(let ((xi (##sys#slot xs 0)))
          (unless (exact-integer? xi) (bad-exact 'bitwise-ior xi))
	  (loop
	   ((##core#primitive "C_u_2_integer_bitwise_ior") x xi)
	   (##sys#slot xs 1) ) ) ) ) )

(define (bitwise-xor . xs)
  (let loop ((x 0) (xs xs))
    (if (null? xs)
	x
	(let ((xi (##sys#slot xs 0)))
          (unless (exact-integer? xi) (bad-exact 'bitwise-xor xi))
	  (loop
	   ((##core#primitive "C_u_2_integer_bitwise_xor") x xi)
	   (##sys#slot xs 1) ) ) ) ) )

(define (bitwise-not n)
  (unless (exact-integer? n) (bad-exact 'bitwise-not n))
  ((##core#primitive "C_u_2_integer_minus") -1 n) )

(define (arithmetic-shift n m)
  (unless (exact-integer? n) (bad-exact 'arithmetic-shift n))
  ;; Strictly speaking, shifting *right* is okay for any number
  ;; (ie, shifting by a negative bignum would just result in 0 or -1)...
  (unless (##core#inline "C_fixnump" m)
    (##sys#signal-hook #:type-error 'arithmetic-shift
                       "can only shift by fixnum amounts" n m))
  ((##core#primitive "C_u_integer_shift") n m))

(define @extended-number->string
  (let ((string-append string-append))
    (lambda (n base)
      (cond
       ((ratnum? n) (string-append (number->string (ratnum-numerator n) base)
                                   "/"
                                   (number->string (ratnum-denominator n) base)))
       ;; What about bases that include an "i"?  That could lead to
       ;; ambiguous results.
       ((cplxnum? n) (let ((r (compnum-real n))
                           (i (compnum-imag n)) )
                       (string-append
                        (number->string r base)
                        ;; The infinities and NaN always print their sign
                        (if (and (finite? i) (%> 'number->string i 0)) "+" "")
                        (number->string i base) "i") ))
       (else (bad-number 'number->string n)))  ) ) )

(define number->string (##core#primitive "C_basic_number_to_string"))
(define ##sys#number->string number->string) ; for printer

;; We try to prevent memory exhaustion attacks by limiting the
;; maximum exponent value.
;; TODO: Make this a parameter?  Would probably slow things down even more...
(define-constant +maximum-allowed-exponent+ 10000)

;; Shorthand
(define-inline (%subchar s i) (##core#inline "C_subchar" s i))
(define (%string->compnum radix str offset exactness)
  (define (go-inexact!)
    ;; Go inexact unless exact was requested (with #e prefix)
    (unless (eq? exactness 'e) (set! exactness 'i)))
  (define (safe-exponent value e)
    (and e (if (not value)
               0
               (cond
                ((%> 'string->number e +maximum-allowed-exponent+)
                 (and (eq? exactness 'i)
                      (cond ((zero? value) 0.0)
                            ((%> 'string->number value 0.0) +inf.0)
                            (else -inf.0))))
                ((%< 'string->number e ((##core#primitive "C_u_integer_negate")
                                        +maximum-allowed-exponent+))
                 (and (eq? exactness 'i) +0.0))
                (else (%* value (expt 10 e)))))))
  (let* ((len (##sys#size str))
         (r..9 (integer->char (fx+ (char->integer #\0) radix)))
         (r..a (integer->char (fx+ (char->integer #\a) (fx- radix 10))))
         (r..A (integer->char (fx+ (char->integer #\A) (fx- radix 10))))
         ;; Ugly flag which we need (note that "exactness" is mutated too!)
         ;; Since there is (almost) no backtracking we can do this.
         (seen-hashes? #f)
         ;; All these procedures return #f or an object consed onto an end
         ;; position.  If the cdr is false, that's the end of the string.
         ;; If just #f is returned, the string contains invalid number syntax.
         (scan-digits
          (lambda (start)
            (let lp ((i start))
              (if (fx= i len)
                  (and (fx> i start) (cons i #f))
                  (let ((c (%subchar str i)))
                    (if (fx<= radix 10)
                        (if (and (char>=? c #\0) (char<=? c r..9))
                            (lp (fx+ i 1))
                            (and (fx> i start) (cons i i)))
                        (if (or (and (char>=? c #\0) (char<=? c #\9))
                                (and (char>=? c #\a) (char<=? c r..a))
                                (and (char>=? c #\A) (char<=? c r..A)))
                            (lp (fx+ i 1))
                            (and (fx> i start) (cons i i)))))))))
         (scan-hashes
          (lambda (start)
            (let lp ((i start))
              (if (fx= i len)
                  (and (fx> i start) (cons i #f))
                  (let ((c (%subchar str i)))
                    (if (eq? c #\#)
                        (lp (fx+ i 1))
                        (and (fx> i start) (cons i i))))))))
         (scan-digits+hashes
          (lambda (start neg? all-hashes-ok?)
            (let* ((digits (and (not seen-hashes?) (scan-digits start)))
                   (hashes (if digits
                               (and (cdr digits) (scan-hashes (cdr digits)))
                               (and all-hashes-ok? (scan-hashes start))))
                   (end (or hashes digits)))
              (and-let* ((end)
                         (num ((##core#primitive "C_digits_to_integer")
                               str start (car end) radix neg?)))
                (when hashes            ; Eeewww. Feeling dirty yet?
                  (set! seen-hashes? #t)
                  (go-inexact!))
                (cons num (cdr end))))))
         (scan-exponent
          (lambda (start)
            (and (fx< start len)
                 (let ((sign (case (%subchar str start)
                               ((#\+) 'pos) ((#\-) 'neg) (else #f))))
                   (and-let* ((start (if sign (fx+ start 1) start))
                              (end (scan-digits start)))
                     (go-inexact!)
                     (cons ((##core#primitive "C_digits_to_integer")
                            str start (car end) radix (eq? sign 'neg))
                           (cdr end)))))))
         (scan-decimal-tail
          (lambda (start neg? decimal-head)
            (and (fx< start len)
                 (let* ((tail (scan-digits+hashes start neg? decimal-head))
                        (next (if tail (cdr tail) start)))
                   (and (or decimal-head (not next)
                            (fx> next start)) ; Don't allow empty "."
                        (case (and next (%subchar str next))
                          ((#\e #\s #\f #\d #\l
                            #\E #\S #\F #\D #\L)
                           (and-let* (((fx> len next))
                                      (ee (scan-exponent (fx+ next 1)))
                                      (e (car ee))
                                      (h (safe-exponent decimal-head e)))
                             (let* ((te (and tail (fx- e (fx- (cdr tail) start))))
                                    (num (and tail (car tail)))
                                    (t (safe-exponent num te)))
                               (cons (if t (%+ h t) h) (cdr ee)))))
                          (else (let* ((last (or next len))
                                       (te (and tail (fx- start last)))
                                       (num (and tail (car tail)))
                                       (t (safe-exponent num te))
                                       (h (or decimal-head 0)))
                                  (cons (if t (%+ h t) h) next)))))))))
         (scan-ureal
          (lambda (start neg?)
            (if (and (fx> len (fx+ start 1)) (eq? radix 10)
                     (eq? (%subchar str start) #\.))
                (begin
                  (go-inexact!)
                  (scan-decimal-tail (fx+ start 1) neg? #f))
                (and-let* ((end (scan-digits+hashes start neg? #f)))
                  (case (and (cdr end) (%subchar str (cdr end)))
                    ((#\.)
                     (go-inexact!)
                     (and (eq? radix 10)
                          (if (fx> len (fx+ (cdr end) 1))
                              (scan-decimal-tail (fx+ (cdr end) 1) neg? (car end))
                              (cons (car end) #f))))
                    ((#\e #\s #\f #\d #\l
                      #\E #\S #\F #\D #\L)
                     (and-let* (((eq? radix 10))
                                ((fx> len (cdr end)))
                                (ee (scan-exponent (fx+ (cdr end) 1)))
                                (num (car end))
                                (val (safe-exponent num (car ee))))
                       (cons val (cdr ee))))
                    ((#\/)
                     (set! seen-hashes? #f) ; Reset flag for denominator
                     (and-let* (((fx> len (cdr end)))
                                (d (scan-digits+hashes (fx+ (cdr end) 1) #f #f))
                                (num (car end))
                                (denom (car d)))
                       (if (not (eq? denom 0))
                           (cons (%/ num denom) (cdr d))
                           ;; Hacky: keep around an inexact until we decide we
                           ;; *really* need exact values, then fail at the end.
                           (and (not (eq? exactness 'e))
                                (case (signum num)
                                  ((-1) (cons -inf.0 (cdr d)))
                                  ((0)  (cons +nan.0 (cdr d)))
                                  ((+1) (cons +inf.0 (cdr d))))))))
                    (else end))))))
         (scan-real
          (lambda (start)
            (and (fx< start len)
                 (let* ((sign (case (%subchar str start)
                                ((#\+) 'pos) ((#\-) 'neg) (else #f)))
                        (next (if sign (fx+ start 1) start)))
                   (and (fx< next len)
                        (case (%subchar str next)
                          ((#\i #\I)
                           (or (and sign
                                    (cond
                                     ((fx= (fx+ next 1) len) ; [+-]i
                                      (cons (if (eq? sign 'neg) -1 1) next))
                                     ((and (fx<= (fx+ next 5) len)
                                           (string-ci=? (substring str next (fx+ next 5)) "inf.0"))
                                      (go-inexact!)
                                      ;; XXX TODO: Make this +inf/-inf
                                      (cons (fp/ (if (eq? sign 'neg) -1.0 1.0) 0.0)
                                            (and (fx< (fx+ next 5) len)
                                                 (fx+ next 5))))
                                     (else #f)))
                               (scan-ureal next (eq? sign 'neg))))
                          ((#\n #\N)
                           (or (and sign
                                    (fx<= (fx+ next 5) len)
                                    (string-ci=? (substring str next (fx+ next 5)) "nan.0")
                                    (begin (go-inexact!)
                                           ;; XXX TODO: Make this +nan.0
                                           (cons (fp/ 0.0 0.0)
                                                 (and (fx< (fx+ next 5) len)
                                                      (fx+ next 5)))))
                               (scan-ureal next (eq? sign 'neg))))
                          (else (scan-ureal next (eq? sign 'neg)))))))))
         (number (and-let* ((r1 (scan-real offset)))
                   (case (and (cdr r1) (%subchar str (cdr r1)))
                     ((#f) (car r1))
                     ((#\i #\I) (and (fx= len (fx+ (cdr r1) 1))
                                     (or (eq? (%subchar str offset) #\+) ; ugh
                                         (eq? (%subchar str offset) #\-))
                                     (make-rectangular 0 (car r1))))
                     ((#\+ #\-)
                      (set! seen-hashes? #f) ; Reset flag for imaginary part
                      (and-let* ((r2 (scan-real (cdr r1)))
                                 ((cdr r2))
                                 ((fx= len (fx+ (cdr r2) 1)))
                                 ((or (eq? (%subchar str (cdr r2)) #\i)
                                      (eq? (%subchar str (cdr r2)) #\I))))
                        (make-rectangular (car r1) (car r2))))
                     ((#\@)
                      (set! seen-hashes? #f) ; Reset flag for angle
                      (and-let* ((r2 (scan-real (fx+ (cdr r1) 1)))
                                 ((not (cdr r2))))
                        (make-polar (car r1) (car r2))))
                     (else #f)))))
    (and number (if (eq? exactness 'i)
                    (exact->inexact number)
                    ;; Ensure we didn't encounter +inf.0 or +nan.0 with #e
                    (and (finite? number) number)))))

(define (string->number str #!optional (base 10))
  (##sys#check-string str 'string->number)
  (unless (and (##core#inline "C_fixnump" base)
               (fx< 1 base) (fx< base 37)) ; We only have 0-9 and the alphabet!
    (bad-base 'string->number base))
  (let scan-prefix ((i 0)
                    (exness #f)
                    (radix #f)
                    (len (##sys#size str)))
    (if (and (fx< (fx+ i 2) len) (eq? (%subchar str i) #\#))
        (case (%subchar str (fx+ i 1))
          ((#\i #\I) (and (not exness) (scan-prefix (fx+ i 2) 'i radix len)))
          ((#\e #\E) (and (not exness) (scan-prefix (fx+ i 2) 'e radix len)))
          ((#\b #\B) (and (not radix) (scan-prefix (fx+ i 2) exness 2 len)))
          ((#\o #\O) (and (not radix) (scan-prefix (fx+ i 2) exness 8 len)))
          ((#\d #\D) (and (not radix) (scan-prefix (fx+ i 2) exness 10 len)))
          ((#\x #\X) (and (not radix) (scan-prefix (fx+ i 2) exness 16 len)))
          (else #f))
        (%string->compnum (or radix base) str i exness))))

;;; Reader hook
(define (##sys#string->number str #!optional (radix 10) exactness)
  (%string->compnum radix str 0 exactness))

(define (randomize . n)
  (let ((nn (if (null? n)
		(##sys#flo2fix (fp/ (current-seconds) 1000.0)) ; wall clock time
		(car n))))
    (unless (exact-integer? nn) (bad-exact 'randomize nn))
    (##core#inline "C_u_i_integer_randomize" nn) ) )

(define (random n)
  (unless (exact-integer? n) (bad-exact 'random n))
  (if (eq? n 0)
      0
      ((##core#primitive "C_u_integer_random") n) ) )


;;; Non-standard type procedures

(define (basic-number? x) (##core#inline "C_i_basic_numberp" x))

(define (extended-number? x) ; This does _not_ "include" basics; see "number?"
  (and (##core#inline "C_blockp" x)
       (##sys#generic-structure? x)
       (or (eq? (##sys#slot x 0) 'ratnum)
           (eq? (##sys#slot x 0) 'compnum))))

(define (bignum? x) (##core#inline "C_i_bignump" x))

(define (nan? x) (##core#inline "C_i_nanp" x))
(define (infinite? x) (##core#inline "C_i_numbers_infinitep" x))
(define (finite? x) (##core#inline "C_i_numbers_finitep" x))

(define (ratnum? x) (##sys#structure? x 'ratnum)) ; rational number

;; XXX THE ONES BELOW ARE EXTREMELY CONFUSING
;; Especially considering the type tag in a complex number is "compnum"!
;; Best to rename cplxnum? to compnum? and ditch the rest.
;; A user can easily define them themselves
(define (cplxnum? x) (##sys#structure? x 'compnum)) ; complex number

(define (rectnum? x)    ; "exact" complex number
  (and (cplxnum? x)
       (integer? (compnum-real x))
       (integer? (compnum-imag x))))

(define (compnum? x)    ; inexact complex number
  (and (cplxnum? x)
       (inexact? (compnum-real x))
       (inexact? (compnum-imag x))))

(define (cintnum? x)    ; integer number
  (or (integer? x) (rectnum? x)) )

(define (cflonum? x)    ; floatingpoint number
  (or (##core#inline "C_i_flonump" x) (compnum? x)) )

;;; What we provide

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

)