std/math: Add ^ overload for float32 and float64 (#20898)

I have added a new overload of `^` for float exponents.
Is two overloads for `float32` and `float64` better than just one
overload with `SomeFloat` type ?
I guess this would not work with `SomeFloat`, as `pow` is not defined
for `float`.

Another remark. Maybe we should catch exponents with 0.5 and call `sqrt`
instead ?

---------

Co-authored-by: Clay Sweetser <Varriount@users.noreply.github.com>
Co-authored-by: metagn <metagngn@gmail.com>
(cherry picked from commit e9a4d096ab)

changelog.md

@@ -10,6 +10,8 @@ rounding guarantees (via the
 
 ## Standard library additions and changes
 
+[//]: # "Changes:"
+- `std/math` The `^` symbol now supports floating-point as exponent in addition to the Natural type.
 
 ## Language changes
 

lib/pure/math.nim

@@ -78,7 +78,7 @@ when not defined(js) and not defined(nimscript): # C
       importc: "frexpf", header: "<math.h>".}
   func c_frexp2(x: cdouble, exponent: var cint): cdouble {.
       importc: "frexp", header: "<math.h>".}
-  
+
   type
     div_t {.importc, header: "<stdlib.h>".} = object
       quot: cint
@@ -89,13 +89,13 @@ when not defined(js) and not defined(nimscript): # C
     lldiv_t {.importc, header: "<stdlib.h>".} = object
       quot: clonglong
       rem: clonglong
-  
+
   when cint isnot clong:
     func divmod_c(x, y: cint): div_t {.importc: "div", header: "<stdlib.h>".}
   when clong isnot clonglong:
     func divmod_c(x, y: clonglong): lldiv_t {.importc: "lldiv", header: "<stdlib.h>".}
   func divmod_c(x, y: clong): ldiv_t {.importc: "ldiv", header: "<stdlib.h>".}
-  func divmod*[T: SomeInteger](x, y: T): (T, T) {.inline.} = 
+  func divmod*[T: SomeInteger](x, y: T): (T, T) {.inline.} =
     ## Specialized instructions for computing both division and modulus.
     ## Return structure is: (quotient, remainder)
     runnableExamples:
@@ -628,11 +628,11 @@ when not defined(js): # C
   func pow*(x, y: float64): float64 {.importc: "pow", header: "<math.h>".} =
     ## Computes `x` raised to the power of `y`.
     ##
-    ## To compute the power between integers (e.g. 2^6),
-    ## use the `^ func <#^,T,Natural>`_.
+    ## You may use the `^ func <#^, T, U>`_ instead.
     ##
     ## **See also:**
-    ## * `^ func <#^,T,Natural>`_
+    ## * `^ (SomeNumber, Natural) func <#^,T,Natural>`_
+    ## * `^ (SomeNumber, SomeFloat) func <#^,T,U>`_
     ## * `sqrt func <#sqrt,float64>`_
     ## * `cbrt func <#cbrt,float64>`_
     runnableExamples:
@@ -828,14 +828,14 @@ else: # JS
       doAssert -6.5 mod  2.5 == -1.5
       doAssert  6.5 mod -2.5 ==  1.5
       doAssert -6.5 mod -2.5 == -1.5
-  
-  func divmod*[T:SomeInteger](num, denom: T): (T, T) = 
+
+  func divmod*[T:SomeInteger](num, denom: T): (T, T) =
     runnableExamples:
       doAssert  divmod(5, 2) ==  (2, 1)
       doAssert divmod(5, -3) == (-1, 2)
     result[0] = num div denom
     result[1] = num mod denom
-  
+
 
 func round*[T: float32|float64](x: T, places: int): T =
   ## Decimal rounding on a binary floating point number.
@@ -1186,8 +1186,8 @@ func `^`*[T: SomeNumber](x: T, y: Natural): T =
   ## `pow <#pow,float64,float64>`_ for negative exponents.
   ##
   ## **See also:**
-  ## * `pow func <#pow,float64,float64>`_ for negative exponent or
-  ##   floats
+  ## * `^ func <#^,T,U>`_ for negative exponent or floats
+  ## * `pow func <#pow,float64,float64>`_ for `float32` or `float64` output
   ## * `sqrt func <#sqrt,float64>`_
   ## * `cbrt func <#cbrt,float64>`_
   runnableExamples:
@@ -1211,6 +1211,48 @@ func `^`*[T: SomeNumber](x: T, y: Natural): T =
         break
       x *= x
 
+func isInteger(y: SomeFloat): bool =
+  ## Determines if a float represents an integer
+  return round(y) == y
+
+func `^`*[T: SomeNumber, U: SomeFloat](x: T, y: U): float =
+  ## Computes `x` to the power of `y`.
+  ##
+  ## Error handling follows the C++ specification even for the JS backend
+  ## https://en.cppreference.com/w/cpp/numeric/math/pow
+  ##
+  ## **See also:**
+  ## * `^ func <#^,T,Natural>`_
+  ## * `pow func <#pow,float64,float64>`_ for `float32` or `float64` output
+  ## * `sqrt func <#sqrt,float64>`_
+  ## * `cbrt func <#cbrt,float64>`_
+  runnableExamples:
+    doAssert almostEqual(5.5 ^ 2.2, 42.540042248725975)
+    doAssert 1.0 ^ Inf == 1.0
+  let
+    isZero_x: bool = (x == 0.0 or x == -0.0)
+    isNegZero: bool = classify(x) == fcNegZero
+    isPosZero: bool = classify(x) == fcZero
+    yIsFinite: bool = (y != Inf and y != -Inf)
+    yIsOddInteger: bool = (isInteger(y) and yIsFinite and (abs(int(y) mod 2) == 1))
+
+  assert not(isPosZero and y < 0 and yIsOddInteger)
+  assert not(isNegZero and y < 0 and yIsOddInteger)
+  assert not(isZero_x and y < 0 and y != -Inf)
+  assert not(isZero_x and y == -Inf)
+  assert not(x < 0 and not isInteger(x) and yIsFinite and not yIsOddInteger)
+  when defined(js):
+    # JS behavior follows an old version of IEEE 754 for compatibility reasons
+    # See https://262.ecma-international.org/#sec-numeric-types-number-exponentiate
+    if (x == 1.0 or x == -1.0) and not yIsFinite:
+      float(1.0)
+    elif x == 1.0 and y.isNan():
+      float(1.0)
+    else:
+      float(pow(x, y))
+  else:
+    float(pow(x, y))
+
 func gcd*[T](x, y: T): T =
   ## Computes the greatest common (positive) divisor of `x` and `y`.
   ##

tests/stdlib/tmath.nim

@@ -186,7 +186,7 @@ template main() =
     when not defined(nimTmathCase2):
       doAssert classify(trunc(f_nan.float32)) == fcNan
     doAssert classify(trunc(0.0'f32)) == fcZero
-  
+
   block: # divmod
     doAssert divmod(int.high, 1) == (int.high, 0)
     doAssert divmod(-1073741823, 17) == (-63161283, -12)
@@ -201,7 +201,7 @@ template main() =
         doAssertRaises(OverflowDefect, (discard divmod(clong.low, -1.clong)))
         doAssertRaises(OverflowDefect, (discard divmod(clonglong.low, -1.clonglong)))
         doAssertRaises(DivByZeroDefect, (discard divmod(1, 0)))
-  
+
   block: # log
     doAssert log(4.0, 3.0) ==~ ln(4.0) / ln(3.0)
     doAssert log2(8.0'f64) == 3.0'f64
@@ -251,7 +251,61 @@ template main() =
     doAssert: compiles(5.5 ^ 2.int8)
     doAssert: compiles(5.5 ^ 2.uint)
     doAssert: compiles(5.5 ^ 2.uint8)
-    doAssert: not compiles(5.5 ^ 2.2)
+
+  block:
+    doAssert: compiles(5.5 ^ 2.2)
+    when not defined(danger):
+      doAssertRaises(AssertionDefect): discard (0.0 ^ -5.0)
+      doAssertRaises(AssertionDefect): discard (-0.0 ^ -5.0)
+      doAssertRaises(AssertionDefect): discard (-0.0 ^ -5.3)
+      doAssertRaises(AssertionDefect): discard (-0.0 ^ -4.0)
+      doAssertRaises(AssertionDefect): discard (-0.0 ^ -5.3)
+      doAssertRaises(AssertionDefect): discard (-0.0 ^ -4.0)
+      doAssertRaises(AssertionDefect): discard (-0.0 ^ -4.0)
+      # Base finite, negative and exponent finite, non-integer returns NaN and raises Error
+      doAssertRaises(AssertionDefect): discard (-5.5 ^ 2.2)
+    doAssert: -1.0 ^ Inf == 1.0
+    doAssert: -1.0 ^ -Inf == 1.0
+    doAssert: 1.0 ^ Inf == 1.0
+    doAssert: 1.0 ^ -Inf == 1.0
+    doAssert: 1.0 ^ -5.4 == 1.0
+    doAssert: 1.0 ^ 1.0 == 1.0
+
+    # ^-Inf or ^Inf returns 0 or Inf depending on base absolute value relative to 1
+    doAssert: 0.5 ^ -Inf == Inf
+    doAssert: -0.5 ^ -Inf == Inf
+    doAssert: 1.01 ^ -Inf == 0.0
+    doAssert: -1.01 ^ -Inf == 0.0
+    doAssert: 0.5 ^ Inf == 0.0
+    doAssert: -0.5 ^ Inf == 0.0
+    doAssert: 1.01 ^ Inf == Inf
+    doAssert: -1.01 ^ Inf == Inf
+
+    # -Inf base (the sign depends on the exponent parity)
+    doAssert: -Inf ^ -5.0 == -0.0
+    doAssert: -Inf ^ -4.0 == 0.0
+    doAssert: -Inf ^ -5.1 == 0.0
+
+    doAssert: -Inf ^ 3.0 == -Inf
+    doAssert: -Inf ^ 3.1 == Inf
+    doAssert: -Inf ^ 2.0 == Inf
+
+    doAssert: Inf ^ -5.0 == 0.0
+    doAssert: Inf ^ -5.1 == 0.0
+
+    doAssert: Inf ^ 4.3 == Inf
+    doAssert: Inf ^ 5.0 == Inf
+
+    doAssert: -Inf ^ -Inf == 0.0
+    doAssert: -Inf ^ Inf == Inf
+    doAssert: Inf ^ -Inf == 0.0
+    doAssert: Inf ^ Inf == Inf
+
+  block:
+    # doAssert: 1.0 ^ NaN == 1.0
+    # doAssert: NaN ^ 0.0 == 1.0
+    doAssert: (0.0 ^ NaN).isNaN
+    doAssert: (NaN ^ 1.0).isNaN
 
   block: # isNaN
     doAssert NaN.isNaN