Skip to content

math.pow: add f16, f80, f128, c_longdouble support #23631

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
+37 −15
Open

math.pow: add f16, f80, f128, c_longdouble support #23631

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 37 additions & 15 deletions lib/std/math/pow.zig
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@ pub fn pow(comptime T: type, x: T, y: T) T {
return math.powi(T, x, y) catch unreachable;
}

if (T != f32 and T != f64) {
@compileError("pow not implemented for " ++ @typeName(T));
if (@typeInfo(T) != .float) {
@compileError("pow not implemented for comptime_float");
}

// pow(x, +-0) = 1 for all x
Expand All @@ -60,15 +60,15 @@ pub fn pow(comptime T: type, x: T, y: T) T {
if (x == 0) {
if (y < 0) {
// pow(+-0, y) = +-inf for y an odd integer
if (isOddInteger(y)) {
if (isOddInteger(T, y)) {
return math.copysign(math.inf(T), x);
}
// pow(+-0, y) = +inf for y an even integer
else {
return math.inf(T);
}
} else {
if (isOddInteger(y)) {
if (isOddInteger(T, y)) {
return x;
} else {
return 0;
Expand Down Expand Up @@ -178,29 +178,37 @@ pub fn pow(comptime T: type, x: T, y: T) T {
return math.scalbn(a1, ae);
}

fn isOddInteger(x: f64) bool {
if (@abs(x) >= 1 << 53) {
// From https://golang.org/src/math/pow.go
// 1 << 53 is the largest exact integer in the float64 format.
fn isOddInteger(comptime T: type, x: T) bool {
if (@abs(x) >= 1 << (math.floatFractionalBits(T) + 1)) {
// From https://golang.org/src/math/pow.go (adapted for generic type T)
// 1 << (floatFractionalBits(T) + 1) is the largest exact integer in the floating point format.
// Any number outside this range will be truncated before the decimal point and therefore will always be
// an even integer.
// Without this check and if x overflows i64 the @intFromFloat(r.ipart) conversion below will panic
// Without this check and if x overflows i128 the @intFromFloat(r.ipart) conversion below will panic
return false;
}
const r = math.modf(x);
return r.fpart == 0.0 and @as(i64, @intFromFloat(r.ipart)) & 1 == 1;
return r.fpart == 0.0 and @as(i128, @intFromFloat(r.ipart)) & 1 == 1;
}

test isOddInteger {
try expect(isOddInteger(math.maxInt(i64) * 2) == false);
try expect(isOddInteger(math.maxInt(i64) * 2 + 1) == false);
try expect(isOddInteger(1 << 53) == false);
try expect(isOddInteger(12.0) == false);
try expect(isOddInteger(15.0) == true);
try expect(isOddInteger(f64, math.maxInt(i64) * 2) == false);
try expect(isOddInteger(f64, math.maxInt(i64) * 2 + 1) == false);
try expect(isOddInteger(f64, 1 << 53) == false);
try expect(isOddInteger(f64, 12.0) == false);
try expect(isOddInteger(f64, 15.0) == true);
}

test pow {
const epsilon = 0.000001;
const epsilon16 = 0.01;

try expect(math.approxEqAbs(f16, pow(f16, 0.0, 3.3), 0.0, epsilon16));
try expect(math.approxEqAbs(f16, pow(f16, 0.8923, 3.3), 0.686572, epsilon16));
try expect(math.approxEqAbs(f16, pow(f16, 0.2, 3.3), 0.004936, epsilon16));
try expect(math.approxEqAbs(f16, pow(f16, 1.5, 3.3), 3.811546, epsilon16));
try expect(math.approxEqAbs(f16, pow(f16, 37.45, 3.3), 155736.7160616, epsilon16));
try expect(math.approxEqAbs(f16, pow(f16, 89.123, 3.3), 2722490.231436, epsilon16));
Comment on lines +210 to +211
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might be a dumb question, but shouldn't these two result in +inf with the results being way over the maximum value of f16?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch! You're right, the expected result for f16 should be +inf.

The test likely passed because the large literal 155736.7... also overflows to +inf when cast to f16, causing approxEqAbs to hit the x == y fast path.

Maybe it's clearer to test explicitly with == math.inf(f16) in this overflow case?

Copy link
Contributor

@fardragon fardragon Apr 23, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I guess that's true. I've mixed up inf comparison and NaN comparison in my mind, but you're right, that's probably what's happening here. But I agree that explicitly comparing to inf or even using math.isInf would be more clear here. At least in my opinion.


try expect(math.approxEqAbs(f32, pow(f32, 0.0, 3.3), 0.0, epsilon));
try expect(math.approxEqAbs(f32, pow(f32, 0.8923, 3.3), 0.686572, epsilon));
Expand All @@ -215,6 +223,20 @@ test pow {
try expect(math.approxEqAbs(f64, pow(f64, 1.5, 3.3), 3.811546, epsilon));
try expect(math.approxEqAbs(f64, pow(f64, 37.45, 3.3), 155736.7160616, epsilon));
try expect(math.approxEqAbs(f64, pow(f64, 89.123, 3.3), 2722490.231436, epsilon));

try expect(math.approxEqAbs(f80, pow(f80, 0.0, 3.3), 0.0, epsilon));
try expect(math.approxEqAbs(f80, pow(f80, 0.8923, 3.3), 0.686572, epsilon));
try expect(math.approxEqAbs(f80, pow(f80, 0.2, 3.3), 0.004936, epsilon));
try expect(math.approxEqAbs(f80, pow(f80, 1.5, 3.3), 3.811546, epsilon));
try expect(math.approxEqAbs(f80, pow(f80, 37.45, 3.3), 155736.7160616, epsilon));
try expect(math.approxEqAbs(f80, pow(f80, 89.123, 3.3), 2722490.231436, epsilon));

try expect(math.approxEqAbs(f128, pow(f128, 0.0, 3.3), 0.0, epsilon));
try expect(math.approxEqAbs(f128, pow(f128, 0.8923, 3.3), 0.686572, epsilon));
try expect(math.approxEqAbs(f128, pow(f128, 0.2, 3.3), 0.004936, epsilon));
try expect(math.approxEqAbs(f128, pow(f128, 1.5, 3.3), 3.811546, epsilon));
try expect(math.approxEqAbs(f128, pow(f128, 37.45, 3.3), 155736.7160616, epsilon));
try expect(math.approxEqAbs(f128, pow(f128, 89.123, 3.3), 2722490.231436, epsilon));
}

test "special" {
Expand Down
Loading