mathematical.lisp 6.38 KB
Newer Older
1 2
;; Mathematical functions
;; Liam Healy, Wed Mar  8 2006 - 22:09
3
;; Time-stamp: <2008-12-31 18:44:26EST mathematical.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

8
;;; Disable floating point traps for each CL implementation. 
9

10 11 12 13 14 15
;;; Not ported because they are all C macros or inline functions:
;;;   Mathematical Constants
;;;   Testing the Sign of Numbers
;;;   Testing for Odd and Even Numbers
;;;   Maximum and Minimum functions
;;; Does CL need the small integer powers?
16

17 18
#+clisp (setf sys::*inhibit-floating-point-underflow* t)

19 20 21 22 23 24 25 26 27 28
;;;;****************************************************************************
;;; Infinities and Not-a-number
;;;;****************************************************************************

(defmacro pmnil (x)
  "+1, -1, or nil"
  `(let ((v ,x))
     (when (or (= 1 v) (= -1 v))
       v)))

29 30 31 32 33 34
#+sbcl
(eval-when (:compile-toplevel :load-toplevel :execute)
  (sb-int:set-floating-point-modes :traps nil)
  (import
   '(sb-ext:double-float-negative-infinity
     sb-ext:double-float-positive-infinity)))
35

36 37
(export '(+nan+ +positive-infinity+ +negative-infinity+))

38
(defconstant +nan+
39 40 41
  (ignore-errors
    (cffi:foreign-funcall "gsl_nan" :double)))

42
(defconstant +positive-infinity+
43 44 45
  (ignore-errors
    (cffi:foreign-funcall "gsl_posinf" :double)))

46
(defconstant +negative-infinity+
47 48 49
  (ignore-errors
    (cffi:foreign-funcall "gsl_neginf" :double)))

liam's avatar
liam committed
50
(defmfun nanp (x)
51 52
  "gsl_isnan" ((x :double))
  :c-return (cr :int)
liam's avatar
liam committed
53 54 55
  :return ((= 1 cr))
  :documentation			; FDL
  "Return T if x is a double-float NaN.")
56

liam's avatar
liam committed
57
(defmfun infinityp (x)
58 59
  "gsl_isinf" ((x :double))
  :c-return (cr :int)
liam's avatar
liam committed
60 61 62 63
  :return ((pmnil cr))
  :documentation			; FDL
  "Return +1 if x is positive infinity, -1 if negative infinity
   nil if finite.")
64

liam's avatar
liam committed
65 66
(defmfun finitep (x)
  "gsl_finite" ((x :double))
67
  :c-return (cr :int)
liam's avatar
liam committed
68 69 70
  :return ((= 1 cr))
  :documentation			; FDL
  "Return T if x is finite.")
71 72 73 74 75

;;;;****************************************************************************
;;; Elementary functions
;;;;****************************************************************************

liam's avatar
liam committed
76
(defmfun log+1 (x)
77 78
  "gsl_log1p" ((x :double))
  :c-return :double
liam's avatar
liam committed
79
  :documentation			; FDL
80 81
  "log(1+x), computed in a way that is accurate for small x.")

liam's avatar
liam committed
82
(defmfun exp-1 (x)
83 84
  "gsl_expm1" ((x :double))
  :c-return :double
liam's avatar
liam committed
85
  :documentation			; FDL
86 87
  "exp(x)-1, computed in a way that is accurate for small x.")

liam's avatar
liam committed
88
(defmfun hypotenuse* (x y)
89 90
  ;; This is redundant; there is "gsl_sf_hypot_e" defined as
  ;; #'hyptoenuse.
91 92
  "gsl_hypot" ((x :double) (y :double))
  :c-return :double
liam's avatar
liam committed
93
  :documentation			;FDL
94
  "The hypotenuse sqrt{x^2 + y^2} computed in a way that avoids overflow.")
95 96

;; Not clear why this function exists
liam's avatar
liam committed
97
(defmfun gsl-asinh (x)
98 99
   "gsl_asinh" ((x :double))
  :c-return :double
liam's avatar
liam committed
100 101
  :documentation			; FDL
  "Arc hyperbolic sine.")
102 103

;; Not clear why this function exists
liam's avatar
liam committed
104
(defmfun gsl-atanh (x)
105 106
   "gsl_atanh" ((x :double))
  :c-return :double
liam's avatar
liam committed
107 108
  :documentation			; FDL
  "Arc hyperbolic tangent.")
109 110 111 112 113 114 115 116 117 118 119 120

;;; gsl_ldexp
;;; gsl_frexp
;;; not mapped because CL has equivalents.

;;;;****************************************************************************
;;; Small integer powers
;;;;****************************************************************************

;;; Does CL need these?

#|
liam's avatar
liam committed
121
;; FDL
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
A common complaint about the standard C library is its lack of a
function for calculating (small) integer powers. GSL provides a
simple functions to fill this gap. For reasons of efficiency,
these functions do not check for overflow or underflow
conditions.

Function: double gsl_pow_int (double x, int n)
    This routine computes the power x^n for integer n. The power is computed efficiently--for example, x^8 is computed as ((x^2)^2)^2, requiring only 3 multiplications. A version of this function which also computes the numerical error in the result is available as gsl_sf_pow_int_e. 

Function: double gsl_pow_2 (const double x)
Function: double gsl_pow_3 (const double x)
Function: double gsl_pow_4 (const double x)
Function: double gsl_pow_5 (const double x)
Function: double gsl_pow_6 (const double x)
Function: double gsl_pow_7 (const double x)
Function: double gsl_pow_8 (const double x)
Function: double gsl_pow_9 (const double x)
    These functions can be used to compute small integer powers x^2, x^3, etc. efficiently. The functions will be inlined when possible so that use of these functions should be as efficient as explicitly writing the corresponding product expression.

|#

143 144 145 146 147 148 149 150 151
;;;; Testing the Sign of Numbers
;;; is all macros

;;;; Testing for Odd and Even Numbers
;;; is all macros

;;;; Maximum and Minimum functions
;;; is all macros and inline functions that have CL equivalents

152 153 154 155
;;;;****************************************************************************
;;; Approximate Comparison of Floating Point Numbers
;;;;****************************************************************************

liam's avatar
liam committed
156
;;; FDL
157 158 159 160 161 162
;;; It is sometimes useful to be able to compare two floating point
;;; numbers approximately, to allow for rounding and truncation
;;; errors. This function implements the approximate
;;; floating-point comparison algorithm proposed by D.E. Knuth in
;;; Section 4.2.2 of Seminumerical Algorithms (3rd edition).

liam's avatar
liam committed
163
(defmfun double-float-unequal (x y epsilon)
164 165
  "gsl_fcmp" ((x :double) (y :double) (epsilon :double))
  :c-return (cr :int)
liam's avatar
liam committed
166
  :documentation			; FDL
167 168 169 170
  "This function determines whether x and y are approximately equal
    to a relative accuracy epsilon.

    The relative accuracy is measured using an interval of size 2
liam's avatar
liam committed
171
    delta, where delta = 2^k \epsilon and k is the maximum
172 173 174 175
    base-2 exponent of x and y as computed by the function
    frexp().

    If x and y lie within this interval, they are considered
176
    approximately equal and the function returns nil. Otherwise if
177
    x < y, the function returns -1, or if x > y, the function
178 179 180 181 182 183 184
    returns +1."
  :return ((pmnil cr)))

;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
185
(save-test mathematical
liam's avatar
liam committed
186 187 188 189
  (log+1 0.001d0)
  (exp-1 0.001d0)
  (hypotenuse 3.0d0 4.0d0))

Liam Healy's avatar
Liam Healy committed
190 191
#| 
;;; I would like to add
192 193 194 195 196 197
(lisp-unit:define-test mathematical
  (lisp-unit:assert-true (nanp +nan+))
  (lisp-unit:assert-false (nanp 1.0d0))
  (lisp-unit:assert-true (finitep 1.0d0))
  (lisp-unit:assert-false (infinityp 1.0d0))
  (lisp-unit:assert-eq 1 (infinityp +positive-infinity+))
Liam Healy's avatar
Liam Healy committed
198 199
  (lisp-unit:assert-false (finitep +positive-infinity+)))
|#