This is the mail archive of the
newlib@sourceware.org
mailing list for the newlib project.
erf(x)+erfc(x)=1 precision problem
- From: "Reini Urban" <rurban at x-ray dot at>
- To: newlib at sourceware dot org
- Date: Sun, 1 Jun 2008 11:19:37 +0200
- Subject: erf(x)+erfc(x)=1 precision problem
While testing the new clisp-cvs on cygwin I came to the following
erf/erfc precision problems with the double definitions.
There are some minor deviations from the glibc2 (linux) results, but a
bigger problem is that erf(-1) + erfc(-1) is unequal to (double) 1.0.
The error is 1.1102230246251565e-16
This looks like just one digit, but it is the only significant digit
deviation in the erf(x) + erfc(x)
test between [-10, 10].
erfc(27) leads to floating point underflow too, but this is probably acceptable.
Sorry, I'm not qualified at all to offer a patch, just the link to sun's
free erf/erfc implementation [1], which lcc also uses.
Results in detail:
CORRECT is glibc, CLISP is cygwin newlib gcc-3.4.4 on a i686 with HW FPU.
Form: (ERFC -1)
CORRECT: 1.842700792949715d0
CLISP : 1.8427007929497148d0
Form: (ERFC 3)
CORRECT: 2.2090496998585438d-5
CLISP : 2.209049699858544d-5
Form: (ERFC 5)
CORRECT: 1.5374597944280351d-12
CLISP : 1.537459794428035d-12
Form: (ERFC 20)
CORRECT: 5.3958656116079005d-176
CLISP : 5.395865611607901d-176
Form: (ERFC 22)
CORRECT: 1.6219058609334724d-212
CLISP : 1.6219058609334726d-212
Form: (ERFC 23)
CORRECT: 4.441265948088058d-232
CLISP : 4.441265948088057d-232
Form: (ERFC 25)
CORRECT: 8.300172571196522d-274
CLISP : 8.300172571196523d-274
Form: (LOOP :FOR I :FROM -10 :TO 10 :ALWAYS (= 1 (+ (ERF I) (ERFC I))))
CORRECT: T
CLISP : NIL
Form: (Y0 1.0)
CORRECT: 0.08825696421567698d0
CLISP : 0.08825696421567696d0
Form: (Y0 10.0)
CORRECT: 0.055671167283599395d0
CLISP : 0.05567116728359938d0
Form: (Y1 10.0)
CORRECT: 0.2490154242069538d0
CLISP : 0.24901542420695386d0
Form: (YN 2 10.0)
CORRECT: -0.0058680824422086345d0
CLISP : -0.005868082442208612d0
[6]> (loop :for i :from -10 :to 10 :always (print (abs (- (+ (os:erf
i) (os:erfc i)) 1))))
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
1.1102230246251565d-16
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
0.0d0
That means that for double precision between -10 and +10 only -1 is imprecise.
glibc is correct here.
(loop :for i :from -100 :to 100 :always (print (list i (abs (- (+
(os:erf i) (os:erfc i)) 1)))))
says that also only -1 has this problem, -100 to +26 is acceptable,
and that erfc(27) leads to floating point underflow.
[1] http://www.velocityreviews.com/forums/t317358-erf-function-in-c.html
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
--
Reini Urban
http://phpwiki.org/ http://murbreak.at/