#floatingpoint

After having some trouble with importing double precision strings into Mathematica, as described in previous few posts, I decided to write something myself: StringToDouble.wl

It is simple to use, just import the package and StringToDouble[str] will return the correctly rounded double precision value.

The implementation takes great advantage of the arbitrary precision integer arithmetic Mathematica provides. Arbitrary precision generally can't be avoided in this kind of tasks as input string can be of any length, containing any number of decimal digits, so it is the only way to get the correctly rounded result.

Algorithm is based on AlgorithmM from this article. This reference and many other details from Rick Regan's blog were of great help during the process.

As string to double conversion is historically proved to be a tricky problem, all compilers or systems had or still have considerable problems with it (e.g. MinGW's atof is heavily broken), rigorous testing was needed. For this purpose I used a few of my tests, but also collected many cases from the most respected strtod implementation development (David Gay's strtod available in atof.c file), Rick Regan's blog, GCC test suite, and other sources. This Mathematica implementation passes them all. Complete test collection is available in TOML format.

It's not only surprising tricky to import valid C-style double precision literal as Mathematica double precision value, but it is still not known how to do it correctly. One could think there are quite a few solutions available using various Mathematica API. The problem resulted in many forum posts (StackOverflow, StackExchange...) with many answers available, all of them containing serious drawbacks. In this post I will present a few different approaches, all of which fail for one reason or another.

Methods tested:

They were checked against the following tests:

test IDinput stringC-style %.17g formatexact value11e-4000022.4703282292062326e-3240032.4703282292062327e-3240042.4703282292062328e-3244.9406564584124654e-3241*2^-107454.9406564584124654e-3244.9406564584124654e-3241*2^-107461e-3209.9998886718268301e-3212024*2^-107472.2250738585072009e-3082.2250738585072009e-3084503599627370495*2^-107482.2250738585072014e-3082.2250738585072014e-3084503599627370496*2^-107491e-251e-258711228593176025*2^-136101e-239.9999999999999996e-246805647338418769*2^-129111.23456789012345678901.23456789012345675559999489923579*2^-52121.13e111.3000000000000016361334473660826*2^-49131e+191e+194882812500000000*2^11141.9e+221.9000000000000002e+224529953002929688*2^22152.7e222.7000000000000002e+226437301635742188*2^22163.1e223.1000000000000002e+227390975952148438*2^22177e227.0000000000000004e+228344650268554688*2^23181.7976931348623158e+3081.7976931348623157e+3089007199254740991*2^971191.7976931348623159e+308inf$\infty$201e+400inf$\infty$

Example of checking method A against test 14:

Results:

method IDfailed tests ID'scommentA1, 2, 3, 4, 5, 6, 9, 10, 12, 14, 19, 20B1, 2, 3, 4, 6, 9, 10, 19, 20C1, 2, 3, 4, 5, 6, 9, 10, 13, 19, 20returns 0 for test 13 due to a bug in Read[#, Number]D1, 2, 3, 4, 5, 6, 7, 11, 19, 20E2, 3, 19, 20tests 2-3 fail due to a bug in 64-bit .NETF4, 14, 15, 16, 17, 19, 20G1, 2, 3, 4, 6, 9, 10, 12, 14, 19, 20H1, 2, 3, 4, 6, 19, 20I1, 2, 3, 4, 6, 19, 20J1, 2, 3, 4, 6, 19, 20K4, 5, 7

Most methods are close, but often last bit is not correct, which nevertheless makes the results unusable. None of the native methods handle underflow, overflow and subnormal (denormalized) cases at all. Method D detects these cases but raises various errors. In other cases literals of denormalized values result tends to contain considerable error as they are deserialized using complete 53-bit precision, no matter IEEE-754 double precision format limits. There is a hack to overcome this problems using ImportString/ExportString with "Real64" type argument exploited in method K.

The biggest disappointments are methods D, E, F. They should be able to import any valid double literal. In the case of F, JSON number literals really can't end with the decimal point (RFC 4627, Section 2.4), but it also contains a bug preventing it to correctly deserialize numbers in some rather normal cases (tests 4, 14-17).

Motivated by the discovery of high error level of Wolfram Mathematica's implementation of gamma function for positive real arguments, as described in the previous post, I investigated the quality of following three notable implementations: Matlab R2015b, Maple 2015 and Netlib Cephes Math Library.

Matlab and Maple are well known commercial products. Netlib is collection of free numerical library sources in C and Fortran available without restrictions, widely used in projects such as CERN's ROOT framework. Gamma function is available in Cephes cprob.tgz archive.

Again, results are a bit surprising. Matlab exhibits high level of inaccuracy with error signature similar to that of Mathematica. Looks like it uses good approximation for $x \lt 12$, but then most probably switches to loggamma technique which soon generates errors in a range of few hundred ULPs.

Maple is an interesting case. I have never seen such error signature. Also, error is very low for $0.75 \lt x \lt 2.25$ with value mostly below 0.5 ULPs (theoretical limit) and always below 0.6 ULPs. This is very impressive! It escalates from that point forward in an interesting steps pattern and soon reaches hundreds ULPs level. It is generaly better then Matlab and Mathematica.

Netlib results are excellent with error always below 6 ULPs. It is clearly visible that Netlib error is equivalent to Scipy error suggesting it uses this implementation as well. I was not aware of this fact before starting to write this post. According to the documentation, recurrence reduction is used for $x \le 34$ and Stirling's formula for $x \gt 34$. Judging from the error plot, it would probably be better to put this domain cut at $x = 15$, it could potentially reduce the maximal error level.

Floating point calculation is rarely correct to all significant bits. Operating in some fixed precision format, most often there are no means to take care of all bits needed to get perfectly correct result. Some of just a few exceptions are multiplication by $0.5$ and division by $2$, which can be executed without loosing a single significant bit for most representable values. But for instance this is not the case for division by $3$, and in a lot of cases you will loose last significant bit so after multiplying with $3$ you will not get the original value back.

Combining more numerical operations into some complex method it is possible to loose a lot of precision (bits) in the process. This is the main distinction between great, good, bad or catastrophic implementations. It is said only 5% of developers are even aware of this aspect of programming, and only a few pursue enhancements in this aspect.

To be able to analyse the quality of any floating point result, it is essential to obtain the more precise one, even infinitely precise if possible. This is perfect task for Wolfram Mathematica. Default floating point format in WM is IEEE-754 compliant double precision format, obtained with N[x] without additional arguments, i.e. do not use something like N[x,5] for this purpose. Other way around, exact (infinitely precise) value of some double precision number $x$ is obtained with SetPrecision[x,\[Infinity]]. Now lets check the error of WM Sqrt[x] function for $x=11.13$ double precision value.

The result is $-0.47$ ulps. On most architectures sqrt(x) method is implemented to maximal accuracy, so error should never be larger than $0.5$ ulps in absolute value, meaning all bits are correct. It is very important to note that exact value of double precision value $11.13$ is not $1113/100$, as it is not representable by the format, but it equals to: $$xExact=11.1300000000000007815970093361102044582366943359375,$$ or more compactly: $$xExact=6265632981579203/562949953421312.$$ This value, calculated at line 8, has to be used to calculate the exact result, not $1113/100$! There are some more details in the workflow so it should be followed precisely.

The following WM code computes and plots ULP error of WM's Sqrt[x] for $i/100,\, i\in[1,10000]$ double precision values. This approach is suitable for every function which can be evaluated to high precision in WM.

Using the methods and criteria established in the previous posts, I will compare the accuracy of Wolfram Mathematica 10, Python 3 and SciPy 0.16 implementations of gamma function for double precision arguments.

As presented in the plots below, both Python math.gamma(x) and SciPy scipy.special.gamma(x) implementations are very accurate with the average error below $1\, ulp$ and maximal error below $6\, ulps$ on most of the domain.

On the other hand, Mathematica implementation excels for $x \lt 10$, but something goes very wrong starting from $x \approx 34$, and gets totally out of control for $x \gt 50$. Error becomes larger than $500\, ulps$, reaching $1800\, ulps$ at some points. It also gets seriously disbalanced with error being spread asymmetrically about x axis.

After a short review what could go wrong, I recalled some tend to use loggamma(x) function to evaluate gamma(x) for large arguments. After checking Mathematica LogGamma[x] 'error signature' the proof was in that this really is the case. As illustrated in the last plot LogGamma[x] has the same domain cut-points and disbalances. This most certainly means Mathematica exponentiates LogGamma[x] to evaluate Gamma[x] for $x\gt 50$. This approach is classical numerical pitfall which generates inaccuracy in last 10 mantissa bits and catapults the error into $1000\, ulps$ range. Log-Exp interplay was the very reason why extended precision formats were introduced in the first place! Wolfram should seriously rework the implementation by either removing the dependency on LogGamma[x] or by computing it in extended precision format.

ULP (unit in the last place) is the most convenient method for measuring numerical error. It is common to see absolute or relative errors being used, but these approaches does not come naturally in most cases. Using ULP one can intuitively see how bad the result is, could it be a little or much better, is there a need for some additional improvements... The most elementary numerical operations (e.g. addition, division) should be calculated within 0.5 ulps, simple transcendental methods (e.g. sin, atan2) usually guarantee an error below 5 ulps, while the most complicated functions can have an error of about 20 ulps and more.

ULP was invented by William M. Kahan in 1960 who originaly defined it as: ulp (x) is the gap between the two floating-point numbers nearest x, even if x is one of them. Although this definition is simple and has many advantages over others, more modern developments have been made since [Muller, 2005].

The following image visualizes many basic concepts regarding this topic. Vertical lines represent consecutive representable numbers. Their distribution is not uniform. As a consequence of the format defined by IEEE-754 standard, density changes at every perfect power of 2. In this example 1 is used for reference, although any perfect power of 2 could have been used. Connection between ulp(1) and machine $\epsilon$ is clearly visible. Value of $z$ is not representable within the format, but ulp(z) is simply the distance between two representable numbers enclosing $z$.

Error of some result $x$ is estimated relative to $correct$ value scaled to $ulp$ 'units': $$error(x) = \frac{x - correct}{ulp(correct)}.$$ Always use the correct value to calculate ULP!

One example, for two double precision numbers $x=2.5445347236749\color{red}{499}$ and $correct=2.5445347236749\color[rgb]{0.2,0.7,0.2}{384}$, error is: $$error(x) = \frac{x - correct}{ulp(correct)} = \frac{1.154\mbox{e-}14}{4.44\mbox{e-}16} = 26 \, ulps$$

In this post I will focus on IEEE-754 double precision numbers. Following battery of tests should be used to check any implementation of double precision ULP calculation. Although it is possible to construct many others, these are carefully selected to detect the most common pitfalls.

$$\begin{align*}\mathcal{min} &= 2.2250738585072014\mbox{e-}308 &\text{(minimal normal double)}\cr \mathcal{max} &= 1.7976931348623157\mbox{e+}308 &\text{(maximal double)}\cr \epsilon &= 2.220446049250313\mbox{e-}16 &\text{(double precision epsilon)} \end{align*}$$

$x$$ulp(x)$1$0.0$$2^{-1074} = 5\mbox{e-}324$2$\epsilon \, \mathcal{min} = 4.9406564584124654\mbox{e-}324$$2^{-1074} = 5\mbox{e-}324$3$\mathcal{min} = 2.2250738585072014\mbox{e-}308$$2^{-1074} = 5\mbox{e-}324$4$2 \, \, \mathcal{min} * (1 - \epsilon) = 4.4501477170144018\mbox{e-}308$$2^{-1074} = 5\mbox{e-}324$5$2^{i}, \quad i\in[-1021, 1023]$$2^{i - 53}$6$2^{i} * (1 + \epsilon), \quad i\in[-1021, 1023]$$2^{i - 52}$7$\mathcal{max} = 1.7976931348623157\mbox{e+}308$$2^{971} = 1.99584030953472\mbox{e+}292$

It is not an easy task to find implementations of ulp(x) containing no errors and passing all the tests. This topic can be considered simple, but there are lot of misconceptions and bugs out there. The most common are messing up the sign, using next representable number to get ulp(x) for positive x and not handling denormalized numbers correctly. For positive x, previous representable number should be used, as clearly demonstrated in the illustration. Otherwise you could end up with underestimated error or overflow for ulp(max). It is advisible to use nextafter in C or struct.pack in Python to obtain the most robust result. If not available, ulp(x) can be calculated via explicit calculation, an approach most convenient for porting to other languages and architectures, but make sure implementation passes all the tests!

Far too often, software implementations of various numerical algorithms use the following relation for triangle law of cosines:

\begin{equation} c^2 = a^2 + b^2 - 2\,a\,b\,cos(\gamma). \end{equation} Although mathematically correct, this relation is numerically unstable. In case of $a\approx b$ and $\gamma\approx 0$, triangle side $c$ will be calculated with a considerable numerical error. The reason is that in this case the difference of two large quantities, $a^2 + b^2$ and $2\,a\,b\,cos(\gamma)$, is very small, so most of the significant bits are canceled out.

This situation can be avoided using numerically stable and always correct formula. It seems this can't be googled out, so I hope this post will make a difference. According to William M. Kahan's Miscalculating Area and Angles of a Needle-like Triangle one should use:

\begin{equation} c^2 = (a - b)^2 + 4\,a\,b\,sin(\gamma / 2)^2. \end{equation}

In this approach difference $a-b$ is calculated early, and for the critical case the result is completely contained within the second term, which is calculated precisely.

Comparison of both relations in case $a=b=1$ computed in double precision and displayed to 15 significant digits is shown in the following table.