This implementation of mod2pi in Julia, for example, uses double doubles ("TwicePrecision") here and there:
https://github.com/pkofod/RemPiO2.jl/blob/master/src/RemPiO2...
And, pertinent to the article, here's a bug in mod2pi that was probably triggered by 80 bit extended floats:
===Update===
Couldn't find the one I was looking for, but I did dig up Eckert, Punched Card Methods in Scientific Computation (1940).
From Chapter VII "The Multiplication of Series":
p.68 "By a continuation of this process we obtain about 100 groups of cards each with its rate card. Since the capacity of the multiplier is eight digits, those groups which have 9 and 10 digit multipliers and eight digit multiplicand are done in two groups. The few terms with large multiplier and multiplicand are done by hand. In a series with reasonable convergence there may be about 20,000 cards altogether. The reproducer is of course used to check the punching."
10 decimal digits is what, 33 bits?
p.74 "The machine time for 100 series is as follows:
Original punching and verifying . . . . 2 days
Listing and summary punching checksums. 1/4 day
Sorting and gang punching duplicates. . 1 day
Multiplying . . . . . . . . . . . . . . 1 day
Sorting, tabulating, and summary punch. 1 day
In this example the checks would have been exact had we included a card for each product regardless of size, but this would have required about forty thousand cards instead of four or five thousand, and the multiplying time would have been over a week instead of one day."I guess in principle a punched card could have contained 80 decimal places, but it seemed like the ALU equivalents were much narrower.
(this was published by the Thomas J. Watson Astronomical Computing Bureau)
Interestingly, you don't need multiprecision floating point for every pixel: you can do precise calculations for a few tentpole pixels and calculate most of the rest at adouble precision based on the differences between their location and a tentpole. But you also need some significant cleverness to accurately fill in the gap betwen "most" and "all":
http://www.fractalforums.com/announcements-and-news/pertubat...
Algorithms stability. As precision decreases when magnitude rises, many tricks are used in sci. comp. to scale everything into ranges where computations won't diverge due to a lack of precision.
https://www.jpl.nasa.gov/edu/news/2016/3/16/how-many-decimal...
* JPL's highest accuracy calculations, which are for interplanetary navigation, we use 3.141592653589793.
* How many digits of pi would we need to calculate the circumference of a circle with a radius of 46 billion light years to an accuracy equal to the diameter of a hydrogen atom (the simplest atom)? The answer is that you would need 39 or 40 decimal places.
By 1596 Ludolph van Ceulen had 20 digits, and even more by the time he designed his tombstone: https://upload.wikimedia.org/wikipedia/commons/e/e8/Monument...
39 or 40 decimal places would have to wait until 1699.
https://en.wikipedia.org/wiki/Chronology_of_computation_of_π
If you are using a 64-bit float, then you have 53 bits of precision in 1D, 53 bits of precision in 2D, and 53 bits of precision in 3D. The precision doesn’t change as you add or remove dimensions.
[1] http://web.mit.edu/tabbott/Public/quaddouble-debian/qd-2.3.4...
(i mean this is a legitimate reason, it's just amusing the extremes you'd have to go to to break 64bits).
Microseconds gets you to the year 2255; milliseconds to heat death.
1. https://www.wolframalpha.com/input/?i=January+1%2C+1970+%2B+...
It is easy to run out of bits when trying to use a single number as both an absolute date/time, and to compute relative durations between timestamps with reasonable precision. I've lost track of the number of bugs I've fixed where someone assumed a double would be big enough without doing the math.
The NTP nanokernel uses a 64.64 integer represention https://papers.freebsd.org/2000/phk-nanokernel/ because 32 bits isn't enough above the binary point because it runs out in 2038, and below the binary point you need 64 bits to accurately represent the nanosecond-scale frequency of the CPU's cycle counter to translate it to wall-clock time.
Addition has problems where x is close to negative y, sure.
So of course you want to avoid error states but a double already has an exponent range well beyond beyond astronomical. The width of the universe is around 10^27, a planck length is around 10^-35, and the limit of the format is ±10^308.
If 64 bit IEEE-754 floating point numbers do not support your needs, adding more bits will almost certainly not fix the problem. The problem isn't lack of bits, it's going to be something fundamental about the nature of IEEE-754. You'll probably need something like a rational data type, or a symbolic algebra system that can abstractly represent everything that shows up in your domain without losing any precision. Generally you'll be required to do this in software.
32 bit is the sweet spot between precision and performance. There are only a handful of things that require 64 bits, (notably latitude and longitude) but these things are common enough that hardware support is (IMO) valuable.
Basically, they're summing up a bunch of intermediate calculations. The intermediate calculations are 64 bit floats, and are precise to 17 decimal digits or whatever. They sum all the intermediate components, and get a result that's only precise to 13 digits or whatever. They use double doubles to perform the summation, and they get a result that's back to being precise to 17 decimal digits. Great.
Now imagine doing this with 128 bit floats, which are precise to 33 digits. So you sum the intermediate results, now you're precise to 29 digits. So you've lost 4 digits of precision again. So you add a 192 bit floating point type...
IEEE-754 floating points are always going to have that problem. If you add two numbers of differing magnitude, you're going to lose precision. Finding the sum of a large sequence is generally going to result in the addition of a very large running total with small individual members.
(I perused the other examples, but they seem to be a variation on the same theme)
> This permitted benchmark results to be accurately reproduced for a significantly longer time, with virtually no change in total run time
They clearly understand that bigger floats do not magically change the fundamental behavior, but the additional headroom makes a significant difference
> A larger example of this sort arose in an atmospheric model (a component of large climate model). While such computations are by their fundamental nature “chaotic,” so that computations will eventually depart from any benchmark standard case, nonetheless it is essential to distinguish avoidable numerical error from fundamental chaos.
> Researchers working with this atmospheric model were perplexed by the difficulty of reproducing benchmark results. Even when their code was ported from one system to another, or when the number of processors used was changed, the computed data diverged from a benchmark run after just a few days of simulated time. As a result, they could never be sure that in the process of porting their code or changing the number of processors that they did not introduce a bug into their code.
> After an in-depth analysis of this code, He and Ding found that merely by employing double-double arithmetic in two critical global summations, almost all of this numerical variability was eliminated. This permitted benchmark results to be accurately reproduced for a significantly longer time, with virtually no change in total run time [3].
The keyword is "numerical variability." Here's the referenced article: https://link.springer.com/article/10.1023/A:1008153532043 Choice quote from the article:
> In climate model simulations, for example, the initial conditions and boundary forcings can seldomly be measured more accurately than a few percent. Thus in most situations, we only require 2 decimal digits accuracy in final results. But this does not imply that 2 decimal digits accuracy arithmetic (or 6-7 bits mantissa plus exponents) can be employed during the internal intermediate calculations. In fact, double precision arithmetic is usually required.
The problem isn't lack of precision. The problem is numerical instability when adding up a bunch of numbers with high absolute values but since they were roughly evenly positive/negative, their sum was approximately 1. IEEE-754 floats, as useful as they are, are just bad at this, and adding more bits isn't a solution, it's a punt. They used Kahan summation or Bailey summation and the problem went away. No 128 bit hardware floats required. (Kahan summation is very well known, Bailey summation is new to me)
Here's my point: if double precision floating point doesn't satisfy your needs, you should dig into the problem and understand why. Understand first, write code second. 999/1000 the solution isn't "we need 128 bit floats", and for that .1%, we're waaaay better off telling those people "Sorry, do it in software and take the performance hit."