;;; Copyright (C) William D Clinger (2017).
;;;
;;; Permission is hereby granted, free of charge, to any person
;;; obtaining a copy of this software and associated documentation
;;; files (the "Software"), to deal in the Software without
;;; restriction, including without limitation the rights to use,
;;; copy, modify, merge, publish, distribute, sublicense, and/or
;;; sell copies of the Software, and to permit persons to whom the
;;; Software is furnished to do so, subject to the following
;;; conditions:
;;;
;;; The above copyright notice and this permission notice shall be
;;; included in all copies or substantial portions of the Software.
;;;
;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;; OTHER DEALINGS IN THE SOFTWARE.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Tests these (srfi 144) constants and procedures:
;;;
;;;     fl-e
;;;     fl-1/e
;;;     fl-e-2
;;;     fl-e-pi/4
;;;     fl-log2-e
;;;     fl-log10-e
;;;     fl-log-2
;;;     fl-1/log-2
;;;     fl-log-3
;;;     fl-log-pi
;;;     fl-log-10
;;;     fl-1/log-10
;;;     fl-pi
;;;     fl-1/pi
;;;     fl-2pi
;;;     fl-pi/2
;;;     fl-pi/4
;;;     fl-2/sqrt-pi
;;;     fl-pi-squared
;;;     fl-degree
;;;     fl-2/pi
;;;     fl-sqrt-2
;;;     fl-sqrt-3
;;;     fl-sqrt-5
;;;     fl-sqrt-10
;;;     fl-1/sqrt-2
;;;     fl-cbrt-2
;;;     fl-cbrt-3
;;;     fl-4thrt-2
;;;     fl-phi
;;;     fl-log-phi
;;;     fl-1/log-phi
;;;     fl-euler
;;;     fl-e-euler
;;;     fl-sin-1
;;;     fl-cos-1
;;;     fl-gamma-1/2
;;;     fl-gamma-1/3
;;;     fl-gamma-2/3
;;;
;;;     fl-greatest
;;;     fl-least
;;;     fl-epsilon
;;;     fl-fast-fl+*
;;;     fl-integer-exponent-zero
;;;     fl-integer-exponent-nan
;;;
;;;     flonum
;;;     fladjacent
;;;     flcopysign
;;;     make-flonum
;;;
;;;     flinteger-fraction
;;;     flexponent
;;;     flinteger-exponent
;;;     flnormalized-fraction-exponent
;;;     flsign-bit
;;;
;;;     flonum?
;;;     fl=?
;;;     fl<?
;;;     fl>?
;;;     fl<=?
;;;     fl>=?
;;;     flunordered?
;;;     flmax
;;;     flmin
;;;     flinteger?
;;;     flzero?
;;;     flpositive?
;;;     flnegative?
;;;     flodd?
;;;     fleven?
;;;     flfinite?
;;;     flinfinite?
;;;     flnan?
;;;     flnormalized?
;;;     fldenormalized?
;;;
;;;     fl+
;;;     fl*
;;;     fl+*
;;;     fl-
;;;     fl/
;;;     flabs
;;;     flabsdiff
;;;     flposdiff
;;;     flsgn
;;;     flnumerator
;;;     fldenominator
;;;     flfloor
;;;     flceiling
;;;     flround
;;;     fltruncate
;;;
;;;     flexp
;;;     flexp2
;;;     flexp-1
;;;     flsquare
;;;     flsqrt
;;;     flcbrt
;;;     flhypot
;;;     flexpt
;;;     fllog
;;;     fllog1+
;;;     fllog2
;;;     fllog10
;;;     make-fllog-base
;;;
;;;     flsin
;;;     flcos
;;;     fltan
;;;     flasin
;;;     flacos
;;;     flatan
;;;     flsinh
;;;     flcosh
;;;     fltanh
;;;     flasinh
;;;     flacosh
;;;     flatanh
;;;
;;;     flquotient
;;;     flremainder
;;;     flremquo
;;;
;;;     flgamma
;;;     flloggamma
;;;     flfirst-bessel
;;;     flsecond-bessel
;;;     flerf
;;;     flerfc


(import r7rs
        (scheme base)
        (srfi 144)
        test
        (scheme inexact)
        (only (chicken string) ->string)
        (only (srfi 1) filter iota))


(define nnan? (conjoin number? nan?))

(define old-comparator (current-test-comparator))
(current-test-comparator
 (lambda (a b)
   (or (and (nnan? a) (nnan? b) #t)
       (old-comparator
        (if (list? a)
            (map (lambda (x) (if (nnan? x) 'NAN x)) a)
            a)
        (if (list? b)
            (map (lambda (x) (if (nnan? x) 'NAN x)) b)
            b)))))

;; convenient values for test cases

(define posints (map flonum '(1 2 3 4 5 10 65536 1e23)))
(define nats (cons (flonum 0) posints))
(define ints (append (map flonum '(-20 -8 -2 -1)) nats))
(define posfracs (map flonum '(1/1000 1/10 1/3 1/2)))
(define extremes
  (list (fl- fl-greatest) (fl- fl-least) fl-least fl-greatest))
(define infinities (map flonum (list -inf.0 +inf.0)))
(define weird (append infinities (list (flonum +nan.0))))

(define somereals (append (map flonum
                               (list (fl- fl-greatest)
                                     -10
                                     (fl- fl-least)
                                     0))
                          posfracs
                          posints))
(define somereals+weird
  (append somereals weird))

(define negzero (flonum -0.0))
(define zero (flonum 0))
(define one (flonum 1))
(define two (flonum 2))

(define neginf (flonum -inf.0))
(define posinf (flonum +inf.0))
(define nan (flonum +nan.0))

(test-begin "srfi-144")
;; Mathematical constants

(test 2.718281828459045235360287                    fl-e)
(test .3678794411714423215955238                    fl-1/e)
(test 7.389056098930650227230427                    fl-e-2)
(test 2.1932800507380154566                         fl-e-pi/4)
(test 1.4426950408889634073599246810018921374266    fl-log2-e)
(test .4342944819032518276511289                    fl-log10-e)
(test .6931471805599453094172321                    fl-log-2)
(test 1.4426950408889634073599246810018921374266    fl-1/log-2)
(test 1.0986122886681096913952452                   fl-log-3)
(test 1.144729885849400174143427                    fl-log-pi)
(test 2.3025850929940456840179915                   fl-log-10)
(test 0.4342944819032518276511289189166050822944    fl-1/log-10)

(test 3.1415926535897932384626433832795028841972    fl-pi)
(test 0.3183098861837906715377675267450287240689    fl-1/pi)
(test 6.283185307179586476925287                    fl-2pi)
(test 1.570796326794896619231322                    fl-pi/2)
(test .7853981633974483096156608                    fl-pi/4)
(test .5641895835477562869480795                    (/ fl-2/sqrt-pi 2))
(test 9.869604401089358618834491                    fl-pi-squared)
(test 0.0174532925199432957692369076848861271344    fl-degree)
(test .3183098861837906715377675                    (/ fl-2/pi 2))

(test 1.4142135623730950488016887242096980785697    fl-sqrt-2)
(test 1.7320508075688772935274463415058723669428    fl-sqrt-3)
(test 2.2360679774997896964091736687311762354406    fl-sqrt-5)
(test 3.1622776601683793319988935444327185337196    fl-sqrt-10)
(test 1.4142135623730950488016887242096980785697    (* 2 fl-1/sqrt-2))
(test 1.2599210498948731647672106072782283505703    fl-cbrt-2)
(test 1.4422495703074083823216383107801095883919    fl-cbrt-3)
(test 1.1892071150027210667174999705604759152930    fl-4thrt-2)

(test 1.6180339887498948482045868343656381177203    fl-phi)
(test 0.4812118250596034474977589134243684231352    fl-log-phi)
(test 2.0780869212350275376013226061177957677422    fl-1/log-phi)
(test 0.5772156649015328606065120900824024310422    fl-euler)
(test 1.7810724179901979852365041031071795491696    fl-e-euler)

(test 0.8414709848078965066525023216302989996226    fl-sin-1)
(test 0.5403023058681397174009366074420766037323    fl-cos-1)

(test 1.7724538509055160272981674833411451827975    fl-gamma-1/2)
(test 2.6789385347077476336556929409746776441287    fl-gamma-1/3)
(test 1.3541179394264004169452880281545137855193    fl-gamma-2/3)

;; Implementation Constants

(test-assert (inexact? fl-greatest))
(test-assert (inexact? fl-least))
(test-assert (inexact? fl-epsilon))

(test-assert (real? fl-greatest))
(test-assert (real? fl-least))
(test-assert (real? fl-epsilon))

(test-assert (flonum? fl-greatest))
(test-assert (flonum? fl-least))
(test-assert (flonum? fl-epsilon))

(test-assert (< 0.0
                fl-least
                fl-epsilon
                1.0
                (+ 1.0 fl-epsilon)
                fl-greatest
                posinf))
(test-assert (= (* 2 fl-greatest) posinf))
(test-assert (= 1 (/ (+ 1 (+ 1.0 fl-epsilon)) 2)))
(test-assert (= 0.0 (/ fl-least 2)))

(test-assert (boolean? fl-fast-fl+*))
(test-assert (exact-integer? fl-integer-exponent-zero))
(test-assert (exact-integer? fl-integer-exponent-nan))

;; Constructors

(test (flonum 3) (flonum 3.0))
(test somereals (map flonum somereals))
(test weird (map flonum weird))

(test somereals (map fladjacent somereals somereals))
(test weird (map fladjacent weird weird))

(test fl-least (fladjacent zero posinf))
(test (fl- fl-least) (fladjacent zero neginf))
(test (fl+ fl-least fl-least) (fladjacent fl-least posinf))
(test zero (fladjacent fl-least neginf))
(test negzero (fladjacent (fl- fl-least) posinf))
(test (fl* -2.0 fl-least) (fladjacent (fl- fl-least) neginf))

(test fl-least (fladjacent zero one))
(test (fl- fl-least) (fladjacent zero (fl- one)))
(test (fl+ fl-least fl-least) (fladjacent fl-least one))
(test zero (fladjacent fl-least (fl- one)))
(test negzero (fladjacent (fl- fl-least) one))
(test (fl* -2.0 fl-least) (fladjacent (fl- fl-least) (fl- one)))

(test fl-epsilon (fl- (fladjacent one fl-greatest) one))
(test (fl/ fl-epsilon 2.0) (fl- one (fladjacent one zero)))

(test fl-greatest (fladjacent posinf zero))
(test (fl- fl-greatest) (fladjacent neginf zero))

(test zero (flcopysign zero posinf))
(test negzero (flcopysign zero neginf))
(test zero (flcopysign zero one))
(test negzero (flcopysign zero (fl- one)))
(test one (flcopysign one fl-least))
(test (fl- one) (flcopysign one (fl- fl-greatest)))
(test one (flcopysign (fl- one) zero))
(test somereals (map flcopysign somereals somereals))
(test (map fl- somereals)
      (map flcopysign somereals (map fl- somereals)))
(test infinities (map flcopysign infinities infinities))
(test (reverse infinities)
      (map flcopysign infinities (reverse infinities)))

(test zero (make-flonum zero 12))
(test zero (make-flonum zero -24))
(test zero (make-flonum zero 0))
(test somereals
      (map make-flonum somereals (map (lambda (x) 0) somereals)))
(test (map (lambda (x) (fl* (flonum 4) x)) somereals)
      (map make-flonum somereals (map (lambda (x) 2) somereals)))
(test (map (lambda (x) (fl/ x (flonum 16))) somereals)
      (map make-flonum somereals (map (lambda (x) -4) somereals)))
(test posinf (make-flonum fl-greatest 1))
(test neginf (make-flonum (fl- fl-greatest) 1))
(test (fl/ fl-greatest two) (make-flonum fl-greatest -1))
(test (fl- (fl/ fl-greatest two)) (make-flonum (fl- fl-greatest) -1))
(test (fl* two fl-least) (make-flonum fl-least 1))
(test (fl- (fl* two fl-least)) (make-flonum (fl- fl-least) 1))
(test zero (make-flonum fl-least -1))
(test negzero (make-flonum (fl- fl-least) -1))

;; Accessors

(call-with-values
    (lambda () (flinteger-fraction 3.75))
  (lambda (q r)
    (test q (flonum 3))
    (test r (flonum .75))))

(call-with-values
    (lambda () (flinteger-fraction -3.75))
  (lambda (q r)
    (test q (flonum -3))
    (test r (flonum -.75))))

(test (flexponent (flexpt two (flonum 12)))
      (flonum 12.0))
(test (flonum 12)
      (flexponent (flexpt two (flonum 12.5))))
(test (flexponent (flexpt two (flonum -5)))
      (flonum -5.0))
(test (flonum +4)
      (flexponent (flexpt two (flonum +4.5))))
(test (flonum -5)
      (flexponent (flexpt two (flonum -4.5))))

(test 12 (flinteger-exponent (flexpt two (flonum 12))))
(test 12 (flinteger-exponent (flexpt two (flonum 12.5))))
(test -5 (flinteger-exponent (flexpt two (flonum -5))))
(test -5 (flinteger-exponent (flexpt two (flonum -4.5))))

(let* ((correct?
        (lambda (x y n)
          (or (fl=? x (* y (expt two n)))
              (fl=? x (* 4.00 y (expt two (- n 2))))
              (fl=? x (* 0.25 y (expt two (+ n 2)))))))
       (test-flnormalized-fraction-exponent
        (lambda (x)
          (call-with-values
              (lambda () (flnormalized-fraction-exponent x))
            (lambda (y n)
              (list (flonum? y)
                    (exact-integer? n)
                    (fl<=? (flonum 0.5) (flabs y))
                    (fl<? (flabs y) one)
                    (correct? x y n)))))))
  (test '(#t #t #f #t #t)
        (test-flnormalized-fraction-exponent zero))
  (test '(#t #t #f #t #t)
        (test-flnormalized-fraction-exponent negzero))
  (test '(#t #t #t #t #t)
        (test-flnormalized-fraction-exponent one))
  (test '(#t #t #t #t #t)
        (test-flnormalized-fraction-exponent two))
  (test '(#t #t #t #t #t)
        (test-flnormalized-fraction-exponent fl-least))
  (test '(#t #t #t #t #t)
        (test-flnormalized-fraction-exponent fl-greatest))
  (test '(#t #t #t #t #t)
        (test-flnormalized-fraction-exponent (fl- fl-least)))
  (test '(#t #t #t #t #t)
        (test-flnormalized-fraction-exponent (fl- fl-greatest)))
  (test '(#t #t #t #t #t)
        (test-flnormalized-fraction-exponent posinf))
  (test '(#t #t #t #t #t)
        (test-flnormalized-fraction-exponent neginf))
  (test '(#t #t #f #f #f)
        (test-flnormalized-fraction-exponent nan)))

(test 0 (flsign-bit one))
(test 0 (flsign-bit zero))
(test 1 (flsign-bit negzero))
(test 1 (flsign-bit (flonum -2)))
(test 0 (flsign-bit posinf))
(test 1 (flsign-bit neginf))

;; Predicates

(let ((alltrue  (map (lambda (x) #t) somereals))
      (allfalse (map (lambda (x) #f) somereals)))

  (test alltrue (map flonum? somereals))
  (test '(#t #t #t) (map flonum? weird))

  (test-assert (not (fl=? zero fl-least)))
  (test-assert (fl=? fl-least fl-least))
  (test-assert (not (fl=? one fl-least)))
  (test-assert (fl=? neginf neginf))
  (test-assert (not (fl=? neginf posinf)))
  (test-assert (not (fl=? posinf neginf)))
  (test-assert (fl=? posinf posinf))
  (test-assert (not (fl=? zero nan)))
  (test-assert (not (fl=? nan one)))
  (test alltrue (map fl=? somereals somereals))
  (test (cdr allfalse) (map fl=? somereals (cdr somereals)))
  (test (cdr allfalse) (map fl=? (cdr somereals) somereals))

  (test-assert (fl<? zero fl-least))
  (test-assert (not (fl<? fl-least fl-least)))
  (test-assert (not (fl<? one fl-least)))
  (test-assert (not (fl<? neginf neginf)))
  (test-assert (fl<? neginf posinf))
  (test-assert (not (fl<? posinf neginf)))
  (test-assert (not (fl<? posinf posinf)))
  (test-assert (not (fl<? zero nan)))
  (test-assert (not (fl<? nan one)))
  (test allfalse (map fl<? somereals somereals))
  (test (cdr alltrue) (map fl<? somereals (cdr somereals)))
  (test (cdr allfalse) (map fl<? (cdr somereals) somereals))

  (test-assert (not (fl>? zero fl-least)))
  (test-assert (not (fl>? fl-least fl-least)))
  (test-assert (fl>? one fl-least))
  (test-assert (not (fl>? neginf neginf)))
  (test-assert (not (fl>? neginf posinf)))
  (test-assert (fl>? posinf neginf))
  (test-assert (not (fl>? posinf posinf)))
  (test-assert (not (fl>? zero nan)))
  (test-assert (not (fl>? nan one)))
  (test allfalse (map fl>? somereals somereals))
  (test (cdr allfalse) (map fl>? somereals (cdr somereals)))
  (test (cdr alltrue) (map fl>? (cdr somereals) somereals))

  (test-assert (fl<=? zero fl-least))
  (test-assert (fl<=? fl-least fl-least))
  (test-assert (not (fl<=? one fl-least)))
  (test-assert (fl<=? neginf neginf))
  (test-assert (fl<=? neginf posinf))
  (test-assert (not (fl<=? posinf neginf)))
  (test-assert (fl<=? posinf posinf))
  (test-assert (not(fl<=? zero nan)))
  (test-assert (not (fl<=? nan one)))
  (test alltrue (map fl<=? somereals somereals))
  (test (cdr alltrue) (map fl<=? somereals (cdr somereals)))
  (test (cdr allfalse) (map fl<=? (cdr somereals) somereals))

  (test-assert   (not (fl>=? zero fl-least)))
  (test-assert (fl>=? fl-least fl-least))
  (test-assert (fl>=? one fl-least))
  (test-assert (fl>=? neginf neginf))
  (test-assert   (not (fl>=? neginf posinf)))
  (test-assert (fl>=? posinf neginf))
  (test-assert (fl>=? posinf posinf))
  (test-assert   (not (fl>=? zero nan)))
  (test-assert   (not (fl>=? nan one)))
  (test alltrue (map fl>=? somereals somereals))
  (test (cdr allfalse) (map fl>=? somereals (cdr somereals)))
  (test (cdr alltrue) (map fl>=? (cdr somereals) somereals))

  (test-assert (not (flunordered? zero fl-least)))
  (test-assert (not (flunordered? fl-least fl-least)))
  (test-assert (not (flunordered? one fl-least)))
  (test-assert (not (flunordered? neginf neginf)))
  (test-assert (not (flunordered? neginf posinf)))
  (test-assert (not (flunordered? posinf neginf)))
  (test-assert (not (flunordered? posinf posinf)))
  (test-assert (flunordered? zero nan))
  (test-assert (flunordered? nan one))
  (test allfalse (map flunordered? somereals somereals))
  (test (cdr allfalse) (map flunordered? somereals (cdr somereals)))
  (test (cdr allfalse) (map flunordered? (cdr somereals) somereals))

  )

(test neginf (flmax))
(test zero (flmax zero))
(test one (flmax zero one))
(test one (flmax one zero))
(test (car (reverse somereals)) (apply flmax somereals))

(test posinf (flmin))
(test one (flmin one))
(test zero (flmin zero one))
(test zero (flmin one zero))
(test (car somereals) (apply flmin somereals))

(test (map fl=?
           somereals
           (map flround somereals))
      (map flinteger? somereals))

(test-assert (not (flzero? neginf)))
(test-assert (not (flzero? (fl- fl-least))))
(test-assert (flzero? negzero))
(test-assert (flzero? zero))
(test-assert (not (flzero? fl-least)))
(test-assert (not (flzero? posinf)))

(test-assert (not (flpositive? neginf)))
(test-assert (not (flpositive? (fl- fl-least))))
(test-assert (not (flpositive? negzero)))
(test-assert (not (flpositive? zero)))
(test-assert (flpositive? fl-least))
(test-assert (flpositive? posinf))

(test-assert (flnegative? neginf))
(test-assert (flnegative? (fl- fl-least)))
(test-assert (not (flnegative? negzero)))    ; explicit in SRFI 144
(test-assert (not (flnegative? zero)))
(test-assert (not (flnegative? fl-least)))
(test-assert (not (flnegative? posinf)))

(test-assert (not (flodd? zero)))
(test-assert (flodd? one))

(test-assert (fleven? zero))
(test-assert   (not (fleven? one)))

(test (map (lambda (x) #t) somereals)
      (map flfinite? somereals))
(test (map (lambda (x) #f) weird)
      (map flfinite? weird))

(test-assert (flinfinite? neginf))
(test-assert (flinfinite? posinf))
(test-assert   (not (flinfinite? nan)))
(test (map (lambda (x) #f) somereals)
      (map flinfinite? somereals))

(test-assert (not (flnan? neginf)))
(test-assert (not (flnan? posinf)))
(test-assert (flnan? nan))
(test (map (lambda (x) #f) somereals)
      (map flnan? somereals))

(test-assert (flnormalized? fl-greatest))
(test-assert (not (flnormalized? fl-least)))
(test-assert (not (fldenormalized? fl-greatest)))
(test-assert (fldenormalized? fl-least))

;; Arithmetic

(test zero (fl+))
(test zero (fl+ zero))
(test #t (flzero? (fl+ negzero)))
(test one (fl+ one))
(test two (fl+ one one))
(test-assert (flnan? (fl+ nan one)))
(test-assert (flnan? (fl+ one nan)))
(test (map (lambda (x) (fl* (flonum 3) x))
           somereals)
      (map fl+ somereals somereals somereals))
(test infinities (map fl+ infinities infinities))
(test (map (lambda (x) #t) infinities)
      (map flnan?
           (map fl+ infinities (reverse infinities))))

(test one (fl*))
(test zero (fl* zero))
(test #t (flzero? (fl* negzero)))
(test one (fl* one))
(test one (fl* one one))
(test-assert (flnan? (fl* nan one)))
(test-assert (flnan? (fl* one nan)))
(test-group "(map (lambda (x) (flonum (expt x 3))) somereals)"
  (for-each
   (lambda (a b)
     (test a b))
   (map fl* somereals somereals somereals)
   (map (lambda (x) (flonum (expt x 3)))
        somereals)))
(test (map (lambda (x) posinf) infinities)
      (map fl* infinities infinities))
(test (map (lambda (x) neginf) infinities)
      (map fl* infinities (reverse infinities)))

(let ((three (flonum 3))
      (four  (flonum 4))
      (five  (flonum 5))
      (x23   (flonum 23))
      (ten11 (flonum (expt 10 11)))
      (ten12 (flonum (expt 10 12))))
  (test x23 (fl+* four five three))
  (test (flonum (+ (* (exact ten11) (exact ten12)) (exact one)))
        (fl+* ten11 ten12 one))
  (test (flonum (+ (* (exact ten11) (exact ten12)) (exact (fl- one))))
        (fl+* ten11 ten12 (fl- one)))

  ;; FIXME: the following test assumes IEEE double precision,
  ;; in which (expt 10 23) lies exactly halfway between the
  ;; two nearest flonums.

  (test-assert (not (fl=? (fl+* ten11 ten12 one)
                          (fl+* ten11 ten12 (fl- one))))))

(test-assert (flnan? (fl+* zero posinf one)))
(test-assert (flnan? (fl+* zero neginf one)))
(test-assert (flnan? (fl+* posinf zero one)))
(test-assert (flnan? (fl+* neginf zero one)))
(test-assert (flnan? (fl+* zero posinf nan)))
(test-assert (flnan? (fl+* zero neginf nan)))
(test-assert (flnan? (fl+* posinf zero nan)))
(test-assert (flnan? (fl+* neginf zero nan)))
(test neginf (fl+* fl-greatest fl-greatest neginf))
(test posinf (fl+* fl-greatest (fl- fl-greatest) posinf))
(test-assert (flnan? (fl+* nan one one)))
(test-assert (flnan? (fl+* one nan one)))
(test-assert (flnan? (fl+* one one nan)))

(test negzero (fl- zero))
(test zero (fl- negzero))
(test (flonum -1) (fl- one))
(test zero (fl- one one))
(test-assert (flnan? (fl- nan one)))
(test-assert (flnan? (fl- one nan)))
(test (map (lambda (x) (if (eqv? x zero) zero (fl- x)))
           somereals)
      (map fl- somereals somereals somereals))
(test '(#t #t)
      (map flnan? (map fl- infinities infinities)))
(test infinities
      (map fl- infinities (reverse infinities)))

(test posinf (fl/ zero))
(test neginf (fl/ negzero))
(test one (fl/ one))
(test one (fl/ one one))
(test-assert (flnan? (fl/ nan one)))
(test-assert (flnan? (fl/ one nan)))
(test (map (lambda (x) (if (flzero? x) (fl/ zero zero) (fl/ x)))
           somereals)
      (map fl/ somereals somereals somereals))
(test '(#t #t)
      (map flnan? (map fl/ infinities infinities)))
(test '(#t #t)
      (map flnan? (map fl/ infinities (reverse infinities))))

(test zero (flabs zero))
(test zero (flabs negzero))
(test one (flabs one))
(test (flonum 5.25) (flabs (flonum -5.25)))

(test one (flabsdiff zero one))
(test one (flabsdiff one zero))
(test zero (flabsdiff one one))
(test posinf (flabsdiff posinf neginf))
(test posinf (flabsdiff neginf posinf))
#|
(test (flposdiff zero one) zero)
(test (flposdiff one zero) one)
(test (flposdiff one one) zero)
(test (flposdiff posinf neginf) posinf)
(test (flposdiff neginf posinf) posinf)
|#
(test one (flsgn posinf))
(test (fl- one) (flsgn neginf))
(test one (flsgn zero))
(test (fl- one) (flsgn negzero))
(test one (flsgn two))
(test (fl- one) (flsgn (fl- two)))

(test (flonum 9) (flnumerator (flonum 2.25)))
(test (flonum 4) (fldenominator (flonum 2.25)))
(test (flonum -9) (flnumerator (flonum -2.25)))
(test (flonum 4) (fldenominator (flonum -2.25)))
(test ints (map flnumerator ints))
(test (map (lambda (x) one) ints)
      (map fldenominator ints))
;; CHICKEN dosen't work this way
;; (test weird (map flnumerator weird))
;; (test (map fldenominator infinities) (list one one))
;; (test-assert (flnan? (flnumerator nan)))
(test-assert (flnan? (fldenominator nan)))

(test (flonum -4) (flfloor    (flonum -3.125)))
(test (flonum -3) (flceiling  (flonum -3.125)))
(test (flonum -3) (flround    (flonum -3.125)))
(test (flonum -3) (fltruncate (flonum -3.125)))

(test (flonum -4) (flfloor    (flonum -3.75)))
(test (flonum -3) (flceiling  (flonum -3.75)))
(test (flonum -4) (flround    (flonum -3.75)))
(test (flonum -3) (fltruncate (flonum -3.75)))

(test (flonum -4) (flfloor    (flonum -3.5)))
(test (flonum -3) (flceiling  (flonum -3.5)))
(test (flonum -4) (flround    (flonum -3.5)))
(test (flonum -3) (fltruncate (flonum -3.5)))

(test ints (map flfloor    ints))
(test ints (map flceiling  ints))
(test ints (map flround    ints))
(test ints (map fltruncate ints))

(test (map (lambda (x) zero) posfracs) (map flfloor    posfracs))
(test (map (lambda (x) one) posfracs) (map flceiling  posfracs))
(test (map (lambda (x) zero) posfracs) (map flround    posfracs))
(test (map (lambda (x) zero) posfracs) (map fltruncate posfracs))

(test weird (map flfloor    weird))
(test weird (map flceiling  weird))
(test weird (map flround    weird))
(test weird (map fltruncate weird))

;; Exponents and logarithms

(test one (flexp negzero))
(test one (flexp zero))
(test fl-e (flexp one))
(test fl-1/e (flexp (fl- one)))
(test fl-e-2 (flexp two))
(test fl-e-pi/4 (flexp fl-pi/4))
(test posinf (flexp posinf))
(test posinf (flexp fl-greatest))
(test one (flexp fl-least))
(test zero (flexp (fl- fl-greatest)))
(test zero (flexp neginf))

(test one (fl+ one (flexp-1 negzero)))
(test one (fl+ one (flexp-1 zero)))
(test fl-e (fl+ one (flexp-1 one)))
(test fl-1/e (fl+ one (flexp-1 (fl- one))))
(test fl-e-2 (fl+ one (flexp-1 two)))
(test fl-e-pi/4 (fl+ one (flexp-1 fl-pi/4)))
(test posinf (fl+ one (flexp-1 posinf)))
(test posinf (fl+ one (flexp-1 fl-greatest)))
(test one (fl+ one (flexp-1 fl-least)))
(test zero (fl+ one (flexp-1 (fl- fl-greatest))))
(test zero (fl+ one (flexp-1 neginf)))

(test one (flexp2 negzero))
(test one (flexp2 zero))
(test two (flexp2 one))
(test (fl/ two) (flexp2 (fl- one)))
(test (fl* two two) (flexp2 two))
(test fl-e (flexp2 fl-log2-e))
(test fl-e (flexp2 fl-log2-e))
(test posinf (flexp2 posinf))
(test posinf (flexp2 fl-greatest))
(test one (flexp2 fl-least))
(test zero (flexp2 (fl- fl-greatest)))
(test zero (flexp2 neginf))

(test zero (flsquare zero))
(test one (flsquare one))
(test (fl+ two two) (flsquare two))
(test two (flsquare fl-sqrt-2))
(test (flonum 3) (flsquare fl-sqrt-3))
(test (flonum 5) (flsquare fl-sqrt-5))
(test (flonum 10) (flsquare fl-sqrt-10))
(test (flonum 25) (flsquare (flonum -5)))
(test posinf (flsquare neginf))
(test posinf (flsquare posinf))

(test zero (flsqrt zero))
(test one (flsqrt one))
(test fl-sqrt-2 (flsqrt two))
(test fl-sqrt-3 (flsqrt (flonum 3)))
(test fl-sqrt-5 (flsqrt (flonum 5)))
(test fl-sqrt-10 (flsqrt (flonum 10)))
(test (flonum 26.419689627245813) (flsqrt (flonum 698)))
(test posinf (flsqrt posinf))

(test zero (flcbrt zero))
(test one (flcbrt one))
(test fl-cbrt-2 (flcbrt two))
(test fl-cbrt-3 (flcbrt (flonum 3)))
(test (flonum 8.8705757224791313) (flcbrt (flonum 698)))
(test (flonum 2.25) (flcbrt (flonum 11.390625)))
(test (flonum -2.25) (flcbrt (flonum -11.390625)))
(test posinf (flcbrt posinf))
(test neginf (flcbrt neginf))

(test zero (flhypot zero zero))
(test one (flhypot zero one))
(test fl-sqrt-5 (flhypot two one))
(test fl-sqrt-5 (flhypot (fl- two) one))
(test fl-sqrt-5 (flhypot two (fl- one)))
(test fl-sqrt-5 (flhypot (fl- two) (fl- one)))
(test (fl/ fl-greatest fl-sqrt-2)
      (flhypot (fl/ fl-greatest two) (fl/ fl-greatest two)))
(test posinf (flhypot zero posinf))
(test posinf (flhypot neginf zero))

(test one (flexpt two zero))
(test two (flexpt two one))
(test (flonum 4) (flexpt two two))
(test fl-sqrt-2 (flexpt two (fl/ two)))
(test (flonum 2.7821842944695155e26)
      (flexpt (flonum 441) (flonum 10)))
(test (flonum 3.37977444523542851)
      (flexpt (flonum 441) (fl/ (flonum 5))))
(for-each (lambda (x)
            (for-each (lambda (frac)
                        (test x
                              (flexpt (flexpt x frac) (fl/ frac))))
                      posfracs))
          (filter flpositive? somereals))

(test neginf (fllog zero))
(test zero (fllog one))
(test fl-log-2 (fllog two))
(test fl-log-3 (fllog (flonum 3)))
(test fl-log-pi (fllog fl-pi))
(test fl-log-10 (fllog (flonum 10)))
(test posinf (fllog posinf))
(for-each (lambda (x)
            (test x (flexp (fllog x))))
          (filter flpositive? somereals))

(test neginf (fllog2 zero))
(test zero (fllog2 one))
(test one (fllog2 two))
(test fl-log2-e (fllog2 fl-e))
(test posinf (fllog2 posinf))
(for-each (lambda (x)
            (test x (flexpt two (fllog2 x))))
          (filter flpositive? somereals))

(test neginf (fllog10 zero))
(test zero (fllog10 one))
(test fl-log10-e (fllog10 fl-e))
(test one (fllog10 (flonum 10)))
(test posinf (fllog10 posinf))
(for-each (lambda (x)
            (test x (flexpt (flonum 10) (fllog10 x))))
          (filter flpositive? somereals))

(test-assert (flpositive? (fllog1+ fl-least)))
(test neginf (fllog1+ (fl- zero one)))
(test zero (fllog1+ (fl- one one)))
(test fl-log-2 (fllog1+ (fl- two one)))
(test fl-log-3 (fllog1+ (fl- (flonum 3) one)))
(test fl-log-pi (fllog1+ (fl- fl-pi one)))
(test fl-log-10 (fllog1+ (fl- (flonum 10) one)))
(test posinf (fllog1+ (fl- posinf one)))
(for-each (lambda (x)
            (test x (flexp (fllog1+ (fl- x one)))))
          (filter flpositive? somereals))

(test fl-log2-e ((make-fllog-base two) fl-e))
(test fl-log10-e ((make-fllog-base (flonum 10)) fl-e))
(test fl-log-2 ((make-fllog-base fl-e) two))
(for-each (lambda (base)
            (let ((f (make-fllog-base base)))
              (for-each (lambda (x)
                          (test x (flexpt (flonum base) (f x))))
                        (filter flpositive? somereals))))
          (map flonum '(3 7 19)))

;; Trigonometric functions

(test zero                       (flsin zero))
(test one                        (flcos zero))
(test zero                       (fltan zero))
(test 0.19866933079506121545941  (flsin (flonum 0.2)))
(test 0.98006657784124163112420  (flcos (flonum 0.2)))
(test 0.47942553860420300027329  (flsin (flonum 0.5)))
(test 0.87758256189037271611628  (flcos (flonum 0.5)))
(test 0.64421768723769105367261  (flsin (flonum 0.7)))
(test 0.76484218728448842625586  (flcos (flonum 0.7)))
(test fl-1/sqrt-2                (flsin fl-pi/4))
(test fl-1/sqrt-2                (flcos fl-pi/4))
(test 0.84147098480789651665250  (flsin one))
(test 0.54030230586813971740094  (flcos one))
(test one                        (flsin fl-pi/2))
(test zero                       (flcos fl-pi/2))
(test 0.90929742682568169539602  (flsin two))
(test -0.41614683654714238699757 (flcos two))
(test 0.14112000805986722210074  (flsin (flonum 3)))
(test -0.98999249660044545727157 (flcos (flonum 3)))
(test zero                       (flsin fl-pi))
(test (fl- one)                  (flcos fl-pi))
(test zero                       (flsin fl-2pi))
(test one                        (flcos fl-2pi))
(test -0.42818266949615100440675 (flsin (flonum 35)))
(test -0.90369220509150675984730 (flcos (flonum 35)))

(for-each (lambda (x)
            (test (fl- (flsin (fl- x))) (flsin x))
            (test (flcos (fl- fl-2pi x)) (flcos x))
            (test (fl/ (flsin x) (flcos x)) (fltan x))
            (test one (flhypot (flsin x) (flcos x))))
          (filter (lambda (x)
                    (and (flnormalized? x)
                         (fl<? (flabs x) (flonum 10000))))
                  (append somereals posfracs)))

(for-each (lambda (x)
            (test x (flsin (flasin x)))
            (test x (flcos (flacos x))))
          (filter (lambda (x) (fl<=? (fl- one) x one))
                  ;; Skip fl-greatest and fl-least, since CHICKEN gives us nan
                  ;; for acos of those
                  (let ((lst (cdr (append somereals posfracs))))
                    (cons (car lst)
                          (cddr lst)))))


(let ((xs (filter (lambda (x)
                    (and (flnormalized? x)
                         (fl<? (flabs x) (flonum 10000))))
                  (append somereals posfracs))))
  (for-each (lambda (x)
              (let ((theta (flatan x)))
                (test x (fl/ (flsin theta) (flcos theta))))
              (for-each (lambda (y)
                          (let ((theta (flatan x)))
                            (test theta
                                  (flatan (flsin theta)
                                          (flcos theta)))))
                        xs))
            xs))

(test zero              (flsinh zero))
(test one               (flcosh zero))
(test 0.201336002541094 (flsinh (flonum 0.2)))
(test 1.020066755619076 (flcosh (flonum 0.2)))
(test 0.521095305493747 (flsinh (flonum 0.5)))
(test 1.127625965206381 (flcosh (flonum 0.5)))
(test 0.758583701839534 (flsinh (flonum 0.7)))
(test 1.255169005630943 (flcosh (flonum 0.7)))
(test 1.175201193643801 (flsinh one))
(test 1.543080634815244 (flcosh one))
(test 3.626860407847019 (flsinh two))
(test 3.762195691083631 (flcosh two))
(test 10.01787492740990 (flsinh (flonum 3)))
(test 10.06766199577777 (flcosh (flonum 3)))
(test 11013.23287470339 (flsinh (flonum 10)))
(test 11013.23292010332 (flcosh (flonum 10)))

(for-each (lambda (x)
            (test x (flasinh (flsinh x)))
            (test (flabs x) (flacosh (flcosh x)))
            (test x (flatanh (fltanh x))))
          (filter (lambda (x)
                    (fl<? (flabs x) (flonum 100)))
                  (append somereals posfracs (map fl- posfracs))))

;; Integer division

(test (flonum 19)    (flquotient  (flonum 9.75) (flonum 0.5)))
(test (flonum 0.25)  (flremainder (flonum 9.75) (flonum 0.5)))
(test (flonum -19)   (flquotient  (flonum -9.75) (flonum 0.5)))
(test (flonum -0.25) (flremainder (flonum -9.75) (flonum 0.5)))
(test (flonum -19)   (flquotient  (flonum 9.75) (flonum -0.5)))
(test (flonum 0.25)  (flremainder (flonum 9.75) (flonum -0.5)))
(test (flonum 19)    (flquotient  (flonum -9.75) (flonum -0.5)))
(test (flonum -0.25) (flremainder (flonum -9.75) (flonum -0.5)))

(let ((f (lambda (x y)
           (call-with-values
               (lambda () (flremquo x y))
             (lambda (r q)
               (list r (modulo q 8)))))))

  (test (list (flonum -0.125) 0) (f (flonum 15.875) (flonum 0.5)))
  (test (list (flonum +0.000) 0) (f (flonum 16.000) (flonum 0.5)))
  (test (list (flonum +0.125) 0) (f (flonum 16.125) (flonum 0.5)))
  (test (list (flonum +0.250) 0) (f (flonum 16.250) (flonum 0.5)))

  (test (list (flonum -0.125) 2) (f (flonum 16.875) (flonum 0.5)))
  (test (list (flonum +0.000) 2) (f (flonum 17.000) (flonum 0.5)))
  (test (list (flonum +0.125) 2) (f (flonum 17.125) (flonum 0.5)))
  (test (list (flonum +0.250) 2) (f (flonum 17.250) (flonum 0.5)))

  (test (list (flonum -0.125) 4) (f (flonum 17.875) (flonum 0.5)))
  (test (list (flonum +0.000) 4) (f (flonum 18.000) (flonum 0.5)))
  (test (list (flonum +0.125) 4) (f (flonum 18.125) (flonum 0.5)))
  (test (list (flonum +0.250) 4) (f (flonum 18.250) (flonum 0.5)))

  (test (list (flonum -0.125) 6) (f (flonum 18.875) (flonum 0.5)))
  (test (list (flonum +0.000) 6) (f (flonum 19.000) (flonum 0.5)))
  (test (list (flonum +0.125) 6) (f (flonum 19.125) (flonum 0.5)))
  (test (list (flonum +0.250) 6) (f (flonum 19.250) (flonum 0.5)))

  (test (list (flonum -0.125) 7) (f (flonum 15.375) (flonum 0.5)))
  (test (list (flonum +0.000) 1) (f (flonum 16.500) (flonum 0.5)))
  (test (list (flonum +0.125) 1) (f (flonum 16.625) (flonum 0.5)))
  (test (list (flonum -0.250) 2) (f (flonum 16.750) (flonum 0.5)))

  (test (list (flonum -0.125) 1) (f (flonum 16.375) (flonum 0.5)))
  (test (list (flonum +0.000) 3) (f (flonum 17.500) (flonum 0.5)))
  (test (list (flonum +0.125) 3) (f (flonum 17.625) (flonum 0.5)))
  (test (list (flonum -0.250) 4) (f (flonum 17.750) (flonum 0.5)))

  (test (list (flonum -0.125) 3) (f (flonum 17.375) (flonum 0.5)))
  (test (list (flonum +0.000) 5) (f (flonum 18.500) (flonum 0.5)))
  (test (list (flonum +0.125) 5) (f (flonum 18.625) (flonum 0.5)))
  (test (list (flonum -0.250) 6) (f (flonum 38.750) (flonum 0.5)))

  (test (list (flonum -0.125) 5) (f (flonum 18.375) (flonum 0.5)))
  (test (list (flonum +0.000) 7) (f (flonum 19.500) (flonum 0.5)))
  (test (list (flonum +0.125) 7) (f (flonum 19.625) (flonum 0.5)))
  (test (list (flonum -0.250) 0) (f (flonum 19.750) (flonum 0.5)))

  )

;; Special functions

(test fl-gamma-1/2 (flgamma (flonum 0.5)))
(test fl-gamma-1/3 (flgamma (flonum #i1/3)))
(test fl-gamma-2/3 (flgamma (flonum #i2/3)))
(test (flonum 1.225416702465178) (flgamma (flonum 0.75)))
(test one (flgamma one))
(test (flonum 0.906402477055477) (flgamma (flonum 1.25)))
(test (flonum 0.888537149431357) (flgamma (flonum 1.38)))
(test (fl/ (flsqrt fl-pi) two) (flgamma (flonum 1.50)))
(test one (flgamma two))
(test (flonum 3628800) (flgamma (flonum 11)))
(test (* 10.5 9.5 8.5 7.5 6.5 5.5 4.5 3.5 2.5 1.5 0.5
         fl-gamma-1/2)
      (flgamma (flonum 11.5)))
(test (* 0.3 -0.7 -1.7 -2.7 -3.7 -4.7 -5.7 -6.7
         (flgamma (flonum -6.7)))
      (flgamma (flonum 1.3)))
(test (* 0.3 -0.7 -1.7 -2.7 -3.7 -4.7 -5.7 -6.7 -7.7
         (flgamma (flonum -7.7)))
      (flgamma (flonum 1.3)))
(test (flonum 2.159256448318447e110)
      (flgamma (flonum 76.5)))
(test (flonum 9.332621544439442e155)
      (flgamma (flonum 100)))

(let* ((f (lambda (x)
            (call-with-values
                (lambda () (flloggamma x))
              list)))
       (g (lambda (x) (car (f x))))
       (s (lambda (x) (cadr (f x)))))
  (test (log fl-gamma-1/2) (g (flonum 0.5)))
  (test (log fl-gamma-1/3) (g (flonum #i1/3)))
  (test (log fl-gamma-2/3) (g (flonum #i2/3)))
  (test (log (flonum 1.225416702465178)) (g (flonum 0.75)))
  (test (log one) (g one))
  (test (log (flonum 0.906402477055477)) (g (flonum 1.25)))
  (test (log (flonum 0.888537149431357)) (g (flonum 1.38)))
  (test (log (fl/ (flsqrt fl-pi) two)) (g (flonum 1.50)))
  (test (log one) (g two))
  (test (log (flonum 3628800)) (g (flonum 11)))
  (test (apply + (map log (cdr (iota 1000))))
        (g (flonum 1000)))
  (test (log (flabs (flgamma (flonum -0.5))))
        (g (flonum -0.5)))
  (test (log (flabs (flgamma (flonum -49.2))))
        (g (flonum -49.2)))
  (test (log (flabs (flgamma (flonum -50.3))))
        (g (flonum -50.3)))
  (test (log (flabs (flgamma (flonum -100.5))))
        (g (flonum -100.5)))
  (test posinf (g posinf))
  ;; Apparently this is fine in CHICKEN
  ;; (test-error (g neginf))
  ;; (test-error (g zero))
  ;; (test-error (g nan))

  (test one (s (flonum 0.5)))
  (test one (s (flonum #i1/3)))
  (test one (s (flonum #i2/3)))
  (test one (s (flonum 0.75)))
  (test one (s one))
  (test one (s (flonum 1.25)))
  (test one (s (flonum 1.38)))
  (test one (s (flonum 1.50)))
  (test one (s two))
  (test one (s (flonum 11)))
  (test one (s (flonum 1000)))
  (test (fl- one)  (s (flonum -0.5)))
  (test one (s (flonum -49.2)))
  (test (fl- one) (s (flonum -50.3)))
  (test one (s (flonum -999.5)))
  (test (fl- one) (s (flonum -1000.5)))
  )

(test one (flfirst-bessel 0 zero))
(test zero (flfirst-bessel 1 zero))
(test zero (flfirst-bessel 2 zero))
(test zero (flfirst-bessel 3 zero))
(test zero (flfirst-bessel 4 zero))
(test zero (flfirst-bessel 5 zero))
(test zero (flfirst-bessel 6 zero))
(test zero (flfirst-bessel 7 zero))
(test zero (flfirst-bessel 8 zero))
(test zero (flfirst-bessel 9 zero))

(test (flonum 0.960398226659563) (flfirst-bessel 0 (flonum 0.4)))
(test (flonum 0.196026577955319) (flfirst-bessel 1 (flonum 0.4)))
(test (flonum 0.01973466311703027) (flfirst-bessel 2 (flonum 0.4)))
(test (flonum 1.320053214983958e-3) (flfirst-bessel 3 (flonum 0.4)))
(test (flonum 6.613510772909677e-5) (flfirst-bessel 4 (flonum 0.4)))
(test (flonum 2.648939597977585e-6) (flfirst-bessel 5 (flonum 0.4)))
(test (flonum 8.838222034285578e-8) (flfirst-bessel 6 (flonum 0.4)))
(test (flonum 2.527012308088268e-9) (flfirst-bessel 7 (flonum 0.4)))
(test (flonum 6.321044023358773e-11) (flfirst-bessel 8 (flonum 0.4)))
(test (flonum 1.405301255241872e-12) (flfirst-bessel 9 (flonum 0.4)))

(test (flonum 0.765197686557967) (flfirst-bessel 0 one))
(test (flonum 0.440050585744934) (flfirst-bessel 1 one))
(test (flonum 0.114903484931900) (flfirst-bessel 2 one))
(test (flonum 1.956335398266841e-2) (flfirst-bessel 3 one))
(test (flonum 2.476638964109955e-3) (flfirst-bessel 4 one))
(test (flonum 2.497577302112344e-4) (flfirst-bessel 5 one))
(test (flonum 2.093833800238927e-5) (flfirst-bessel 6 one))
(test (flonum 1.502325817436808e-6) (flfirst-bessel 7 one))
(test (flonum 9.422344172604501e-8) (flfirst-bessel 8 one))
(test (flonum 5.249250179911875e-9) (flfirst-bessel 9 one))

(test (flonum -0.1698542521511835) (flfirst-bessel 0 (flonum 17)))
(test (flonum -0.09766849275778065) (flfirst-bessel 1 (flonum 17)))
(test (flonum 0.1583638412385035) (flfirst-bessel 2 (flonum 17)))
(test (flonum 0.1349305730591932) (flfirst-bessel 3 (flonum 17)))
(test (flonum -0.1107412860446706) (flfirst-bessel 4 (flonum 17)))
(test (flonum -0.1870441194231559) (flfirst-bessel 5 (flonum 17)))
(test (flonum 0.0007153334428141831) (flfirst-bessel 6 (flonum 17)))
(test (flonum 0.1875490606769070) (flfirst-bessel 7 (flonum 17)))
(test (flonum 0.1537368341734622) (flfirst-bessel 8 (flonum 17)))
(test (flonum -0.04285556969011908) (flfirst-bessel 9 (flonum 17)))

(test (flonum +5.581232766925182e-2) (flfirst-bessel 0 (flonum 50)))
(test (flonum -9.751182812517514e-2) (flfirst-bessel 1 (flonum 50)))
(test (flonum -5.971280079425882e-2) (flfirst-bessel 2 (flonum 50)))
(test (flonum +9.273480406163443e-2) (flfirst-bessel 3 (flonum 50)))
(test (flonum +7.084097728165495e-2) (flfirst-bessel 4 (flonum 50)))
(test (flonum -8.140024769656964e-2) (flfirst-bessel 5 (flonum 50)))
(test (flonum -8.712102682096888e-2) (flfirst-bessel 6 (flonum 50)))
(test (flonum +6.049120125953711e-2) (flfirst-bessel 7 (flonum 50)))
(test (flonum +1.040585631736393e-1) (flfirst-bessel 8 (flonum 50)))
(test (flonum -2.719246104397254e-2) (flfirst-bessel 9 (flonum 50)))

(let ((f (lambda (n x y)
           (test (flexpt x (flonum n))
                 (fl/ (flfirst-bessel n x) y)))))
  (f 10 (flonum 8) 5.659370382712330e-11)
  (f 11 (flonum 8) 2.979844833403859e-12)
  (f 20 (flonum 8) 1.804618055659560e-25)
  (f 21 (flonum 8) 4.456237531615702e-27))

(test 2.074861066333589e-1 (flfirst-bessel 10 (flonum 10)))
(test 4.507973143721253e-3 (flfirst-bessel 15 (flonum 10)))
(test 1.151336924781340e-5 (flfirst-bessel 20 (flonum 10)))
(test 1.784513607871595e-30 (flfirst-bessel 50 (flonum 10)))
(test 6.597316064155381e-89 (flfirst-bessel 100 (flonum 10)))

(test -5.473217693547201e-2 (flfirst-bessel 10 (flonum 100)))
(test +1.519812122392732e-2 (flfirst-bessel 15 (flonum 100)))
(test +6.221745849833875e-2 (flfirst-bessel 20 (flonum 100)))
(test -3.869833972852538e-2 (flfirst-bessel 50 (flonum 100)))
(test +9.636667329586156e-2 (flfirst-bessel 100 (flonum 100)))


(test neginf (flsecond-bessel 0 zero))
(test neginf (flsecond-bessel 1 zero))
(test neginf (flsecond-bessel 2 zero))
(test neginf (flsecond-bessel 3 zero))
(test neginf (flsecond-bessel 4 zero))
(test neginf (flsecond-bessel 5 zero))
(test neginf (flsecond-bessel 6 zero))
(test neginf (flsecond-bessel 7 zero))
(test neginf (flsecond-bessel 8 zero))
(test neginf (flsecond-bessel 9 zero))
(test neginf (flsecond-bessel 10 zero))
(test neginf (flsecond-bessel 20 zero))

(test (flonum -0.6060245684270096) (flsecond-bessel 0 (flonum 0.4)))
(test (flonum -1.780872044270051) (flsecond-bessel 1 (flonum 0.4)))
(test (flonum -8.298335652923247) (flsecond-bessel 2 (flonum 0.4)))
(test (flonum -8.120248448496242e1) (flsecond-bessel 3 (flonum 0.4)))
(test (flonum -1.209738931621513e3) (flsecond-bessel 4 (flonum 0.4)))
(test (flonum -2.411357614794530e4) (flsecond-bessel 5 (flonum 0.4)))
(test (flonum -6.016296647670110e5) (flsecond-bessel 6 (flonum 0.4)))
(test (flonum -1.802477636686238e7) (flsecond-bessel 7 (flonum 0.4)))
(test (flonum -6.302655431754164e8) (flsecond-bessel 8 (flonum 0.4)))
(test (flonum -2.519259695064979e10) (flsecond-bessel 9 (flonum 0.4)))

(test (flonum +0.08825696421567696) (flsecond-bessel 0 one))
(test (flonum -0.7812128213002887) (flsecond-bessel 1 one))
(test (flonum -1.650682606816254) (flsecond-bessel 2 one))
(test (flonum -5.821517605964729) (flsecond-bessel 3 one))
(test (flonum -3.327842302897212e1) (flsecond-bessel 4 one))
(test (flonum -2.604058666258122e2) (flsecond-bessel 5 one))
(test (flonum -2.570780243229150e3) (flsecond-bessel 6 one))
(test (flonum -3.058895705212399e4) (flsecond-bessel 7 one))
(test (flonum -4.256746184865067e5) (flsecond-bessel 8 one))
(test (flonum -6.780204938731983e6) (flsecond-bessel 9 one))

(test (flonum -0.09263719844232369) (flsecond-bessel 0 (flonum 17)))
(test (flonum +0.1672050360772337) (flsecond-bessel 1 (flonum 17)))
(test (flonum +0.1123083791572924) (flsecond-bessel 2 (flonum 17)))
(test (flonum -0.1407795350990472) (flsecond-bessel 3 (flonum 17)))
(test (flonum -0.1619952738981326) (flsecond-bessel 4 (flonum 17)))
(test (flonum +0.06454646502933781) (flsecond-bessel 5 (flonum 17)))
(test (flonum +0.1999637827389195) (flsecond-bessel 6 (flonum 17)))
(test (flonum +0.07660444043342891) (flsecond-bessel 7 (flonum 17)))
(test (flonum -0.1368777729702134) (flsecond-bessel 8 (flonum 17)))
(test (flonum -0.2054305796995121) (flsecond-bessel 9 (flonum 17)))

(test (flonum -9.806499547007708e-2) (flsecond-bessel 0 (flonum 50)))
(test (flonum -5.679566856201477e-2) (flsecond-bessel 1 (flonum 50)))
(test (flonum +9.579316872759649e-2) (flsecond-bessel 2 (flonum 50)))
(test (flonum +6.445912206022249e-2) (flsecond-bessel 3 (flonum 50)))
(test (flonum -8.805807408036979e-2) (flsecond-bessel 4 (flonum 50)))
(test (flonum -7.854841391308165e-2) (flsecond-bessel 5 (flonum 50)))
(test (flonum +7.234839129775346e-2) (flsecond-bessel 6 (flonum 50)))
(test (flonum +9.591202782454248e-2) (flsecond-bessel 7 (flonum 50)))
(test (flonum -4.549302350688156e-2) (flsecond-bessel 8 (flonum 50)))
(test (flonum -1.104697953467446e-1) (flsecond-bessel 9 (flonum 50)))

(let ((f (lambda (n x y)
           (test (flexpt x (flonum (- n)))
                 (fl/ (flsecond-bessel n x) y)))))
  (f 10 (flonum 8) -0.97356279e9)
  (f 20 (flonum 8) -0.962608e23)
  )

(test -3.598141521834027e-1 (flsecond-bessel 10 (flonum 10)))
(test -6.364745876939129e0  (flsecond-bessel 15 (flonum 10)))
(test -1.597483848269626e3  (flsecond-bessel 20 (flonum 10)))
(test -3.641066501800740e27 (flsecond-bessel 50 (flonum 10)))
(test -4.849148271180607e85 (flsecond-bessel 100 (flonum 10)))

(test +5.833157423641493e-2 (flsecond-bessel 10 (flonum 100)))
(test +7.879068694670284e-2 (flsecond-bessel 15 (flonum 100)))
(test +5.124797307618842e-2 (flsecond-bessel 20 (flonum 100)))
(test +7.650526394480304e-2 (flsecond-bessel 50 (flonum 100)))
(test -1.669214114175765e-1 (flsecond-bessel 100 (flonum 100)))


(test zero (flerf zero))
(test one (flerf posinf))
(test (fl- one) (flerf neginf))

(test -0.5204998778130465 (flerf (flonum -0.5)))
(test -0.8427007929497149 (flerf (flonum -1.0)))
(test -0.9661051464753107 (flerf (flonum -1.5)))
(test -0.9953222650189527 (flerf (flonum -2.0)))

(test 0.2763263901682369  (flerf (flonum 0.25)))
(test 0.5204998778130465  (flerf (flonum 0.50)))
(test 0.7111556336535151  (flerf (flonum 0.75)))
(test 0.8427007929497149  (flerf (flonum 1.00)))
(test 0.9229001282564582  (flerf (flonum 1.25)))
(test 0.9661051464753107  (flerf (flonum 1.50)))
(test 0.9866716712191824  (flerf (flonum 1.75)))
(test 0.9953222650189527  (flerf (flonum 2.00)))

(test one (flerfc zero))
(test zero (flerfc posinf))
(test two (flerfc neginf))

(test 1.520499877813047 (flerfc (flonum -0.5)))
(test 1.842700792949715 (flerfc (flonum -1.0)))
(test 1.966105146475310 (flerfc (flonum -1.5)))
(test 1.995322265018953 (flerfc (flonum -2.0)))

(test (- 1.0 0.2763263901682369) (flerfc (flonum 0.25)))
(test (- 1.0 0.5204998778130465) (flerfc (flonum 0.50)))
(test (- 1.0 0.7111556336535151) (flerfc (flonum 0.75)))
(test (- 1.0 0.8427007929497149) (flerfc (flonum 1.00)))
(test (- 1.0 0.9229001282564582) (flerfc (flonum 1.25)))
(test (- 1.0 0.9661051464753107) (flerfc (flonum 1.50)))
(test (- 1.0 0.9866716712191824) (flerfc (flonum 1.75)))
(test (- 1.0 0.9953222650189527) (flerfc (flonum 2.00)))

(let* ((f (lambda (x y)
            (test y
                  (fl* x (flexp (flsquare x)) (flerfc x)))))
       (g (lambda (x^-2 y)
            (f (flsqrt (fl/ x^-2)) y))))
  (g (flonum 0.250) 0.5107913526210115)
  (g (flonum 0.125) 0.5340672374463438)
  (g (flonum 0.100) 0.5394141079847177)
  (g (flonum 0.075) 0.5450551142044333)
  (g (flonum 0.050) 0.5510294916242269)
  (g (flonum 0.025) 0.5573864515094235)
  (g (flonum 0.020) 0.5587090260270524)
  (g (flonum 0.015) 0.5600499741400112)
  (g (flonum 0.010) 0.5614099274382259)
  (g (flonum 0.005) 0.5627895581750243))

(test-end)
(test-exit)