467,154 Members | 919 Online

# When (32-bit) double precision isn't precise enough

• viewed: 5693
Share:
5 Replies
 DAVID SCHULMAN wrote: .... ================================================== ======================== [code] ================================================== ======================== And here's its output: P(8.5 s) = 0.000000000000000000E+000 1.#INF00000000000000E+000 .... Can a calculation like this be performed with some kind of arbitrary- precision numeric package (something like Michael Ring's MAPM or PHP BCMath), or evaluated directly by Mathematica or the java.math package? Mathematica, MathCAD, MathLab, Maple all can do such stuff. Here is what Maple spits out: Digits := 190: 1 - erf(20.); 1 / %; -175 .539586561160790 10 .1853270766878889323178572243436799325655920174831 38279103445269\ 28082079337097653296424608926588679681171314780140 46442587\ 58970925582679708132243363506261583082363313448399 25105294\ 176 06688155039 10 Jirka Nov 13 '05 #2
 "DAVID SCHULMAN" writes: [Long explanation snipped. The poster wants to calculate the exclusion probability corresponding to 20 standard deviations of a normal distribution.] In C99, you can use the `erfc' function. Many C89 implementations provide it as an extension. #include #include int main (void) { printf ("%e\n", erfc (20.0 / sqrt (2.0))); return 0; } This program prints `5.507248e-89'. Note that although the /mathematical/ definition of `erfc(x)' is `1.0 - erf(x)', this is not a good way to actually calculate the function. The `erf' function would have to return a value which is 5.507248e-89 less than 1.0, which would require at least 294 significant bits to represent. Usually, floating point numbers are much less precise in C, so `erf' just returns the closest value it can represent, 1.0. So I decided to try to tackle this in C. I downloaded a package called the GNU Scientific Library (http://www.gnu.org/software/gsl/) - there's a nice precompiled binary for Microsoft Visual C++ at http://www.network-theory.co.uk/gsl/freedownloads.html. I wrote a little piece of code to try it out: ================================================== ======================== #include #include #include int main(void) { double odds; double index; double sigma; double result; for (index = 1.0; index <= 8.5; index += 0.5) { sigma = index / M_SQRT2; result = 1 - (gsl_sf_erf (sigma)); odds = 1 / result; printf("P(%2.1f \345) = %.18LE \t %-1.18lE\n", index, result, odds); } return 0; } ================================================== ======================== While I'm not familiar with the GNU Scientific Library, I would guess that it also provides an `erfc' function if it provides an `erf' function. My Visual C++ compiler (v. 6.0) specifies that both double and long double use an 8-byte representation: "The long double contains 80 bits: 1 for sign, 15 for exponent, and 64 for mantissa. Its range is +/- 1.2E4932 with at least 19 digits of precision." As explained above `erf(20.0/sqrt(2.0))' can only be represented as a value different from 1.0 with at least 294 bits for the mantissa. Martin Nov 13 '05 #3