mathematical.lisp 7.5 KB
Newer Older
1 2
;; Mathematical functions
;; Liam Healy, Wed Mar  8 2006 - 22:09
3
;; Time-stamp: <2010-08-07 21:35:54EDT mathematical.lisp>
4
;;
5
;; Copyright 2006, 2007, 2008, 2009, 2010 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program 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 General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20 21 22

(in-package :gsl)

23
;;; Disable floating point traps for each CL implementation. 
24

25 26 27 28 29 30
;;; 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?
31

32 33
#+clisp (setf sys::*inhibit-floating-point-underflow* t)

34 35 36 37 38 39 40 41 42 43
;;;;****************************************************************************
;;; Infinities and Not-a-number
;;;;****************************************************************************

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

44 45 46 47 48 49
#+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)))
50

51
(export '(+nan+ +positive-infinity+ +negative-infinity+))
52 53 54
(map-name '+positive-infinity+ "gsl_posinf")
(map-name '+negative-infinity+ "gsl_neginf")
(map-name '+nan+ "gsl_nan")
55

56
(defconstant +nan+
57 58 59
  (ignore-errors
    (cffi:foreign-funcall "gsl_nan" :double)))

60
(defconstant +positive-infinity+
61 62 63
  (ignore-errors
    (cffi:foreign-funcall "gsl_posinf" :double)))

64
(defconstant +negative-infinity+
65 66 67
  (ignore-errors
    (cffi:foreign-funcall "gsl_neginf" :double)))

liam's avatar
liam committed
68
(defmfun nanp (x)
69 70
  "gsl_isnan" ((x :double))
  :c-return (cr :int)
liam's avatar
liam committed
71 72 73
  :return ((= 1 cr))
  :documentation			; FDL
  "Return T if x is a double-float NaN.")
74

liam's avatar
liam committed
75
(defmfun infinityp (x)
76 77
  "gsl_isinf" ((x :double))
  :c-return (cr :int)
liam's avatar
liam committed
78 79 80
  :return ((pmnil cr))
  :documentation			; FDL
  "Return +1 if x is positive infinity, -1 if negative infinity
81
   nil if finite.  Some platforms will return only +1 for either sign.")
82

liam's avatar
liam committed
83 84
(defmfun finitep (x)
  "gsl_finite" ((x :double))
85
  :c-return (cr :int)
liam's avatar
liam committed
86 87 88
  :return ((= 1 cr))
  :documentation			; FDL
  "Return T if x is finite.")
89 90 91 92 93

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

liam's avatar
liam committed
94
(defmfun log+1 (x)
95 96
  "gsl_log1p" ((x :double))
  :c-return :double
liam's avatar
liam committed
97
  :documentation			; FDL
98 99
  "log(1+x), computed in a way that is accurate for small x.")

liam's avatar
liam committed
100
(defmfun exp-1 (x)
101 102
  "gsl_expm1" ((x :double))
  :c-return :double
liam's avatar
liam committed
103
  :documentation			; FDL
104 105
  "exp(x)-1, computed in a way that is accurate for small x.")

liam's avatar
liam committed
106
(defmfun hypotenuse* (x y)
107
  ;; This is redundant; there is "gsl_sf_hypot_e" defined as
108
  ;; #'hypotenuse.
109 110
  "gsl_hypot" ((x :double) (y :double))
  :c-return :double
liam's avatar
liam committed
111
  :documentation			;FDL
112
  "The hypotenuse sqrt{x^2 + y^2} computed in a way that avoids overflow.")
113 114

;; Not clear why this function exists
liam's avatar
liam committed
115
(defmfun gsl-asinh (x)
116 117
   "gsl_asinh" ((x :double))
  :c-return :double
liam's avatar
liam committed
118 119
  :documentation			; FDL
  "Arc hyperbolic sine.")
120 121

;; Not clear why this function exists
liam's avatar
liam committed
122
(defmfun gsl-atanh (x)
123 124
   "gsl_atanh" ((x :double))
  :c-return :double
liam's avatar
liam committed
125 126
  :documentation			; FDL
  "Arc hyperbolic tangent.")
127 128 129 130 131 132 133 134 135 136 137 138

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

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

;;; Does CL need these?

#|
liam's avatar
liam committed
139
;; FDL
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
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.

|#

161 162 163 164 165 166 167 168 169
;;;; 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

170 171 172 173
;;;;****************************************************************************
;;; Approximate Comparison of Floating Point Numbers
;;;;****************************************************************************

liam's avatar
liam committed
174
;;; FDL
175 176 177 178 179 180
;;; 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
181
(defmfun double-float-unequal (x y epsilon)
182 183
  "gsl_fcmp" ((x :double) (y :double) (epsilon :double))
  :c-return (cr :int)
liam's avatar
liam committed
184
  :documentation			; FDL
185 186 187 188
  "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
189
    delta, where delta = 2^k \epsilon and k is the maximum
190 191 192 193
    base-2 exponent of x and y as computed by the function
    frexp().

    If x and y lie within this interval, they are considered
194
    approximately equal and the function returns nil. Otherwise if
195
    x < y, the function returns -1, or if x > y, the function
196 197 198
    returns +1."
  :return ((pmnil cr)))

199 200 201 202 203 204
;;;;****************************************************************************
;;;; pi
;;;;****************************************************************************

(defconstant dpi (coerce pi 'double-float))

205 206 207 208
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
209
(save-test mathematical
liam's avatar
liam committed
210 211 212 213
  (log+1 0.001d0)
  (exp-1 0.001d0)
  (hypotenuse 3.0d0 4.0d0))

Liam Healy's avatar
Liam Healy committed
214 215
#| 
;;; I would like to add
216 217 218 219 220 221
(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
222 223
  (lisp-unit:assert-false (finitep +positive-infinity+)))
|#