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/__rem_pio2f.c
 5
 6const std = @import("std");
 7const rem_pio2_large = @import("rem_pio2_large.zig").rem_pio2_large;
 8const math = std.math;
 9
10const toint = 1.5 / math.floatEps(f64);
11// pi/4
12const pio4 = 0x1.921fb6p-1;
13// invpio2:  53 bits of 2/pi
14const invpio2 = 6.36619772367581382433e-01; // 0x3FE45F30, 0x6DC9C883
15// pio2_1:   first 25 bits of pi/2
16const pio2_1 = 1.57079631090164184570e+00; // 0x3FF921FB, 0x50000000
17// pio2_1t:  pi/2 - pio2_1
18const pio2_1t = 1.58932547735281966916e-08; // 0x3E5110b4, 0x611A6263
19
20// Returns the remainder of x rem pi/2 in *y
21// use double precision for everything except passing x
22// use rem_pio2_large() for large x
23pub fn rem_pio2f(x: f32, y: *f64) i32 {
24    var tx: [1]f64 = undefined;
25    var ty: [1]f64 = undefined;
26    var @"fn": f64 = undefined;
27    var ix: u32 = undefined;
28    var n: i32 = undefined;
29    var sign: bool = undefined;
30    var e0: u32 = undefined;
31    var ui: u32 = undefined;
32
33    ui = @bitCast(x);
34    ix = ui & 0x7fffffff;
35
36    // 25+53 bit pi is good enough for medium size
37    if (ix < 0x4dc90fdb) { // |x| ~< 2^28*(pi/2), medium size
38        // Use a specialized rint() to get fn.
39        @"fn" = @as(f64, @floatCast(x)) * invpio2 + toint - toint;
40        n = @intFromFloat(@"fn");
41        y.* = x - @"fn" * pio2_1 - @"fn" * pio2_1t;
42        // Matters with directed rounding.
43        if (y.* < -pio4) {
44            n -= 1;
45            @"fn" -= 1;
46            y.* = x - @"fn" * pio2_1 - @"fn" * pio2_1t;
47        } else if (y.* > pio4) {
48            n += 1;
49            @"fn" += 1;
50            y.* = x - @"fn" * pio2_1 - @"fn" * pio2_1t;
51        }
52        return n;
53    }
54    if (ix >= 0x7f800000) { // x is inf or NaN
55        y.* = x - x;
56        return 0;
57    }
58    // scale x into [2^23, 2^24-1]
59    sign = ui >> 31 != 0;
60    e0 = (ix >> 23) - (0x7f + 23); // e0 = ilogb(|x|)-23, positive
61    ui = ix - (e0 << 23);
62    tx[0] = @as(f32, @bitCast(ui));
63    n = rem_pio2_large(&tx, &ty, @as(i32, @intCast(e0)), 1, 0);
64    if (sign) {
65        y.* = -ty[0];
66        return -n;
67    }
68    y.* = ty[0];
69    return n;
70}