Replaces the 2017 Go port of std.math.pow (the last one in std.math) with the musl algorithm (ARM optimized-routines, also used by glibc). The Go algorithm exponentiates by squaring per exponent bit, which roughly doubles the relative error with every significant bit of the exponent (golang/go#25270). This is the same problem that caused the bundled libc pow to be reverted to musl C in #31794 (see also #31207).
f64: faithful port of musl'spow, with both the fma and non-fma code paths, selected at comptime like musl's__FP_FAST_FMA.f16(new) andf32: computed asexp2(y * log2(|x|))inf64with a single rounding to the result type.- C99 special-case layer unchanged;
f80/f128remainTODO.
Validation, against a 128-bit MPFR reference using the methodology of Gladman/Innocente/Zimmermann, "Accuracy of Mathematical Functions":
- over 8M bit-space samples,
f64max error improves from 2777.9 ulp to 0.505 ulp andf32from 860.6 ulp to 0.500 ulp, with zero results off by more than 1 ulp (previously 62737). Targeted probing of the old algorithm shows the error passing tens of millions of ulps as the exponent gains bits; the new maximum stays within musl's analytical bound of 0.54 ulp, and the known worst case of this code base measures 0.522908 ulp, matching the published 0.523. - bit-identical to the vendored musl
pow.ccompiled withzig cc, over 8M bit-space samples on each of the fma and non-fma paths. - special-case truth table (
+-0,+-1,+-inf,nan, subnormals, odd/even/non-integer/huge/tiny exponents) passes forf16/f32/f64and matches musl bit for bit onf64;nanresults keep Zig's canonical NaN, matching existingstd.mathbehavior, rather than musl's payload propagation. - one test expectation corrected:
pow(f32, 89.123, 3.3)is 2722489.75 when correctly rounded; the old value encoded the Go algorithm's error.
Replaces the 2017 Go port of `std.math.pow` (the last one in `std.math`) with the musl algorithm (ARM optimized-routines, also used by glibc). The Go algorithm exponentiates by squaring per exponent bit, which roughly doubles the relative error with every significant bit of the exponent ([golang/go#25270](https://github.com/golang/go/issues/25270)). This is the same problem that caused the bundled libc `pow` to be reverted to musl C in #31794 (see also #31207).
* `f64`: faithful port of musl's `pow`, with both the fma and non-fma code paths, selected at comptime like musl's `__FP_FAST_FMA`.
* `f16` (new) and `f32`: computed as `exp2(y * log2(|x|))` in `f64` with a single rounding to the result type.
* C99 special-case layer unchanged; `f80`/`f128` remain `TODO`.
Validation, against a 128-bit MPFR reference using the methodology of Gladman/Innocente/Zimmermann, ["Accuracy of Mathematical Functions"](https://homepages.loria.fr/pzimmermann/papers/glibc238-20230921.pdf):
* over 8M bit-space samples, `f64` max error improves from 2777.9 ulp to 0.505 ulp and `f32` from 860.6 ulp to 0.500 ulp, with zero results off by more than 1 ulp (previously 62737). Targeted probing of the old algorithm shows the error passing tens of millions of ulps as the exponent gains bits; the new maximum stays within musl's analytical bound of 0.54 ulp, and the known worst case of this code base measures 0.522908 ulp, matching the published 0.523.
* bit-identical to the vendored musl `pow.c` compiled with `zig cc`, over 8M bit-space samples on each of the fma and non-fma paths.
* special-case truth table (`+-0`, `+-1`, `+-inf`, `nan`, subnormals, odd/even/non-integer/huge/tiny exponents) passes for `f16`/`f32`/`f64` and matches musl bit for bit on `f64`; `nan` results keep Zig's canonical NaN, matching existing `std.math` behavior, rather than musl's payload propagation.
* one test expectation corrected: `pow(f32, 89.123, 3.3)` is 2722489.75 when correctly rounded; the old value encoded the Go algorithm's error.