Commit a26ad0bb authored by Raymond Toy's avatar Raymond Toy
Browse files

Merge branch 'issue-18-more-tests-incomplete-gamma' into 'master'

Fix #18: Add more tests incomplete gamma

Closes #18

See merge request !19
parents 0f3c4ea4 b9c5af9a
Pipeline #2433 passed with stage
in 4 minutes and 54 seconds
......@@ -400,14 +400,44 @@
(setf ans (add-qd-d ans p)))))
(defun make-float (sign int-part frac-part scale exp)
"Create a quad-float given the SIGN, INT-PART, FRAC-PART, SCALE, and EXP.
For example, if we're reading #q-123.45678q100, MAKE-FLOAT is called like
(MAKE-FLOAT -1 123 45678 5 100)"
(declare (type (member -1 1) sign)
(type unsigned-byte int-part frac-part)
(fixnum scale exp))
(let ((qd (rational-to-qd (* (+ int-part (/ frac-part (expt 10 scale)))
(expt 10 exp)))))
(if (minusp sign)
(neg-qd qd)
qd)))
;; Quick exit when the int part and frac part are zero.
(when (and (zerop int-part) (zerop frac-part))
(return-from make-float (make-qd-d (* sign 0d0))))
;; Let i, f, s, e be int-part, frac-part, scale and exp, respectively. Compute the
;; rational number
;;
;; r = (i+(f/10^s))*10^e = (10^s*i+f)/10^s*10^e
;;
;; Take base 10 log:
;;
;; log10(r) = log10(10^s*i+f) + e - s
;;
;; This will underflow if log10(r) <=
;; log10(least-positive-double-float) and overflow if log10(r) >=
;; log10(most-positive-double-float). Let's catch these cases early
;; so that we don't end up computing huge rationals for 10^(e-s).
(let* ((log-res (+ (log (+ frac-part (* int-part (expt 10 scale))) 10)
(- exp scale)))
(qd
(cond ((< log-res (log least-positive-double-float 10))
;; Underflow. Silently set to 0.
(rational-to-qd 0))
((> log-res (log most-positive-double-float 10))
;; Overflow.
(error 'floating-point-overflow))
(t
(rational-to-qd (* (+ int-part (/ frac-part (expt 10 scale)))
(expt 10 exp)))))))
(if (minusp sign)
(neg-qd qd)
qd)))
;; This seems to work, but really needs to be rewritten!
(defun read-qd (stream)
......
......@@ -188,6 +188,7 @@
add-dd-qd
add-qd add-qd-t
add-qd-d-t
neg-qd-t
neg-qd
sub-qd sub-qd-dd sub-qd-d sub-d-qd
mul-qd-d mul-qd-dd
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -91,7 +91,13 @@
bessel-j.neg-order.d.4
bessel-j-complex-arg.d.5
bessel-j-1/2-complex.d.1
bessel-j-1/2-complex.q.1))
bessel-j-1/2-complex.q.1
incomplete-gamma-small-upper-normalized
incomplete-gamma-small-lower-normalized
incomplete-gamma-large-upper
incomplete-gamma-large-upper-normalized
incomplete-gamma-large-lower
incomplete-gamma-large-lower-normalized))
;;; Some simple tests from the Yozo Hida's qd package.
......@@ -1844,3 +1850,242 @@
(result (incomplete-gamma-upper-normalized #q200 #q100)))
(check-accuracy 211.99 result true))
nil)
;; More tests using data from BOOST.
(load "tests/data/incomplete-gamma-medium-data.lisp")
(rt:deftest incomplete-gamma-medium-upper
(let* ((test-data incomplete-gamma-medium-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for inc-upper = (aref test-data k 2)
for test-inc-upper =
(check-accuracy 196.15 (incomplete-gamma-upper a z) inc-upper)
when test-inc-upper append (list (list (list k a z) test-inc-upper))))
nil)
(rt:deftest incomplete-gamma-medium-lower
(let* ((test-data incomplete-gamma-medium-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for inc-upper = (aref test-data k 4)
for test-inc-upper =
(check-accuracy 195.83 (incomplete-gamma-lower a z) inc-upper)
when test-inc-upper append (list (list (list k a z) test-inc-upper))))
nil)
(rt:deftest incomplete-gamma-medium-lower-normalized
(let* ((test-data incomplete-gamma-medium-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for inc-upper = (aref test-data k 5)
for test-inc-upper =
(check-accuracy 196.84 (incomplete-gamma-lower a z :normalized-p t) inc-upper)
when test-inc-upper append (list (list (list k a z) test-inc-upper))))
nil)
(load "tests/data/incomplete-gamma-integer-data.lisp")
(rt:deftest incomplete-gamma-integer-upper
(let* ((test-data incomplete-gamma-integer-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 2)
for result = (check-accuracy 196.08 (incomplete-gamma-upper a z) true)
when result
append (list (list (list k a z) result))))
nil)
(rt:deftest incomplete-gamma-integer-upper-normalized
(let* ((test-data incomplete-gamma-integer-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 3)
for result = (check-accuracy 197.62 (incomplete-gamma-upper a z :normalized-p t) true)
when result
append (list (list (list k a z) result))))
nil)
(rt:deftest incomplete-gamma-integer-lower
(let* ((test-data incomplete-gamma-integer-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 4)
for result = (check-accuracy 196.51 (incomplete-gamma-lower a z) true)
when result
append (list (list (list k a z) result))))
nil)
(rt:deftest incomplete-gamma-integer-lower-normalized
(let* ((test-data incomplete-gamma-integer-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 5)
for result = (check-accuracy 197.66 (incomplete-gamma-lower a z :normalized-p t) true)
when result
append (list (list (list k a z) result))))
nil)
(load "tests/data/incomplete-gamma-small-data.lisp")
(rt:deftest incomplete-gamma-small-upper
(let* ((test-data incomplete-gamma-small-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 2)
for result = (check-accuracy 162.37 (incomplete-gamma-upper a z) true)
when result
append (list (list (list k a z) result))))
nil)
;; This test fails with errors because the continued fraction fails to
;; converge.
(rt:deftest incomplete-gamma-small-upper-normalized
(let* ((test-data incomplete-gamma-small-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 3)
for result = (multiple-value-bind (ok error)
(ignore-errors (incomplete-gamma-upper a z :normalized-p t))
(if ok
(check-accuracy 300 ok true)
error))
when result
append (list (list (list k a z) result))))
nil)
(rt:deftest incomplete-gamma-small-lower
(let* ((test-data incomplete-gamma-small-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 4)
for result = (check-accuracy 196.68 (incomplete-gamma-lower a z) true)
when result
append (list (list (list k a z) result))))
nil)
;; This test fails with errors because the continued fraction fails to
;; converge.
(rt:deftest incomplete-gamma-small-lower-normalized
(let* ((test-data incomplete-gamma-small-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 5)
for actual = (multiple-value-list (ignore-errors (incomplete-gamma-lower a z :normalized-p t)))
for result = (multiple-value-bind (ok error)
(ignore-errors (incomplete-gamma-lower a z :normalized-p t))
(if ok
(check-accuracy 300 ok true)
error))
when result
append (list (list (list k a z) result))))
nil)
(load "tests/data/incomplete-gamma-large-data.lisp")
(rt:deftest incomplete-gamma-large-upper
(let* ((test-data incomplete-gamma-large-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 2)
for result = (progn
(multiple-value-bind (ok error)
(ignore-errors (incomplete-gamma-upper a z))
(if ok
(check-accuracy 188.07 ok true)
error)))
when result
append (list (list (list k a z) result))))
nil)
(rt:deftest incomplete-gamma-large-upper-normalized
(let* ((test-data incomplete-gamma-large-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 3)
for result = (progn
(multiple-value-bind (ok error)
(ignore-errors (incomplete-gamma-upper a z :normalized-p t))
(if ok
(check-accuracy 500 ok true)
error)))
when result
append (list (list (list k a z) result))))
nil)
(rt:deftest incomplete-gamma-large-lower
(let* ((test-data incomplete-gamma-large-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 4)
for result = (progn
(multiple-value-bind (ok error)
(ignore-errors (incomplete-gamma-lower a z))
(if ok
(check-accuracy 188.07 ok true)
error)))
when result
append (list (list (list k a z) result))))
nil)
(rt:deftest incomplete-gamma-large-lower-normalized
(let* ((test-data incomplete-gamma-large-data)
(n-tests (array-dimension test-data 0)))
;;
(loop for k from 0 below n-tests
for a = (aref test-data k 0)
for z = (aref test-data k 1)
for true = (aref test-data k 5)
for result = (progn
(multiple-value-bind (ok error)
(ignore-errors (incomplete-gamma-lower a z :normalized-p t))
(if ok
(check-accuracy 500 ok true)
error)))
when result
append (list (list (list k a z) result))))
nil)
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