#lang scheme
;;; PLT Scheme Science Collection
;;; math.ss
;;; Copyright (c) 2004-2008 M. Douglas Williams
;;;
;;; This library is free software; you can redistribute it and/or
;;; modify it under the terms of the GNU Lesser General Public
;;; License as published by the Free Software Foundation; either
;;; version 2.1 of the License, or (at your option) any later version.
;;;
;;; This library is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;;; Lesser General Public License for more details.
;;;
;;; You should have received a copy of the GNU Lesser General Public
;;; License along with this library; if not, write to the Free
;;; Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
;;; 02111-1307 USA.
;;;
;;; ----------------------------------------------------------------------------
;;;
;;; This code is based on the Mathematical Functions in the GNU
;;; Scientific Library (GSL).
;;;
;;; Version  Date      Description
;;; 0.1.0    08/13/04  This is the initial realease of the mathematical function
;;;                    routines ported from GSL. (Doug Williams)
;;; 1.0.0    09/28/04  Removed the code for small integer powers. Marked as ready
;;;                    for Release 1.0. Added contracts for functions (Doug
;;;                    Williams)
;;; 2.0.0    11/17/07  Added unchecked versions and preparing for PLT Scheme
;;;                    V4.0. (Doug Williams)
;;; 3.0.0    06/09/08  Changes required for V4.0. (Doug Williams)
;;; 4.0.0    10/04/09  Added use of unsafe operations. Added real->float to
;;;                    protect unsafe operations. Added log10 and dB. Removed
;;;                    the unchecked operations. (Doug Williams)

(require scheme/unsafe/ops)
(require "machine.ss")

;;; -----------------------------------------------------------------
;;;
;;; Mathematical Constants
;;;
;;; The library ensures that the standard BSD mathematical constants
;;; are defined.

;;; The base of exponentials, e
(define e 2.71828182845904523536028747135)

;;; The base-2 logarithm of e, log_2(e)
(define log2e 1.44269504088896340735992468100)

;;; The base-10 logarithm of e, log_10(e)
(define log10e 0.43429448190325182765112891892)

;;; The square root of two, sqrt(2)
(define sqrt2 1.41421356237309504880168872421)

;;; The square root of one half, sqrt(1/2)
(define sqrt1/2 0.70710678118654752440084436210)

;;; The square root of three, sqrt(3)
(define sqrt3 1.73205080756887729352744634151)

;;; The constants pi, 2*pi, and 4*pi
(define pi 3.14159265358979323846264338328)
(define 2*pi (* 2.0 pi))
(define 4*pi (* 4.0 pi))

;;; Pi divided by two, pi/2
(define pi/2 1.57079632679489661923132169164)

;;; Pi divided by four, pi/4
(define pi/4 0.78539816339744830966156608458)

;;; The square root of pi, sqrt(pi)
(define sqrtpi 1.77245385090551602729816748334)

;;; Two divided by the square root of pi, 2/sqrt(pi)
(define 2/sqrtpi 1.12837916709551257389615890312)

;;; The reciprocal of pi, 1/pi
(define 1/pi 0.31830988618379067153776752675)

;;; Twice the reciprocal of pi, 2/pi
(define 2/pi 0.63661977236758134307553505349)

;;; The natural logarithm of ten, ln(10), and its inverse
(define ln10 2.30258509299404568401799145468)
(define 1/ln10 (/ ln10))

;;; The natural logarithm of two, ln(2), and its inverse
(define ln2 0.69314718055994530941723212146)
(define 1/ln2 (/ ln2))

;;; The natural logarithm of pi, ln(pi)
(define lnpi 1.14472988584940017414342735135)

;;; Euler's constant, gamma
(define euler 0.57721566490153286060651209008)

;;; -----------------------------------------------------------------
;;;
;;; Infinities and Not-a-Number
;;;
;;; PLT Scheme provides +inf.0 (infinity), -inf.0 (negative
;;; infinity), +nan.0 (not a number) and -nan.0 (same as +nan.0) as
;;; inexact numerical constants.

;;; (nan? x) -> boolean?
;;;  any/c
;;; Returns true if x is equivalent to +nan.0.
(define (nan? x)
  (eqv? x +nan.0))

;;; (infinity? x) -> (or/c -1 1 #f)
;;;   x : any/c
;;; Returns 1 if x is +inf.0, -1 if x is -inf.0, and false otherwise.
(define (infinite? x)
  (if (eqv? x +inf.0)
      1
      (if (eqv? x -inf.0)
          -1
          #f)))

;;; (finite? x) -> boolean?
;;;   x : real?
;;; Returns true of x is a finite, real number.
(define (finite? x)
  (and (real? x)
       (not (nan? x))
       (not (infinite? x))))

;;; -----------------------------------------------------------------------------
;;;
;;; Elementary Functions

;;; (log1p x) -> number?
;;;   x : real?
;;; Computes the value of log(1+x) is a way that is accurate for small x.
(define (log1p x)
  (with-float (x)
    (let ((y (unsafe-fl+ 1.0 x)))
      (- (log y)
         (unsafe-fl/ (unsafe-fl- (unsafe-fl- y 1.0) x) y)))))

;;; (expm1 x) -> inexact-real?
;;;   x : real?
;;; Computes the value of exp(x)-1 in a way that is accurate for
;;; small x.
(define (expm1 x)
  (with-float (x)
    (if (unsafe-fl< (unsafe-flabs x) ln2)
        ;; Compute the taylor series S = x + (1/2!) x^2 + (1/3!) x^3 +
        ;; ...
        (let ((i 1.0)
              (sum x)
              (term (unsafe-fl/ x 1.0)))
          (let loop ()
            (set! i (unsafe-fl+ i 1.0))
            (set! term (unsafe-fl* term (unsafe-fl/ x i)))
            (set! sum (unsafe-fl+ sum term))
            (when (unsafe-fl> (unsafe-flabs term)
                              (unsafe-fl* (unsafe-flabs sum) double-epsilon))
              (loop)))
          sum)
        (unsafe-fl- (exp x) 1.0))))

;;; (hypot x y) -> inexact-real?
;;;   x : real?
;;;   y : real?
;;; Computes the value of sqrt(x^2 + y^2) in a way that avoids overflow.
(define (hypot x y)
  (with-float (x y)
    (let ((xabs (abs x))
          (yabs (abs y))
          (min 0.0)
          (max 0.0))
      (if (unsafe-fl< xabs yabs)
          (begin
            (set! min xabs)
            (set! max yabs))
          (begin
            (set! min yabs)
            (set! max xabs)))
      (when (unsafe-fl= min 0.0)
        max)
      (let ((u (unsafe-fl/ min max)))
        (unsafe-fl* max (sqrt (unsafe-fl+ 1.0 (unsafe-fl* u u))))))))

;;; (acosh x) -> inexact-real?
;;;   x : real?
;;; Computes the value of arccosh(x).
(define (acosh x)
  (with-float (x)
    (cond ((unsafe-fl> x (unsafe-fl/ 1.0 sqrt-double-epsilon))
           (unsafe-fl+ (log x) ln2))
          ((unsafe-fl> x 2.0)
           (log (unsafe-fl-
                 (unsafe-fl* 2.0 x)
                 (unsafe-fl/ 1.0
                             (unsafe-fl+ (sqrt (unsafe-fl- (unsafe-fl* x x) 1.0))
                                         x)))))
          ((unsafe-fl> x 1.0)
           (let ((t (unsafe-fl- x 1.0)))
             (log1p (unsafe-fl+ (unsafe-fl* 2.0 t) (unsafe-fl* t t)))))
          ((unsafe-fl= x 1.0)
           0.0)
          (else
           +nan.0))))

;;; (asinh x) -> inexact-real?
;;;   x : real?
;;; Computes the value of arcsinh(x).
(define (asinh x)
  (with-float (x)
    (let ((a (abs x))
          (s (if (unsafe-fl< x 0.0) -1.0 1.0)))
      (cond ((unsafe-fl> a (unsafe-fl/ 1.0 sqrt-double-epsilon))
             (unsafe-fl* s (unsafe-fl+ (log a) ln2)))
            ((unsafe-fl> a 2.0)
             (unsafe-fl* s (log (unsafe-fl+
                                 (unsafe-fl* 2.0 a)
                                 (unsafe-fl/
                                  a
                                  (sqrt (unsafe-fl+ (unsafe-fl* a a) 1.0)))))))
            ((unsafe-fl> a sqrt-double-epsilon)
             (let ((a2 (unsafe-fl* a a)))
               (unsafe-fl*
                s
                (log1p (unsafe-fl+
                        a
                        (unsafe-fl/
                         a2
                         (unsafe-fl+ 1.0 (sqrt (unsafe-fl+ 1.0 a2)))))))))
            (else
             x)))))

;;; (atanh x) -> inexact-real?
;;;   x : real?
;;; Computes the value of arctanh(x).
(define (atanh x)
  (with-float (x)
    (let ((a (abs x))
          (s (if (unsafe-fl< x 0.0) -1.0 1.0)))
      (cond ((unsafe-fl> a 1.0)
             +nan.0)
            ((unsafe-fl= a 1.0)
             (if (unsafe-fl< x 0.0) -inf.0 +inf.0))
            ((unsafe-fl>= x 0.5)
             (unsafe-fl* s
                         (unsafe-fl* 0.5
                                     (log1p (unsafe-fl/
                                             (unsafe-fl* 2.0 a)
                                             (unsafe-fl- 1.0 a))))))
            ((unsafe-fl> a double-epsilon)
             (unsafe-fl* s
                         (unsafe-fl*
                          0.5
                          (log1p (unsafe-fl+ (unsafe-fl* 2.0 a) 
                                             (unsafe-fl/
                                              (unsafe-fl* 2.0
                                                          (unsafe-fl* a a))
                                              (unsafe-fl- 1.0 a)))))))
            (else
             x)))))

;;; (ldexp x e) -> inexact-real?
;;;   x : real?
;;;   r : integer?
;;; Computes the value of x * 2^e.
(define (ldexp x e)
  (with-float (x)
    (let ((p2 (expt 2.0 e)))
      (unsafe-fl* x p2))))

;;; (frexp x) -> inexact-real? integer?
;;; Splits the number x into its normalized fraction f and exponent
;;; e, such that x = f * 2^e and 0.5 <= f < 1.
(define (frexp x)
  (with-float (x)
    (if (unsafe-fl= x 0.0)
        (values
         0.0
         0)
        (let* ((ex (ceiling (log (unsafe-fl/ (unsafe-flabs x) ln2))))
               (ei (inexact->exact ex))
               (f (ldexp x (- ei))))
          (let loop ()
            (when (unsafe-fl>= (unsafe-flabs f) 1.0)
              (set! ei (unsafe-fx+ ei 1))
              (set! f (unsafe-fl/ f 2.0))
              (loop)))
          (let loop ()
            (when (unsafe-fl< (abs f) 0.5)
              (set! ei (unsafe-fx- ei 1))
              (set! f (unsafe-fl* f 2.0))
              (loop)))
          (values
           f
           ei)))))

;;; -----------------------------------------------------------------
;;; Testing the Sign of Numbers

;;; (sign x) -> (or/c -1 1)
;;; Return the sign of a real number, 1 for positive (or zero),
;;; -1 for negative.
(define (sign x)
  (if (>= x 0) 1 -1))

;;; -----------------------------------------------------------------
;;; Testing form Odd and Even Numbers
;;;
;;; Scheme provides even? and odd?

;;; -----------------------------------------------------------------
;;; Maximimum and Minimum Functions
;;;
;;; Scheme provides min and max

;;; -----------------------------------------------------------------
;;; Approximate Comparison of Floating Point Numbers

;;; (fcmp x1 x2) -> (or/c -1 0 1)
;;;   x1 : real?
;;;   x2 : real?
;;; Compare two real for difference less than the specified epsilon.
;;; Returns -1 if x1 < x2; 0 if x1 ~=  x2; 1 if x1 > x2.
(define (fcmp x1 x2 epsilon)
  (with-float (x1 x2)
    (let ((difference (unsafe-fl- x1 x2))
          (max (if (unsafe-fl> (unsafe-flabs x1) (unsafe-flabs x2)) x1 x2)))
      (let*-values (((dummy exponent) (frexp max))
                    ((delta) (ldexp epsilon exponent)))
        (if (unsafe-fl> difference delta)
            1                             ; x1 > x2
            (if (unsafe-fl< difference (unsafe-fl- 0.0 delta))
                -1                        ; x1 < x2
                0))))))                   ; x1 ~= x2

;;; log10 and dB

;;; (log10 x) -> inexact-real?
;;;   x : real?
;;; Returns the log base 10 of x.
(define (log10 x)
  (with-float (x)
    (unsafe-fl* 1/ln10 (log x))))

;;; (dB x) -> inexact-real?
;;;   x : real?
;;; Returns the value of x in decibels, 10*log10(x).
(define (dB x)
  (with-float (x)
    (unsafe-fl* 10.0 (unsafe-fl* 1/ln10 (log x)))))

;;; Assure floats for unsafe code

;;; (real->float x) -> inexact-real?
;;;   x : real?
;;; Returns an inexact real (i.e., a float) given real x. Raises an error if x
;;; is not a real. This can be used to ensure a real value is a float, even in
;;; unsafe code.
(define (real->float x)
  (if (real? x)
      (exact->inexact x)
      (error "expected real, given" x)))

;;; (real-vector->float-vector v) -> (vectorof inexact-real?)
;;;   v : (vectorof real?)
;;; Returns a vector of inexact reals (i.e., floats) given a vector of reals, v.
;;; Raises an error if an element of v is not a real.
(define (real-vector->float-vector v)
  (build-vector
   (vector-length v)
   (lambda (i)
     (real->float (vector-ref v i)))))

;;; (with-float (x ...)
;;;   expr ...)
;;; Executes the expr's with the x's guaranteed to be floats. All of the x's
;;; must be identifiers.
(define-syntax (with-float stx)
  (syntax-case stx ()
    ((_ (x ...) expr ...)
     (for ((id (in-list (syntax->list #'(x ...)))))
       (unless (identifier? id)
         (raise-syntax-error #f
                             "not an identifier"
                             stx
                             id)))
     #`(let ((x (if (real? x)
                    (exact->inexact x)
                    (error "expected real, given" x)))
             ...)
         expr ...))))

;;; Module Contracts

(provide
 e
 log2e
 log10e
 sqrt2
 sqrt1/2
 sqrt3
 pi
 2*pi
 4*pi
 pi/2
 pi/4
 sqrtpi
 2/sqrtpi
 1/pi
 2/pi
 ln10
 1/ln10
 ln2
 1/ln2
 lnpi
 euler
 real->float
 with-float)

(provide/contract
 (nan?
  (-> any/c boolean?))
 (infinite?
  (-> any/c (or/c -1 #f 1)))
 (finite?
  (-> any/c boolean?))
 (log1p
  (-> real? number?))
 (expm1
  (-> real? inexact-real?))
 (hypot
  (-> real? real? inexact-real?))
 (acosh
  (-> real? inexact-real?))
 (asinh
  (-> real? inexact-real?))
 (atanh
  (-> real? inexact-real?))
 (ldexp
  (-> real? integer? inexact-real?))
 (frexp
  (-> real? (values inexact-real? integer?)))
 (sign
  (-> real? (or/c -1 1)))
 (fcmp
  (-> real? real? real? (or/c -1 0 1)))
 (log10
  (-> real? inexact-real?))
 (dB
  (-> real? inexact-real?))
 (real-vector->float-vector
  (-> (vectorof real?) (vectorof inexact-real?))))