;;; Copyright (C) William D Clinger (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
;;;
;;; Milton Abramowitz and Irene A Stegun [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 Bessel and Struve Functions.
;;; 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 http://vrici.lojban.org/~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

(foreign-declare "#include \"float.h\"")

;; Implementation Constants

(define fl-greatest maximum-flonum)

(define fl-least (foreign-value "DBL_TRUE_MIN" double))

(define fl-epsilon flonum-epsilon)

(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

; Implements post-finalization note 1
(define (flonum x)
  (if (real? x)
      (inexact x)
      +nan.0))

(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)))))
                 (else    ; 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)
  (cond ((flnan? 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))))
        ((fl=? x 0.0)    ; 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
  (case-lambda
    (() -inf.0)
    ((x) x)
    ((x y) (fpmax x y))
    ((x . rest)
     (fpmax (fpmax x (car rest))
            (apply flmax (cdr rest))))))

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

(define flmin
  (case-lambda
    (() +inf.0)
    ((x) x)
    ((x y) (fpmin x y))
    ((x . rest)
     (fpmin (fpmin x (car rest))
            (apply flmin (cdr rest))))))

;(define flinteger? R6RS)            ; defined by (rnrs arithmetic flonums)
;(define flzero? R6RS)               ; defined by (rnrs arithmetic flonums)
;(define flpositive? 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+* fma))

;(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)
             (cond ((fl<? (flabs x) 0.5)    ; 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)

;;; Returns log(x+1), as in C99 log1p.

;;; See

;;; https://stat.ethz.ch/pipermail/r-devel/2003-August/027396.html
;;; https://books.google.com/books?id=OjUyDwAAQBAJ&pg=PA290&lpg=PA290&dq=beebe+log1p&source=bl&ots=VLxmiSk1fA&sig=ACfU3U0_8tqKemomSjKW73iJ0zUO1u3p3Q&hl=en&sa=X&ved=2ahUKEwjfxZbE8LvhAhVNm-AKHWScB7w4ChDoATAAegQICRAB#v=onepage&q=beebe%20log1p&f=false

;;; for justification


(define fllog1+
  (flop1 'fllog1+
         (lambda (x)
           (let ((u (fl+ 1.0 x)))
             (cond ((fl=? u 1.0)
                    x) ;; gets sign of zero result correct
                   ((fl=? u x)
                    (fllog u)) ;; large arguments and infinities
                   (else
                    (fl* (fllog u) (fl/ x (fl- u 1.0)))))))))


(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)
           (cond ((or (flinfinite? x)
                      (flnan? x))
                  x)
                 ((flnegative? x)
                  (fl- (flasinh (fl- x))))
                 ((fl<? x 3.725290298461914e-9)   ;; (flexpt 2. -28.)
                  x)
                 ((fl<? x 2.)
                  ;; the naive formula is
                  ;; (log (+ x (sqrt (+ (square x) 1))))

                  ;; We want to
                  ;; 1.  Use exact operations when possible (no roundoff)
                  ;; 2.  Add or subtract things of differing magnitudes,
                  ;;     so for most arguments at most one roundoff error.

                  ;; Biggeset possible problem near x=0, so we write
                  ;; (sqrt (+ 1 (square x)))
                  ;; as
                  ;; (+ 1 (- (sqrt (+ 1 (square x))) 1))
                  ;; and then multiply the second part in
                  ;; numerator and denominator by
                  ;; (+ (sqrt (+ 1 (square x))) 1)

                  (let ((x^2 (flsquare x)))
                    (fllog1+ (fl+ x
                                  (fl/ x^2
                                       (fl+ 1.
                                            (flsqrt (fl+ 1.0 x^2))))))))
                 ((fl<? x 268435456.) ;; (flexpt 2. 28.)
                  (let ((x^2 (flsquare x)))
                    (fllog (fl+ (fl* 2. x) ;; exact
                                ;; the rest is small
                                (fl/ 1.
                                     (fl+ x
                                          (flsqrt (fl+ 1.0 x^2))))))))
                 (else
                  (fl+ (fllog x) fl-log-2))))))

(define flacosh
  (flop1 'flacosh
         (lambda (x)
           (cond ((flnan? x) x)
                 ((fl<? x 1.0) +nan.0)
                 ((fl<? x 2.0)
                  ;; the naive formula is
                  ;; (log (+ x (sqrt (- (square x) 1))))

                  ;; We want to
                  ;; 1.  Use exact operations when possible (no roundoff)
                  ;; 2.  Add or subtract things of differing magnitudes,
                  ;;     so for most arguments at most one roundoff error.

                  (let ((x-1 (fl- x 1.))) ;; exact
                    (fllog1+ (fl+ x-1 ;; smaller than next expression
                                  (flsqrt (fl+ (fl* 2. x-1) ;; exact
                                               ;; relatively small
                                               (flsquare x-1)))))))
                 ((fl<? x 268435456.) ;; (flexpt 2. 28.)
                  (fllog (fl- (fl* 2. x) ;; exact
                              ;; next is smaller
                              (fl/ (fl+ x (flsqrt (fl* (fl- x 1.) ;; exact
                                                       (fl+ x 1.) ;; exact
                                                       )))))))
                 (else
                  (fl+ (fllog x) fl-log-2))))))

(define flatanh
  (flop1 'flatanh
         (lambda (x)
           (cond ((fl<? x 0.)
                  (fl- (flatanh (fl- x))))
                 ;; we rewrite
                 ;; (/ (+ 1 x) (- 1 x))
                 ;; as
                 ;; (+ 1 (* 2 (/ x (- 1 x))))
                 ;; and call fllog1+ instead of fllog
                 (else
                  (fl* +0.5                                    ;; exact
                       (fllog1+ (fl* +2.0                      ;; exact
                                     (fl/ x
                                          (fl- 1.0 x)))))))))) ;; exact

;;; 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