  Copyright ( C )    ( 2016 ) .
;;; 
;;; 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. 

;;; References
;;;
    and    [ editors ] .
;;; Handbook of Mathematical Functions With Formulas, Graphs, and
  Mathematical Tables .   United States Department of Commerce .
  National Bureau of Standards Applied Mathematics Series , 55 ,
  June 1964 .   Fifth Printing , August 1966 , with corrections .
;;;
  R W Hamming .   Numerical Methods for Scientists and Engineers .
  McGraw - Hill , 1962 .
;;;
  Donald E Knuth .   The Art of Computer Programming , Volume 2 ,
  Seminumerical Algorithms , Second Edition .   Addison - Wesley , 1981 .
;;;
  J N Newman .   Approximations for the  and   .
  Mathematics of Computation , 43(168 ) , October 1984 , pages 551 - 556 .

;;; I have deliberately avoided recent references, and have also
;;; avoided looking at any code or pseudocode for these or similar
;;; functions.

  Quick - and - dirty implementation of a draft of SRFI 144 ( flonums ) ,
  as specified at /~cowan/temp/srfi-144.html
  as of 4 June 2017 .
;;;
;;; FIXME: not as accurate as it should be
;;; FIXME: not as fast as it should be
  FIXME : assumes IEEE arithmetic or similar
;;; FIXME: assumes all inexact reals are flonums
;;; FIXME: assumes (scheme inexact)

  Mathematical Constants
;;;
;;; The mathematical constants are now defined in 144.constants.scm


;; Implementation Constants

(define fl-greatest
  (let loop ((x (- (expt 2.0 precision-bits) 1.0)))
    (if (finite? (* 2.0 x))
        (loop (* 2.0 x))
        x)))

(define fl-least
  (let loop ((x 1.0))
    (if (> (/ x 2.0) 0.0)
        (loop (/ x 2.0))
        x)))

(define fl-epsilon
  (let loop ((eps 1.0))
    (if (= 1.0 (+ 1.0 eps))
        (* 2.0 eps)
        (loop (/ eps 2.0)))))

(define fl-integer-exponent-zero                ; arbitrary
  (exact (- (log fl-least 2.0) 1.0)))

(define fl-integer-exponent-nan                 ; arbitrary
  (- fl-integer-exponent-zero 1))

;;; Constructors

(define (flonum x)
  (if (real? x)
      (inexact x)
      (error "bad argument passed to flonum" x)))

(define fladjacent
  (flop2 'fladjacent
         (lambda (x y)
           (define (loop y)
             (let* ((y3 (fl+ (fl* 0.999755859375 x) (fl* 0.000244140625 y))))
               (cond ((fl<? x y3 y)
                      (loop y3))
                     ((fl<? y y3 x)
                      (loop y3))
                     (else
                      (loop2 y)))))
           (define (loop2 y)
             (let* ((y2 (fl/ (fl+ x y) 2.0))
                    (y2 (if (flinfinite? y2)
                            (fl+ (fl* 0.5 x) (fl* 0.5 y))
                            y2)))
               (cond ((fl=? x y2)
                      y)
                     ((fl=? y y2)
                      y)
                     (else
                      (loop2 y2)))))
           (cond ((flinfinite? x)
                  (cond ((fl<? x y) (fl- fl-greatest))
                        ((fl>? x y) fl-greatest)
                        (else x)))
                 ((fl=? x y)
                  x)
                 ((flzero? x)
                  (if (flpositive? y)
                      fl-least
                      (fl- fl-least)))
                 ((fl<? x y)
                  (loop (flmin y
                               fl-greatest
                               (flmax (* 2.0 x)
                                      (* 0.5 x)))))
                 ((fl>? x y)
                  (loop (flmax y
                               (fl- fl-greatest)
                               (flmin (* 2.0 x)
                                      (* 0.5 x)))))
  x or y is a NaN
                  x)))))

(define flcopysign
  (flop2 'flcopysign
         (lambda (x y)
           (if (= (flsign-bit x) (flsign-bit y))
               x
               (fl- x)))))

(define (make-flonum x n)
  (let ((y (expt 2.0 n)))
    (cond ((or (not (flonum? x))
               (not (exact-integer? n)))
           (error "bad arguments to make-flonum" x n))
          ((finite? y)
           (* x y))
          (else
           (inexact (* (exact x) (expt 2 n)))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  Accessors

(define (flinteger-fraction x)
  (check-flonum! 'flinteger-fraction x)
  (let* ((result1 (fltruncate x))
         (result2 (fl- x result1)))
    (values result1 result2)))

(define (flexponent x)
  (floor (fllog2 (flabs x))))

(define (flinteger-exponent x)
  (exact (flexponent x)))

(define (flnormalized-fraction-exponent x)
  (define (return result1 result2)
    (cond ((fl<? result1 0.5)
           (values (fl* 2.0 result1) (- result2 1)))
          ((fl>=? result1 1.0)
           (values (fl* 0.5 result1) (+ result2 1)))
          (else
           (values result1 result2))))
  (check-flonum! 'flnormalized-fraction-exponent x)
  unspecified for NaN
         (values x 0))
        ((fl<? x 0.0)
         (call-with-values
          (lambda () (flnormalized-fraction-exponent (fl- x)))
          (lambda (y n)
            (values (fl- y) n))))
  unspecified for 0.0
         (values 0.0 0))
        ((flinfinite? x)
         (values 0.5 (+ 3 (exact (round (fllog2 fl-greatest))))))
        ((flnormalized? x)
         (let* ((result2 (exact (flround (fllog2 x))))
                (result2 (if (integer? result2)
                             result2
                             (round result2)))
                (two^result2 (inexact (expt 2.0 result2))))
           (if (flinfinite? two^result2)
               (call-with-values
                (lambda () (flnormalized-fraction-exponent (fl/ x 4.0)))
                (lambda (y n)
                  (values y (+ n 2))))
               (return (fl/ x two^result2) result2))))
        (else
         (let* ((k (+ 2 precision-bits))
                (two^k (expt 2 k)))
           (call-with-values
            (lambda ()
              (flnormalized-fraction-exponent (fl* x (inexact two^k))))
            (lambda (y n)
              (return y (- n k))))))))

(define (flsign-bit x)
  (check-flonum! 'flsign-bit x)
  (cond ((fl<? x 0.0)
         1)
        ((eqv? x -0.0)
         1)
        (else
         0)))


;;; Predicates

;(define flonum? R6RS)               ; defined by (rnrs arithmetic flonums)
;(define fl=? R6RS)                  ; defined by (rnrs arithmetic flonums)
;(define fl<? R6RS)                  ; defined by (rnrs arithmetic flonums)
;(define fl>? R6RS)                  ; defined by (rnrs arithmetic flonums)
;(define fl<=? R6RS)                 ; defined by (rnrs arithmetic flonums)
;(define fl>=? R6RS)                 ; defined by (rnrs arithmetic flonums)

(define (flunordered? x y)
  (or (flnan? x) (flnan? y)))

  incompatible with ( rnrs arithmetic flonums ) in zero - argument case

(define flmax
  (let ((flmax2 (flop2 'flmax max)))
    (lambda args
      (cond ((null? args)
             -inf.0)
            ((null? (cdr args))
             (car args))
            ((null? (cddr args))
             (flmax2 (car args) (cadr args)))
            (else
             (flmax2 (flmax2 (car args) (cadr args))
                     (apply flmax (cddr args))))))))

  incompatible with ( rnrs arithmetic flonums ) in zero - argument case

(define flmin
  (let ((flmin2 (flop2 'flmin min)))
    (lambda args
      (cond ((null? args)
             +inf.0)                 ; spec says fl-least, but that's wrong
            ((null? (cdr args))
             (car args))
            ((null? (cddr args))
             (flmin2 (car args) (cadr args)))
            (else
             (flmin2 (flmin2 (car args) (cadr args))
                     (apply flmin (cddr args))))))))

;(define flinteger? R6RS)            ; defined by (rnrs arithmetic flonums)
;(define flzero? R6RS)               ; defined by (rnrs arithmetic flonums)
( define  ? R6RS )            ; defined by ( rnrs arithmetic flonums )
;(define flnegative? R6RS)           ; defined by (rnrs arithmetic flonums)
( define flodd ? R6RS )                 ; defined by ( rnrs arithmetic flonums )
;(define fleven? R6RS)               ; defined by (rnrs arithmetic flonums)
( define flfinite ? R6RS )              ; defined by ( rnrs arithmetic flonums )
;(define flinfinite? R6RS)           ; defined by (rnrs arithmetic flonums)
( define flnan ? R6RS )                 ; defined by ( rnrs arithmetic flonums )

(define flnormalized?
  (lambda (x)
    (check-flonum! 'flnormalized? x)
    (let ((x (flabs x)))
      (and (flfinite? x)
           (fl<? (fl/ fl-greatest) x)))))

(define fldenormalized?
  (lambda (x)
    (check-flonum! 'fldenormalized? x)
    (let ((x (flabs x)))
      (and (flfinite? x)
           (fl<? 0.0 x)
           (fl<=? x (fl/ fl-greatest))))))

;;; Arithmetic

;(define fl+ R6RS)                   ; defined by (rnrs arithmetic flonums)
;(define fl* R6RS)                   ; defined by (rnrs arithmetic flonums)

  Spec says " as if to infinite precision and rounded only once " .

(define fl+*
  (flop3 'fl+*
         (lambda (x y z)
           (cond (c-functions-are-available
                  (fma x y z))
                 ((and (flfinite? x) (flfinite? y))
                  (if (flfinite? z)
                      (let ((x (exact x))
                            (y (exact y))
                            (z (exact z)))
                        (flonum (+ (* x y) z)))
                      z))
                 (else
                  (fl+ (fl* x y) z))))))

;(define fl- R6RS)                   ; defined by (rnrs arithmetic flonums)
;(define fl/ R6RS)                   ; defined by (rnrs arithmetic flonums)
;(define flabs R6RS)                 ; defined by (rnrs arithmetic flonums)

(define (flabsdiff x y)
  (flabs (fl- x y)))

(define (flposdiff x y)
  (let ((diff (fl- x y)))
    (if (flnegative? diff)
        0.0
        diff)))

(define (flsgn x)
  (flcopysign 1.0 x))

  ( flnumerator + nan.0 ) and ( fldenominator + nan.0 ) must be NaNs , which
  is not required by the R6RS specification of ( rnrs arithmetic flonums ) .

(define flnumerator
  (flop1 'flnumerator
         (lambda (x)
           (if (flnan? x)
               x
               (r6rs:flnumerator x)))))

(define fldenominator
  (flop1 'fldenominator
         (lambda (x)
           (if (flnan? x)
               x
               (r6rs:fldenominator x)))))

;(define flfloor R6RS)               ; defined by (rnrs arithmetic flonums)
;(define flceiling R6RS)             ; defined by (rnrs arithmetic flonums)
;(define flround R6RS)               ; defined by (rnrs arithmetic flonums)
;(define fltruncate R6RS)            ; defined by (rnrs arithmetic flonums)

;;; Exponents and logarithms

;(define flexp R6RS)                 ; defined by (rnrs arithmetic flonums)

(define flexp2 (flop1 'flexp2 (lambda (x) (flexpt 2.0 x))))

  e^x = \sum_n ( z^n / ( n ! ) )
;;;
  FIXME : the number of terms and the constant 0.5 seem reasonable
  for IEEE double precision , but the number of terms might need
;;; to be increased for higher precisions.

(define flexp-1
  (flop1 'flexp-1
         (let ((coefs (cons 0.0
                            (map fl/
                                 (map factorial
                                      '(1.0 2.0 3.0 4.0 5.0
                                        6.0 7.0 8.0 9.0 10.0
                                        11.0 12.0 13.0 14.0 15.0))))))
           (lambda (x)
  FIXME
                    (polynomial-at x coefs))
                   (else
                    (fl- (flexp x) 1.0)))))))

(define flsquare (flop1 'flsquare (lambda (x) (fl* x x))))

( define flsqrt R6RS )                 ; defined by ( rnrs arithmetic flonums )

(define flcbrt
  (flop1 'flcbrt
         (lambda (x)
           (cond ((flnegative? x)
                  (fl- (flcbrt (fl- x))))
                 (else
                  (flexpt x (fl/ 3.0)))))))

(define flhypot
  (flop2 'flhypot
         (lambda (x y)
           (cond ((flzero? x) (flabs y))
                 ((flzero? y) (flabs x))
                 ((or (flinfinite? x) (flinfinite? y)) +inf.0)
                 ((flnan? x) x)
                 ((flnan? y) y)
                 ((fl>? y x) (flhypot y x))
                 (else
                  (let* ((y/x (fl/ y x))
                         (root (flsqrt (fl+ 1.0 (fl* y/x y/x)))))
                    (fl* (flabs x) root)))))))

;(define flexpt R6RS)                ; defined by (rnrs arithmetic flonums)
;(define fllog R6RS)                 ; defined by (rnrs arithmetic flonums)

    ) , as in C99 log1p .
;;;
;;; log (x + 1) = \sum_{n=1}^\infty - (-1)^n x^n/n

(define fllog1+
  (flop1 'fllog1+
         (let ((coefs (cons 0.0
                            (map fl/
                                 '(1.0 -2.0 3.0 -4.0 5.0
                                   -6.0 7.0 -8.0 9.0 -10.0)))))
           (lambda (x)
  FIXME
                    (polynomial-at x coefs))
                   (else
                    (fllog (fl+ 1.0 x))))))))
           

(define fllog2 (flop1 'fllog2 (lambda (x) (log x 2.0))))

(define fllog10 (flop1 'fllog10 (lambda (x) (log x 10.0))))

(define (make-fllog-base base)
  (check-flonum! 'make-fllog-base base)
  (if (fl>? base 1.0)
      (flop1 'procedure-created-by-make-fllog-base
             (lambda (x) (log x base)))
      (error "argument to make-fllog-base must be greater than 1.0" base)))

;;; Trigonometric functions

;(define flsin R6RS)                 ; defined by (rnrs arithmetic flonums)
;(define flcos R6RS)                 ; defined by (rnrs arithmetic flonums)
;(define fltan R6RS)                 ; defined by (rnrs arithmetic flonums)
;(define flasin R6RS)                ; defined by (rnrs arithmetic flonums)
( define flacos R6RS )                 ; defined by ( rnrs arithmetic flonums )
;(define flatan R6RS)                ; defined by (rnrs arithmetic flonums)

(define flsinh
  (flop1 'flsinh
         (lambda (x)
           (cond ((not (flfinite? x)) x)
                 ((fl<? (flabs x) 0.75)
                  (fl/ (fl- (flexp-1 x) (flexp-1 (fl- x))) 2.0))
                 (else
                  (fl/ (fl- (flexp x) (flexp (fl- x))) 2.0))))))

(define flcosh
  (flop1 'flcosh
         (lambda (x)
           (cond ((not (flfinite? x)) (flabs x))
                 ((fl<? (flabs x) 0.75)
                  (fl+ 1.0 (fl/ (fl+ (flexp-1 x) (flexp-1 (fl- x))) 2.0)))
                 (else
                  (fl/ (fl+ (flexp x) (flexp (fl- x))) 2.0))))))

(define fltanh
  (flop1 'fltanh
         (lambda (x)
           (cond ((flinfinite? x) (flcopysign 1.0 x))
                 ((flnan? x) x)
                 (else
                  (let ((a (flsinh x))
                        (b (flcosh x)))
                    (cond ((fl=? a b)
                           1.0)
                          ((fl=? a (fl- b))
                           -1.0)
                          (else
                           (fl/ (flsinh x) (flcosh x))))))))))

;;; inverse hyperbolic functions

(define flasinh
  (flop1 'flasinh
         (lambda (x)
           (define (eqn4.6.31 x^2 c k k0)
  FIXME
                 c
                 (fl+ c
                      (fl* x^2
                           (eqn4.6.31 x^2
                                      (fl- (fl/ (fl* c k0 k0)
                                                (fl* (fl+ k0 1.0)
                                                     (fl+ k0 2.0))))
                                      (+ k 2)
                                      (fl+ k0 2.0))))))
           (cond ((flzero? x) x)
                 ((not (flfinite? x)) x)
                 ((fl<? x 0.0)
                  (fl- (flasinh (fl- x))))
                 ((fl<? x 0.5)
                  (fl* x (eqn4.6.31 (fl* x x) 1.0 1 1.0)))
                 (else
                  (let* ((x^2+1 (fl+ (fl* x x) 1.0))
                         (root (if (flfinite? x^2+1)
                                   (flsqrt x^2+1)
                                   x))
                         (a (fl+ x root)))
                    (if (flfinite? a)
                        (fllog a)
                        (fl+ fl-log-2 (fllog x)))))))))

(define flacosh
  (flop1 'flacosh
         (lambda (x)
           (if (fl<=? x 0.5)
               +nan.0
               (let* ((x^2-1 (fl- (fl* x x) 1.0))
                      (root (if (flfinite? x^2-1)
                                (flsqrt x^2-1)
                                x))
                      (a (fl+ x root)))
                 (if (flfinite? a)
                     (fllog a)
                     (fl+ fl-log-2 (fllog x))))))))

(define flatanh
  (flop1 'flatanh
         (lambda (x)
           (define (eqn4.6.33 x^2 k k0)
  FIXME
                 (fl/ k0)
                 (fl+ (fl/ k0)
                      (fl* x^2
                           (eqn4.6.33 x^2 (+ k 2) (+ k0 2.0))))))
           (cond ((flzero? x) x)
                 ((not (flfinite? x)) x)
                 ((fl<? x 0.0)
                  (fl- (flatanh (fl- x))))
                 ((fl<? x 0.5)
                  (fl* x (eqn4.6.33 (fl* x x) 1 1.0)))
                 (else
                  (fl* 0.5 (fllog (fl/ (fl+ 1.0 x) (fl- 1.0 x)))))))))

  Integer division

(define flquotient
  (flop2 'flquotient
         (lambda (x y)
           (fltruncate (fl/ x y)))))

;;; FIXME: should probably implement the following part of the C spec:
;;; "If the returned value is 0, it will have the same sign as x."

(define flremainder
  (flop2 'flremainder
         (lambda (x y)
           (fl- x (fl* y (flquotient x y))))))

(define (flremquo x y)
  (check-flonum! 'flremquo x)
  (check-flonum! 'flremquo y)
  (let* ((quo (flround (fl/ x y)))
         (rem (fl- x (fl* y quo))))
    (values rem
            (exact quo))))

  Special functions are defined in 144.special.scm

  eof

 
 Permission is hereby granted, free of charge, to any person
 obtaining a copy of this software and associated documentation
 restriction, including without limitation the rights to use,
 copy, modify, merge, publish, distribute, sublicense, and/or
 Software is furnished to do so, subject to the following
 conditions:
 
 The above copyright notice and this permission notice shall be
 
 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. 
 References

 Handbook of Mathematical Functions With Formulas, Graphs, and



 I have deliberately avoided recent references, and have also
 avoided looking at any code or pseudocode for these or similar
 functions.

 FIXME: not as accurate as it should be
 FIXME: not as fast as it should be
 FIXME: assumes all inexact reals are flonums
 FIXME: assumes (scheme inexact)

 The mathematical constants are now defined in 144.constants.scm
 Implementation Constants
 arbitrary
 arbitrary
 Constructors

 Predicates
(define flonum? R6RS)               ; defined by (rnrs arithmetic flonums)
(define fl=? R6RS)                  ; defined by (rnrs arithmetic flonums)
(define fl<? R6RS)                  ; defined by (rnrs arithmetic flonums)
(define fl>? R6RS)                  ; defined by (rnrs arithmetic flonums)
(define fl<=? R6RS)                 ; defined by (rnrs arithmetic flonums)
(define fl>=? R6RS)                 ; defined by (rnrs arithmetic flonums)
 spec says fl-least, but that's wrong
(define flinteger? R6RS)            ; defined by (rnrs arithmetic flonums)
(define flzero? R6RS)               ; defined by (rnrs arithmetic flonums)
 defined by ( rnrs arithmetic flonums )
(define flnegative? R6RS)           ; defined by (rnrs arithmetic flonums)
 defined by ( rnrs arithmetic flonums )
(define fleven? R6RS)               ; defined by (rnrs arithmetic flonums)
 defined by ( rnrs arithmetic flonums )
(define flinfinite? R6RS)           ; defined by (rnrs arithmetic flonums)
 defined by ( rnrs arithmetic flonums )
 Arithmetic
(define fl+ R6RS)                   ; defined by (rnrs arithmetic flonums)
(define fl* R6RS)                   ; defined by (rnrs arithmetic flonums)
(define fl- R6RS)                   ; defined by (rnrs arithmetic flonums)
(define fl/ R6RS)                   ; defined by (rnrs arithmetic flonums)
(define flabs R6RS)                 ; defined by (rnrs arithmetic flonums)
(define flfloor R6RS)               ; defined by (rnrs arithmetic flonums)
(define flceiling R6RS)             ; defined by (rnrs arithmetic flonums)
(define flround R6RS)               ; defined by (rnrs arithmetic flonums)
(define fltruncate R6RS)            ; defined by (rnrs arithmetic flonums)
 Exponents and logarithms
(define flexp R6RS)                 ; defined by (rnrs arithmetic flonums)

 to be increased for higher precisions.
 defined by ( rnrs arithmetic flonums )
(define flexpt R6RS)                ; defined by (rnrs arithmetic flonums)
(define fllog R6RS)                 ; defined by (rnrs arithmetic flonums)

 log (x + 1) = \sum_{n=1}^\infty - (-1)^n x^n/n
 Trigonometric functions
(define flsin R6RS)                 ; defined by (rnrs arithmetic flonums)
(define flcos R6RS)                 ; defined by (rnrs arithmetic flonums)
(define fltan R6RS)                 ; defined by (rnrs arithmetic flonums)
(define flasin R6RS)                ; defined by (rnrs arithmetic flonums)
 defined by ( rnrs arithmetic flonums )
(define flatan R6RS)                ; defined by (rnrs arithmetic flonums)
 inverse hyperbolic functions
 FIXME: should probably implement the following part of the C spec:
 "If the returned value is 0, it will have the same sign as x."
  Copyright ( C )    ( 2016 ) .
  files ( the " Software " ) , to deal in the Software without
  sell copies of the Software , and to permit persons to whom the
  included in all copies or substantial portions of the Software .
  THE SOFTWARE IS PROVIDED " AS IS " , WITHOUT WARRANTY OF ANY KIND ,

    and    [ editors ] .
  Mathematical Tables .   United States Department of Commerce .
  National Bureau of Standards Applied Mathematics Series , 55 ,
  June 1964 .   Fifth Printing , August 1966 , with corrections .
  R W Hamming .   Numerical Methods for Scientists and Engineers .
  McGraw - Hill , 1962 .
  Donald E Knuth .   The Art of Computer Programming , Volume 2 ,
  Seminumerical Algorithms , Second Edition .   Addison - Wesley , 1981 .
  J N Newman .   Approximations for the  and   .
  Mathematics of Computation , 43(168 ) , October 1984 , pages 551 - 556 .


  Quick - and - dirty implementation of a draft of SRFI 144 ( flonums ) ,
  as specified at /~cowan/temp/srfi-144.html
  as of 4 June 2017 .
  FIXME : assumes IEEE arithmetic or similar

  Mathematical Constants



(define fl-greatest
  (let loop ((x (- (expt 2.0 precision-bits) 1.0)))
    (if (finite? (* 2.0 x))
        (loop (* 2.0 x))
        x)))

(define fl-least
  (let loop ((x 1.0))
    (if (> (/ x 2.0) 0.0)
        (loop (/ x 2.0))
        x)))

(define fl-epsilon
  (let loop ((eps 1.0))
    (if (= 1.0 (+ 1.0 eps))
        (* 2.0 eps)
        (loop (/ eps 2.0)))))

  (exact (- (log fl-least 2.0) 1.0)))

  (- fl-integer-exponent-zero 1))


(define (flonum x)
  (if (real? x)
      (inexact x)
      (error "bad argument passed to flonum" x)))

(define fladjacent
  (flop2 'fladjacent
         (lambda (x y)
           (define (loop y)
             (let* ((y3 (fl+ (fl* 0.999755859375 x) (fl* 0.000244140625 y))))
               (cond ((fl<? x y3 y)
                      (loop y3))
                     ((fl<? y y3 x)
                      (loop y3))
                     (else
                      (loop2 y)))))
           (define (loop2 y)
             (let* ((y2 (fl/ (fl+ x y) 2.0))
                    (y2 (if (flinfinite? y2)
                            (fl+ (fl* 0.5 x) (fl* 0.5 y))
                            y2)))
               (cond ((fl=? x y2)
                      y)
                     ((fl=? y y2)
                      y)
                     (else
                      (loop2 y2)))))
           (cond ((flinfinite? x)
                  (cond ((fl<? x y) (fl- fl-greatest))
                        ((fl>? x y) fl-greatest)
                        (else x)))
                 ((fl=? x y)
                  x)
                 ((flzero? x)
                  (if (flpositive? y)
                      fl-least
                      (fl- fl-least)))
                 ((fl<? x y)
                  (loop (flmin y
                               fl-greatest
                               (flmax (* 2.0 x)
                                      (* 0.5 x)))))
                 ((fl>? x y)
                  (loop (flmax y
                               (fl- fl-greatest)
                               (flmin (* 2.0 x)
                                      (* 0.5 x)))))
  x or y is a NaN
                  x)))))

(define flcopysign
  (flop2 'flcopysign
         (lambda (x y)
           (if (= (flsign-bit x) (flsign-bit y))
               x
               (fl- x)))))

(define (make-flonum x n)
  (let ((y (expt 2.0 n)))
    (cond ((or (not (flonum? x))
               (not (exact-integer? n)))
           (error "bad arguments to make-flonum" x n))
          ((finite? y)
           (* x y))
          (else
           (inexact (* (exact x) (expt 2 n)))))))


  Accessors

(define (flinteger-fraction x)
  (check-flonum! 'flinteger-fraction x)
  (let* ((result1 (fltruncate x))
         (result2 (fl- x result1)))
    (values result1 result2)))

(define (flexponent x)
  (floor (fllog2 (flabs x))))

(define (flinteger-exponent x)
  (exact (flexponent x)))

(define (flnormalized-fraction-exponent x)
  (define (return result1 result2)
    (cond ((fl<? result1 0.5)
           (values (fl* 2.0 result1) (- result2 1)))
          ((fl>=? result1 1.0)
           (values (fl* 0.5 result1) (+ result2 1)))
          (else
           (values result1 result2))))
  (check-flonum! 'flnormalized-fraction-exponent x)
  unspecified for NaN
         (values x 0))
        ((fl<? x 0.0)
         (call-with-values
          (lambda () (flnormalized-fraction-exponent (fl- x)))
          (lambda (y n)
            (values (fl- y) n))))
  unspecified for 0.0
         (values 0.0 0))
        ((flinfinite? x)
         (values 0.5 (+ 3 (exact (round (fllog2 fl-greatest))))))
        ((flnormalized? x)
         (let* ((result2 (exact (flround (fllog2 x))))
                (result2 (if (integer? result2)
                             result2
                             (round result2)))
                (two^result2 (inexact (expt 2.0 result2))))
           (if (flinfinite? two^result2)
               (call-with-values
                (lambda () (flnormalized-fraction-exponent (fl/ x 4.0)))
                (lambda (y n)
                  (values y (+ n 2))))
               (return (fl/ x two^result2) result2))))
        (else
         (let* ((k (+ 2 precision-bits))
                (two^k (expt 2 k)))
           (call-with-values
            (lambda ()
              (flnormalized-fraction-exponent (fl* x (inexact two^k))))
            (lambda (y n)
              (return y (- n k))))))))

(define (flsign-bit x)
  (check-flonum! 'flsign-bit x)
  (cond ((fl<? x 0.0)
         1)
        ((eqv? x -0.0)
         1)
        (else
         0)))




(define (flunordered? x y)
  (or (flnan? x) (flnan? y)))

  incompatible with ( rnrs arithmetic flonums ) in zero - argument case

(define flmax
  (let ((flmax2 (flop2 'flmax max)))
    (lambda args
      (cond ((null? args)
             -inf.0)
            ((null? (cdr args))
             (car args))
            ((null? (cddr args))
             (flmax2 (car args) (cadr args)))
            (else
             (flmax2 (flmax2 (car args) (cadr args))
                     (apply flmax (cddr args))))))))

  incompatible with ( rnrs arithmetic flonums ) in zero - argument case

(define flmin
  (let ((flmin2 (flop2 'flmin min)))
    (lambda args
      (cond ((null? args)
            ((null? (cdr args))
             (car args))
            ((null? (cddr args))
             (flmin2 (car args) (cadr args)))
            (else
             (flmin2 (flmin2 (car args) (cadr args))
                     (apply flmin (cddr args))))))))


(define flnormalized?
  (lambda (x)
    (check-flonum! 'flnormalized? x)
    (let ((x (flabs x)))
      (and (flfinite? x)
           (fl<? (fl/ fl-greatest) x)))))

(define fldenormalized?
  (lambda (x)
    (check-flonum! 'fldenormalized? x)
    (let ((x (flabs x)))
      (and (flfinite? x)
           (fl<? 0.0 x)
           (fl<=? x (fl/ fl-greatest))))))



  Spec says " as if to infinite precision and rounded only once " .

(define fl+*
  (flop3 'fl+*
         (lambda (x y z)
           (cond (c-functions-are-available
                  (fma x y z))
                 ((and (flfinite? x) (flfinite? y))
                  (if (flfinite? z)
                      (let ((x (exact x))
                            (y (exact y))
                            (z (exact z)))
                        (flonum (+ (* x y) z)))
                      z))
                 (else
                  (fl+ (fl* x y) z))))))


(define (flabsdiff x y)
  (flabs (fl- x y)))

(define (flposdiff x y)
  (let ((diff (fl- x y)))
    (if (flnegative? diff)
        0.0
        diff)))

(define (flsgn x)
  (flcopysign 1.0 x))

  ( flnumerator + nan.0 ) and ( fldenominator + nan.0 ) must be NaNs , which
  is not required by the R6RS specification of ( rnrs arithmetic flonums ) .

(define flnumerator
  (flop1 'flnumerator
         (lambda (x)
           (if (flnan? x)
               x
               (r6rs:flnumerator x)))))

(define fldenominator
  (flop1 'fldenominator
         (lambda (x)
           (if (flnan? x)
               x
               (r6rs:fldenominator x)))))




(define flexp2 (flop1 'flexp2 (lambda (x) (flexpt 2.0 x))))

  e^x = \sum_n ( z^n / ( n ! ) )
  FIXME : the number of terms and the constant 0.5 seem reasonable
  for IEEE double precision , but the number of terms might need

(define flexp-1
  (flop1 'flexp-1
         (let ((coefs (cons 0.0
                            (map fl/
                                 (map factorial
                                      '(1.0 2.0 3.0 4.0 5.0
                                        6.0 7.0 8.0 9.0 10.0
                                        11.0 12.0 13.0 14.0 15.0))))))
           (lambda (x)
  FIXME
                    (polynomial-at x coefs))
                   (else
                    (fl- (flexp x) 1.0)))))))

(define flsquare (flop1 'flsquare (lambda (x) (fl* x x))))


(define flcbrt
  (flop1 'flcbrt
         (lambda (x)
           (cond ((flnegative? x)
                  (fl- (flcbrt (fl- x))))
                 (else
                  (flexpt x (fl/ 3.0)))))))

(define flhypot
  (flop2 'flhypot
         (lambda (x y)
           (cond ((flzero? x) (flabs y))
                 ((flzero? y) (flabs x))
                 ((or (flinfinite? x) (flinfinite? y)) +inf.0)
                 ((flnan? x) x)
                 ((flnan? y) y)
                 ((fl>? y x) (flhypot y x))
                 (else
                  (let* ((y/x (fl/ y x))
                         (root (flsqrt (fl+ 1.0 (fl* y/x y/x)))))
                    (fl* (flabs x) root)))))))


    ) , as in C99 log1p .

(define fllog1+
  (flop1 'fllog1+
         (let ((coefs (cons 0.0
                            (map fl/
                                 '(1.0 -2.0 3.0 -4.0 5.0
                                   -6.0 7.0 -8.0 9.0 -10.0)))))
           (lambda (x)
  FIXME
                    (polynomial-at x coefs))
                   (else
                    (fllog (fl+ 1.0 x))))))))
           

(define fllog2 (flop1 'fllog2 (lambda (x) (log x 2.0))))

(define fllog10 (flop1 'fllog10 (lambda (x) (log x 10.0))))

(define (make-fllog-base base)
  (check-flonum! 'make-fllog-base base)
  (if (fl>? base 1.0)
      (flop1 'procedure-created-by-make-fllog-base
             (lambda (x) (log x base)))
      (error "argument to make-fllog-base must be greater than 1.0" base)))



(define flsinh
  (flop1 'flsinh
         (lambda (x)
           (cond ((not (flfinite? x)) x)
                 ((fl<? (flabs x) 0.75)
                  (fl/ (fl- (flexp-1 x) (flexp-1 (fl- x))) 2.0))
                 (else
                  (fl/ (fl- (flexp x) (flexp (fl- x))) 2.0))))))

(define flcosh
  (flop1 'flcosh
         (lambda (x)
           (cond ((not (flfinite? x)) (flabs x))
                 ((fl<? (flabs x) 0.75)
                  (fl+ 1.0 (fl/ (fl+ (flexp-1 x) (flexp-1 (fl- x))) 2.0)))
                 (else
                  (fl/ (fl+ (flexp x) (flexp (fl- x))) 2.0))))))

(define fltanh
  (flop1 'fltanh
         (lambda (x)
           (cond ((flinfinite? x) (flcopysign 1.0 x))
                 ((flnan? x) x)
                 (else
                  (let ((a (flsinh x))
                        (b (flcosh x)))
                    (cond ((fl=? a b)
                           1.0)
                          ((fl=? a (fl- b))
                           -1.0)
                          (else
                           (fl/ (flsinh x) (flcosh x))))))))))


(define flasinh
  (flop1 'flasinh
         (lambda (x)
           (define (eqn4.6.31 x^2 c k k0)
  FIXME
                 c
                 (fl+ c
                      (fl* x^2
                           (eqn4.6.31 x^2
                                      (fl- (fl/ (fl* c k0 k0)
                                                (fl* (fl+ k0 1.0)
                                                     (fl+ k0 2.0))))
                                      (+ k 2)
                                      (fl+ k0 2.0))))))
           (cond ((flzero? x) x)
                 ((not (flfinite? x)) x)
                 ((fl<? x 0.0)
                  (fl- (flasinh (fl- x))))
                 ((fl<? x 0.5)
                  (fl* x (eqn4.6.31 (fl* x x) 1.0 1 1.0)))
                 (else
                  (let* ((x^2+1 (fl+ (fl* x x) 1.0))
                         (root (if (flfinite? x^2+1)
                                   (flsqrt x^2+1)
                                   x))
                         (a (fl+ x root)))
                    (if (flfinite? a)
                        (fllog a)
                        (fl+ fl-log-2 (fllog x)))))))))

(define flacosh
  (flop1 'flacosh
         (lambda (x)
           (if (fl<=? x 0.5)
               +nan.0
               (let* ((x^2-1 (fl- (fl* x x) 1.0))
                      (root (if (flfinite? x^2-1)
                                (flsqrt x^2-1)
                                x))
                      (a (fl+ x root)))
                 (if (flfinite? a)
                     (fllog a)
                     (fl+ fl-log-2 (fllog x))))))))

(define flatanh
  (flop1 'flatanh
         (lambda (x)
           (define (eqn4.6.33 x^2 k k0)
  FIXME
                 (fl/ k0)
                 (fl+ (fl/ k0)
                      (fl* x^2
                           (eqn4.6.33 x^2 (+ k 2) (+ k0 2.0))))))
           (cond ((flzero? x) x)
                 ((not (flfinite? x)) x)
                 ((fl<? x 0.0)
                  (fl- (flatanh (fl- x))))
                 ((fl<? x 0.5)
                  (fl* x (eqn4.6.33 (fl* x x) 1 1.0)))
                 (else
                  (fl* 0.5 (fllog (fl/ (fl+ 1.0 x) (fl- 1.0 x)))))))))

  Integer division

(define flquotient
  (flop2 'flquotient
         (lambda (x y)
           (fltruncate (fl/ x y)))))


(define flremainder
  (flop2 'flremainder
         (lambda (x y)
           (fl- x (fl* y (flquotient x y))))))

(define (flremquo x y)
  (check-flonum! 'flremquo x)
  (check-flonum! 'flremquo y)
  (let* ((quo (flround (fl/ x y)))
         (rem (fl- x (fl* y quo))))
    (values rem
            (exact quo))))

  Special functions are defined in 144.special.scm

  eof
