Fix exponent wrapping in Math.frexp(BigFloat) for very large values (…

…#14971)

`BigFloat`s represent their base-`256 ** sizeof(LibGMP::MpLimb)` exponent with a `LibGMP::MpExp` field, but `LibGMP.mpf_get_d_2exp` only returns the base-2 exponent as a `LibC::Long`, so values outside `(2.0.to_big_f ** -0x80000001)...(2.0.to_big_f ** 0x7FFFFFFF)` lead to an exponent overflow on Windows or 32-bit platforms:

```crystal
require "big"

Math.frexp(2.0.to_big_f ** 0xFFFFFFF5)  # => {1.55164027193164307015e+1292913986, -10}
Math.frexp(2.0.to_big_f ** -0xFFFFFFF4) # => {1.61119819150333097422e-1292913987, 13}
Math.frexp(2.0.to_big_f ** 0x7FFFFFFF)  # raises OverflowError
```

This patch fixes it by computing the exponent ourselves.

spec/std/big/big_float_spec.cr

@@ -548,7 +548,23 @@ end
 
 describe "BigFloat Math" do
   it ".frexp" do
+    Math.frexp(0.to_big_f).should eq({0.0, 0})
+    Math.frexp(1.to_big_f).should eq({0.5, 1})
     Math.frexp(0.2.to_big_f).should eq({0.8, -2})
+    Math.frexp(2.to_big_f ** 63).should eq({0.5, 64})
+    Math.frexp(2.to_big_f ** 64).should eq({0.5, 65})
+    Math.frexp(2.to_big_f ** 200).should eq({0.5, 201})
+    Math.frexp(2.to_big_f ** -200).should eq({0.5, -199})
+    Math.frexp(2.to_big_f ** 0x7FFFFFFF).should eq({0.5, 0x80000000})
+    Math.frexp(2.to_big_f ** 0x80000000).should eq({0.5, 0x80000001})
+    Math.frexp(2.to_big_f ** 0xFFFFFFFF).should eq({0.5, 0x100000000})
+    Math.frexp(1.75 * 2.to_big_f ** 0x123456789).should eq({0.875, 0x12345678A})
+    Math.frexp(2.to_big_f ** -0x80000000).should eq({0.5, -0x7FFFFFFF})
+    Math.frexp(2.to_big_f ** -0x80000001).should eq({0.5, -0x80000000})
+    Math.frexp(2.to_big_f ** -0x100000000).should eq({0.5, -0xFFFFFFFF})
+    Math.frexp(1.75 * 2.to_big_f ** -0x123456789).should eq({0.875, -0x123456788})
+    Math.frexp(-(2.to_big_f ** 0x7FFFFFFF)).should eq({-0.5, 0x80000000})
+    Math.frexp(-(2.to_big_f ** -0x100000000)).should eq({-0.5, -0xFFFFFFFF})
   end
 
   it ".sqrt" do

src/big/big_float.cr

@@ -537,15 +537,24 @@ end
 module Math
   # Decomposes the given floating-point *value* into a normalized fraction and an integral power of two.
   def frexp(value : BigFloat) : {BigFloat, Int64}
-    LibGMP.mpf_get_d_2exp(out exp, value) # we need BigFloat frac, so will skip Float64 one.
+    return {BigFloat.zero, 0_i64} if value.zero?
+
+    # We compute this ourselves since `LibGMP.mpf_get_d_2exp` only returns a
+    # `LibC::Long` exponent, which is not sufficient for 32-bit `LibC::Long` and
+    # 32-bit `LibGMP::MpExp`, e.g. on 64-bit Windows.
+    # Since `0.5 <= frac.abs < 1.0`, the radix point should be just above the
+    # most significant limb, and there should be no leading zeros in that limb.
+    leading_zeros = value.@mpf._mp_d[value.@mpf._mp_size.abs - 1].leading_zeros_count
+    exp = 8_i64 * sizeof(LibGMP::MpLimb) * value.@mpf._mp_exp - leading_zeros
+
     frac = BigFloat.new do |mpf|
-      if exp >= 0
-        LibGMP.mpf_div_2exp(mpf, value, exp)
-      else
-        LibGMP.mpf_mul_2exp(mpf, value, -exp)
-      end
+      # remove leading zeros in the most significant limb
+      LibGMP.mpf_mul_2exp(mpf, value, leading_zeros)
+      # reset the exponent manually
+      mpf.value._mp_exp = 0
     end
-    {frac, exp.to_i64}
+
+    {frac, exp}
   end
 
   # Calculates the square root of *value*.
@@ -559,21 +568,3 @@ module Math
     BigFloat.new { |mpf| LibGMP.mpf_sqrt(mpf, value) }
   end
 end
-
-# :nodoc:
-struct Crystal::Hasher
-  def self.reduce_num(value : BigFloat)
-    float_normalize_wrap(value) do |value|
-      # more exact version of `Math.frexp`
-      LibGMP.mpf_get_d_2exp(out exp, value)
-      frac = BigFloat.new do |mpf|
-        if exp >= 0
-          LibGMP.mpf_div_2exp(mpf, value, exp)
-        else
-          LibGMP.mpf_mul_2exp(mpf, value, -exp)
-        end
-      end
-      float_normalize_reference(value, frac, exp)
-    end
-  end
-end