Commit cbcc15dd authored by Liam Healy's avatar Liam Healy

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.
parent 04f30601
;; 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
......
;; 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)
......
......@@ -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)))))))
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