From 2fa23a009a8cf64b194d1ce9f313f4dc8f355b0d Mon Sep 17 00:00:00 2001 From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30> Date: Thu, 23 Mar 2006 04:23:46 +0000 Subject: [PATCH] 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 --- general/mathematical.lisp | 50 +++++++++++++++++++++++++++++---------- interface.lisp | 24 +++++++++++-------- polynomial.lisp | 43 ++++++++++++++++++++++++++++++++- 3 files changed, 94 insertions(+), 23 deletions(-) diff --git a/general/mathematical.lisp b/general/mathematical.lisp index 161bc98d..3c496373 100644 --- a/general/mathematical.lisp +++ b/general/mathematical.lisp @@ -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))) - diff --git a/interface.lisp b/interface.lisp index acec769e..de56d463 100644 --- a/interface.lisp +++ b/interface.lisp @@ -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 diff --git a/polynomial.lisp b/polynomial.lisp index b96de5a4..6327f859 100644 --- a/polynomial.lisp +++ b/polynomial.lisp @@ -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))) -- GitLab