Commit aa68fe67 authored by Liam M. Healy's avatar Liam M. Healy
Browse files

Correct return values from jacobian-elliptic-functions; add tests

The GSL function gsl_sf_elljac_e writes sn, cn, dn values to doubles
for which the pointer has been passed in, as shown in
/usr/include/gsl/gsl_sf_elljac.h, and not to sf-result as one might be
led to conclude from the "_e" at the end of the name.  This is now
fixed in GSLL, so that the tests, which sometimes passed and sometimes
didn't depending on what was run before it, now should always work.
The single non-error test been replaced by the 35 GSL tests in
gsl/specfunc/test_sf.c, and they all pass in SBCL Linux AMD64.
parent 49cb0c32
;; Jacobian elliptic functions
;; Liam Healy, Mon Mar 20 2006 - 22:21
;; Time-stamp: <2009-12-27 10:10:05EST elliptic-functions.lisp>
;; Time-stamp: <2010-04-14 13:06:10EDT elliptic-functions.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -20,25 +20,59 @@
(in-package :gsl)
;; /usr/include/gsl/gsl_sf_elljac.h
(defmfun jacobian-elliptic-functions (u m)
"gsl_sf_elljac_e"
((u :double) (m :double) (sn sf-result) (cn sf-result) (dn sf-result))
:return ((val sn) (val cn) (val dn) (err sn) (err cn) (err dn))
((u :double) (m :double)
(sn (:pointer :double)) (cn (:pointer :double)) (dn (:pointer :double)))
:documentation ; FDL
"The Jacobian elliptic functions sn(u|m),
cn(u|m), dn(u|m) computed by descending Landen transformations.")
;;; > (jacobian-elliptic-functions 0.61802d0 0.5d0)
;;; 0.564575752943391
;;; 0.8253812568676386
;;; 0.916857191493965
;;; > (jacobian-elliptic-functions 0.2d0 0.81d0)
;;; 0.19762082367187697
;;; 0.9802785369736752
;;; 0.9840560289645665
;;; > (jacobian-elliptic-functions 0.61802d0 1.5d0)
;;; ;;;error
(defvar *elljac-K* (/ pi (* 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))
(defvar *elljac-C2* (sqrt 0.9d0))
(save-test elliptic-functions
(jacobian-elliptic-functions 0.2d0 0.81d0)
(jacobian-elliptic-functions 0.61802d0 1.5d0))
(save-test
elliptic-functions
;; Tests from gsl/specfunc/test_sf.c
(jacobian-elliptic-functions 0.5d0 0.5d0)
(jacobian-elliptic-functions 1.0d0 0.3d0)
(jacobian-elliptic-functions 1.0d0 0.6d0)
(jacobian-elliptic-functions 2.0d0 0.999999d0)
(jacobian-elliptic-functions 1.69695970624443d0 0.270378013104138d0)
(jacobian-elliptic-functions 0.0d0 0.1d0)
(jacobian-elliptic-functions -1.0d-10 0.1d0)
(jacobian-elliptic-functions 1.0d-10 0.1d0)
(jacobian-elliptic-functions 1.0d-30 0.1d0)
(jacobian-elliptic-functions (- (/ *elljac-K* 2.0d0) 1.0d-10) 0.1d0)
(jacobian-elliptic-functions (/ *elljac-K* 2.0d0) 0.1d0)
(jacobian-elliptic-functions (+ (/ *elljac-K* 2.0d0) 1.0d-10) 0.1d0)
(jacobian-elliptic-functions (- *elljac-K* 1.0d-10) 0.1d0)
(jacobian-elliptic-functions *elljac-K* 0.1d0)
(jacobian-elliptic-functions (+ *elljac-K* 1.0d-10) 0.1d0)
(jacobian-elliptic-functions (* 3/2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (- (* 2 *elljac-K*) 1.0d-10) 0.1d0)
(jacobian-elliptic-functions (* 2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (+ (* 2 *elljac-K*) 1.0d-10) 0.1d0)
(jacobian-elliptic-functions (* 5/2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (- (* 3 *elljac-K*) 1.0d-10) 0.1d0)
(jacobian-elliptic-functions (* 3 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (+ (* 3 *elljac-K*) 1.0d-10) 0.1d0)
(jacobian-elliptic-functions (* 7/2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (- (* 4 *elljac-K*) 1.0d-10) 0.1d0)
(jacobian-elliptic-functions (* 4 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (* 9/2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (* -1/2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (* -1 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (* -3/2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (* -2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (* -5/2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (* -3 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (* -7/2 *elljac-K*) 0.1d0)
(jacobian-elliptic-functions (* -4 *elljac-K*) 0.1d0)
;; Signal an error
(jacobian-elliptic-functions 0.61802d0 1.5d0))
;; Regression test ELLIPTIC-FUNCTIONS for GSLL, automatically generated
;; Regression test ELLIPTIC-FUNCTIONS for GSLL
;;
;; 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
......@@ -18,13 +18,152 @@
(in-package :gsl)
(LISP-UNIT:DEFINE-TEST ELLIPTIC-FUNCTIONS
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.19762082367187703d0 0.9802785369736752d0
0.9840560289645665d0 0.0d0 0.0d0 0.0d0)
(MULTIPLE-VALUE-LIST
(JACOBIAN-ELLIPTIC-FUNCTIONS 0.2d0 0.81d0)))
(LISP-UNIT:ASSERT-ERROR 'INPUT-DOMAIN
(JACOBIAN-ELLIPTIC-FUNCTIONS
0.61802d0 1.5d0)))
(LISP-UNIT:DEFINE-TEST
ELLIPTIC-FUNCTIONS
(let ((lisp-unit:*epsilon* (* 2 2 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 0.4707504736556572833d0 0.8822663948904402865d0 0.9429724257773856873d0)
(multiple-value-list (jacobian-elliptic-functions 0.5d0 0.5d0))))
(let ((lisp-unit:*epsilon* (* 2 2 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 0.8187707145344889190d0 0.5741206467465548795d0 0.8938033089590823040d0)
(multiple-value-list (jacobian-elliptic-functions 1.0d0 0.3d0))))
(let ((lisp-unit:*epsilon* (* 2 2 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 0.7949388393365780943d0 0.6066895760718277578d0 0.7879361300438814425d0)
(multiple-value-list (jacobian-elliptic-functions 1.0d0 0.6d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 0.96402778575700186570d0 0.26580148285600686381d0 0.26580323105264131136d0)
(multiple-value-list (jacobian-elliptic-functions 2.0d0 0.999999d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
;; The first one should actually be TOL0
(list 1.0d0 0.0d0 0.8541791304497336d0)
(multiple-value-list (jacobian-elliptic-functions 1.69695970624443d0 0.270378013104138d0))))
(let ((lisp-unit:*epsilon* (* 2 2 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 0.0d0 1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions 0.0d0 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 2 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list -1.0d-10 1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions -1.0d-10 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 2 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 1.0d-10 1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions 1.0d-10 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 2 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 1.0d-30 1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions 1.0d-30 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 256 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list (- *elljac-A* (* 1.0d-10 *elljac-B* *elljac-C*))
(+ *elljac-B* (* 1.0d-10 *elljac-A* *elljac-C*))
(+ *elljac-C* (* 0.1d0 1.0d-10 *elljac-A* *elljac-B*)))
(multiple-value-list (jacobian-elliptic-functions (- (/ *elljac-K* 2.0d0) 1.0d-10) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 256 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list *elljac-A* *elljac-B* *elljac-C*)
(multiple-value-list (jacobian-elliptic-functions (/ *elljac-K* 2.0d0) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list (+ *elljac-A* (* 1.0d-10 *elljac-B* *elljac-C*))
(- *elljac-B* (* 1.0d-10 *elljac-A* *elljac-C*))
(- *elljac-C* (* 0.1d0 1.0d-10 *elljac-A* *elljac-B*)))
(multiple-value-list (jacobian-elliptic-functions (+ (/ *elljac-K* 2.0d0) 1.0d-10) 0.1d0))))
(let ((lisp-unit:*epsilon* 1.0d-5)) ; Varying tolerances, max = 10*TEST_SNGL
(lisp-unit::assert-numerical-equal
(list 1.0d0 (* 1.0d-10 *elljac-C2*) *elljac-C2* )
(multiple-value-list (jacobian-elliptic-functions (- *elljac-K* 1.0d-10) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 256 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 1.0d0 0.0d0 *elljac-C2*)
(multiple-value-list (jacobian-elliptic-functions *elljac-K* 0.1d0))))
(let ((lisp-unit:*epsilon* 1.0d-5)) ; Varying tolerances, max = 10*TEST_SNGL
(lisp-unit::assert-numerical-equal
(list 1.0d0 (* -1.0d-10 *elljac-C2*) *elljac-C2*)
(multiple-value-list (jacobian-elliptic-functions (+ *elljac-K* 1.0d-10) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 256 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list *elljac-A* (- *elljac-B*) *elljac-C*)
(multiple-value-list (jacobian-elliptic-functions (* 3/2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* 1.0d-5)) ; Varying tolerances, max = 10*TEST_SNGL
(lisp-unit::assert-numerical-equal
(list 1.0d-10 -1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions (- (* 2 *elljac-K*) 1.0d-10) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 0.0d0 -1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions (* 2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* 1.0d-5)) ; Varying tolerances, max = 10*TEST_SNGL
(lisp-unit::assert-numerical-equal
(list -1.0d-10 -1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions (+ (* 2 *elljac-K*) 1.0d-10) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list (- *elljac-A*) (- *elljac-B*) *elljac-C*)
(multiple-value-list (jacobian-elliptic-functions (* 5/2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* 1.0d-5)) ; Varying tolerances, max = 10*TEST_SNGL
(lisp-unit::assert-numerical-equal
(list -1.0d0 (* -1.0d-10 *elljac-C2*) *elljac-C2*)
(multiple-value-list (jacobian-elliptic-functions (- (* 3 *elljac-K*) 1.0d-10) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list -1.0d0 0.0d0 *elljac-C2*)
(multiple-value-list (jacobian-elliptic-functions (* 3 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* 1.0d-5)) ; Varying tolerances, max = 10*TEST_SNGL
(lisp-unit::assert-numerical-equal
(list -1.0d0 (* 1.0d-10 *elljac-C2*) *elljac-C2*)
(multiple-value-list (jacobian-elliptic-functions (+ (* 3 *elljac-K*) 1.0d-10) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list (- *elljac-A*) *elljac-B* *elljac-C*)
(multiple-value-list (jacobian-elliptic-functions (* 7/2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* 1.0d-5)) ; Varying tolerances, max = 10*TEST_SNGL
(lisp-unit::assert-numerical-equal
(list -1.0d-10 1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions (- (* 4 *elljac-K*) 1.0d-10) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 0.0d0 1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions (* 4 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list *elljac-A* *elljac-B* *elljac-C*)
(multiple-value-list (jacobian-elliptic-functions (* 9/2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list (- *elljac-A*) *elljac-B* *elljac-C*)
(multiple-value-list (jacobian-elliptic-functions (* -1/2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list -1.0d0 0.0d0 *elljac-C2*)
(multiple-value-list (jacobian-elliptic-functions (* -1 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list (- *elljac-A*) (- *elljac-B*) *elljac-C*)
(multiple-value-list (jacobian-elliptic-functions (* -3/2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 2 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 0.0d0 -1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions (* -2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list *elljac-A* (- *elljac-B*) *elljac-C*)
(multiple-value-list (jacobian-elliptic-functions (* -5/2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 1.0d0 0.0d0 *elljac-C2*)
(multiple-value-list (jacobian-elliptic-functions (* -3 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list *elljac-A* *elljac-B* *elljac-C*)
(multiple-value-list (jacobian-elliptic-functions (* -7/2 *elljac-K*) 0.1d0))))
(let ((lisp-unit:*epsilon* (* 2 2 double-float-epsilon)))
(lisp-unit::assert-numerical-equal
(list 0.0d0 1.0d0 1.0d0)
(multiple-value-list (jacobian-elliptic-functions (* -4 *elljac-K*) 0.1d0))))
(LISP-UNIT:ASSERT-ERROR 'INPUT-DOMAIN (JACOBIAN-ELLIPTIC-FUNCTIONS 0.61802d0 1.5d0)))
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