zig

fork of https://codeberg.org/ziglang/zig
Log | Files | Refs | README | LICENSE

blob cd3b71ca (4211B) - Raw

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