Menu
▾
▴
[Maxima-commits] Maxima, A Computer Algebra System branch master updated. 5.28.0-120-g400cb79
|
From: Raymond T. <rt...@us...> - 2012-11-15 03:32:16
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima, A Computer Algebra System".
The branch, master has been updated
via 400cb7904ee5ecb1f7b3c790c8c66e453dcd0f95 (commit)
via fee6d7400a72bf2809a014d64e60004d1c05d4b0 (commit)
via 2284ca4017b07e9c5551e86b7f17d9ae65f233af (commit)
via cd29013f2f1eb364493ac547365317c82c6618d8 (commit)
from e7f991bca3b4e287b164f84a702bb29e69cae393 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 400cb7904ee5ecb1f7b3c790c8c66e453dcd0f95
Merge: fee6d74 e7f991b
Author: Raymond Toy <toy...@gm...>
Date: Wed Nov 14 19:31:59 2012 -0800
Merge branch 'master' of ssh://maxima.git.sourceforge.net/gitroot/maxima/maxima
commit fee6d7400a72bf2809a014d64e60004d1c05d4b0
Author: Raymond Toy <toy...@gm...>
Date: Wed Nov 14 19:31:24 2012 -0800
Fix Bug 3587304 erfc(x) for x > 6 is wrong
src/gamma.lisp:
o Add bf-erfc to evaluate erfc more accurately for large x.
o Use it in simp-erfc.
tests/rtest_gamma.mac
o Add tests.
diff --git a/src/gamma.lisp b/src/gamma.lisp
index db5d140..128e888 100755
--- a/src/gamma.lisp
+++ b/src/gamma.lisp
@@ -2032,6 +2032,40 @@
(bigfloat (bigfloat (maxima::bfloat-erf (maxima::to z))))
(complex-bigfloat (bigfloat (maxima::complex-bfloat-erf (maxima::to z))))))))
+(defun bf-erfc (z)
+ ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is
+ ;; near 1. Wolfram says
+ ;;
+ ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2))
+ ;;
+ ;; For real(z) > 0, sqrt(z^2)/z is 1 so
+ ;;
+ ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
+ ;; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)
+ ;;
+ ;; For real(z) < 0, sqrt(z^2)/z is -1 so
+ ;;
+ ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
+ ;; = 2 - 1/sqrt(pi)*gamma_incomplete(1/2,z^2)
+ (flet ((gamma-inc (z)
+ (etypecase z
+ (cl:number
+ (maxima::gamma-incomplete 1/2 z))
+ (bigfloat
+ (bigfloat:to (maxima::$bfloat
+ (maxima::bfloat-gamma-incomplete (maxima::$bfloat maxima::1//2)
+ (maxima::to z)))))
+ (complex-bigfloat
+ (bigfloat:to (maxima::$bfloat
+ (maxima::complex-bfloat-gamma-incomplete (maxima::$bfloat maxima::1//2)
+ (maxima::to z))))))))
+ (if (>= (realpart z) 0)
+ (/ (gamma-inc (* z z))
+ (sqrt (%pi z)))
+ (- 2
+ (/ (gamma-inc (* z z))
+ (sqrt (%pi z)))))))
+
(in-package :maxima)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
@@ -2253,21 +2287,13 @@
((eq z '$inf) 0)
((eq z '$minf) 2)
- ;; Check for numerical evaluation. Use erfc(z) = 1-erf(z).
+ ;; Check for numerical evaluation.
- ((float-numerical-eval-p z)
- (- 1.0 (erf ($float z))))
- ((complex-float-numerical-eval-p z)
- (complexify
- (- 1.0
- (complex-erf
- (complex ($float ($realpart z)) ($float ($imagpart z)))))))
- ((bigfloat-numerical-eval-p z)
- (sub 1.0 (bfloat-erf ($bfloat z))))
- ((complex-bigfloat-numerical-eval-p z)
- (sub 1.0
- (complex-bfloat-erf
- (add ($bfloat ($realpart z)) (mul '$%i ($bfloat ($imagpart z)))))))
+ ((or (float-numerical-eval-p z)
+ (complex-float-numerical-eval-p z)
+ (bigfloat-numerical-eval-p z)
+ (complex-bigfloat-numerical-eval-p z))
+ (to (bigfloat::bf-erfc (bigfloat:to z))))
;; Argument simplification
diff --git a/tests/rtest_gamma.mac b/tests/rtest_gamma.mac
index 9891e98..3ac95cb 100755
--- a/tests/rtest_gamma.mac
+++ b/tests/rtest_gamma.mac
@@ -2579,6 +2579,31 @@ relerror(
1b-65);
true;
+/* Bug 3587304 erfc(x) for x > 6 is wrong */
+relerror(
+ erfc(6.0),
+ 2.15197367124989131659335039918738463047751406e-17,
+ 1e-15);
+true;
+
+relerror(
+ erfc(-4.0),
+ 2-erfc(4.0),
+ 1e-15);
+true;
+
+relerror(
+ erfc(6b0),
+ 2.1519736712498913116593350399187384630477514061688542100527892051056337238484927860b-17,
+ 1b-64);
+true;
+
+relerror(
+ erfc(-6b0),
+ 1.9999999999999999784802632875010868834066496008126153695224859383114578994b0,
+ 1b-64);
+true;
+
/* We have done a check for the numerical algorithm of the Erf function which
calls the Incomplete Gamma function.
We do not do further numerical tests for the other Error functions,
commit 2284ca4017b07e9c5551e86b7f17d9ae65f233af
Author: Raymond Toy <toy...@gm...>
Date: Wed Nov 14 17:23:16 2012 -0800
Correct the threshold in bf-erf.
diff --git a/src/gamma.lisp b/src/gamma.lisp
index ba4a4b8..db5d140 100755
--- a/src/gamma.lisp
+++ b/src/gamma.lisp
@@ -2004,7 +2004,7 @@
(cond ((typep z 'cl:real)
;; Use Slatec derf, which should be faster than the series.
(maxima::erf z))
- ((<= (abs z) 0.47329)
+ ((<= (abs z) 0.476936)
;; Use the series A&S 7.1.5 for small x:
;;
;; erf(z) = 2*z/sqrt(%pi) * sum((-1)^n*z^(2*n)/n!/(2*n+1), n, 0, inf)
commit cd29013f2f1eb364493ac547365317c82c6618d8
Author: Raymond Toy <toy...@gm...>
Date: Wed Nov 14 15:43:51 2012 -0800
Fix Bug 3587191: erf inaccurate for small bigfloat values
src/gamma.lisp:
o For small complex, bigfloat, or complex bigfloat values, update
bf-erf to use a series instead of using gamma-incomplete to compute
the value.
o Add comment that complex-erf, bfloat-erf, and complex-bfloat-erf
have round-off problems for small arguments.
tests/rtest_gamma.mac:
o Add test.
diff --git a/src/gamma.lisp b/src/gamma.lisp
index 00bd194..ba4a4b8 100755
--- a/src/gamma.lisp
+++ b/src/gamma.lisp
@@ -1005,7 +1005,7 @@
(when (eq ($sign (sub (simplify (list '(mabs) del))
(mul (simplify (list '(mabs) sum)) gm-eps)))
'$neg)
- (when *debug-gamma* (format t "~&Series converged.~%"))
+ (when *debug-gamma* (format t "~&Series converged to ~A.~%" sum))
(return
(sub (simplify (list '(%gamma) a))
($rectform
@@ -1885,15 +1885,15 @@
;; Check for numerical evaluation
((float-numerical-eval-p z)
- (erf ($float z)))
+ (bigfloat::bf-erf ($float z)))
((complex-float-numerical-eval-p z)
(complexify
- (complex-erf (complex ($float ($realpart z)) ($float ($imagpart z))))))
+ (bigfloat::bf-erf (complex ($float ($realpart z)) ($float ($imagpart z))))))
((bigfloat-numerical-eval-p z)
- (bfloat-erf ($bfloat z)))
+ (to (bigfloat::bf-erf (bigfloat:to ($bfloat z)))))
((complex-bigfloat-numerical-eval-p z)
- (complex-bfloat-erf
- (add ($bfloat ($realpart z)) (mul '$%i ($bfloat ($imagpart z))))))
+ (to (bigfloat::bf-erf
+ (bigfloat:to (add ($bfloat ($realpart z)) (mul '$%i ($bfloat ($imagpart z))))))))
;; Argument simplification
@@ -1944,6 +1944,7 @@
(defun complex-erf (z)
(let ((result
+ ;; Warning! This has round-off problems when abs(z) is very small.
(*
(/ (sqrt (expt z 2)) z)
(- 1.0
@@ -1959,6 +1960,7 @@
result))))
(defun bfloat-erf (z)
+ ;; Warning! This has round-off problems when abs(z) is very small.
(let ((1//2 ($bfloat '((rat simp) 1 2))))
;; The argument is real, the result is real too
($realpart
@@ -1970,6 +1972,7 @@
(bfloat-gamma-incomplete 1//2 ($bfloat (power z 2)))))))))
(defun complex-bfloat-erf (z)
+ ;; Warning! This has round-off problems when abs(z) is very small.
(let* (($ratprint nil)
(1//2 ($bfloat '((rat simp) 1 2)))
(result
@@ -1992,6 +1995,45 @@
;; A general complex result
result))))
+(in-package :bigfloat)
+
+;; Erf(z) for all z. Z must be a CL real or complex number or a
+;; BIGFLOAT or COMPLEX-BIGFLOAT object. The result will be of the
+;; same type as Z.
+(defun bf-erf (z)
+ (cond ((typep z 'cl:real)
+ ;; Use Slatec derf, which should be faster than the series.
+ (maxima::erf z))
+ ((<= (abs z) 0.47329)
+ ;; Use the series A&S 7.1.5 for small x:
+ ;;
+ ;; erf(z) = 2*z/sqrt(%pi) * sum((-1)^n*z^(2*n)/n!/(2*n+1), n, 0, inf)
+ ;;
+ ;; The threshold is approximately erf(x) = 0.5. (Doesn't
+ ;; have to be super accurate.) This gives max accuracy when
+ ;; using the identity erf(x) = 1 - erfc(x).
+ (let ((z^2 (* z z))
+ (eps (epsilon z)))
+ (do* ((n 0 (+ n 1))
+ (sum 1 (+ sum term))
+ (factor 1/3
+ (/ (+ n n 1) (+ n 1) (+ n n 3)))
+ (term (* (- z^2) factor)
+ (* term (* (- z^2) factor))))
+ ((< (abs term) (* eps (abs sum)))
+ (* 2 z sum (/ (sqrt (%pi z)))))
+ #+nil
+ (format t "n = ~S: sum = ~S term = ~S factor = ~S~%" n sum term factor))))
+ (t
+ ;; The general case.
+ (etypecase z
+ (cl:real (maxima::erf z))
+ (cl:complex (maxima::complex-erf z))
+ (bigfloat (bigfloat (maxima::bfloat-erf (maxima::to z))))
+ (complex-bigfloat (bigfloat (maxima::complex-bfloat-erf (maxima::to z))))))))
+
+(in-package :maxima)
+
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |