fixes #6353 (#6951)

2 files changed, 57 insertions, 4 deletions

diff --git a/lib/pure/math.nim b/lib/pure/math.nim
index a9dabfa48..cbd04a145 100644
--- a/lib/pure/math.nim
+++ b/lib/pure/math.nim
@@ -351,15 +351,19 @@ proc round*[T: float32|float64](x: T, places: int = 0): T =
     result = round0(x*mult)/mult
 
 when not defined(JS):
-  proc frexp*(x: float32, exponent: var int): float32 {.
+  proc c_frexp*(x: float32, exponent: var int32): float32 {.
     importc: "frexp", header: "<math.h>".}
-  proc frexp*(x: float64, exponent: var int): float64 {.
+  proc c_frexp*(x: float64, exponent: var int32): float64 {.
     importc: "frexp", header: "<math.h>".}
+  proc frexp*[T, U](x: T, exponent: var U): T =
     ## Split a number into mantissa and exponent.
     ## `frexp` calculates the mantissa m (a float greater than or equal to 0.5
     ## and less than 1) and the integer value n such that `x` (the original
     ## float value) equals m * 2**n. frexp stores n in `exponent` and returns
     ## m.
+    var exp: int32
+    result = c_frexp(x, exp)
+    exponent = exp
 else:
   proc frexp*[T: float32|float64](x: T, exponent: var int): T =
     if x == 0.0:
@@ -368,9 +372,14 @@ else:
     elif x < 0.0:
       result = -frexp(-x, exponent)
     else:
-      var ex = floor(log2(x))
-      exponent = round(ex)
+      var ex = trunc(log2(x))
+      exponent = int(ex)
       result = x / pow(2.0, ex)
+      if abs(result) >= 1:
+        inc(exponent)
+        result = result / 2
+      if exponent == 1024 and result == 0.0:
+        result = 0.99999999999999988898
 
 proc splitDecimal*[T: float32|float64](x: T): tuple[intpart: T, floatpart: T] =
   ## Breaks `x` into an integral and a fractional part.
diff --git a/tests/stdlib/tfrexp1.nim b/tests/stdlib/tfrexp1.nim
new file mode 100644
index 000000000..c6bb2b38c
--- /dev/null
+++ b/tests/stdlib/tfrexp1.nim
@@ -0,0 +1,44 @@
+discard """
+  targets: "js c c++"
+  output: '''ok'''
+"""
+
+import math
+import strformat
+
+const manualTest = false
+
+proc frexp_test(lo, hi, step: float64) =
+  var exp: int
+  var frac: float64
+
+  var eps = 1e-15.float64
+
+  var x:float64 = lo
+  while x <= hi:
+    frac = frexp(x.float, exp)
+    let rslt = pow(2.0, float(exp)) * frac
+
+    doAssert(abs(rslt - x) < eps)
+
+    when manualTest:
+      echo fmt("x: {x:10.3f} exp: {exp:4d} frac: {frac:24.20f} check: {$(abs(rslt - x) < eps):-5s} {rslt: 9.3f}")
+    x += step
+
+when manualTest:
+  var exp: int
+  var frac: float64
+
+  for flval in [1.7976931348623157e+308, -1.7976931348623157e+308, # max, min float64
+                3.4028234663852886e+38, -3.4028234663852886e+38,   # max, min float32
+                4.9406564584124654e-324, -4.9406564584124654e-324, # smallest/largest positive/negative float64
+                1.4012984643248171e-45, -1.4012984643248171e-45,   # smallest/largest positive/negative float32
+                2.2250738585072014e-308, 1.1754943508222875e-38]:  # smallest normal float64/float32
+    frac = frexp(flval, exp)
+    echo fmt("{flval:25.16e}, {exp: 6d}, {frac: .20f} {frac * pow(2.0, float(exp)): .20e}")
+
+  frexp_test(-1000.0, 1000.0, 0.0125)
+else:
+  frexp_test(-1000000.0, 1000000.0, 0.125)
+
+echo "ok"