mathematical.lisp 7.5 KB
 liam committed Jan 21, 2008 1 2 ;; Mathematical functions ;; Liam Healy, Wed Mar 8 2006 - 22:09  Liam Healy committed Aug 08, 2010 3 ;; Time-stamp: <2010-08-07 21:35:54EDT mathematical.lisp>  Liam Healy committed Dec 27, 2009 4 ;;  Liam M. Healy committed Jul 19, 2010 5 ;; Copyright 2006, 2007, 2008, 2009, 2010 Liam M. Healy  Liam Healy committed Dec 27, 2009 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 .  liam committed Mar 15, 2006 20 21 22  (in-package :gsl)  liam committed May 30, 2006 23 ;;; Disable floating point traps for each CL implementation.  liam committed Mar 23, 2006 24   liam committed May 30, 2006 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?  liam committed Mar 23, 2006 31   Liam Healy committed Jan 01, 2009 32 33 #+clisp (setf sys::*inhibit-floating-point-underflow* t)  liam committed Mar 15, 2006 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)))  liam committed Mar 23, 2006 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)))  liam committed Mar 23, 2006 50   liam committed May 30, 2006 51 (export '(+nan+ +positive-infinity+ +negative-infinity+))  Liam Healy committed May 22, 2010 52 53 54 (map-name '+positive-infinity+ "gsl_posinf") (map-name '+negative-infinity+ "gsl_neginf") (map-name '+nan+ "gsl_nan")  liam committed May 30, 2006 55   liam committed Mar 23, 2006 56 (defconstant +nan+  liam committed Mar 23, 2006 57 58 59  (ignore-errors (cffi:foreign-funcall "gsl_nan" :double)))  liam committed Mar 23, 2006 60 (defconstant +positive-infinity+  liam committed Mar 23, 2006 61 62 63  (ignore-errors (cffi:foreign-funcall "gsl_posinf" :double)))  liam committed Mar 23, 2006 64 (defconstant +negative-infinity+  liam committed Mar 23, 2006 65 66 67  (ignore-errors (cffi:foreign-funcall "gsl_neginf" :double)))  liam committed Feb 18, 2008 68 (defmfun nanp (x)  liam committed May 30, 2006 69 70  "gsl_isnan" ((x :double)) :c-return (cr :int)  liam committed Feb 18, 2008 71 72 73  :return ((= 1 cr)) :documentation ; FDL "Return T if x is a double-float NaN.")  liam committed Mar 15, 2006 74   liam committed Feb 18, 2008 75 (defmfun infinityp (x)  liam committed May 30, 2006 76 77  "gsl_isinf" ((x :double)) :c-return (cr :int)  liam committed Feb 18, 2008 78 79 80  :return ((pmnil cr)) :documentation ; FDL "Return +1 if x is positive infinity, -1 if negative infinity  Liam Healy committed Apr 18, 2010 81  nil if finite. Some platforms will return only +1 for either sign.")  liam committed May 30, 2006 82   liam committed Feb 18, 2008 83 84 (defmfun finitep (x) "gsl_finite" ((x :double))  liam committed May 30, 2006 85  :c-return (cr :int)  liam committed Feb 18, 2008 86 87 88  :return ((= 1 cr)) :documentation ; FDL "Return T if x is finite.")  liam committed Mar 15, 2006 89 90 91 92 93  ;;;;**************************************************************************** ;;; Elementary functions ;;;;****************************************************************************  liam committed Feb 18, 2008 94 (defmfun log+1 (x)  liam committed May 30, 2006 95 96  "gsl_log1p" ((x :double)) :c-return :double  liam committed Feb 18, 2008 97  :documentation ; FDL  liam committed May 30, 2006 98 99  "log(1+x), computed in a way that is accurate for small x.")  liam committed Feb 18, 2008 100 (defmfun exp-1 (x)  liam committed May 30, 2006 101 102  "gsl_expm1" ((x :double)) :c-return :double  liam committed Feb 18, 2008 103  :documentation ; FDL  liam committed May 30, 2006 104 105  "exp(x)-1, computed in a way that is accurate for small x.")  liam committed Feb 18, 2008 106 (defmfun hypotenuse* (x y)  liam committed Jan 21, 2008 107  ;; This is redundant; there is "gsl_sf_hypot_e" defined as  Liam M. Healy committed Jul 19, 2010 108  ;; #'hypotenuse.  liam committed May 30, 2006 109 110  "gsl_hypot" ((x :double) (y :double)) :c-return :double  liam committed Feb 18, 2008 111  :documentation ;FDL  liam committed Jan 21, 2008 112  "The hypotenuse sqrt{x^2 + y^2} computed in a way that avoids overflow.")  liam committed Mar 15, 2006 113 114  ;; Not clear why this function exists  liam committed Feb 18, 2008 115 (defmfun gsl-asinh (x)  liam committed May 30, 2006 116 117  "gsl_asinh" ((x :double)) :c-return :double  liam committed Feb 18, 2008 118 119  :documentation ; FDL "Arc hyperbolic sine.")  liam committed May 30, 2006 120 121  ;; Not clear why this function exists  liam committed Feb 18, 2008 122 (defmfun gsl-atanh (x)  liam committed May 30, 2006 123 124  "gsl_atanh" ((x :double)) :c-return :double  liam committed Feb 18, 2008 125 126  :documentation ; FDL "Arc hyperbolic tangent.")  liam committed Mar 15, 2006 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 committed Feb 18, 2008 139 ;; FDL  liam committed Mar 15, 2006 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. |#  liam committed Mar 23, 2006 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  liam committed Mar 15, 2006 170 171 172 173 ;;;;**************************************************************************** ;;; Approximate Comparison of Floating Point Numbers ;;;;****************************************************************************  liam committed Feb 18, 2008 174 ;;; FDL  liam committed Mar 15, 2006 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 committed Feb 18, 2008 181 (defmfun double-float-unequal (x y epsilon)  liam committed May 30, 2006 182 183  "gsl_fcmp" ((x :double) (y :double) (epsilon :double)) :c-return (cr :int)  liam committed Feb 18, 2008 184  :documentation ; FDL  liam committed Mar 15, 2006 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 committed Feb 18, 2008 189  delta, where delta = 2^k \epsilon and k is the maximum  liam committed Mar 15, 2006 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  liam committed May 30, 2006 194  approximately equal and the function returns nil. Otherwise if  liam committed Mar 15, 2006 195  x < y, the function returns -1, or if x > y, the function  liam committed May 30, 2006 196 197 198  returns +1." :return ((pmnil cr)))  Liam Healy committed Aug 08, 2010 199 200 201 202 203 204 ;;;;**************************************************************************** ;;;; pi ;;;;**************************************************************************** (defconstant dpi (coerce pi 'double-float))  liam committed May 30, 2006 205 206 207 208 ;;;;**************************************************************************** ;;;; Examples and unit test ;;;;****************************************************************************  Liam Healy committed Oct 25, 2008 209 (save-test mathematical  liam committed Feb 18, 2008 210 211 212 213  (log+1 0.001d0) (exp-1 0.001d0) (hypotenuse 3.0d0 4.0d0))  Liam Healy committed Oct 25, 2008 214 215 #| ;;; I would like to add  liam committed May 30, 2006 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 committed Oct 25, 2008 222 223  (lisp-unit:assert-false (finitep +positive-infinity+))) |#`