Skip to content
Snippets Groups Projects
Commit 6438e06d authored by liam's avatar liam
Browse files

Quadratic and cubic solvers returned integer is the number of real

roots, not a GSL error code.  Add another argument :return-code to
defun-sf to handle the number of return values properly.  Break out 
#'check-gsl-status from #'defun-sf.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@2991 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 0e3b5615
No related branches found
No related tags found
No related merge requests found
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
; description: Macros to interface GSL functions. ; description: Macros to interface GSL functions.
; date: Mon Mar 6 2006 - 22:35 ; date: Mon Mar 6 2006 - 22:35
; author: Liam M. Healy ; author: Liam M. Healy
; modified: Tue Mar 21 2006 - 23:10 ; modified: Wed Mar 22 2006 - 10:08
;******************************************************** ;********************************************************
(in-package :gsl) (in-package :gsl)
...@@ -78,7 +78,7 @@ and a scaling exponent e10, such that the value is val*10^e10." ...@@ -78,7 +78,7 @@ and a scaling exponent e10, such that the value is val*10^e10."
(defun pick-result (decl) (defun pick-result (decl)
(if (listp (second decl)) (if (listp (second decl))
`((loop for i from 0 below ,(second (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) (case (second decl)
(sf-result (sf-result
`((cffi:foreign-slot-value ,(first decl) 'sf-result 'val) `((cffi:foreign-slot-value ,(first decl) 'sf-result 'val)
...@@ -101,10 +101,18 @@ and a scaling exponent e10, such that the value is val*10^e10." ...@@ -101,10 +101,18 @@ and a scaling exponent e10, such that the value is val*10^e10."
`(',(second d)) `(',(second d))
`(',(first (second d)) ,(second (second d)))))) `(',(first (second d)) ,(second (second d))))))
(defun check-gsl-status (status-code context)
"Check the return status code from a GSL function and signal a warning
if it is not :SUCCESS."
(unless (eql :success (cffi:foreign-enum-keyword 'gsl-errorno status-code))
(warn 'gsl-warning
:gsl-errno status :gsl-context context)))
;;; Warning isn't quite right for lambdas. ;;; Warning isn't quite right for lambdas.
;;; New name? ;;; New name?
(defmacro defun-sf (defmacro defun-sf
(cl-name arguments gsl-name &key documentation return mode) (cl-name arguments gsl-name
&key documentation return mode (return-code :check-status))
"Define a mathematical special function from GSL using the _e form "Define a mathematical special function from GSL using the _e form
GSL function definition. If cl-name is :lambda, make a lambda." GSL function definition. If cl-name is :lambda, make a lambda."
(let ((args (mapcar #'first arguments)) (let ((args (mapcar #'first arguments))
...@@ -122,26 +130,30 @@ and a scaling exponent e10, such that the value is val*10^e10." ...@@ -122,26 +130,30 @@ and a scaling exponent e10, such that the value is val*10^e10."
`(,@(if (eq cl-name :lambda) `(,@(if (eq cl-name :lambda)
'(lambda) '(lambda)
`(defunx-map ,cl-name ,gsl-name)) `(defunx-map ,cl-name ,gsl-name))
,(if mode ,(if mode
`(,@args &optional (mode :double-prec)) `(,@args &optional (mode :double-prec))
`(,@args)) `(,@args))
,@(when documentation (list documentation)) ,@(when documentation (list documentation))
(cffi:with-foreign-objects (cffi:with-foreign-objects
,(mapcar #'wfo-declare return-symb-type) ,(mapcar #'wfo-declare return-symb-type)
(let ((status (let ((status
(cffi:foreign-funcall (cffi:foreign-funcall
,gsl-name ,gsl-name
,@(mapcan (lambda (ar) (list (second ar) (first ar))) ,@(mapcan (lambda (ar) (list (second ar) (first ar)))
arguments) arguments)
,@(when mode '(sf-mode mode)) ,@(when mode '(sf-mode mode))
,@(mapcan (lambda (r) `(:pointer ,(first r))) return-symb-type) ,@(mapcan (lambda (r) `(:pointer ,(first r))) return-symb-type)
:int))) :int)))
(unless (eql :success (cffi:foreign-enum-keyword 'gsl-errorno status)) ,@(if (eq return-code :check-status)
(warn 'gsl-warning `((check-gsl-status status `(,',cl-name ,,@args))))
:gsl-errno status :gsl-context `(,',cl-name ,,@args))) (values
(values ,@(if (eq return-code :number-of-answers)
,@(mapcan #'pick-result return-symb-type))))))) (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))))))))
;;; arguments: list like ((x :double) (y :double)) ;;; arguments: list like ((x :double) (y :double))
;;; mode: t or nil ;;; mode: t or nil
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
; description: Polynomials ; description: Polynomials
; date: Tue Mar 21 2006 - 18:33 ; date: Tue Mar 21 2006 - 18:33
; author: Liam M. Healy ; author: Liam M. Healy
; modified: Tue Mar 21 2006 - 23:07 ; modified: Wed Mar 22 2006 - 10:06
;******************************************************** ;********************************************************
;;; $Id: $ ;;; $Id: $
...@@ -12,9 +12,6 @@ ...@@ -12,9 +12,6 @@
;;; To do: ;;; To do:
;;; Awaits incorporation of arrays into CFFI. ;;; Awaits incorporation of arrays into CFFI.
;;; Divided difference needs to be checked. ;;; Divided difference needs to be checked.
;;; Quadratic solver always gives output range warning, except when it needs to.
;;; Cubic solver always gives output range warning.
;;; finish ;;; finish
;;;;**************************************************************************** ;;;;****************************************************************************
...@@ -39,7 +36,7 @@ ...@@ -39,7 +36,7 @@
;;; Use with-divided-difference to compute the divided difference, an ;;; Use with-divided-difference to compute the divided difference, an
;;; opaque object that is then passed to ;;; opaque object that is then passed to
;;; get-divided-difference, eval-divided-difference ;;; get-divided-difference, eval-divided-difference or taylor-divided-difference
;;; in the body. ;;; in the body.
(export (export
...@@ -116,20 +113,14 @@ ...@@ -116,20 +113,14 @@
;;;; Quadratic Equations ;;;; Quadratic Equations
;;;;**************************************************************************** ;;;;****************************************************************************
;;; GSL seems to always give a output range error warning:
;;; WARNING:
;;; GSL condition ERANGE (2), Output range error in
;;; (SOLVE-QUADRATIC 1.0d0 2.0d0 1.0d0)
;;; Oddly, there is NO warning if a complex answer is generated
;;; from the real root solver.
(defun-sf solve-quadratic ((a :double) (b :double) (c :double)) (defun-sf solve-quadratic ((a :double) (b :double) (c :double))
"gsl_poly_solve_quadratic" "gsl_poly_solve_quadratic"
:documentation :documentation
"The real roots of the quadratic equation a x^2 + b x + c = 0. "The real roots of the quadratic equation a x^2 + b x + c = 0.
Warning: if roots are not real, nonsense is produced and Two values are always returned; if the roots are not real, these
no warning is given." values are NIL."
:return (:double :double)) :return (:double :double)
:return-code :number-of-answers)
(defun-sf solve-quadratic-complex ((a :double) (b :double) (c :double)) (defun-sf solve-quadratic-complex ((a :double) (b :double) (c :double))
"gsl_poly_complex_solve_quadratic" "gsl_poly_complex_solve_quadratic"
...@@ -141,14 +132,6 @@ ...@@ -141,14 +132,6 @@
;;;; Cubic Equations ;;;; Cubic Equations
;;;;**************************************************************************** ;;;;****************************************************************************
;;; GSL seems to always give warning:
;;; WARNING:
;;; GSL condition EFAULT (3), Invalid pointer in
;;; (SOLVE-CUBIC -3.0d0 3.0d0 -1.0d0)
;;; GSL appears to give a different return code (EDOM) if there is
;;; only one real root.
;;; (solve-cubic -6.0d0 -13.0d0 42.0d0) ;;; (solve-cubic -6.0d0 -13.0d0 42.0d0)
;;; -3.0d0 ;;; -3.0d0
;;; 1.9999999999999996d0 ;;; 1.9999999999999996d0
...@@ -156,16 +139,18 @@ ...@@ -156,16 +139,18 @@
;;; (solve-cubic -1.0d0 1.0d0 -1.0d0) ;;; (solve-cubic -1.0d0 1.0d0 -1.0d0)
;;; 1.0d0 ;;; 1.0d0
;;; 1.9999999999999996d0 ;;; NIL
;;; 7.0d0 ;;; NIL
(defun-sf solve-cubic ((a :double) (b :double) (c :double)) (defun-sf solve-cubic ((a :double) (b :double) (c :double))
"gsl_poly_solve_cubic" "gsl_poly_solve_cubic"
:documentation :documentation
"Find the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0 "Find the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0
with a leading coefficient of unity. The real roots are returned. with a leading coefficient of unity. The roots are given
No special indication is given if there is only one real root." in ascending order. Three values are always returned;
:return (:double :double :double)) if a roots is not real, the value returned for it will be NIL."
:return (:double :double :double)
:return-code :number-of-answers)
(defun-sf solve-cubic-complex ((a :double) (b :double) (c :double)) (defun-sf solve-cubic-complex ((a :double) (b :double) (c :double))
"gsl_poly_complex_solve_cubic" "gsl_poly_complex_solve_cubic"
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment