When Double Precision Is Not Enough

An interesting tool to help preserve floating point precision is Herbie (<a href="https:&#x2F;&#x2F;herbie.uwplse.org&#x2F;" rel="nofollow">https:&#x2F;&#x2F;herbie.uwplse.org&#x2F;</a>), which takes an expression as input and outputs an optimized version of it that will retain precision.<p>The example they give on their site is:<p><pre><code> sqrt(x+1) - sqrt(x) =&gt; 1&#x2F;(sqrt(x+1) + sqrt(x))</code></pre>

There are a few edge cases in geospatial computational geometry that require greater than quad precision to correctly compute a double precision result. This fact is sometimes used as a litmus test for implementation correctness. Many implementations don&#x27;t have a conditional escape for these edge cases that switches to arbitrary precision algorithms in their double precision algorithms, which is easy to verify.

I remember a class I had where we were shown a methods for exact PC arithmetic with integers, rationals (e.g., 1&#x2F;2) and algebraic numbers (e.g., sqrt(2)). Only thing it couldn&#x27;t deal with were transcendental numbers (e.g., pi)<p>I think it worked by representing the numbers by integer matrices, all operations were then matrix operations. Unfortunately, the matrices were gaining dimension when new algebraic numbers got involved, so it probably wasn&#x27;t useful for anything :-).<p>Anyway, it it blew me away anyway at the time, one of the few things that sticked with me from school.

If anyone ever runs out of precision using IIR filters: convert to second-order section (SOS). I stumbled upon this trick a few years ago and it&#x27;s a game changer. No need for FIR.

Came here expecting an article about Kahan summation, left having learned about a multiple-precision Python library...

IIUC most people instead precondition their matrices and stick to double precision instead of quad or deep precision.

This reminded me of Posits[1] and Unums[2]. I dabbled with my own implementation some years ago, but nothing serious.<p>Anyone tried them out? I see there&#x27;s some FPGA implementations and a RISC-V extension idea[3], would be fun to try it on a softcore.<p>[1]: <a href="https:&#x2F;&#x2F;posithub.org&#x2F;docs&#x2F;posit_standard.pdf" rel="nofollow">https:&#x2F;&#x2F;posithub.org&#x2F;docs&#x2F;posit_standard.pdf</a><p>[2]: <a href="http:&#x2F;&#x2F;johngustafson.net&#x2F;unums.html" rel="nofollow">http:&#x2F;&#x2F;johngustafson.net&#x2F;unums.html</a><p>[3]: <a href="https:&#x2F;&#x2F;www.posithub.org&#x2F;khub_doc" rel="nofollow">https:&#x2F;&#x2F;www.posithub.org&#x2F;khub_doc</a>

&gt;&gt; Numerically stable algorithms exist for such problems, but these occasionally fail leaving us to brute force the calculation with higher precision to minimize floating point rounding errors.<p>Do algorithm authors really implement two code paths?<p>One that uses a cost function iterative algorithm to solve the matrix (this is what everybody does), and a second algorithm using the brute force approach the author showed (known not to work).

First thought: the curve in that first plot reflects the pseudospectrum? (<a href="https:&#x2F;&#x2F;en.wikipedia.org&#x2F;wiki&#x2F;Pseudospectrum" rel="nofollow">https:&#x2F;&#x2F;en.wikipedia.org&#x2F;wiki&#x2F;Pseudospectrum</a>) Spectra of nonnormal matrices can be very hard to calculate accurately.

The are algorithms taught in numerical analysis class for re-ordering calculations to maximize precision. Computing in order you wrote on paper or the way you were initially taught may not be the best.

I&#x27;ll admit being lost in some of the math here but I wonder perhaps naively, if a language like Julia, which does fractions as fractions, would handle it better.

Thankfully someone figured out that quad should be in C++ a while back, unfortunately they didnt figure out that the next few dozen should be too... the arbitrary precision stuff is an alternative in theory, and would be if any of the hundreds of half completed projects to convert unintelligeble, and unmaintainable, unclear, etc, C to the in this specific case extremely much better C++ interfaces was complete.

There is a whole field about that, called numerical analysis: <a href="https:&#x2F;&#x2F;en.m.wikipedia.org&#x2F;wiki&#x2F;Numerical_analysis" rel="nofollow">https:&#x2F;&#x2F;en.m.wikipedia.org&#x2F;wiki&#x2F;Numerical_analysis</a>

Sounds like bad algorithm.

Obligatory reference to computable analysis: <a href="https:&#x2F;&#x2F;en.wikipedia.org&#x2F;wiki&#x2F;Computable_analysis" rel="nofollow">https:&#x2F;&#x2F;en.wikipedia.org&#x2F;wiki&#x2F;Computable_analysis</a><p>This shows that you <i>can</i> do exact real arithmetic without worrying about round-off errors. But it&#x27;s expensive.<p>If you want something more practical (but dodgy), you can use arbitrary-precision arithmetic instead of exact real arithmetic. To that end, there&#x27;s the Decimal module in the Python standard lib. There&#x27;s also the BC programming language which is part of POSIX. If you want to do matrix arithmetic, you can use mpmath for Python like they did in this blog post, or Eigen3 for C++ with an arbitrary-precision type. This is dodgy of course because the arithmetic remains inexact, albeit with smaller rounding errors than double-precision floats.

also worth noting that you can ask the fpu to trap to catch some types of precision loss or fp underflows.

Rational arithmetic needs better matrix libraries. No loss of precision then.

use 80 bits? <a href="https:&#x2F;&#x2F;en.wikipedia.org&#x2F;wiki&#x2F;Standard_Apple_Numerics_Environment" rel="nofollow">https:&#x2F;&#x2F;en.wikipedia.org&#x2F;wiki&#x2F;Standard_Apple_Numerics_Enviro...</a>