mathematical.lisp 6.32 KB
Newer Older
1 2
;; Mathematical functions
;; Liam Healy, Wed Mar  8 2006 - 22:09
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2008-10-25 19:09:09EDT 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 19 20 21 22 23 24 25 26
;;;;****************************************************************************
;;; Infinities and Not-a-number
;;;;****************************************************************************

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

27 28 29 30 31 32
#+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)))
33

34 35
(export '(+nan+ +positive-infinity+ +negative-infinity+))

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

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

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

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

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

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

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

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

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

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

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

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

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

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

;;; Does CL need these?

#|
liam's avatar
liam committed
119
;; FDL
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
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.

|#

141 142 143 144 145 146 147 148 149
;;;; 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

150 151 152 153
;;;;****************************************************************************
;;; Approximate Comparison of Floating Point Numbers
;;;;****************************************************************************

liam's avatar
liam committed
154
;;; FDL
155 156 157 158 159 160
;;; 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
161
(defmfun double-float-unequal (x y epsilon)
162 163
  "gsl_fcmp" ((x :double) (y :double) (epsilon :double))
  :c-return (cr :int)
liam's avatar
liam committed
164
  :documentation			; FDL
165 166 167 168
  "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
169
    delta, where delta = 2^k \epsilon and k is the maximum
170 171 172 173
    base-2 exponent of x and y as computed by the function
    frexp().

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

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

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

Liam Healy's avatar
Liam Healy committed
188 189
#| 
;;; I would like to add
190 191 192 193 194 195
(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
196 197
  (lisp-unit:assert-false (finitep +positive-infinity+)))
|#