Commit 2fa23a00 authored by liam's avatar liam

Add ability to return value of C function in defun-sf. Start form for

general polynomial solver.  If on SBCL, generate nans and infinities.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@2994 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 2bf37cb3
......@@ -3,13 +3,17 @@
; description: Mathematical functions
; date: Wed Mar 8 2006 - 22:09
; author: Liam M. Healy
; modified: Sat Mar 18 2006 - 15:13
; modified: Wed Mar 22 2006 - 20:28
;********************************************************
(in-package :gsl)
(export '(nanp infinityp finitep log+1 exp-1 hypotenuse approximately=))
;;;; Mathematical Constants
;;; Is all macros
;;;;****************************************************************************
;;; Infinities and Not-a-number
;;;;****************************************************************************
......@@ -20,7 +24,21 @@
(when (or (= 1 v) (= -1 v))
v)))
(defun nanp (x)
#+sbcl (sb-int:set-floating-point-modes :traps nil)
(defparameter *nan*
(ignore-errors
(cffi:foreign-funcall "gsl_nan" :double)))
(defparameter *positive-infinity*
(ignore-errors
(cffi:foreign-funcall "gsl_posinf" :double)))
(defparameter *negative-infinity*
(ignore-errors
(cffi:foreign-funcall "gsl_neginf" :double)))
(defunx-map nanp "gsl_isnan" (x)
"Return T if x is a double-float NaN."
(= 1
(cffi:foreign-funcall
......@@ -28,7 +46,7 @@
:double x
:int)))
(defun infinityp (x)
(defunx-map infinityp "gsl_isinf" (x)
"Return +1 if x is positive infinity, -1 if negative infinity
nil if finite."
(pmnil
......@@ -37,7 +55,7 @@
:double x
:int)))
(defun finitep (x)
(defunx-map finitep "gsl_finite" (x)
"Return T if finite."
(= 1
(cffi:foreign-funcall
......@@ -49,21 +67,21 @@
;;; Elementary functions
;;;;****************************************************************************
(defun log+1 (x)
(defunx-map log+1 "gsl_log1p" (x)
"log(1+x), computed in a way that is accurate for small x."
(cffi:foreign-funcall
"gsl_log1p"
:double x
:double))
(defun exp-1 (x)
(defunx-map exp-1 "gsl_expm1" (x)
"exp(x)-1, computed in a way that is accurate for small x."
(cffi:foreign-funcall
"gsl_expm1"
:double x
:double))
(defun hypotenuse (x y)
(defunx-map hypotenuse "gsl_hypot" (x y)
"sqrt{x^2 + y^2} computed in a way that avoids overflow."
(cffi:foreign-funcall
"gsl_hypot"
......@@ -72,7 +90,7 @@
:double))
;; Not clear why this function exists
(defun gsl-asinh (x)
(defunx-map gsl-asinh "gsl_asinh" (x)
"arcsinh"
(cffi:foreign-funcall
"gsl_asinh"
......@@ -80,8 +98,8 @@
:double))
;;; Not clear why this function exists
(defun gsl-atanh (x)
"arcsinh"
(defunx-map gsl-atanh "gsl_atanh" (x)
"arctanh"
(cffi:foreign-funcall
"gsl_atanh"
:double x
......@@ -119,6 +137,15 @@ Function: double gsl_pow_9 (const double x)
|#
;;;; 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
;;;;****************************************************************************
;;; Approximate Comparison of Floating Point Numbers
;;;;****************************************************************************
......@@ -129,7 +156,7 @@ Function: double gsl_pow_9 (const double x)
;;; floating-point comparison algorithm proposed by D.E. Knuth in
;;; Section 4.2.2 of Seminumerical Algorithms (3rd edition).
(defun approximately= (x y epsilon)
(defunx-map double-float-equal "gsl_fcmp" (x y epsilon)
"This function determines whether x and y are approximately equal
to a relative accuracy epsilon.
......@@ -152,4 +179,3 @@ Function: double gsl_pow_9 (const double x)
:double y
:double epsilon
:double)))
......@@ -3,7 +3,7 @@
; description: Macros to interface GSL functions.
; date: Mon Mar 6 2006 - 22:35
; author: Liam M. Healy
; modified: Wed Mar 22 2006 - 12:20
; modified: Wed Mar 22 2006 - 23:12
;********************************************************
(in-package :gsl)
......@@ -78,7 +78,8 @@ and a scaling exponent e10, such that the value is val*10^e10."
(defun pick-result (decl)
(if (listp (second decl))
`((loop for i from 0 below ,(second (second decl))
collect (cffi:mem-aref ,(first decl) ,(first (second decl)) i)))
collect
(cffi:mem-aref ,(first decl) ,(first (second decl)) i)))
(case (second decl)
(sf-result
`((cffi:foreign-slot-value ,(first decl) 'sf-result 'val)
......@@ -112,7 +113,7 @@ and a scaling exponent e10, such that the value is val*10^e10."
;;; New name?
(defmacro defun-sf
(cl-name arguments gsl-name
&key documentation return mode (c-return-value :error-code))
&key documentation return mode (c-return-value :error-code))
"Define a CL function that provides an interface to a GSL function.
If cl-name is :lambda, make a lambda. Arguments:
arguments: a list of input arguments (symbol type) to the GSL function
......@@ -153,13 +154,16 @@ and a scaling exponent e10, such that the value is val*10^e10."
,@(if (eq c-return-value :error-code)
`((check-gsl-status status `(,',cl-name ,,@args))))
(values
,@(if (eq c-return-value :number-of-answers)
(mapcan
(lambda (decl seq)
`((when (> status ,seq) ,@(pick-result decl))))
return-symb-type
(loop for i below (length return) collect i))
(mapcan #'pick-result return-symb-type))))))))
,@(case c-return-value
(:number-of-answers
(mapcan
(lambda (decl seq)
`((when (> status ,seq) ,@(pick-result decl))))
return-symb-type
(loop for i below (length return) collect i)))
(:return '(status))
(:error-code
(mapcan #'pick-result return-symb-type)))))))))
;;; arguments: list like ((x :double) (y :double))
;;; mode: t or nil
......@@ -3,7 +3,7 @@
; description: Polynomials
; date: Tue Mar 21 2006 - 18:33
; author: Liam M. Healy
; modified: Wed Mar 22 2006 - 12:10
; modified: Wed Mar 22 2006 - 23:17
;********************************************************
;;; $Id: $
......@@ -161,3 +161,44 @@
if a root does not exist, the value returned for it will be NIL."
:return (gsl-complex gsl-complex gsl-complex)
:c-return-value :number-of-answers)
;;;;****************************************************************************
;;;; General Polynomial Equations
;;;;****************************************************************************
;;; See /usr/include/gsl/gsl_poly.h
(cffi:defcstruct poly-complex-workspace
(nc :uint)
(matrix :pointer))
(export 'with-poly-complex-workspace)
(defmacro with-poly-complex-workspace (workspace &body body)
"Macro to create and cleanup workspace for polynomial root solver."
`(let ((,workspace
(funcall
(defun-sf :lambda ((n :uint))
"gsl_poly_complex_workspace_alloc"
:return
(poly-complex-workspace)
:c-return-value :return))))
(unwind-protect
(progn ,@body)
(funcall
(defun-sf :lambda ((,workspace poly-complex-workspace))
"gsl_poly_complex_workspace_free"
:c-return-value :void)))))
(defun-sf polynomial-solve
((a :pointer) (n :uint) (workspace poly-complex-workspace))
"gsl_poly_complex_solve"
:documentation
"The roots of the general polynomial
@c{$P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1}$}
@math{P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_@{n-1@} x^@{n-1@}} using
balanced-QR reduction of the companion matrix. The parameter @var{n}
specifies the length of the coefficient array. The coefficient of the
highest order term must be non-zero. The function requires a workspace
@var{w} of the appropriate size. The @math{n-1} roots are returned in
the packed complex array @var{z} of length @math{2(n-1)}, alternating
real and imaginary parts."
:return ((gsl-complex n)))
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment