#lang scheme/base
;;; 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)

(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
 (rename-out (nan? unchecked-nan?)
             (infinite? unchecked-infinite?)
             (finite? unchecked-finite?)
             (log1p unchecked-log1p)
             (expm1 unchecked-expm1)
             (hypot unchecked-hypot)
             (acosh unchecked-acosh)
             (asinh unchecked-asinh)
             (atanh unchecked-atanh)
             (ldexp unchecked-ldexp)
             (frexp unchecked-frexp)
             (sign unchecked-sign)
             (fcmp unchecked-fcmp)))

(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))
          (when (> (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)))
    (when (= 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 ()
          (when (>= (abs f) 1.0)
            (set! ei (+ ei 1))
            (set! f (/ f 2.0))
            (loop)))
        (let loop ()
          (when (< (abs f) 0.5)
            (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