From cbcc15ddbf8e24511b8b67ccc41bcd1cd48faec7 Mon Sep 17 00:00:00 2001
From: Liam Healy <liam@thinkpad.local>
Date: Sat, 31 Jan 2009 22:04:01 -0500
Subject: [PATCH] Polynomial uses #'evaluate

Polynomial evaluation is now done with a method of #'evaluate.  We
can't do complex evaluation because the return value would be complex,
and there is no way to handle that in CFFI yet.  Evaluation with
divided differences is a keyword switched part of the same method; it
works but gives an puzzling style warning in SBCL about the keyword.
The use of an marray in an optional/key argument necessitated a change
to defmfun that redefined #'body-expand from #'body-no-optional-arg,
made #'body-no-optional-arg a new function that wraps necessary
array-handling forms around the expanded unswitched body form from
body-expand.
---
 init/defmfun-single.lisp |  33 ++++++----
 polynomial.lisp          |  72 ++++++++++++++--------
 tests/polynomial.lisp    | 130 +++++++++++++++++++--------------------
 3 files changed, 133 insertions(+), 102 deletions(-)

diff --git a/init/defmfun-single.lisp b/init/defmfun-single.lisp
index 7a2d2259..7a78524b 100644
--- a/init/defmfun-single.lisp
+++ b/init/defmfun-single.lisp
@@ -1,6 +1,6 @@
 ;; Helpers that define a single GSL function interface
 ;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp
-;; Time-stamp: <2009-01-25 10:37:13EST defmfun-single.lisp>
+;; Time-stamp: <2009-01-31 21:19:28EST defmfun-single.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -77,10 +77,9 @@
      (mapdown (eq body-maker 'body-optional-arg)))
   "A complete definition form, starting with defun, :method, or defmethod."
   (destructuring-bind
-	(&key documentation inputs outputs before after
+	(&key documentation before after
 	      qualifier gsl-version &allow-other-keys)
       key-args
-    (declare (ignorable inputs outputs)) ; workaround for compiler errors that don't see it's used
     (if (have-at-least-gsl-version gsl-version)
 	`(,definition
 	     ,@(when (and name (not (defgeneric-method-p name)))
@@ -108,12 +107,7 @@
 			      'append
 			      (mapcar 'rest (subseq arglist (1+ auxstart))))))))))))
 	   ,@(when documentation (list documentation))
-	   #-native
-	   ,(funcall body-maker name arglist gsl-name c-arguments key-args)
-	   #+native
-	   ,(native-pointer
-	     (union inputs outputs)
-	     (funcall body-maker name arglist gsl-name c-arguments key-args)))
+	   ,(funcall body-maker name arglist gsl-name c-arguments key-args))
         `(,definition
 	     ,@(when (and name (not (defgeneric-method-p name)))
 		     (list name))
@@ -131,6 +125,23 @@
 	  `(progn ,@body))))
 
 (defun body-no-optional-arg (name arglist gsl-name c-arguments key-args)
+  "Wrap necessary array-handling forms around the expanded unswitched
+  body form."
+  #-native
+  (destructuring-bind (&key inputs &allow-other-keys) key-args
+    `(progn
+       ,@(mapcar (lambda (v) `(copy-cl-to-c ,v))
+		 (intersection
+		  inputs (arglist-plain-and-categories arglist nil)))
+       ,(body-expand name arglist gsl-name c-arguments key-args)))
+  #+native
+  (destructuring-bind (&key inputs outputs &allow-other-keys) key-args
+    (native-pointer
+     (intersection
+      (union inputs outputs) (arglist-plain-and-categories arglist nil))
+     (body-expand name arglist gsl-name c-arguments key-args))))
+
+(defun body-expand (name arglist gsl-name c-arguments key-args)
   "Expand the body (computational part) of the defmfun."
   (with-defmfun-key-args key-args
     (let* ((cret-type (if (member c-return *special-c-return*)
@@ -163,9 +174,7 @@
 	   allocated-decl
 	   (mapcar #'wfo-declare allocated-decl)
 	   'cffi:with-foreign-objects
-	   `(#-native
-	     ,@(mapcar (lambda (v) `(copy-cl-to-c ,v)) inputs)
-	     ,@before
+	   `(,@before
 	     (let ((,cret-name
 		    (cffi:foreign-funcall
 		     ,gsl-name
diff --git a/polynomial.lisp b/polynomial.lisp
index 1063aecc..18db9462 100644
--- a/polynomial.lisp
+++ b/polynomial.lisp
@@ -1,6 +1,6 @@
 ;; Polynomials
 ;; Liam Healy, Tue Mar 21 2006 - 18:33
-;; Time-stamp: <2008-12-30 21:33:37EST polynomial.lisp>
+;; Time-stamp: <2009-01-31 19:34:56EST polynomial.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -14,14 +14,50 @@
 ;;;; Polynomial Evaluation
 ;;;;****************************************************************************
 
-(defmfun polynomial-eval (coefficients x)
-  "gsl_poly_eval"
-  (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet) (x :double))
-  :inputs (coefficients)
+(defmfun evaluate
+    ((coefficients vector-double-float) (x float) &key divided-difference)
+  ("gsl_poly_eval" "gsl_poly_dd_eval")
+  ((((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet)
+    (x :double))
+   (((c-pointer divided-difference) :pointer)
+    ((c-pointer coefficients) :pointer)
+    ((dim0 coefficients) sizet)
+    (x :double)))
+  :definition :method
+  :inputs (coefficients divided-difference)
   :c-return :double
   :documentation			; FDL
   "Evaluate the polyonomial with coefficients at the point x.")
 
+#|
+;;; These won't work until we can take returned complex (structs).
+(defmfun evaluate
+    ((coefficients vector-double-float) (x complex)
+     &key &allow-other-keys)
+  "gsl_poly_complex_eval"
+  (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet)
+   (x complex-double-c))
+  :definition :method
+  :gsl-version (1 11)
+  :inputs (coefficients)
+  :c-return complex-double-c
+  :documentation			; FDL
+  "Evaluate the polyonomial with coefficients at the complex value x.")
+
+(defmfun evaluate
+    ((coefficients vector-complex-double-float) (x complex)
+     &key &allow-other-keys)
+  "gsl_complex_poly_complex_eval"
+  (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet)
+   (x complex-double-c))
+  :definition :method
+  :gsl-version (1 11)
+  :inputs (coefficients)
+  :c-return complex-double-c
+  :documentation			; FDL
+  "Evaluate the polyonomial with coefficients at the complex value x.")
+|#
+
 ;;;;****************************************************************************
 ;;;; Divided Difference Representation of Polynomials
 ;;;;****************************************************************************
@@ -41,18 +77,6 @@
    divided-differences of (xa,ya) are stored in the array
    dd, of the same length.")
 
-(defmfun polynomial-eval-divided-difference (dd xa x)
-  "gsl_poly_dd_eval"
-  (((c-pointer dd) :pointer)
-   ((c-pointer xa) :pointer)
-   ((dim0 xa) sizet)
-   (x :double))
-  :inputs (dd xa)
-  :c-return :double
-  :documentation			; FDL
-  "Evaluate the polynomial stored in divided-difference form
-   in the arrays dd and xa at the point x.")
-
 (defmfun taylor-divided-difference (coefs xp dd xa workspace)
   "gsl_poly_dd_taylor"
   (((c-pointer coefs) :pointer)
@@ -154,16 +178,16 @@
 
 (save-test polynomial
  (let ((xa #m(0.0d0 1.0d0 2.0d0 3.0d0))
-	(ya #m(2.5d0 7.2d0 32.7d0 91.0d0))
-	(dd (make-marray 'double-float :dimensions 4)))
+       (ya #m(2.5d0 7.2d0 32.7d0 91.0d0))
+       (dd (make-marray 'double-float :dimensions 4)))
    (divided-difference dd xa ya)
    (list
-    (polynomial-eval-divided-difference dd xa 0.0d0)
-    (polynomial-eval-divided-difference dd xa 1.0d0)
-    (polynomial-eval-divided-difference dd xa 2.0d0)
-    (polynomial-eval-divided-difference dd xa 3.0d0)))
+    (evaluate xa 0.0d0 :divided-difference dd)
+    (evaluate xa 1.0d0 :divided-difference dd)
+    (evaluate xa 2.0d0 :divided-difference dd)
+    (evaluate xa 3.0d0 :divided-difference dd)))
  (let ((vec #m(1.0d0 2.0d0 3.0d0)))
-   (polynomial-eval vec -1.0d0))
+   (evaluate vec -1.0d0))
  (solve-quadratic 1.0d0 0.0d0 1.0d0)
  (solve-quadratic 1.0d0 -2.0d0 1.0d0)
  (solve-quadratic-complex 1.0d0 -2.0d0 1.0d0)
diff --git a/tests/polynomial.lisp b/tests/polynomial.lisp
index d71a547a..63f3263e 100644
--- a/tests/polynomial.lisp
+++ b/tests/polynomial.lisp
@@ -3,70 +3,68 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST POLYNOMIAL
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST (LIST 2.5d0 7.2d0 32.7d0 91.0d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((XA
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                             '(0.0d0 1.0d0 2.0d0 3.0d0)))
-                               (YA
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                             '(2.5d0 7.2d0 32.7d0 91.0d0)))
-                               (DD (MAKE-MARRAY 'DOUBLE-FLOAT :DIMENSIONS 4)))
-                           (DIVIDED-DIFFERENCE DD XA YA)
-                           (LIST
-                            (POLYNOMIAL-EVAL-DIVIDED-DIFFERENCE DD XA 0.0d0)
-                            (POLYNOMIAL-EVAL-DIVIDED-DIFFERENCE DD XA 1.0d0)
-                            (POLYNOMIAL-EVAL-DIVIDED-DIFFERENCE DD XA 2.0d0)
-                            (POLYNOMIAL-EVAL-DIVIDED-DIFFERENCE DD XA
-                                                                3.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 2.0d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (LET ((VEC
-                                                                  (MAKE-MARRAY
-                                                                   'DOUBLE-FLOAT
-                                                                   :INITIAL-CONTENTS
-                                                                   '(1.0d0
-                                                                     2.0d0
-                                                                     3.0d0))))
-                                                             (POLYNOMIAL-EVAL
-                                                              VEC -1.0d0))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST (LIST) (LIST))
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (SOLVE-QUADRATIC
-                                                            1.0d0 0.0d0
-                                                            1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.0d0 1.0d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (SOLVE-QUADRATIC
-                                                            1.0d0 -2.0d0
-                                                            1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST #C(1.0d0 0.0d0) #C(1.0d0 0.0d0))
-                        (MULTIPLE-VALUE-LIST
-                         (SOLVE-QUADRATIC-COMPLEX 1.0d0 -2.0d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST -3.000000000000001d0 1.9999999999999996d0
-                              7.000000000000001d0)
-                        (MULTIPLE-VALUE-LIST
-                         (SOLVE-CUBIC -6.0d0 -13.0d0 42.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST #C(-5.551115123125783d-17 -0.9999999999999999d0)
-                              #C(-5.551115123125783d-17 0.9999999999999999d0)
-                              #C(1.0d0 0.0d0))
-                        (MULTIPLE-VALUE-LIST
-                         (SOLVE-CUBIC-COMPLEX -1.0d0 1.0d0 -1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         #(#C(-0.8090169943749477d0 0.5877852522924734d0)
-                           #C(-0.8090169943749477d0 -0.5877852522924734d0)
-                           #C(0.3090169943749475d0 0.951056516295153d0)
-                           #C(0.3090169943749475d0 -0.951056516295153d0)
-                           #C(0.9999999999999999d0 0.0d0)))
-                        (MULTIPLE-VALUE-LIST
-                         (CL-ARRAY
-                          (POLYNOMIAL-SOLVE
-                           (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                        '(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0
-                                          1.0d0)))))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST (LIST 2.5d0 7.2d0 32.7d0 91.0d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((XA
+	   (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+			'(0.0d0 1.0d0 2.0d0 3.0d0)))
+	  (YA
+	   (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+			'(2.5d0 7.2d0 32.7d0 91.0d0)))
+	  (DD (MAKE-MARRAY 'DOUBLE-FLOAT :DIMENSIONS 4)))
+      (DIVIDED-DIFFERENCE DD XA YA)
+      (LIST (EVALUATE XA 0.0d0 :DIVIDED-DIFFERENCE DD)
+	    (EVALUATE XA 1.0d0 :DIVIDED-DIFFERENCE DD)
+	    (EVALUATE XA 2.0d0 :DIVIDED-DIFFERENCE DD)
+	    (EVALUATE XA 3.0d0 :DIVIDED-DIFFERENCE DD)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 2.0d0)
+				     (MULTIPLE-VALUE-LIST
+				      (LET ((VEC
+					     (MAKE-MARRAY
+					      'DOUBLE-FLOAT
+					      :INITIAL-CONTENTS
+					      '(1.0d0
+						2.0d0
+						3.0d0))))
+					(EVALUATE
+					 VEC -1.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST (LIST) (LIST))
+				     (MULTIPLE-VALUE-LIST
+				      (SOLVE-QUADRATIC
+				       1.0d0 0.0d0
+				       1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.0d0 1.0d0)
+				     (MULTIPLE-VALUE-LIST
+				      (SOLVE-QUADRATIC
+				       1.0d0 -2.0d0
+				       1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #C(1.0d0 0.0d0) #C(1.0d0 0.0d0))
+   (MULTIPLE-VALUE-LIST
+    (SOLVE-QUADRATIC-COMPLEX 1.0d0 -2.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -3.000000000000001d0 1.9999999999999996d0
+	 7.000000000000001d0)
+   (MULTIPLE-VALUE-LIST
+    (SOLVE-CUBIC -6.0d0 -13.0d0 42.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #C(-5.551115123125783d-17 -0.9999999999999999d0)
+	 #C(-5.551115123125783d-17 0.9999999999999999d0)
+	 #C(1.0d0 0.0d0))
+   (MULTIPLE-VALUE-LIST
+    (SOLVE-CUBIC-COMPLEX -1.0d0 1.0d0 -1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(#C(-0.8090169943749477d0 0.5877852522924734d0)
+      #C(-0.8090169943749477d0 -0.5877852522924734d0)
+      #C(0.3090169943749475d0 0.951056516295153d0)
+      #C(0.3090169943749475d0 -0.951056516295153d0)
+      #C(0.9999999999999999d0 0.0d0)))
+   (MULTIPLE-VALUE-LIST
+    (CL-ARRAY
+     (POLYNOMIAL-SOLVE
+      (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+		   '(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0
+		     1.0d0)))))))
 
-- 
GitLab