rem_pio2

pub fn rem_pio2(x: f64, y: *[2]f64) i32

File

Code

pub fn rem_pio2(x: f64, y: *[2]f64) i32 {
    var z: f64 = undefined;
    var tx: [3]f64 = undefined;
    var ty: [2]f64 = undefined;
    var n: i32 = undefined;
    var ix: u32 = undefined;
    var sign: bool = undefined;
    var i: i32 = undefined;
    var ui: u64 = undefined;

    ui = @bitCast(x);
    sign = ui >> 63 != 0;
    ix = @truncate((ui >> 32) & 0x7fffffff);
    if (ix <= 0x400f6a7a) { // |x| ~<= 5pi/4
        if ((ix & 0xfffff) == 0x921fb) { // |x| ~= pi/2 or 2pi/2
            return medium(ix, x, y);
        }
        if (ix <= 0x4002d97c) { // |x| ~<= 3pi/4
            if (!sign) {
                z = x - pio2_1; // one round good to 85 bits
                y[0] = z - pio2_1t;
                y[1] = (z - y[0]) - pio2_1t;
                return 1;
            } else {
                z = x + pio2_1;
                y[0] = z + pio2_1t;
                y[1] = (z - y[0]) + pio2_1t;
                return -1;
            }
        } else {
            if (!sign) {
                z = x - 2 * pio2_1;
                y[0] = z - 2 * pio2_1t;
                y[1] = (z - y[0]) - 2 * pio2_1t;
                return 2;
            } else {
                z = x + 2 * pio2_1;
                y[0] = z + 2 * pio2_1t;
                y[1] = (z - y[0]) + 2 * pio2_1t;
                return -2;
            }
        }
    }
    if (ix <= 0x401c463b) { // |x| ~<= 9pi/4
        if (ix <= 0x4015fdbc) { // |x| ~<= 7pi/4
            if (ix == 0x4012d97c) { // |x| ~= 3pi/2
                return medium(ix, x, y);
            }
            if (!sign) {
                z = x - 3 * pio2_1;
                y[0] = z - 3 * pio2_1t;
                y[1] = (z - y[0]) - 3 * pio2_1t;
                return 3;
            } else {
                z = x + 3 * pio2_1;
                y[0] = z + 3 * pio2_1t;
                y[1] = (z - y[0]) + 3 * pio2_1t;
                return -3;
            }
        } else {
            if (ix == 0x401921fb) { // |x| ~= 4pi/2 */
                return medium(ix, x, y);
            }
            if (!sign) {
                z = x - 4 * pio2_1;
                y[0] = z - 4 * pio2_1t;
                y[1] = (z - y[0]) - 4 * pio2_1t;
                return 4;
            } else {
                z = x + 4 * pio2_1;
                y[0] = z + 4 * pio2_1t;
                y[1] = (z - y[0]) + 4 * pio2_1t;
                return -4;
            }
        }
    }
    if (ix < 0x413921fb) { // |x| ~< 2^20*(pi/2), medium size
        return medium(ix, x, y);
    }
    // all other (large) arguments
    if (ix >= 0x7ff00000) { // x is inf or NaN
        y[0] = x - x;
        y[1] = y[0];
        return 0;
    }
    // set z = scalbn(|x|,-ilogb(x)+23)
    ui = @bitCast(x);
    ui &= std.math.maxInt(u64) >> 12;
    ui |= @as(u64, 0x3ff + 23) << 52;
    z = @bitCast(ui);

    i = 0;
    while (i < 2) : (i += 1) {
        tx[@intCast(i)] = @floatFromInt(@as(i32, @intFromFloat(z)));
        z = (z - tx[@intCast(i)]) * 0x1p24;
    }
    tx[@intCast(i)] = z;
    // skip zero terms, first term is non-zero
    while (tx[@intCast(i)] == 0.0) {
        i -= 1;
    }
    n = rem_pio2_large(tx[0..], ty[0..], @as(i32, @intCast((ix >> 20))) - (0x3ff + 23), i + 1, 1);
    if (sign) {
        y[0] = -ty[0];
        y[1] = -ty[1];
        return -n;
    }
    y[0] = ty[0];
    y[1] = ty[1];
    return n;
}