Bug in std.math.pow with large powers?

Explain
1

The following code prints nan, even though one would reasonably assume it to print 1.

std.debug.print("{}\n", .{std.math.pow(f32, -1, 1e20)});                 // nan?

std.math.pow seems to return nan whenever the base is negative and the power is very large. Curiously, the following 2 lines print what you would expect:

std.debug.print("{}\n", .{std.math.pow(f32, -1, std.math.inf(f32))});    // 1
std.debug.print("{}\n", .{std.math.pow(f32, -std.math.inf(f32), 1e20)}); // inf

Looking at the doc comment for `std.math.pow`, I found this line interesting: 
pow(x, y)      = nan for finite x < 0 and finite non-integer y

It’s the only special case that almost describes our first case. However, if I understand floats correctly, 1e20 and other similarly large numbers (e.g., 1.00000000001e10) should be considered integers because at that size, the gaps between consecutive f32s are so large that f32 can effectively only represent integers.


Running the code through a debugger reveals the offending code: 
const r1 = math.modf(@abs(y));
var yi = r1.ipart;
var yf = r1.fpart;

if (yf != 0 and x < 0) {
    return math.nan(T);
}
if (yi >= 1 << (@typeInfo(T).float.bits - 1)) {
    return @exp(y * @log(x)); // !!
}

The log of a negative number is nan, which then propagates through and we get a nan returned to us. This looks unintentional but I’m not sure. I’m not even sure if the condition right above should be using @typeInfo(T).float.bits - 1 or math.floatFractionalBits(T).

2

I think you’ve answered your question. std.math.pow(f32, -1, 1e20) == nan is intentional.

Also convered by the doc comment and working as expected.

“Almost” doesn’t count. It either fits the description or not.

If you can make a new pow function that extends the valid inputs, I’m sure it would make a valuable contribution.

3

1e20 is an integer.

Either way, the implementation appears to be wrong. IEEE 754-2008 lists the following special cases:

For the pow function (integral exponents get special treatment):
pow(x, ±0) is 1 for any x (even a zero, quiet NaN, or infinity)
pow(±0, y) is ±∞ and signals the divideByZero exception for y an odd integer <0
pow(±0, −∞) is +∞ with no exception
pow(±0, +∞) is +0 with no exception
pow(±0, y) is +∞ and signals the divideByZero exception for finite y<0 and not an odd integer
pow(±0, y) is ±0 for finite y>0 an odd integer
pow(±0, y) is +0 for finite y>0 and not an odd integer
pow(−1, ±∞) is 1 with no exception
pow(+1, y) is 1 for any y (even a quiet NaN)
pow(x, y) signals the invalid operation exception for finite x<0 and finite non-integer y

1 Like
4

I looked a bit -acutally it took me quite a while- at std.math.pow and found an explanation:

For y values that are greater than or equal to 2^31 a code path in std.math.pow is chosen where the result is calculated using @exp(y) * @log(x). Because x is negative, @loglogloglogloglogloglogloglogloglogloglogloglogloglogloglogloglogloglogloglogloglogloglogloglog returns nan. Here is the code I used to find the issue.

The crucial part is here

https://codeberg.org/ziglang/zig/src/branch/master/lib/std/math/pow.zig#L126

const std = @import("std");

pub fn main() void {
    const T = f32;
    const x: T = -1.0;
    const y_one: T = 2147483648  - 128;  // 2**31 - 2**7
    std.debug.print("y_highest_one: 0b{b} => {}\n", .{@as(u32, @bitCast(y_one)),std.math.pow(T, x, y_one)});
    const y_nan: T = 2147483648 ; // 2**31 => nan
    std.debug.print(" y_lowest_nan: 0b{b} => {}\n", .{@as(u32, @bitCast(y_nan)), std.math.pow(T, x, y_nan)}); // nan

    // depending on the following contidtion, the result is calculated using @exp(y * @log(x))
    // for y_one it is calculated by other means
    std.debug.assert(modf_one.ipart < 1 << (@typeInfo(T).float.bits - 1));
    // for y_nan and higher @exp(y * @log(x)) is used
    std.debug.assert(modf_nan.ipart >= 1 << (@typeInfo(T).float.bits - 1));
    // this returnsn nan
    std.debug.assert(std.math.isNan(@exp(y_nan) * @log(x)));

    // the issue is not y though
    std.debug.assert(!std.math.isNan(@exp(y_nan) ));

    // but log(-1)!!!
    std.debug.assert(std.math.isNan(@log(x)));

}

I guess it is a bug, but I am not sure.

Edit: Fixed codeberg link, removed some noise in code

1 Like
5

Ah, so my interpretation was correct.

Will take a look into std to see how it implements those exceptions (I noticed std.math.raiseDivByZero and co. are… empty functions?) and try to implement a fix.

6

You can reference a specific line by adding a url fragment #Lxxx.
https://codeberg.org/ziglang/zig/src/branch/master/lib/std/math/pow.zig#L126 turns into this box below:

2 Likes
7

Ok, I get it now. They meant pow(x,y) = nan for (finite x < 0 and finite non-integer y). I thought they meant pow(x,y) = nan for the cases: 1) (finite x < 0) and 2) (finite non-integer y).

Signaling an invalid operation would cause this function to have side-effects, which would be both unexpected and hinder optimizations. (at least it will once they implement raiseInvalid)

8

I don’t think most languages actually raise FPU exceptions, like is the case with C# (found this during my research).

I’m going to interpret “signals the invalid operation exception” as “return NaN”

That doesn’t seem to be the intention of the Zig code either, as it doesn’t call std.math.raiseInvalid. I plan on sticking with that and just returning nan.

9

The fix turned out to be a simple one-liner lol.

- return @exp(y * @log(x));
+ return @exp(y * @log(@abs(x)));

I did find another weird edge case to fix, anyways here’s the PR for anyone interested:

4 Likes