;;; PLT Scheme Science Collection
;;; math.ss
;;; Copyright (c) 2004 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)

(module math mzscheme
  
  (require (lib "contract.ss"))
  
  (provide
   e
   log2e
   log10e
   sqrt2
   sqrt1/2
   sqrt3
   pi
   pi/2
   pi/4
   sqrtpi
   2/sqrtpi
   1/pi
   2/pi
   ln10
   ln2
   lnpi
   euler)
  
  (provide/contract
   (nan?
    (-> any/c boolean?))
   (infinite?
    (-> any/c (union (integer-in -1 1)
                    boolean?)))
   (finite?
    (-> any/c boolean?))
   (log1p
    (-> real? number?))
   (expm1
    (-> real? real?))
   (hypot
    (-> real? real? real?))
   (acosh
    (-> real? real?))
   (asinh
    (-> real? real?))
   (atanh
    (-> real? real?))
   (ldexp
    (-> real? integer? real?))
   (frexp
    (-> real? (values real? integer?)))
   (sign
    (-> real? (integer-in -1 1)))
   (fcmp
    (-> real? real? real? (integer-in -1 1))))
  
  (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 constant pi
  (define pi 3.14159265358979323846264338328)
  
  ;; 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)
  (define ln10 2.30258509299404568401799145468)
  
  ;; The natural logarithm of two, ln(2)
  (define ln2 0.69314718055994530941723212146)
  
  ;; 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?: any -> boolean
  ;;
  ;; Returns true if x is equivalent to +nan.0.
  (define (nan? x)
    (eqv? x +nan.0))
  
  ;; infinite?: any -> boolean or (-1, 1)
  ;;
  ;; 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?: real -> boolean
  ;;
  ;; Returns true of x is a finite, real number.
  (define (finite? x)
    (and (real? x)
         (not (nan? x))
         (not (infinite? x))))
  
  ;; -----------------------------------------------------------------
  ;;
  ;; Elementary Functions
  
  ;; log1p: real -> number (real? for x > 1, complex for x < 1)
  ;;
  ;; Computes the value of log(1+x) is a way that is accurate for
  ;; small x.
  (define (log1p x)
    (let ((y (+ 1.0 x)))
      (- (log y) (/ (- (- y 1.0) x) y))))
  
  ;; expm1: real -> real
  ;;
  ;; Computes the value of exp(x)-1 in a way that is accurate for
  ;; small x.
  (define (expm1 x)
    (if (< (abs x) ln2)
        ;; Compute the taylor series S = x + (1/2!) x^2 + (1/3!) x^3 +
        ;; ...
        (let ((i 1.0)
              (sum x)
              (term (/ x 1.0)))
          (let loop ()
            (set! i (+ i 1.0))
            (set! term (* term (/ x i)))
            (set! sum (+ sum term))
            (if (> (abs term) (* (abs sum) double-epsilon))
                (loop)))
          sum)
        (- (exp x) 1.0)))
  
  ;; hypot: real x real -> real
  ;;
  ;; Computes the value of sqrt(x^2 + y^2) in a way that avoids
  ;; overflow.
  (define (hypot x y)
    (let ((xabs (abs x))
          (yabs (abs y))
          (min 0)
          (max 0))
      (if (< xabs yabs)
          (begin
            (set! min xabs)
            (set! max yabs))
          (begin
            (set! min yabs)
            (set! max xabs)))
      (if (= min 0)
          max)
      (let ((u (/ min max)))
        (* max (sqrt (+ 1.0 (* u u)))))))
  
  ;; acosh: real -> real
  ;;
  ;; Computes the value of arccosh(x)
  (define (acosh x)
    (cond ((> x (/ 1.0 sqrt-double-epsilon))
           (+ (log x) ln2))
          ((> x 2.0)
           (log (- (* 2.0 x)
                   (/ 1.0 (+ (sqrt (- (* x x) 1.0)) x)))))
          ((> x 1.0)
           (let ((t (- x 1.0)))
             (log1p (+ (* 2.0 t) (* t t)))))
          ((= x 1.0)
           0.0)
          (else
           +nan.0)))
  
  ;; asinh: real -> real
  ;;
  ;; Computes the value of arcsinh(x)
  (define (asinh x)
    (let ((a (abs x))
          (s (if (< x 0) -1.0 1.0)))
      (cond ((> a (/ 1.0 sqrt-double-epsilon))
             (* s (+ (log a) ln2)))
            ((> a 2.0)
             (* s (log (+ (* 2.0 a)
                          (/ a (sqrt (+ (* a a) 1.0)))))))
            ((> a sqrt-double-epsilon)
             (let ((a2 (* a a)))
               (* s (log1p (+ a (/ a2 (+ 1.0 (sqrt (+ 1.0 a2)))))))))
            (else
             x))))
  
  ;; atanh: real -> real
  ;;
  ;; Computes the value of arctanh(x).
  (define (atanh x)
    (let ((a (abs x))
          (s (if (< x 0) -1.0 1.0)))
      (cond ((> a 1)
             +nan.0)
            ((= a 1)
             (if (< x 0) -inf.0 +inf.0))
            ((>= x 0.5)
             (* s 0.5 (log1p (/ (* 2.0 a) (- 1.0 a)))))
            ((> a double-epsilon)
             (* s 0.5 (log1p (+ (* 2.0 a) 
                                (/ (* 2.0 a a) (- 1.0 a))))))
            (else
             x))))
  
  ;; ldexp: real x integer -> real
  ;;
  ;; Computes the value of x * 2^e.
  (define (ldexp x e)
    (let ((p2 (expt 2.0 e)))
      (* x p2)))
  
  ;; frexp: real -> real x 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)
    (if (= x 0)
        (values
         0.0
         0)
        (let* ((ex (ceiling (log (/ (abs x) ln2))))
               (ei (inexact->exact ex))
               (f (ldexp x (- ei))))
          (let loop ()
            (if (>= (abs f) 1.0)
                (begin
                  (set! ei (+ ei 1))
                  (set! f (/ f 2.0))
                  (loop))))
          (let loop ()
            (if (< (abs f) 0.5)
                (begin
                  (set! ei (- ei 1))
                  (set! f (* f 2.0))
                  (loop))))
          (values
           f
           ei))))
  
  ;; -----------------------------------------------------------------
  ;; Testing the Sign of Numbers
  
  ;; sign: real -> integer
  ;;
  ;; 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: real x real x real -> integer
  ;;
  ;; 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)
    (let* ((exponent 0)
           (delta 0.0)
           (difference (- x1 x2))
           (dummy 0.0)
           (max (if (> (abs x1) (abs x2)) x1 x2)))
      (set!-values (dummy exponent) (frexp max))
      (set! delta (ldexp epsilon exponent))
      (if (> difference delta)
          1                             ; x1 > x2
          (if (< difference (- delta))
              -1                        ; x1 < x2
              0))))                     ; x1 ~= x2
  
)