rem_pio2.zig (6045B) - 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_pio2.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.921fb54442d18p-1; 13 // invpio2: 53 bits of 2/pi 14 const invpio2 = 6.36619772367581382433e-01; // 0x3FE45F30, 0x6DC9C883 15 // pio2_1: first 33 bit of pi/2 16 const pio2_1 = 1.57079632673412561417e+00; // 0x3FF921FB, 0x54400000 17 // pio2_1t: pi/2 - pio2_1 18 const pio2_1t = 6.07710050650619224932e-11; // 0x3DD0B461, 0x1A626331 19 // pio2_2: second 33 bit of pi/2 20 const pio2_2 = 6.07710050630396597660e-11; // 0x3DD0B461, 0x1A600000 21 // pio2_2t: pi/2 - (pio2_1+pio2_2) 22 const pio2_2t = 2.02226624879595063154e-21; // 0x3BA3198A, 0x2E037073 23 // pio2_3: third 33 bit of pi/2 24 const pio2_3 = 2.02226624871116645580e-21; // 0x3BA3198A, 0x2E000000 25 // pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) 26 const pio2_3t = 8.47842766036889956997e-32; // 0x397B839A, 0x252049C1 27 28 fn medium(ix: u32, x: f64, y: *[2]f64) i32 { 29 var w: f64 = undefined; 30 var t: f64 = undefined; 31 var r: f64 = undefined; 32 var @"fn": f64 = undefined; 33 var n: i32 = undefined; 34 var ex: i32 = undefined; 35 var ey: i32 = undefined; 36 var ui: u64 = undefined; 37 38 // rint(x/(pi/2)) 39 @"fn" = x * invpio2 + toint - toint; 40 n = @intFromFloat(@"fn"); 41 r = x - @"fn" * pio2_1; 42 w = @"fn" * pio2_1t; // 1st round, good to 85 bits 43 // Matters with directed rounding. 44 if (r - w < -pio4) { 45 n -= 1; 46 @"fn" -= 1; 47 r = x - @"fn" * pio2_1; 48 w = @"fn" * pio2_1t; 49 } else if (r - w > pio4) { 50 n += 1; 51 @"fn" += 1; 52 r = x - @"fn" * pio2_1; 53 w = @"fn" * pio2_1t; 54 } 55 y[0] = r - w; 56 ui = @bitCast(y[0]); 57 ey = @intCast((ui >> 52) & 0x7ff); 58 ex = @intCast(ix >> 20); 59 if (ex - ey > 16) { // 2nd round, good to 118 bits 60 t = r; 61 w = @"fn" * pio2_2; 62 r = t - w; 63 w = @"fn" * pio2_2t - ((t - r) - w); 64 y[0] = r - w; 65 ui = @bitCast(y[0]); 66 ey = @intCast((ui >> 52) & 0x7ff); 67 if (ex - ey > 49) { // 3rd round, good to 151 bits, covers all cases 68 t = r; 69 w = @"fn" * pio2_3; 70 r = t - w; 71 w = @"fn" * pio2_3t - ((t - r) - w); 72 y[0] = r - w; 73 } 74 } 75 y[1] = (r - y[0]) - w; 76 return n; 77 } 78 79 // Returns the remainder of x rem pi/2 in y[0]+y[1] 80 // 81 // use rem_pio2_large() for large x 82 // 83 // caller must handle the case when reduction is not needed: |x| ~<= pi/4 */ 84 pub fn rem_pio2(x: f64, y: *[2]f64) i32 { 85 var z: f64 = undefined; 86 var tx: [3]f64 = undefined; 87 var ty: [2]f64 = undefined; 88 var n: i32 = undefined; 89 var ix: u32 = undefined; 90 var sign: bool = undefined; 91 var i: i32 = undefined; 92 var ui: u64 = undefined; 93 94 ui = @bitCast(x); 95 sign = ui >> 63 != 0; 96 ix = @truncate((ui >> 32) & 0x7fffffff); 97 if (ix <= 0x400f6a7a) { // |x| ~<= 5pi/4 98 if ((ix & 0xfffff) == 0x921fb) { // |x| ~= pi/2 or 2pi/2 99 return medium(ix, x, y); 100 } 101 if (ix <= 0x4002d97c) { // |x| ~<= 3pi/4 102 if (!sign) { 103 z = x - pio2_1; // one round good to 85 bits 104 y[0] = z - pio2_1t; 105 y[1] = (z - y[0]) - pio2_1t; 106 return 1; 107 } else { 108 z = x + pio2_1; 109 y[0] = z + pio2_1t; 110 y[1] = (z - y[0]) + pio2_1t; 111 return -1; 112 } 113 } else { 114 if (!sign) { 115 z = x - 2 * pio2_1; 116 y[0] = z - 2 * pio2_1t; 117 y[1] = (z - y[0]) - 2 * pio2_1t; 118 return 2; 119 } else { 120 z = x + 2 * pio2_1; 121 y[0] = z + 2 * pio2_1t; 122 y[1] = (z - y[0]) + 2 * pio2_1t; 123 return -2; 124 } 125 } 126 } 127 if (ix <= 0x401c463b) { // |x| ~<= 9pi/4 128 if (ix <= 0x4015fdbc) { // |x| ~<= 7pi/4 129 if (ix == 0x4012d97c) { // |x| ~= 3pi/2 130 return medium(ix, x, y); 131 } 132 if (!sign) { 133 z = x - 3 * pio2_1; 134 y[0] = z - 3 * pio2_1t; 135 y[1] = (z - y[0]) - 3 * pio2_1t; 136 return 3; 137 } else { 138 z = x + 3 * pio2_1; 139 y[0] = z + 3 * pio2_1t; 140 y[1] = (z - y[0]) + 3 * pio2_1t; 141 return -3; 142 } 143 } else { 144 if (ix == 0x401921fb) { // |x| ~= 4pi/2 */ 145 return medium(ix, x, y); 146 } 147 if (!sign) { 148 z = x - 4 * pio2_1; 149 y[0] = z - 4 * pio2_1t; 150 y[1] = (z - y[0]) - 4 * pio2_1t; 151 return 4; 152 } else { 153 z = x + 4 * pio2_1; 154 y[0] = z + 4 * pio2_1t; 155 y[1] = (z - y[0]) + 4 * pio2_1t; 156 return -4; 157 } 158 } 159 } 160 if (ix < 0x413921fb) { // |x| ~< 2^20*(pi/2), medium size 161 return medium(ix, x, y); 162 } 163 // all other (large) arguments 164 if (ix >= 0x7ff00000) { // x is inf or NaN 165 y[0] = x - x; 166 y[1] = y[0]; 167 return 0; 168 } 169 // set z = scalbn(|x|,-ilogb(x)+23) 170 ui = @bitCast(x); 171 ui &= std.math.maxInt(u64) >> 12; 172 ui |= @as(u64, 0x3ff + 23) << 52; 173 z = @bitCast(ui); 174 175 i = 0; 176 while (i < 2) : (i += 1) { 177 tx[@intCast(i)] = @floatFromInt(@as(i32, @intFromFloat(z))); 178 z = (z - tx[@intCast(i)]) * 0x1p24; 179 } 180 tx[@intCast(i)] = z; 181 // skip zero terms, first term is non-zero 182 while (tx[@intCast(i)] == 0.0) { 183 i -= 1; 184 } 185 n = rem_pio2_large(tx[0..], ty[0..], @as(i32, @intCast((ix >> 20))) - (0x3ff + 23), i + 1, 1); 186 if (sign) { 187 y[0] = -ty[0]; 188 y[1] = -ty[1]; 189 return -n; 190 } 191 y[0] = ty[0]; 192 y[1] = ty[1]; 193 return n; 194 }