From b0b4b498c864937ba0c7d718fdbe4f0a1347bc1e Mon Sep 17 00:00:00 2001 From: dlesnoff <54949944+dlesnoff@users.noreply.github.com> Date: Thu, 10 Oct 2024 19:30:40 +0200 Subject: [PATCH] =?UTF-8?q?std/math:=C2=A0Add=20`^`=20overload=20for=20flo?= =?UTF-8?q?at32=20and=20float64=20(#20898)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 Co-authored-by: metagn (cherry picked from commit e9a4d096abe8c4abf5cf5ca736828e44ffd691f8) --- changelog.md | 2 ++ lib/pure/math.nim | 64 ++++++++++++++++++++++++++++++++++-------- tests/stdlib/tmath.nim | 60 +++++++++++++++++++++++++++++++++++++-- 3 files changed, 112 insertions(+), 14 deletions(-) diff --git a/changelog.md b/changelog.md index 39f5292451..eb3f123096 100644 --- a/changelog.md +++ b/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 diff --git a/lib/pure/math.nim b/lib/pure/math.nim index ed7d2382fd..dc8257bfb6 100644 --- a/lib/pure/math.nim +++ b/lib/pure/math.nim @@ -78,7 +78,7 @@ when not defined(js) and not defined(nimscript): # C importc: "frexpf", header: "".} func c_frexp2(x: cdouble, exponent: var cint): cdouble {. importc: "frexp", header: "".} - + type div_t {.importc, header: "".} = object quot: cint @@ -89,13 +89,13 @@ when not defined(js) and not defined(nimscript): # C lldiv_t {.importc, header: "".} = object quot: clonglong rem: clonglong - + when cint isnot clong: func divmod_c(x, y: cint): div_t {.importc: "div", header: "".} when clong isnot clonglong: func divmod_c(x, y: clonglong): lldiv_t {.importc: "lldiv", header: "".} func divmod_c(x, y: clong): ldiv_t {.importc: "ldiv", header: "".} - 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: "".} = ## 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`. ## diff --git a/tests/stdlib/tmath.nim b/tests/stdlib/tmath.nim index 22e5f7d88f..b28ec41a31 100644 --- a/tests/stdlib/tmath.nim +++ b/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