1// Ported from musl, which is licensed under the MIT license:
  2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
  3//
  4// https://git.musl-libc.org/cgit/musl/tree/src/math/cbrtf.c
  5// https://git.musl-libc.org/cgit/musl/tree/src/math/cbrt.c
  6
  7const std = @import("../std.zig");
  8const math = std.math;
  9const expect = std.testing.expect;
 10
 11/// Returns the cube root of x.
 12///
 13/// Special Cases:
 14///  - cbrt(+-0)   = +-0
 15///  - cbrt(+-inf) = +-inf
 16///  - cbrt(nan)   = nan
 17pub fn cbrt(x: anytype) @TypeOf(x) {
 18    const T = @TypeOf(x);
 19    return switch (T) {
 20        f32 => cbrt32(x),
 21        f64 => cbrt64(x),
 22        else => @compileError("cbrt not implemented for " ++ @typeName(T)),
 23    };
 24}
 25
 26fn cbrt32(x: f32) f32 {
 27    const B1: u32 = 709958130; // (127 - 127.0 / 3 - 0.03306235651) * 2^23
 28    const B2: u32 = 642849266; // (127 - 127.0 / 3 - 24 / 3 - 0.03306235651) * 2^23
 29
 30    var u = @as(u32, @bitCast(x));
 31    var hx = u & 0x7FFFFFFF;
 32
 33    // cbrt(nan, inf) = itself
 34    if (hx >= 0x7F800000) {
 35        return x + x;
 36    }
 37
 38    // cbrt to ~5bits
 39    if (hx < 0x00800000) {
 40        // cbrt(+-0) = itself
 41        if (hx == 0) {
 42            return x;
 43        }
 44        u = @as(u32, @bitCast(x * 0x1.0p24));
 45        hx = u & 0x7FFFFFFF;
 46        hx = hx / 3 + B2;
 47    } else {
 48        hx = hx / 3 + B1;
 49    }
 50
 51    u &= 0x80000000;
 52    u |= hx;
 53
 54    // first step newton to 16 bits
 55    var t: f64 = @as(f32, @bitCast(u));
 56    var r: f64 = t * t * t;
 57    t = t * (@as(f64, x) + x + r) / (x + r + r);
 58
 59    // second step newton to 47 bits
 60    r = t * t * t;
 61    t = t * (@as(f64, x) + x + r) / (x + r + r);
 62
 63    return @as(f32, @floatCast(t));
 64}
 65
 66fn cbrt64(x: f64) f64 {
 67    const B1: u32 = 715094163; // (1023 - 1023 / 3 - 0.03306235651 * 2^20
 68    const B2: u32 = 696219795; // (1023 - 1023 / 3 - 54 / 3 - 0.03306235651 * 2^20
 69
 70    // |1 / cbrt(x) - p(x)| < 2^(23.5)
 71    const P0: f64 = 1.87595182427177009643;
 72    const P1: f64 = -1.88497979543377169875;
 73    const P2: f64 = 1.621429720105354466140;
 74    const P3: f64 = -0.758397934778766047437;
 75    const P4: f64 = 0.145996192886612446982;
 76
 77    var u = @as(u64, @bitCast(x));
 78    var hx = @as(u32, @intCast(u >> 32)) & 0x7FFFFFFF;
 79
 80    // cbrt(nan, inf) = itself
 81    if (hx >= 0x7FF00000) {
 82        return x + x;
 83    }
 84
 85    // cbrt to ~5bits
 86    if (hx < 0x00100000) {
 87        u = @as(u64, @bitCast(x * 0x1.0p54));
 88        hx = @as(u32, @intCast(u >> 32)) & 0x7FFFFFFF;
 89
 90        // cbrt(+-0) = itself
 91        if (hx == 0) {
 92            return x;
 93        }
 94        hx = hx / 3 + B2;
 95    } else {
 96        hx = hx / 3 + B1;
 97    }
 98
 99    u &= 1 << 63;
100    u |= @as(u64, hx) << 32;
101    var t = @as(f64, @bitCast(u));
102
103    // cbrt to 23 bits
104    // cbrt(x) = t * cbrt(x / t^3) ~= t * P(t^3 / x)
105    const r = (t * t) * (t / x);
106    t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4));
107
108    // Round t away from 0 to 23 bits
109    u = @as(u64, @bitCast(t));
110    u = (u + 0x80000000) & 0xFFFFFFFFC0000000;
111    t = @as(f64, @bitCast(u));
112
113    // one step newton to 53 bits
114    const s = t * t;
115    var q = x / s;
116    const w = t + t;
117    q = (q - t) / (w + q);
118
119    return t + t * q;
120}
121
122test cbrt {
123    try expect(cbrt(@as(f32, 0.0)) == cbrt32(0.0));
124    try expect(cbrt(@as(f64, 0.0)) == cbrt64(0.0));
125}
126
127test cbrt32 {
128    const epsilon = 0.000001;
129
130    try expect(math.isPositiveZero(cbrt32(0.0)));
131    try expect(math.approxEqAbs(f32, cbrt32(0.2), 0.584804, epsilon));
132    try expect(math.approxEqAbs(f32, cbrt32(0.8923), 0.962728, epsilon));
133    try expect(math.approxEqAbs(f32, cbrt32(1.5), 1.144714, epsilon));
134    try expect(math.approxEqAbs(f32, cbrt32(37.45), 3.345676, epsilon));
135    try expect(math.approxEqAbs(f32, cbrt32(123123.234375), 49.748501, epsilon));
136}
137
138test cbrt64 {
139    const epsilon = 0.000001;
140
141    try expect(math.isPositiveZero(cbrt64(0.0)));
142    try expect(math.approxEqAbs(f64, cbrt64(0.2), 0.584804, epsilon));
143    try expect(math.approxEqAbs(f64, cbrt64(0.8923), 0.962728, epsilon));
144    try expect(math.approxEqAbs(f64, cbrt64(1.5), 1.144714, epsilon));
145    try expect(math.approxEqAbs(f64, cbrt64(37.45), 3.345676, epsilon));
146    try expect(math.approxEqAbs(f64, cbrt64(123123.234375), 49.748501, epsilon));
147}
148
149test "cbrt.special" {
150    try expect(math.isPositiveZero(cbrt32(0.0)));
151    try expect(@as(u32, @bitCast(cbrt32(-0.0))) == @as(u32, 0x80000000));
152    try expect(math.isPositiveInf(cbrt32(math.inf(f32))));
153    try expect(math.isNegativeInf(cbrt32(-math.inf(f32))));
154    try expect(math.isNan(cbrt32(math.nan(f32))));
155}
156
157test "cbrt64.special" {
158    try expect(math.isPositiveZero(cbrt64(0.0)));
159    try expect(math.isNegativeZero(cbrt64(-0.0)));
160    try expect(math.isPositiveInf(cbrt64(math.inf(f64))));
161    try expect(math.isNegativeInf(cbrt64(-math.inf(f64))));
162    try expect(math.isNan(cbrt64(math.nan(f64))));
163}