EFT overflows

[munrocket]

@JeffreySarnoff in one rust package https://docs.rs/safeeft/0.0.5/safeeft/ I found new safe EFT.
They said that M.Kashiwagi in his paper (in japanese) found special cases where twoSum and twoProducts not safe enough with IEEE 747-2008 floats. What do you think about that?

twosum ( 3 . 5 6 3 0 6 2 4 4 4 4 8 7 4 5 3 9 e +307 ,
−1 .7976931348623157 e +308)
( −1 .4413868904135704 e +308 , nan )
twoproduct ( 6 . 9 2 9 0 0 1 7 1 3 8 6 9 9 3 6 e +236 ,
2 . 5 9 4 4 4 7 5 2 5 1 9 5 2 0 0 3 e+71)
( 1 . 7 9 7 6 9 3 1 3 4 8 6 2 3 1 5 7 e +308 , i n f )

[dpsanders]

This is very interesting and exactly what we need for IntervalArithmetic.jl.

What is the license on that package? Would you like to port it to Julia?

[dpsanders]

This seems to be the repo, with MIT license!

https://github.com/GNQG/safeeft

[munrocket]

Actually good old EFT just works in most of the cases, especially for usage in graphics. It is possible that nobody found this problem since 1970? Just curious.

[JeffreySarnoff]

I am aware of these outlier cases. I will take a look at the repo code and opine.

[JeffreySarnoff]

AccurateArithmetic.jl is no longer current (though I try). The best code for eft stuff is ErrorfreeArithmetic.jl in my github space.

[munrocket]

Sorry, probably opened issue in wrong place.

[JeffreySarnoff]

safetwosum_branch from the rust repo is < 20% slower than two_sum, the straight line version is much more time consuming.

• [ ]

[JeffreySarnoff]

safetwoproduct_fma from the rust repo is the same as two_prod.
So I could easily add a safe_two_sum that uses the branch logic.
@dpsanders does this work for you?

[JeffreySarnoff]

oic (in the paper) there are also versions of add, sqrt, mul, div that round up and versions round down [but they are not entirely correct :<]

[JeffreySarnoff]

read too quickly (my Japanese is rusty at best [google translate is difficult to read because there are overlapping sections in the translation]. https://translate.google.com/translate?hl=en&sl=ja&u=http://verifiedby.me/adiary/pub/kashi/image/201406/nas2014.pdf&prev=search

The rounded-up and rounded-down routines are intended to deliver Float64 results using two_xxx routines internally. They do work correctly for that @dpsanders they use pretty much the same approach we have used in FastRounding. I could play with recoding that for you. after New Years

[dpsanders]

That sounds great. I don't read Japanese but may be able to get some help with that.

[munrocket]

In Julia you also have emulated fma instruction with EFT too.

[JeffreySarnoff]

yes .. however hardware fma for Float64 and Float32 is widespread, so much so that all of DoubleFloats is built upon the presumption of its presence [and there have been no issues with that decision]. The additional time required absent hardware fma is so great that to require that implementation is to remove one of the essential qualities that others find worthwhile in the pkg.

[munrocket]

So... This is just an overflow problem. If some value overflows, we can say that result become unexpected, that’s all.

[munrocket]

Also this Kashiwagi paper is cool, but without formal proof.

Here is an excerpt from normalization algorithm proof

A missing part in the formal verification is the handling of overflow, as Flocq does not have an upper bound on the exponent range. It is not a real flaw, as an overflow would produce an infinity or a NaN at the end of the algorithm as only additions are involved.
This work was unexpectedly complicated: the formal verification was tedious due to many gaps, both expected (e.g., underflow) or unexpected (error in [10], handling of intermediate zeros that greatly modifies the proofs). More precisely about the intermediate zeros, let us assume that the result of the renormaliza- tion is (y0, y1, y2, 0, y3, y4, 0, y5, 0, 0, 0) with the yi ̸= 0. Then the pen-and-paper theorem of [10] ensures that OV S P (y0, y1), OV S P (y1, y2), and OV S P (y3, y4). Our theorem formally proved in Coq ensures that OV S P (y0, y1), OV S P (y1, y2), OV S P (y2, y3), OV S P (y3, y4), and OV S P (y4, y5). Interleaving zeros may seem a triviality as zeros could easily be removed, except that this is costly as it in- volved testing each value. For this basic block, it is much more efficient to handle zeros inside the proof than to remove them in the algorithm.
We had some hope to be able to automate that kind of proofs, as many such algorithms exist in the literature, but this experiment has shown that, even from a detailed and reasonable pen-and-paper proof, much remains to be done to get a formal one. This paper shows that complex algorithms in FP arithmetic really need formal proofs: it is otherwise impossible to be certain that special cases have not been forgotten or overlooked in the pen-and-paper proof, as was the case here.

@dpsanders you can trust to EFT, when all values without NaN and Infinity.

[JeffreySarnoff]

@munrocket Now ErrorfreeTransformations.jl#master three_fma (and new, two_fma) should work alright -- let me know if they do not. If you find they work, I will move the changes here too.

[munrocket]

Probably I need to close the issue because it’s become clear what they calculate in rust packages. I think that we need more cites or/and formal proves to trust a new papers. Also we need to know when we need safe overflows and how often they occur.

[munrocket]

Found an error in paper here

def succ (x) :
  if (a < 2.**(−1021)):
    return x + a * 2.**(−1074) //<--

[JeffreySarnoff]

You found an error in what paper? and what should it read?

[munrocket]

I am used successor()/predecessor() from M.Kashiwagi paper it and found an error in this test.

succ(2e-300) === 2.0000000000000004e-300

I am just reporting that this paper have some bugs. I am able to fix it. But anyway it feels little bit strange.

return x + 2.**(−1074) //<--fixed

[JeffreySarnoff]

how is that function given in the paper? do you have a link to this paper?

[JeffreySarnoff]

If you have not seen this paper Computing predecessor and successor in rounding to nearest.

[munrocket]

Yea, this technique is described in “Computing predecessor and successor in rounding to nearest”, but it doesn’t have any source code, so I am used paper from Kashiwagi
https://translate.google.com/translate?hl=en&sl=ja&u=http://verifiedby.me/adiary/pub/kashi/image/201406/nas2014.pdf&prev=search
I am just warned you about bugs. They can be in another functions too, seems that code not properly tested.

[munrocket]

If you know some other paper about safe EFT and correctness of algorithms, I will appreciate.

[JeffreySarnoff]

Yes -- I know that paper and it has glitches.

[Double-Double Building Blocks]
M. Joldes, V. Popescu, and J.M. Muller.
Tight and rigorous error bounds for basic building blocks of double-word arithmetic
2016, working paper.
https://hal.archives-ouvertes.fr/hal-01351529v2/document

[Triple-Double Building Blocks]
Christoph Quirin Lauter.
Basic building blocks for a triple-double intermediate format
2005, research report.
https://hal.inria.fr/inria-00070314/document

https://www.researchgate.net/publication/228568591_Error-free_transformations_in_real_and_complex_floating_point_arithmetic

http://perso.ens-lyon.fr/jean-michel.muller/correct_fma_err.pdf

[munrocket]

thanks!