Commit 7aa0111c authored by Liam Healy's avatar Liam Healy

Replace 'pi with 'dpi, the double-float version

The standard declares that 'pi is an implementation-dependent
long-float, but GSL requires double float, so the constant 'dpi has
been defined as (coerce pi 'double-float).  This should make CLISP at
the least happier.
parent 7591d6dd
;; Monte Carlo Integration
;; Liam Healy Sat Feb 3 2007 - 17:42
;; Time-stamp: <2010-06-29 22:15:24EDT monte-carlo.lisp>
;; Time-stamp: <2010-08-07 21:41:14EDT monte-carlo.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -205,13 +205,14 @@
(defun mcrw (x y z)
"Example function for Monte Carlo used in random walk studies."
(* (/ (expt pi 3))
(* (/ (expt dpi 3))
(/ (- 1 (* (cos x) (cos y) (cos z))))))
(defparameter *mc-lower* #m(0.0d0 0.0d0 0.0d0))
(defparameter *mc-upper*
(grid:make-foreign-array 'double-float :initial-contents (list pi pi pi)))
(grid:make-foreign-array
'double-float :initial-contents (make-list 3 :initial-element dpi)))
(defun random-walk-plain-example (&optional (nsamples 500000))
(monte-carlo-integrate-plain 'mcrw *mc-lower* *mc-upper* nsamples))
......
;; Numerical integration techniques that require tables
;; Liam Healy 2009-04-04 15:24:05EDT
;; Time-stamp: <2009-12-27 09:42:09EST numerical-integration-with-tables.lisp>
;; Time-stamp: <2010-08-07 21:41:16EDT numerical-integration-with-tables.lisp>
;;
;; Copyright 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -256,11 +256,11 @@
;; point overflow in making the table. Setting to 100 works and
;; gives the result given in test.c.
(integration-QAWO
'integration-test-f456 0.0d0 (* 10 pi) 1.0d0 :sine 100 0.0d0 1.0d-7)
'integration-test-f456 0.0d0 (* 10 dpi) 1.0d0 :sine 100 0.0d0 1.0d-7)
;; For the negative interval, it was necessary to have the function
;; takethe absolute value of the argument and then transfer the sign of
;; x to the answer.
(integration-QAWO
'integration-test-f456 0.0d0 (* 10 pi) -1.0d0 :sine 100 0.0d0 1.0d-7)
'integration-test-f456 0.0d0 (* 10 dpi) -1.0d0 :sine 100 0.0d0 1.0d-7)
(integration-QAWF
'integration-test-f457 0.0d0 (/ pi 2) 1.0d0 :cosine 1000 1.0d-7))
'integration-test-f457 0.0d0 (/ dpi 2) 1.0d0 :cosine 1000 1.0d-7))
;; Numerical integration
;; Liam Healy, Wed Jul 5 2006 - 23:14
;; Time-stamp: <2010-07-07 14:21:55EDT numerical-integration.lisp>
;; Time-stamp: <2010-08-07 21:46:21EDT numerical-integration.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Copyright 2006, 2007, 2008, 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -321,9 +321,9 @@
(lambda (x) (exp (* alpha x))))
(save-test numerical-integration
(integration-qng 'sin 0.0d0 pi)
(integration-QAG 'sin 0.0d0 pi :gauss15 20)
(integration-QAG 'sin 0.0d0 pi :gauss21 40)
(integration-qng 'sin 0.0d0 dpi)
(integration-QAG 'sin 0.0d0 dpi :gauss15 20)
(integration-QAG 'sin 0.0d0 dpi :gauss21 40)
;; Tests from gsl-1.11/integration/test.c
;; Functions defined in gsl-1.11/integration/tests.c
(integration-QNG (integration-test-f1 2.6d0) 0.0d0 1.0d0 0.1d0 0.0d0)
......
;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18
;; Time-stamp: <2010-06-30 19:57:28EDT spline-example.lisp>
;; Time-stamp: <2010-08-07 21:37:07EDT spline-example.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -43,7 +43,7 @@
over a 0-2pi interval and interpolated with
+periodic-cubic-spline-interpolation+"
(let* ((xarr
(loop with step = (/ (* 2.0 pi) intervals)
(loop with step = (/ (* 2.0 dpi) intervals)
for i from 0 upto intervals
collect (* i step)))
(xmarr (grid:make-foreign-array 'double-float :initial-contents xarr))
......@@ -51,9 +51,7 @@
(grid:make-foreign-array 'double-float :initial-contents (mapcar 'sin xarr)))
(spline
(make-spline +periodic-cubic-spline-interpolation+ xmarr ymarr)))
(evaluate-integral spline 0d0
#+clisp (coerce pi 'double-float) ;; pi=3.14..L0 on clisp
#-clisp pi)))
(evaluate-integral spline 0d0 dpi)))
(save-test interpolation
(spline-example 0.1d0)
......
;; Mathematical functions
;; Liam Healy, Wed Mar 8 2006 - 22:09
;; Time-stamp: <2010-07-20 16:56:44EDT mathematical.lisp>
;; Time-stamp: <2010-08-07 21:35:54EDT mathematical.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -196,6 +196,12 @@ Function: double gsl_pow_9 (const double x)
returns +1."
:return ((pmnil cr)))
;;;;****************************************************************************
;;;; pi
;;;;****************************************************************************
(defconstant dpi (coerce pi 'double-float))
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
......
;; Series acceleration.
;; Liam Healy, Wed Nov 21 2007 - 18:41
;; Time-stamp: <2010-06-30 19:57:28EDT series-acceleration.lisp>
;; Time-stamp: <2010-08-07 21:41:14EDT series-acceleration.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -91,7 +91,7 @@
(defun acceleration-example (&optional (print-explanation t))
(let ((maxterms 20)
(sum 0.0d0)
(zeta2 (/ (expt pi 2) 6)))
(zeta2 (/ (expt dpi 2) 6)))
(let ((levin (make-levin maxterms))
(array (grid:make-foreign-array 'double-float :dimensions maxterms)))
(dotimes (n maxterms)
......
;; Univariate minimization
;; Liam Healy Tue Jan 8 2008 - 21:02
;; Time-stamp: <2010-07-01 11:55:38EDT minimization-one.lisp>
;; Time-stamp: <2010-08-07 21:41:13EDT minimization-one.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -226,14 +226,14 @@
(when print-steps
(format t "~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g~&"
iter lower upper
min (- min pi)
min (- min dpi)
(- upper lower)))
while (and (< iter max-iter)
;; abs and rel error swapped in example?
(not (min-test-interval lower upper 0.001d0 0.0d0)))
finally
(return (values iter lower upper min (- min pi) (- upper lower)))))))
(return (values iter lower upper min (- min dpi) (- upper lower)))))))
(save-test minimization-one
(minimization-one-example +brent-fminimizer+ nil nil)
......
;; Jacobian elliptic functions
;; Liam Healy, Mon Mar 20 2006 - 22:21
;; Time-stamp: <2010-04-14 13:06:10EDT elliptic-functions.lisp>
;; Time-stamp: <2010-08-07 21:41:12EDT elliptic-functions.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -30,7 +30,7 @@
"The Jacobian elliptic functions sn(u|m),
cn(u|m), dn(u|m) computed by descending Landen transformations.")
(defvar *elljac-K* (/ pi (* 2.0d0 0.9741726903999478375938128316d0)))
(defvar *elljac-K* (/ dpi (* 2.0d0 0.9741726903999478375938128316d0)))
(defvar *elljac-A* (/ (sqrt (1+ (sqrt 0.9d0)))))
(defvar *elljac-B* (/ (expt 0.9d0 1/4) (sqrt (1+ (sqrt 0.9d0)))))
(defvar *elljac-C* (expt 0.9d0 1/4))
......
;; Mathieu functions
;; Liam Healy 2009-02-16 16:30:59EST mathieu.lisp
;; Time-stamp: <2010-07-15 10:35:18EDT mathieu.lisp>
;; Time-stamp: <2010-08-07 21:38:53EDT mathieu.lisp>
;;
;; Copyright 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -219,77 +219,77 @@
;;; Tests from gsl-1.11/specfunc/test_mathieu.c
(save-test mathieu
(mathieu-ce 0 0.0d0 0.0d0)
(mathieu-ce 0 0.0d0 (/ pi 2))
(mathieu-ce 0 0.0d0 (/ dpi 2))
(mathieu-ce 0 5.0d0 0.0d0)
(mathieu-ce 0 5.0d0 (/ pi 2))
(mathieu-ce 0 5.0d0 (/ dpi 2))
(mathieu-ce 0 10.0d0 0.0d0)
(mathieu-ce 0 10.0d0 (/ pi 2))
(mathieu-ce 0 10.0d0 (/ dpi 2))
(mathieu-ce 0 15.0d0 0.0d0)
(mathieu-ce 0 15.0d0 (/ pi 2))
(mathieu-ce 0 15.0d0 (/ dpi 2))
(mathieu-ce 0 20.0d0 0.0d0)
(mathieu-ce 0 20.0d0 (/ pi 2))
(mathieu-ce 0 20.0d0 (/ dpi 2))
(mathieu-ce 0 25.0d0 0.0d0)
(mathieu-ce 0 25.0d0 (/ pi 2))
(mathieu-ce 0 25.0d0 (/ dpi 2))
(mathieu-ce 1 0.0d0 0.0d0)
(mathieu-ce 1 5.0d0 0.0d0)
(mathieu-ce 1 10.0d0 0.0d0)
(mathieu-ce 1 15.0d0 0.0d0)
(mathieu-ce 1 20.0d0 0.0d0)
(mathieu-ce 1 25.0d0 0.0d0)
(mathieu-se 1 0.0d0 (/ pi 2))
(mathieu-se 1 5.0d0 (/ pi 2))
(mathieu-se 1 10.0d0 (/ pi 2))
(mathieu-se 1 15.0d0 (/ pi 2))
(mathieu-se 1 20.0d0 (/ pi 2))
(mathieu-se 1 25.0d0 (/ pi 2))
(mathieu-se 1 0.0d0 (/ dpi 2))
(mathieu-se 1 5.0d0 (/ dpi 2))
(mathieu-se 1 10.0d0 (/ dpi 2))
(mathieu-se 1 15.0d0 (/ dpi 2))
(mathieu-se 1 20.0d0 (/ dpi 2))
(mathieu-se 1 25.0d0 (/ dpi 2))
(mathieu-ce 2 0.0d0 0.0d0)
(mathieu-ce 2 0.0d0 (/ pi 2))
(mathieu-ce 2 0.0d0 (/ dpi 2))
(mathieu-ce 2 5.0d0 0.0d0)
(mathieu-ce 2 5.0d0 (/ pi 2))
(mathieu-ce 2 5.0d0 (/ dpi 2))
(mathieu-ce 2 10.0d0 0.0d0)
(mathieu-ce 2 10.0d0 (/ pi 2))
(mathieu-ce 2 10.0d0 (/ dpi 2))
(mathieu-ce 2 15.0d0 0.0d0)
(mathieu-ce 2 15.0d0 (/ pi 2))
(mathieu-ce 2 15.0d0 (/ dpi 2))
(mathieu-ce 2 20.0d0 0.0d0)
(mathieu-ce 2 20.0d0 (/ pi 2))
(mathieu-ce 2 20.0d0 (/ dpi 2))
(mathieu-ce 2 25.0d0 0.0d0)
(mathieu-ce 2 25.0d0 (/ pi 2))
(mathieu-ce 2 25.0d0 (/ dpi 2))
(mathieu-ce 5 0.0d0 0.0d0)
(mathieu-ce 5 5.0d0 0.0d0)
(mathieu-ce 5 10.0d0 0.0d0)
(mathieu-ce 5 15.0d0 0.0d0)
(mathieu-ce 5 20.0d0 0.0d0)
(mathieu-ce 5 25.0d0 0.0d0)
(mathieu-se 5 0.0d0 (/ pi 2))
(mathieu-se 5 5.0d0 (/ pi 2))
(mathieu-se 5 10.0d0 (/ pi 2))
(mathieu-se 5 15.0d0 (/ pi 2))
(mathieu-se 5 20.0d0 (/ pi 2))
(mathieu-se 5 25.0d0 (/ pi 2))
(mathieu-se 5 0.0d0 (/ dpi 2))
(mathieu-se 5 5.0d0 (/ dpi 2))
(mathieu-se 5 10.0d0 (/ dpi 2))
(mathieu-se 5 15.0d0 (/ dpi 2))
(mathieu-se 5 20.0d0 (/ dpi 2))
(mathieu-se 5 25.0d0 (/ dpi 2))
(mathieu-ce 10 0.0d0 0.0d0)
(mathieu-ce 10 0.0d0 (/ pi 2))
(mathieu-ce 10 0.0d0 (/ dpi 2))
(mathieu-ce 10 5.0d0 0.0d0)
(mathieu-ce 10 5.0d0 (/ pi 2))
(mathieu-ce 10 5.0d0 (/ dpi 2))
(mathieu-ce 10 10.0d0 0.0d0)
(mathieu-ce 10 10.0d0 (/ pi 2))
(mathieu-ce 10 10.0d0 (/ dpi 2))
(mathieu-ce 10 15.0d0 0.0d0)
(mathieu-ce 10 15.0d0 (/ pi 2))
(mathieu-ce 10 15.0d0 (/ dpi 2))
(mathieu-ce 10 20.0d0 0.0d0)
(mathieu-ce 10 20.0d0 (/ pi 2))
(mathieu-ce 10 20.0d0 (/ dpi 2))
(mathieu-ce 10 25.0d0 0.0d0)
(mathieu-ce 10 25.0d0 (/ pi 2))
(mathieu-ce 10 25.0d0 (/ dpi 2))
(mathieu-ce 15 0.0d0 0.0d0)
(mathieu-ce 15 5.0d0 0.0d0)
(mathieu-ce 15 10.0d0 0.0d0)
(mathieu-ce 15 15.0d0 0.0d0)
(mathieu-ce 15 20.0d0 0.0d0)
(mathieu-ce 15 25.0d0 0.0d0)
(mathieu-se 15 0.0d0 (/ pi 2))
(mathieu-se 15 5.0d0 (/ pi 2))
(mathieu-se 15 10.0d0 (/ pi 2))
(mathieu-se 15 15.0d0 (/ pi 2))
(mathieu-se 15 20.0d0 (/ pi 2))
(mathieu-se 15 25.0d0 (/ pi 2))
(grid:copy-to (mathieu-ce-array 0.0d0 (/ pi 2) 6))
(mathieu-se 15 0.0d0 (/ dpi 2))
(mathieu-se 15 5.0d0 (/ dpi 2))
(mathieu-se 15 10.0d0 (/ dpi 2))
(mathieu-se 15 15.0d0 (/ dpi 2))
(mathieu-se 15 20.0d0 (/ dpi 2))
(mathieu-se 15 25.0d0 (/ dpi 2))
(grid:copy-to (mathieu-ce-array 0.0d0 (/ dpi 2) 6))
(grid:copy-to (mathieu-ce-array 20.0d0 0.0d0 16))
(grid:copy-to (mathieu-se-array 20.0d0 (/ pi 2) 15 1)))
(grid:copy-to (mathieu-se-array 20.0d0 (/ dpi 2) 15 1)))
This diff is collapsed.
;; Regression test NUMERICAL-INTEGRATION for GSLL, automatically generated
;;
;; Copyright 2009 Liam M. Healy
;; Copyright 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -54,27 +54,27 @@
(LIST -0.12813684839916742d0 6.875028324415666d-12)
(MULTIPLE-VALUE-LIST
(INTEGRATION-QAWO 'INTEGRATION-TEST-F456 0.0d0
(* 10 PI) 1.0d0 :SINE 100 0.0d0
(* 10 dPI) 1.0d0 :SINE 100 0.0d0
1.d-7)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST 0.12813684839916742d0 6.875028324415666d-12)
(MULTIPLE-VALUE-LIST
(INTEGRATION-QAWO 'INTEGRATION-TEST-F456 0.0d0
(* 10 PI) -1.0d0 :SINE 100 0.0d0
(* 10 dPI) -1.0d0 :SINE 100 0.0d0
1.d-7)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST 0.9999999999279803d0 1.556289974669056d-8)
(MULTIPLE-VALUE-LIST
(INTEGRATION-QAWF 'INTEGRATION-TEST-F457 0.0d0
(/ PI 2) 1.0d0 :COSINE 1000 1.d-7)))
(/ dPI 2) 1.0d0 :COSINE 1000 1.d-7)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST 2.0d0 2.220446049250313d-14 21)
(MULTIPLE-VALUE-LIST (INTEGRATION-QNG 'SIN 0.0d0 PI)))
(MULTIPLE-VALUE-LIST (INTEGRATION-QNG 'SIN 0.0d0 dPI)))
(LISP-UNIT:ASSERT-ERROR 'TYPE-ERROR
(INTEGRATION-QAG 'SIN 0.0d0 PI
(INTEGRATION-QAG 'SIN 0.0d0 dPI
:GAUSS15 20))
(LISP-UNIT:ASSERT-ERROR 'TYPE-ERROR
(INTEGRATION-QAG 'SIN 0.0d0 PI
(INTEGRATION-QAG 'SIN 0.0d0 dPI
:GAUSS21 40))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST 0.07716049379303085d0 9.424302194248481d-8 21)
......
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