rem_pio2f.zig (2247B) - Raw
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 6 const std = @import("std"); 7 const rem_pio2_large = @import("rem_pio2_large.zig").rem_pio2_large; 8 const math = std.math; 9 10 const toint = 1.5 / math.floatEps(f64); 11 // pi/4 12 const pio4 = 0x1.921fb6p-1; 13 // invpio2: 53 bits of 2/pi 14 const invpio2 = 6.36619772367581382433e-01; // 0x3FE45F30, 0x6DC9C883 15 // pio2_1: first 25 bits of pi/2 16 const pio2_1 = 1.57079631090164184570e+00; // 0x3FF921FB, 0x50000000 17 // pio2_1t: pi/2 - pio2_1 18 const 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 23 pub 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 }