From 6438e06ddda69d5f60219ac165519b2109aae379 Mon Sep 17 00:00:00 2001
From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30>
Date: Wed, 22 Mar 2006 15:10:14 +0000
Subject: [PATCH] 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
---
 interface.lisp  | 58 +++++++++++++++++++++++++++++--------------------
 polynomial.lisp | 41 +++++++++++-----------------------
 2 files changed, 48 insertions(+), 51 deletions(-)

diff --git a/interface.lisp b/interface.lisp
index 80c2e680..feb5b998 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:    Tue Mar 21 2006 - 23:10
+; modified:    Wed Mar 22 2006 - 10:08
 ;********************************************************
 
 (in-package :gsl)
@@ -78,7 +78,7 @@ 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)
@@ -101,10 +101,18 @@ and a scaling exponent e10, such that the value is val*10^e10."
 	   `(',(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.
 ;;; New name?
 (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
    GSL function definition.  If cl-name is :lambda, make a lambda."
   (let ((args (mapcar #'first arguments))
@@ -122,26 +130,30 @@ and a scaling exponent e10, such that the value is val*10^e10."
     `(,@(if (eq cl-name :lambda)
 	    '(lambda)
 	    `(defunx-map ,cl-name ,gsl-name))
-	,(if mode 
-	     `(,@args &optional (mode :double-prec))
-	     `(,@args))
-	,@(when documentation (list documentation))
-	(cffi:with-foreign-objects
-	    ,(mapcar #'wfo-declare return-symb-type)
-	  (let ((status
-		 (cffi:foreign-funcall
-		  ,gsl-name
-		  ,@(mapcan (lambda (ar) (list (second ar) (first ar)))
-			    arguments)
-		  ,@(when mode '(sf-mode mode))
-		  ,@(mapcan (lambda (r) `(:pointer ,(first r))) return-symb-type)
-		  :int)))
-	    (unless (eql :success (cffi:foreign-enum-keyword 'gsl-errorno status))
-	      (warn 'gsl-warning
-		    :gsl-errno status :gsl-context `(,',cl-name ,,@args)))
-	    (values
-	     ,@(mapcan #'pick-result return-symb-type)))))))
-
+      ,(if mode 
+	   `(,@args &optional (mode :double-prec))
+	   `(,@args))
+      ,@(when documentation (list documentation))
+      (cffi:with-foreign-objects
+	  ,(mapcar #'wfo-declare return-symb-type)
+	(let ((status
+	       (cffi:foreign-funcall
+		,gsl-name
+		,@(mapcan (lambda (ar) (list (second ar) (first ar)))
+			  arguments)
+		,@(when mode '(sf-mode mode))
+		,@(mapcan (lambda (r) `(:pointer ,(first r))) return-symb-type)
+		:int)))
+	  ,@(if (eq return-code :check-status)
+		`((check-gsl-status status `(,',cl-name ,,@args))))
+	  (values
+	   ,@(if (eq return-code :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))))))))
 
 ;;; arguments: list like ((x :double) (y :double))
 ;;; mode: t or nil
diff --git a/polynomial.lisp b/polynomial.lisp
index 91c2f907..b826f0fb 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:    Tue Mar 21 2006 - 23:07
+; modified:    Wed Mar 22 2006 - 10:06
 ;********************************************************
 ;;; $Id: $
 
@@ -12,9 +12,6 @@
 ;;; To do:
 ;;; Awaits incorporation of arrays into CFFI.
 ;;; 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
 
 ;;;;****************************************************************************
@@ -39,7 +36,7 @@
 
 ;;; Use with-divided-difference to compute the divided difference, an
 ;;; 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.
 
 (export
@@ -116,20 +113,14 @@
 ;;;; 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))
   "gsl_poly_solve_quadratic"
   :documentation
   "The real roots of the quadratic equation a x^2 + b x + c = 0.
-   Warning: if roots are not real, nonsense is produced and
-   no warning is given."
-  :return (:double :double))
+   Two values are always returned; if the roots are not real, these
+   values are NIL."
+  :return (:double :double)
+  :return-code :number-of-answers)
 
 (defun-sf solve-quadratic-complex ((a :double) (b :double) (c :double))
   "gsl_poly_complex_solve_quadratic"
@@ -141,14 +132,6 @@
 ;;;; 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)
 ;;; -3.0d0
 ;;; 1.9999999999999996d0
@@ -156,16 +139,18 @@
 
 ;;; (solve-cubic -1.0d0 1.0d0 -1.0d0)
 ;;; 1.0d0
-;;; 1.9999999999999996d0
-;;; 7.0d0
+;;; NIL
+;;; NIL
 
 (defun-sf solve-cubic ((a :double) (b :double) (c :double))
   "gsl_poly_solve_cubic"
   :documentation
   "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.
-   No special indication is given if there is only one real root."
-  :return (:double :double :double))
+   with a leading coefficient of unity.  The roots are given
+   in ascending order.  Three values are always returned;
+   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))
   "gsl_poly_complex_solve_cubic"
-- 
GitLab