util.impala

/*
 * iterators
 */

static VECTOR_LENGTH = 8;

fn @each(a: i32, b: i32, body: fn(i32) -> ()) -> () {
    for i in range_step(a, b, VECTOR_LENGTH) {
        for lane in vectorize(VECTOR_LENGTH) {
            @@body(i + lane);
        }
    }
}


/*
 * misc
 */

fn @lerp(a: f32, b: f32, k: f32) -> f32 {
    (1.0f - k) * a + k * b
}

fn @lerp2(a: f32, b: f32, c: f32, k1: f32, k2: f32) -> f32 {
    (1.0f - k1 - k2) * a + k1 * b + k2 * c
}

fn floatbits(bits: u32) -> f32 { bitcast(bits) }
fn intbits(f: f32) -> u32 { bitcast(f) }

fn @cif(c: bool, t: fn () -> ()) -> () {
    if rv_any(c) {
        if c {
            @@t()
        }
    }
}


/*
 * Vec2, Vec3, Vec4
 */

struct Vec2 {
    x: f32,
    y: f32
}

struct Vec3 {
    x: f32,
    y: f32,
    z: f32
}

struct Vec4 {
    x: f32,
    y: f32,
    z: f32,
    w: f32
}

fn @make_vec2(x: f32, y: f32) -> Vec2 { Vec2 { x: x, y: y } }
fn @make_vec3(x: f32, y: f32, z: f32) -> Vec3 { Vec3 { x: x, y: y, z: z } }
fn @make_vec4(x: f32, y: f32, z: f32, w: f32) -> Vec4 { Vec4 { x: x, y: y, z: z, w: w } }

fn @vec2_to_3(v: Vec2, z: f32) -> Vec3 { make_vec3(v.x, v.y, z) }
fn @vec3_to_2(v: Vec3) -> Vec2 { make_vec2(v.x, v.y) }
fn @vec3_to_4(v: Vec3, w: f32) -> Vec4 { make_vec4(v.x, v.y, v.z, w) }
fn @vec4_to_3(v: Vec4) -> Vec3 { make_vec3(v.x, v.y, v.z) }

fn @vec2_map(v: Vec2, f: fn (f32) -> f32) -> Vec2 { make_vec2(@@f(v.x), @@f(v.y)) }
fn @vec3_map(v: Vec3, f: fn (f32) -> f32) -> Vec3 { make_vec3(@@f(v.x), @@f(v.y), @@f(v.z)) }
fn @vec4_map(v: Vec4, f: fn (f32) -> f32) -> Vec4 { make_vec4(@@f(v.x), @@f(v.y), @@f(v.z), @@f(v.w)) }

fn @vec2_zip(a: Vec2, b: Vec2, f: fn (f32, f32) -> f32) -> Vec2 { make_vec2(@@f(a.x, b.x), @@f(a.y, b.y)) }
fn @vec3_zip(a: Vec3, b: Vec3, f: fn (f32, f32) -> f32) -> Vec3 { make_vec3(@@f(a.x, b.x), @@f(a.y, b.y), @@f(a.z, b.z)) }
fn @vec4_zip(a: Vec4, b: Vec4, f: fn (f32, f32) -> f32) -> Vec4 { make_vec4(@@f(a.x, b.x), @@f(a.y, b.y), @@f(a.z, b.z), @@f(a.w, b.w)) }

fn @vec2_add(a: Vec2, b: Vec2) -> Vec2 { vec2_zip(a, b, |x, y| x + y) }
fn @vec3_add(a: Vec3, b: Vec3) -> Vec3 { vec3_zip(a, b, |x, y| x + y) }
fn @vec4_add(a: Vec4, b: Vec4) -> Vec4 { vec4_zip(a, b, |x, y| x + y) }
fn @vec2_sub(a: Vec2, b: Vec2) -> Vec2 { vec2_zip(a, b, |x, y| x - y) }
fn @vec3_sub(a: Vec3, b: Vec3) -> Vec3 { vec3_zip(a, b, |x, y| x - y) }
fn @vec4_sub(a: Vec4, b: Vec4) -> Vec4 { vec4_zip(a, b, |x, y| x - y) }
fn @vec2_mul(a: Vec2, b: Vec2) -> Vec2 { vec2_zip(a, b, |x, y| x * y) }
fn @vec3_mul(a: Vec3, b: Vec3) -> Vec3 { vec3_zip(a, b, |x, y| x * y) }
fn @vec4_mul(a: Vec4, b: Vec4) -> Vec4 { vec4_zip(a, b, |x, y| x * y) }
fn @vec2_div(a: Vec2, b: Vec2) -> Vec2 { vec2_zip(a, b, |x, y| x / y) }
fn @vec3_div(a: Vec3, b: Vec3) -> Vec3 { vec3_zip(a, b, |x, y| x / y) }
fn @vec4_div(a: Vec4, b: Vec4) -> Vec4 { vec4_zip(a, b, |x, y| x / y) }

fn @vec2_neg(v: Vec2) -> Vec2 { vec2_map(v, |x| -x) }
fn @vec3_neg(v: Vec3) -> Vec3 { vec3_map(v, |x| -x) }
fn @vec4_neg(v: Vec4) -> Vec4 { vec4_map(v, |x| -x) }

fn @vec2_mulf(v: Vec2, t: f32) -> Vec2 { vec2_mul(v, make_vec2(t, t)) }
fn @vec3_mulf(v: Vec3, t: f32) -> Vec3 { vec3_mul(v, make_vec3(t, t, t)) }
fn @vec4_mulf(v: Vec4, t: f32) -> Vec4 { vec4_mul(v, make_vec4(t, t, t, t)) }

fn @vec2_dot(a: Vec2, b: Vec2) -> f32 { a.x * b.x + a.y * b.y }
fn @vec3_dot(a: Vec3, b: Vec3) -> f32 { a.x * b.x + a.y * b.y + a.z * b.z }
fn @vec4_dot(a: Vec4, b: Vec4) -> f32 { a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w }

fn @vec3_cross(a: Vec3, b: Vec3) -> Vec3 {
    make_vec3(a.y * b.z - a.z * b.y,
              a.z * b.x - a.x * b.z,
              a.x * b.y - a.y * b.x)
}

fn @vec2_len2(v: Vec2) -> f32 { vec2_dot(v, v) }
fn @vec3_len2(v: Vec3) -> f32 { vec3_dot(v, v) }
fn @vec4_len2(v: Vec4) -> f32 { vec4_dot(v, v) }

fn @vec2_reflect(v: Vec2, n: Vec2) -> Vec2 { vec2_sub(vec2_mulf(n, 2.0f * vec2_dot(n, v)), v) }
fn @vec3_reflect(v: Vec3, n: Vec3) -> Vec3 { vec3_sub(vec3_mulf(n, 2.0f * vec3_dot(n, v)), v) }
fn @vec4_reflect(v: Vec4, n: Vec4) -> Vec4 { vec4_sub(vec4_mulf(n, 2.0f * vec4_dot(n, v)), v) }

fn @vec2_len(math: Intrinsics, v: Vec2) -> f32 { math.sqrtf(vec2_len2(v)) }
fn @vec3_len(math: Intrinsics, v: Vec3) -> f32 { math.sqrtf(vec3_len2(v)) }
fn @vec4_len(math: Intrinsics, v: Vec4) -> f32 { math.sqrtf(vec4_len2(v)) }

fn @vec2_normalize(math: Intrinsics, v: Vec2) -> Vec2 { vec2_mulf(v, 1.0f / vec2_len(math, v)) }
fn @vec3_normalize(math: Intrinsics, v: Vec3) -> Vec3 { vec3_mulf(v, 1.0f / vec3_len(math, v)) }
fn @vec4_normalize(math: Intrinsics, v: Vec4) -> Vec4 { vec4_mulf(v, 1.0f / vec4_len(math, v)) }

fn @vec2_lerp(a: Vec2, b: Vec2, k: f32) -> Vec2 { vec2_zip(a, b, |x, y| lerp(x, y, k)) }
fn @vec3_lerp(a: Vec3, b: Vec3, k: f32) -> Vec3 { vec3_zip(a, b, |x, y| lerp(x, y, k)) }
fn @vec4_lerp(a: Vec4, b: Vec4, k: f32) -> Vec4 { vec4_zip(a, b, |x, y| lerp(x, y, k)) }

fn @vec2_lerp2(a: Vec2, b: Vec2, c: Vec2, u: f32, v: f32) -> Vec2 {
    Vec2 {
        x: lerp2(a.x, b.x, c.x, u, v),
        y: lerp2(a.y, b.y, c.y, u, v)
    }
}
fn @vec3_lerp2(a: Vec3, b: Vec3, c: Vec3, u: f32, v: f32) -> Vec3 {
    Vec3 {
        x: lerp2(a.x, b.x, c.x, u, v),
        y: lerp2(a.y, b.y, c.y, u, v),
        z: lerp2(a.z, b.z, c.z, u, v)
    }
}
fn @vec4_lerp2(a: Vec4, b: Vec4, c: Vec4, u: f32, v: f32) -> Vec4 {
    Vec4 {
        x: lerp2(a.x, b.x, c.x, u, v),
        y: lerp2(a.y, b.y, c.y, u, v),
        z: lerp2(a.z, b.z, c.z, u, v),
        w: lerp2(a.w, b.w, c.w, u, v)
    }
}

/*
 * RNGState
 */

struct RNGState {
    z1: u32,
    z2: u32,
    z3: u32,
    z4: u32,
}

fn random(state: &mut RNGState) -> u32 {
    let mut b  = ((state.z1 << 6u) ^ state.z1) >> 13u;
    state.z1 = ((state.z1 & 4294967294u) << 18u) ^ b;
    b  = ((state.z2 << 2u) ^ state.z2) >> 27u;
    state.z2 = ((state.z2 & 4294967288u) << 2u) ^ b;
    b  = ((state.z3 << 13u) ^ state.z3) >> 21u;
    state.z3 = ((state.z3 & 4294967280u) << 7u) ^ b;
    b  = ((state.z4 << 3u) ^ state.z4) >> 12u;
    state.z4 = ((state.z4 & 4294967168u) << 13u) ^ b;
    state.z1 ^ state.z2 ^ state.z3 ^ state.z4
}

fn frandom(state: &mut RNGState) -> f32 {
    let mut irand = random(state);
    irand &= (1u<<23u)-1u;
    floatbits(0x3F800000u | irand)-1.0f
}

fn seed_rng(state: &mut RNGState, seed: u32) -> () {
    state.z1 = seed;
    state.z2 = seed ^ 0xbeeff00du;
    state.z3 = ((seed & 0xffffu) << 16u) | (seed >> 16u);
    state.z4 = (((seed & 0xffu) << 24u) | ((seed & 0xff00u)  << 8u) |
                 ((seed & 0xff0000u) >> 8u) | (seed & 0xff000000u) >> 24u);
}

/*
 * math: misc
 */

fn @floor(f: f32) -> f32 { cpu_intrinsics.floorf(f) }
fn @clamp(i: i32, a: i32, b: i32) -> i32 { math.min(math.max(i, a), b) }
fn @clampf(i: f32, a: f32, b: f32) -> f32 { math.fminf(math.fmaxf(i, a), b) }

/*
 * math: exp, pow, log
 */

fn ldexp(x: f32, mut n: u32) -> f32 {
    let mut ex = 0x7F800000u;
    let mut ix = intbits(x);
    ex &= ix;              // extract old exponent;
    ix = ix & !0x7F800000u;  // clear exponent
    n = (n << 23u) + ex;
    ix |= n; // insert new exponent
    floatbits(ix)
}

fn frexp(x: f32) -> (f32, u32) {
    let mut ex = 0x7F800000u;              // exponent mask
    let mut ix = intbits(x);
    ex &= ix;
    ix &= !0x7F800000u;  // clear exponent
    let pw2 = (ex >> 23u) - 126u; // compute exponent
    ix |= 0x3F000000u;         // insert exponent +1 in x
    (floatbits(ix), pw2)
}

fn exp(mut x_full: f32) -> f32 {
    let fast = true;
    if fast {
        let mut z = floor(1.44269504088896341f * x_full + 0.5f);
        x_full -= z * 0.693359375f;
        x_full -= z * -2.12194440e-4f;
        let n = z as int;

        z = x_full * x_full;
        z = (((((1.9875691500E-4f  * x_full + 1.3981999507E-3f) * x_full +
                8.3334519073E-3f) * x_full + 4.1665795894E-2f) * x_full +
                1.6666665459E-1f) * x_full + 5.0000001201E-1f) * z + x_full + 1.f;
        x_full = ldexp(z, n as u32);
        x_full
    } else {
        let ln2_part1 = 0.6931457519f;
        let ln2_part2 = 1.4286067653e-6f;
        let one_over_ln2 = 1.44269502162933349609375f;

        let scaled = x_full * one_over_ln2;
        let k_real = floor(scaled);
        let k = k_real as i32;

        // Reduced range version of x
        let mut x = x_full - k_real * ln2_part1;
        x -= k_real * ln2_part2;

        // These coefficients are for e^x in [0, ln(2)]
        let one = 1.f;
        let c2 = 0.4999999105930328369140625f;
        let c3 = 0.166668415069580078125f;
        let c4 = 4.16539050638675689697265625e-2f;
        let c5 = 8.378830738365650177001953125e-3f;
        let c6 = 1.304379315115511417388916015625e-3f;
        let c7 = 2.7555381529964506626129150390625e-4f;

        let mut result = x * c7 + c6;
        result = x * result + c5;
        result = x * result + c4;
        result = x * result + c3;
        result = x * result + c2;
        result = x * result + one;
        result = x * result + one;

        // Compute 2^k (should differ for uniform float and double, but I'll avoid
        // it for now and just do uniform floats)
        let fpbias = 127;
        let mut biased_n = k + fpbias;
        let overflow = k > fpbias;
        // Minimum exponent is -126, so if k is <= -127 (k + 127 <= 0)
        // we've got underflow. -127 * ln(2) -> -88.02. So the most
        // negative uniform float input that doesn't result in zero is like -88.
        let underflow = (biased_n <= 0);
        let InfBits = 0x7f800000;
        biased_n <<= 23;
        // Reuniform interpret this thing as uniform float
        let two_to_the_n = floatbits(biased_n as u32);
        // Handle both doubles and uniform floats (hopefully eliding the copy for uniform float)
        let elemtype_2n = two_to_the_n;
        result *= elemtype_2n;
        result = select(overflow, floatbits(InfBits as u32), result);
        result = select(underflow, 0.f, result);
        result
    }
}

fn @log(x: f32) -> f32 {
    fn @range_reduce(input: f32) -> (f32, u32) {
        let int_version = intbits(input);
        let nonexponent_mask = 0x807FFFFFu;

        let exponent_neg1 = (126u << 23u);
        let biased_exponent = int_version >> 23u;
        let offset_exponent = biased_exponent + 1u;
        let exponent = offset_exponent - 127u; // get the real value
        let blended = (int_version & nonexponent_mask) | (exponent_neg1);
        let reduced = floatbits(blended);
        (reduced, exponent)
    }

    let fast = true;
    if fast {
        let (x_, e_) = frexp(x);
        let mut x_full = x_;
        let mut e = e_;

        let x_smaller_SQRTHF = select(0.707106781186547524f > x_full, 0xffffffffu, 0u);
        e += x_smaller_SQRTHF;
        let mut ix_add = intbits(x_full);
        ix_add &= x_smaller_SQRTHF;
        x_full += floatbits(ix_add) - 1.f;

        let mut z = x_full * x_full;
        let mut y =
            ((((((((7.0376836292E-2f * x_full
                    + -1.1514610310E-1f) * x_full
                   + 1.1676998740E-1f) * x_full
                  + -1.2420140846E-1f) * x_full
                 + 1.4249322787E-1f) * x_full
                + -1.6668057665E-1f) * x_full
               + 2.0000714765E-1f) * x_full
              + -2.4999993993E-1f) * x_full
             + 3.3333331174E-1f) * x_full * z;

        let fe = e as f32;
        y += fe * -2.12194440e-4f;
        y -= 0.5f * z;
        z  = x_full + y;
        z + 0.693359375f * fe
    } else {
        let NaN_bits = 0x7fc00000u;
        let Neg_Inf_bits = 0xFF800000u;
        let NaN = floatbits(NaN_bits);
        let neg_inf = floatbits(Neg_Inf_bits);
        let use_nan = x < 0.f;
        let use_inf = x == 0.f;
        let exceptional = use_nan || use_inf;
        let one = 1.0f;

        let patched = select(exceptional, one, x);
        let (reduced, exponent) = range_reduce(patched);

        let ln2 = 0.693147182464599609375f;

        let x1 = one - reduced;
        let c1 = 0.50000095367431640625f;
        let c2 = 0.33326041698455810546875f;
        let c3 = 0.2519190013408660888671875f;
        let c4 = 0.17541764676570892333984375f;
        let c5 = 0.3424419462680816650390625f;
        let c6 = -0.599632322788238525390625f;
        let c7 = +1.98442304134368896484375f;
        let c8 = -2.4899270534515380859375f;
        let c9 = +1.7491014003753662109375f;

        let mut result = x1 * c9 + c8;
        result = x1 * result + c7;
        result = x1 * result + c6;
        result = x1 * result + c5;
        result = x1 * result + c4;
        result = x1 * result + c3;
        result = x1 * result + c2;
        result = x1 * result + c1;
        result = x1 * result + one;

        // Equation was for -(ln(red)/(1-red))
        result *= -x1;
        result += (exponent as f32) * ln2;

        select(exceptional, select(use_nan, NaN, neg_inf), result)
    }
}

/*
 * math: trigonometry
 */

fn @sin(x_full: f32) -> f32 {
    static pi_over_two_vec = 1.57079637050628662109375f;
    static two_over_pi_vec = 0.636619746685028076171875f;
    let scaled = x_full * two_over_pi_vec;
    let k_real = floor(scaled);
    let k = k_real as i32;

    // Reduced range version of x
    let x = x_full - k_real * pi_over_two_vec;
    let k_mod4 = k & 3;
    let sin_usecos = (k_mod4 == 1) | (k_mod4 == 3);
    let flip_sign = (k_mod4 > 1);

    // These coefficients are from sollya with fpminimax(sin(x)/x, [|0, 2,
    // 4, 6, 8, 10|], [|single...|], [0;Pi/2]);
    static sin_c2 = -0.16666667163372039794921875f;
    static sin_c4 = 8.333347737789154052734375e-3f;
    static sin_c6 = -1.9842604524455964565277099609375e-4f;
    static sin_c8 = 2.760012648650445044040679931640625e-6f;
    static sin_c10 = -2.50293279435709337121807038784027099609375e-8f;

    static cos_c2 = -0.5f;
    static cos_c4 = 4.166664183139801025390625e-2f;
    static cos_c6 = -1.388833043165504932403564453125e-3f;
    static cos_c8 = 2.47562347794882953166961669921875e-5f;
    static cos_c10 = -2.59630184018533327616751194000244140625e-7f;

    let outside = select(sin_usecos, 1.f    , x);
    let c2      = select(sin_usecos, cos_c2 , sin_c2);
    let c4      = select(sin_usecos, cos_c4 , sin_c4);
    let c6      = select(sin_usecos, cos_c6 , sin_c6);
    let c8      = select(sin_usecos, cos_c8 , sin_c8);
    let c10     = select(sin_usecos, cos_c10, sin_c10);

    let x2 = x * x;
    let mut formula = x2 * c10 + c8;
    formula = x2 * formula + c6;
    formula = x2 * formula + c4;
    formula = x2 * formula + c2;
    formula = x2 * formula + 1.f;
    formula *= outside;

    select(flip_sign, -formula, formula)
}

fn @cos(x_full: f32) -> f32 {
    static pi_over_two_vec = 1.57079637050628662109375f;
    static two_over_pi_vec = 0.636619746685028076171875f;
    let scaled = x_full * two_over_pi_vec;
    let k_real = floor(scaled);
    let k = k_real as i32;

    // Reduced range version of x
    let x = x_full - k_real * pi_over_two_vec;

    let k_mod4 = k & 3;
    let cos_usecos = (k_mod4 == 0) | (k_mod4 == 2);
    let flip_sign = (k_mod4 == 1) | (k_mod4 == 2);

    let sin_c2 = -0.16666667163372039794921875f;
    let sin_c4 = 8.333347737789154052734375e-3f;
    let sin_c6 = -1.9842604524455964565277099609375e-4f;
    let sin_c8 = 2.760012648650445044040679931640625e-6f;
    let sin_c10 = -2.50293279435709337121807038784027099609375e-8f;

    let cos_c2 = -0.5f;
    let cos_c4 = 4.166664183139801025390625e-2f;
    let cos_c6 = -1.388833043165504932403564453125e-3f;
    let cos_c8 = 2.47562347794882953166961669921875e-5f;
    let cos_c10 = -2.59630184018533327616751194000244140625e-7f;

    let outside = select(cos_usecos, 1.0f, x);
    let c2      = select(cos_usecos, cos_c2,  sin_c2 );
    let c4      = select(cos_usecos, cos_c4,  sin_c4 );
    let c6      = select(cos_usecos, cos_c6,  sin_c6 );
    let c8      = select(cos_usecos, cos_c8,  sin_c8 );
    let c10     = select(cos_usecos, cos_c10, sin_c10);

    let x2 = x * x;
    let mut formula = x2 * c10 + c8;
    formula = x2 * formula + c6;
    formula = x2 * formula + c4;
    formula = x2 * formula + c2;
    formula = x2 * formula + 1.f;
    formula *= outside;

    select(flip_sign, -formula, formula)
}

fn @sincos(x_full: f32) -> (f32, f32) {
    let pi_over_two_vec = 1.57079637050628662109375f;
    let two_over_pi_vec = 0.636619746685028076171875f;
    let scaled = x_full * two_over_pi_vec;
    let k_real = floor(scaled);
    let k = k_real as i32;

    // Reduced range version of x
    let x = x_full - k_real * pi_over_two_vec;
    let k_mod4 = k & 3;
    let cos_usecos = (k_mod4 == 0) | (k_mod4 == 2);
    let sin_usecos = (k_mod4 == 1) | (k_mod4 == 3);
    let sin_flipsign = (k_mod4 > 1);
    let cos_flipsign = (k_mod4 == 1) | (k_mod4 == 2);

    let one_vec = 1.f;
    let sin_c2 = -0.16666667163372039794921875f;
    let sin_c4 = 8.333347737789154052734375e-3f;
    let sin_c6 = -1.9842604524455964565277099609375e-4f;
    let sin_c8 = 2.760012648650445044040679931640625e-6f;
    let sin_c10 = -2.50293279435709337121807038784027099609375e-8f;

    let cos_c2 = -0.5f;
    let cos_c4 = 4.166664183139801025390625e-2f;
    let cos_c6 = -1.388833043165504932403564453125e-3f;
    let cos_c8 = 2.47562347794882953166961669921875e-5f;
    let cos_c10 = -2.59630184018533327616751194000244140625e-7f;

    let x2 = x * x;

    let mut sin_formula = x2 * sin_c10 + sin_c8;
    let mut cos_formula = x2 * cos_c10 + cos_c8;
    sin_formula = x2 * sin_formula + sin_c6;
    cos_formula = x2 * cos_formula + cos_c6;

    sin_formula = x2 * sin_formula + sin_c4;
    cos_formula = x2 * cos_formula + cos_c4;

    sin_formula = x2 * sin_formula + sin_c2;
    cos_formula = x2 * cos_formula + cos_c2;

    sin_formula = x2 * sin_formula + one_vec;
    cos_formula = x2 * cos_formula + one_vec;

    sin_formula *= x;

    let mut sin_result = select(sin_usecos, cos_formula, sin_formula);
    let mut cos_result = select(cos_usecos, cos_formula, sin_formula);

    sin_result = select(sin_flipsign, -sin_result, sin_result);
    cos_result = select(cos_flipsign, -cos_result, cos_result);

    (sin_result, cos_result)
}

fn tan(x_full: f32) -> f32 {
    let pi_over_four_vec = 0.785398185253143310546875f;
    let four_over_pi_vec = 1.27323949337005615234375f;

    let x_lt_0 = x_full < 0.f;
    let y = if x_lt_0 { -x_full } else { x_full };
    let scaled = y * four_over_pi_vec;

    let k_real = floor(scaled);
    let k = k_real as i32;

    let mut x = y - k_real * pi_over_four_vec;

    // if k & 1, x -= Pi/4
    let need_offset = (k & 1) != 0;
    x = if need_offset { x - pi_over_four_vec } else { x };

    // if k & 3 == (0 or 3) let z = tan_In...(y) otherwise z = -cot_In0To...
    let k_mod4 = k & 3;
    let use_cotan = (k_mod4 == 1) || (k_mod4 == 2);

    let  one_vec = 1.0f;

    let  tan_c2 = 0.33333075046539306640625f;
    let  tan_c4 = 0.13339905440807342529296875f;
    let  tan_c6 = 5.3348250687122344970703125e-2f;
    let  tan_c8 = 2.46033705770969390869140625e-2f;
    let  tan_c10 = 2.892402000725269317626953125e-3f;
    let  tan_c12 = 9.500005282461643218994140625e-3f;

    let  cot_c2 = -0.3333333432674407958984375f;
    let  cot_c4 = -2.222204394638538360595703125e-2f;
    let  cot_c6 = -2.11752182804048061370849609375e-3f;
    let  cot_c8 = -2.0846328698098659515380859375e-4f;
    let  cot_c10 = -2.548247357481159269809722900390625e-5f;
    let  cot_c12 = -3.5257363606433500535786151885986328125e-7f;

    let x2 = x * x;

    let z =
        if (use_cotan) {
            let mut cot_val = x2 * cot_c12 + cot_c10;
            cot_val = x2 * cot_val + cot_c8;
            cot_val = x2 * cot_val + cot_c6;
            cot_val = x2 * cot_val + cot_c4;
            cot_val = x2 * cot_val + cot_c2;
            cot_val = x2 * cot_val + one_vec;
            // The equation is for x * cot(x) but we need -x * cot(x) for the tan part.
            cot_val /= -x;
            cot_val
        } else {
            let mut tan_val = x2 * tan_c12 + tan_c10;
            tan_val = x2 * tan_val + tan_c8;
            tan_val = x2 * tan_val + tan_c6;
            tan_val = x2 * tan_val + tan_c4;
            tan_val = x2 * tan_val + tan_c2;
            tan_val = x2 * tan_val + one_vec;
            // Equation was for tan(x)/x
            tan_val *= x;
            tan_val
        };

    if x_lt_0 { -z } else { z }
}