Avoid overflows when converting rationals with large denominator and/or numerator to floating point

[fingolfin]

Fixes #2734

[coveralls]

Coverage increased (+3.0e-05%) to 93.73% when pulling 7d8f948 on fingolfin:mh/fix-float-overflow into e62d395 on gap-system:master.

[fingolfin]

Unfortunately, this still isn't good enough if the denominator is still "large" after reduction:

gap> Float((10^309 - 1) / 10^309);
nan

but of course it should give 1.

One idea would be to multiply the "fractional" part by something like 2^53, turn that into an integer, and finally divide by 2.^53, but I'd bet that also has its pitfalls. Perhaps there is a more clever way to do all this? @laurentbartholdi might have an idea?

[wilfwilson]

The change you've made is obviously mathematically correct, and improves the situation.

How do you want to proceed with this PR? If you're happy for it to be merged, then go for it. But if you want this PR to not be merged until it contains a complete fix, then I suggest adding the 'Do not merge' label.

[fingolfin]

I was kinda hoping that @laurentbartholdi or @stevelinton who wrote the float support in GAP would have a comment on this...

[stevelinton]

The very early float support I did, if I recall correctly, did not provide a conversion from Rationals -- my plan was to keep rationals as separate from floats as possible. So I don't think I have any special insight here.

I would still argue that such conversions should not be part of inner loops of calculations -- they are useful for UI purposes and a few "hybrid" algorithms, but performance is not really a concern.

A more accurate conversion is necessarily going to have to know details of the currently selected FP implementation, essentially the precision, so that it can find the nearest FP value to the rational. For IEEE FP that will means first identifying which k so that 2^k < Abs(x) < 2^{k+1} (k may be negative) and then finding the the closest multiple of 2^{k-53} (or something similar) to x. I don't think there is a short and slick way to do it.

[laurentbartholdi]

I'm fine with all solutions, but would prefer just throwing an error when the conversion cannot be done easily.

…

On Fri, Nov 20, 2020 at 12:48 PM Steve Linton ***@***.***> wrote: The very early float support I did, if I recall correctly, did not provide a conversion from Rationals -- my plan was to keep rationals as separate from floats as possible. So I don't think I have any special insight here. I would still argue that such conversions should not be part of inner loops of calculations -- they are useful for UI purposes and a few "hybrid" algorithms, but performance is not really a concern. A more accurate conversion is necessarily going to have to know details of the currently selected FP implementation, essentially the precision, so that it can find the nearest FP value to the rational. For IEEE FP that will means first identifying which k so that 2^k < Abs(x) < 2^{k+1} (k may be negative) and then finding the the closest multiple of 2^{k-53} (or something similar) to x. I don't think there is a short and slick way to do it. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#4144 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AARAQUA3GRJACSBMCEBEK43SQZJSTANCNFSM4STGRFHQ> .

-- Laurent Bartholdi laurent.bartholdi<at>gmail<dot>com Mathematisches Institut, Georg-August Universität zu Göttingen Bunsenstrasse 3-5, D-37073 Göttingen, Germany

[fingolfin]

@laurentbartholdi I think the conversion can always be done "easily", the problem is how accurate it is (right now it can have zero accurate digits...). I am not sure whether "throwing an error when the conversion cannot be done easily" is substantially easier than doing the conversion?

Perhaps it would make more sense to always raise an error in the generic conversion methods, something like "not implemented for this floating point type", and then instead provide methods only for the machine integers, for which we can implement the conversion fairly accurately (I think)?

[laurentbartholdi]

I don't have strong feelings about this; but looking at the PR I think it's fixing the wrong problem: it would be better to ask MakeFloat(filter,10^309+1) to trigger an error, rather than to attempt to force this huge integer into IEEE754. The code for converting rationals to fp should be as clean as possible: convert numerator and denominator, and divide.

…

On Fri, Nov 20, 2020 at 5:27 PM Max Horn ***@***.***> wrote: @laurentbartholdi <https://github.com/laurentbartholdi> I think the conversion can always be done "easily", the problem is how accurate it is (right now it can have zero accurate digits...). I am not sure whether "throwing an error when the conversion cannot be done easily" is substantially easier than doing the conversion? Perhaps it would make more sense to always raise an error in the generic conversion methods, something like "not implemented for this floating point type", and then instead provide methods only for the machine integers, for which we can implement the conversion fairly accurately (I think)? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#4144 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AARAQUDIH54RMNA45QI4UGDSQ2KG7ANCNFSM4STGRFHQ> .

-- Laurent Bartholdi laurent.bartholdi<at>gmail<dot>com Mathematisches Institut, Georg-August Universität zu Göttingen Bunsenstrasse 3-5, D-37073 Göttingen, Germany

[fingolfin]

The question is, what is the goal? I would have thought that it is "produce a useful conversion which minimizes the numerical error". I don't see the value in producing one that is "as clean as possible" in the vein @laurentbartholdi mentions; this is how the code this PR is replacing looks like -- and as documented here, this needlessly produces atrocious numerical errors, and overall seems mostly useless. I think not providing any conversion would be preferable over that, as coding this "clean" conversion is trivial anyway.

@stevelinton has a point in that we can implement a good quality conversion for machine floats, as we know their precision. So, I am tempted to do that. The "generic" conversion then can either be implemented by first converting to machine floats and then to the actual target (assuming that all reasonable float implementations will allow conversion from/to machine floats). Or we replace the "generic" conversion with a "not implemented" errors: I don't think providing a bad default is good, indeed I think it's actively harmful.

So, overall, it seems we all agree that an error would be best in the general case :-)

[codecov]

Codecov Report

Merging #4144 (bd0104a) into master (b1d1416) will increase coverage by 11.17%.
The diff coverage is 100.00%.

❗ Current head bd0104a differs from pull request most recent head ae974ef. Consider uploading reports for the commit ae974ef to get more accurate results

@@             Coverage Diff             @@
##           master    #4144       +/-   ##
===========================================
+ Coverage   82.23%   93.40%   +11.17%     
===========================================
  Files         685      716       +31     
  Lines      290627   812429   +521802     
  Branches     8626        0     -8626     
===========================================
+ Hits       238994   758888   +519894     
- Misses      49743    53541     +3798     
+ Partials     1890        0     -1890     

Impacted Files	Coverage Δ
lib/float.gi	47.81% <ø> (-0.47%)	⬇️
lib/ieee754.g	100.00% <100.00%> (ø)
src/gapw95.c	0.00% <0.00%> (-100.00%)	⬇️
lib/memory.gi	43.12% <0.00%> (-39.92%)	⬇️
lib/methsel.g	51.16% <0.00%> (-24.70%)	⬇️
lib/colorprompt.g	46.15% <0.00%> (-24.16%)	⬇️
lib/matobjplist.gi	42.79% <0.00%> (-20.13%)	⬇️
src/io.c	58.09% <0.00%> (-18.99%)	⬇️
lib/mgmfree.gi	58.53% <0.00%> (-16.47%)	⬇️
lib/matobj.gi	58.95% <0.00%> (-13.95%)	⬇️
... and 397 more