port cbrt from musl

`dmath_cbrt` was `dmath_pow(x, 1.0 / 3.0)`, which returned nan for every
negative input and rode on the spline-table log2/exp2 path, so it was
inexact (cbrt(27) == 2.9999999999999996, 1800 of the first 2000 perfect
cubes wrong) and mishandled subnormals, because dmath_log2 does not
normalize them: cbrt(5e-324) came out 5 orders of magnitude off.

Replace it with a line-by-line port of cbrt64 from the Zig standard
library (musl's s_cbrt.c), next to the asin/acos/atan ports in
dmath_zig.c. cbrt is not required by IEEE 754 to be correctly rounded,
so unlike sqrt there is no hardware instruction to lean on; this is
integer bit manipulation plus basic arithmetic, no libm and no FMA.

Checked against correctly rounded results computed with exact rational
arithmetic over 40000 random inputs: correctly rounded 91.6% of the
time, max error 1 ulp. CPython on Windows (MSVC UCRT) manages 69.6% and
2 ulp, so a few of the new expectations in 930_deterministic_float.py
differ from math.cbrt by 1 ulp on purpose; those are annotated inline.

Also expose math.cbrt, since dmath_cbrt otherwise has no caller
reachable from Python: colorcvt has been disabled since b07db073.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

docs/modules/math.md

@@ -87,6 +87,10 @@ Return `x` raised to the power `y`.
 
 Return the square root of `x`.
 
+### `math.cbrt(x)`
+
+Return the cube root of `x`. Unlike `sqrt`, `x` may be negative.
+
 ### `math.acos(x)`
 
 Return the arc cosine of `x`, in radians.

src/common/dmath.c

@@ -218,10 +218,6 @@ double dmath_sqrt(double x) {
 #endif
 }
 
-double dmath_cbrt(double x) {
-    return dmath_pow(x, 1.0 / 3.0);
-}
-
 // https://github.com/kraj/musl/blob/kraj/master/src/math/sincos.c
 static double __sin(double x, double y, int iy)
 {

src/common/dmath_zig.c

@@ -1,5 +1,5 @@
-// Inverse trigonometric functions ported from the Zig standard library, which
-// in turn ports them from musl (MIT licensed):
+// Math functions ported from the Zig standard library, which in turn ports
+// them from musl (MIT licensed):
 // https://github.com/ziglang/zig/blob/master/lib/std/math/
 //
 // These live in their own file so the provenance stays obvious: each function
@@ -304,3 +304,63 @@ double dmath_atan2(double y, double x) {
         default: return (z - pi_lo) - pi; /* atan(-, -) */
     }
 }
+
+// https://github.com/ziglang/zig/blob/master/lib/std/math/cbrt.zig
+//
+// cbrt is not required by IEEE 754 to be correctly rounded, so there is no
+// hardware instruction to lean on like `dmath_sqrt` does; this software
+// version is what gives the same bits on every platform.
+double dmath_cbrt(double x) {
+    const uint32_t B1 = 715094163; /* (1023 - 1023 / 3 - 0.03306235651) * 2^20 */
+    const uint32_t B2 = 696219795; /* (1023 - 1023 / 3 - 54 / 3 - 0.03306235651) * 2^20 */
+
+    /* |1 / cbrt(x) - p(x)| < 2^-23.5 */
+    const double P0 = 1.87595182427177009643;
+    const double P1 = -1.88497979543377169875;
+    const double P2 = 1.621429720105354466140;
+    const double P3 = -0.758397934778766047437;
+    const double P4 = 0.145996192886612446982;
+
+    union ZigF64 ux = {.f = x};
+    uint64_t u = ux.i;
+    uint32_t hx = (uint32_t)(u >> 32) & 0x7FFFFFFF;
+
+    /* cbrt(nan, inf) = itself */
+    if(hx >= 0x7FF00000) return x + x;
+
+    /* cbrt to ~5bits */
+    if(hx < 0x00100000) {
+        union ZigF64 us = {.f = x * 0x1.0p54};
+        u = us.i;
+        hx = (uint32_t)(u >> 32) & 0x7FFFFFFF;
+
+        /* cbrt(+-0) = itself */
+        if(hx == 0) return x;
+        hx = hx / 3 + B2;
+    } else {
+        hx = hx / 3 + B1;
+    }
+
+    u &= 0x8000000000000000ULL;
+    u |= (uint64_t)hx << 32;
+    union ZigF64 ut = {.i = u};
+    double t = ut.f;
+
+    /* cbrt to 23 bits
+     * cbrt(x) = t * cbrt(x / t^3) ~= t * P(t^3 / x) */
+    double r = (t * t) * (t / x);
+    t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4));
+
+    /* Round t away from 0 to 23 bits */
+    ut.f = t;
+    ut.i = (ut.i + 0x80000000) & 0xFFFFFFFFC0000000ULL;
+    t = ut.f;
+
+    /* one step newton to 53 bits */
+    double s = t * t;
+    double q = x / s;
+    double w = t + t;
+    q = (q - t) / (w + q);
+
+    return t + t * q;
+}

src/modules/math.c

@@ -118,6 +118,7 @@ ONE_ARG_FUNC(log2, dmath_log2)
 ONE_ARG_FUNC(log10, dmath_log10)
 TWO_ARG_FUNC(pow, dmath_pow)
 ONE_ARG_FUNC(sqrt, dmath_sqrt)
+ONE_ARG_FUNC(cbrt, dmath_cbrt)
 
 ONE_ARG_FUNC(acos, dmath_acos)
 ONE_ARG_FUNC(asin, dmath_asin)
@@ -198,6 +199,7 @@ void pk__add_module_math() {
 
     py_bindfunc(mod, "pow", math_pow);
     py_bindfunc(mod, "sqrt", math_sqrt);
+    py_bindfunc(mod, "cbrt", math_cbrt);
 
     py_bindfunc(mod, "acos", math_acos);
     py_bindfunc(mod, "asin", math_asin);

tests/930_deterministic_float.py

@@ -102,6 +102,42 @@ assert math.isnan(math.sqrt(-1))
 assert math.isnan(math.sqrt(-math.inf))
 assert math.isnan(math.sqrt(math.nan))
 
+# test cbrt
+# cbrt is not required by IEEE 754 to be correctly rounded, so unlike sqrt there is
+# no hardware instruction to lean on; these come from the musl port in dmath_zig.c.
+# The cases marked below differ from CPython's math.cbrt by 1 ulp on Windows, where
+# it is the musl port that is correctly rounded and the platform CRT that is not.
+assertEqual(math.cbrt(2), 1.2599210498948732)  # CPython/msvc gives 1.259921049894873
+assertEqual(math.cbrt(3), 1.4422495703074083)
+assertEqual(math.cbrt(10), 2.154434690031884)
+assertEqual(math.cbrt(0.1), 0.4641588833612779)  # CPython/msvc gives 0.464158883361278
+assertEqual(math.cbrt(125.872509), 5.011606490362725)
+assertEqual(math.cbrt(1225.296280), 10.700737364388822)
+assertEqual(math.cbrt(1e300), 1e100)
+assertEqual(math.cbrt(1.7976931348623157e308), 5.643803094122362e102)
+# negative inputs, which the old dmath_pow based cbrt returned nan for
+assertEqual(math.cbrt(-8.0), -2.0)
+assertEqual(math.cbrt(-0.5), -0.7937005259840998)  # CPython/msvc gives -0.7937005259840997
+assertEqual(math.cbrt(-1000.0), -10.0)
+assertEqual(math.cbrt(-1e300), -1e100)
+# subnormal inputs, which the old cbrt got wrong by orders of magnitude
+# because dmath_log2 does not normalize them: cbrt(5e-324) came out as 2.2323972485981933e-103
+assertEqual(math.cbrt(5e-324), 1.7031839360032603e-108)
+assertEqual(math.cbrt(-5e-324), -1.7031839360032603e-108)
+assertEqual(math.cbrt(1e-320), 2.1544266950262728e-107)  # CPython/msvc gives 2.154426695026273e-107
+assertEqual(math.cbrt(2.2250738585072014e-308), 2.812644285236262e-103)  # CPython/msvc gives 2.8126442852362615e-103
+# perfect cubes must be exact; the old cbrt got 1800 of these 2000 wrong,
+# e.g. cbrt(27) == 2.9999999999999996 and cbrt(8) == 1.9999999999999998
+for i in range(2000):
+    assertEqual(math.cbrt(i * i * i), float(i))
+    assertEqual(math.cbrt(-(i * i * i)), -float(i))
+# special values
+assertEqual(math.cbrt(math.inf), math.inf)
+assertEqual(math.cbrt(-math.inf), -math.inf)
+assertEqual(math.copysign(1.0, math.cbrt(0.0)), 1.0)
+assertEqual(math.copysign(1.0, math.cbrt(-0.0)), -1.0)
+assert math.isnan(math.cbrt(math.nan))
+
 # test cos, sin, tan
 assertEqual(math.cos(0), 1.0)
 assertEqual(math.cos(math.pi/2), 6.123233995736766e-17)