Incorrect rounding mode used for fromIntegral :: Integer -> Double/Float

Summary

When converting from Integer to Double (or Float), GHC uses round-toward-zero rounding mode. It should use the default mode round-toward-nearest-ties-to-even instead.

It appears up till GHC 7.8.3, GHC did indeed use this default mode, but something changed after that at some point, and at least in 8.6.5, the incorrect mode is used.

Steps to reproduce

As discussed in https://stackoverflow.com/questions/58035248/a-haskell-floating-point-calculation-anomaly/58037815#58037815

This expression:

*Main> (fromIntegral $ 4141414141414141*4141414141414141) :: Double
1.7151311090705025e31

Produces an unexpected result since the resulting integer is converted using round-toward-zero. While the report is rather silent on rounding-modes, it's general accepted practice that round-nearest-ties-to-even is the default rounding mode to use. As discussed in the stack-overflow question linked, this results in bad results in diverges from other languages. For instance, the equivalent Python produces:

>>> float(long(4141414141414141)*long(4141414141414141))
1.7151311090705027e+31

which is a much better result for this expression. (Note the final digit before the exponent, 5 vs 7.)

Expected behavior

Produce the better conversion result by using round-nearest-ties-to-even rounding mode of IEEE754.

Analysis

The issue goes to this function: https://hackage.haskell.org/package/integer-gmp-1.0.2.0/docs/src/GHC.Integer.Type.html#doubleFromInteger

which reads:

doubleFromInteger :: Integer -> Double#
doubleFromInteger (S# m#) = int2Double# m#
doubleFromInteger (Jp# bn@(BN# bn#))
    = c_mpn_get_d bn# (sizeofBigNat# bn) 0#
doubleFromInteger (Jn# bn@(BN# bn#))
    = c_mpn_get_d bn# (negateInt# (sizeofBigNat# bn)) 0#

which in turn calls: https://github.com/ghc/ghc/blob/master/libraries/integer-gmp/cbits/wrappers.c#L183-L190:

/* Convert bignum to a `double`, truncating if necessary
 * (i.e. rounding towards zero).
 *
 * sign of mp_size_t argument controls sign of converted double
 */
HsDouble
integer_gmp_mpn_get_d (const mp_limb_t sp[], const mp_size_t sn,
                       const HsInt exponent)
{

This function decidedly uses round-towards-zero. Instead, a version that uses round-nearest-ties-to-even should be called.

Integer -> Float conversion suffers from the same problem, as it goes through the Double value. So, fixing the first would fix the second as well.

Environment

• GHC version used: 8.6.5 and newer. On the stack-overflow discussion, people noted that GHC 7.8.3 returns the correct value. So, the behavior somehow changed in between these releases.

Optional:

• Operating System: N/A
• System Architecture: N/A