fmod.zig (12218B) - Raw
1 const builtin = @import("builtin"); 2 const std = @import("std"); 3 const math = std.math; 4 const assert = std.debug.assert; 5 const arch = builtin.cpu.arch; 6 const common = @import("common.zig"); 7 const normalize = common.normalize; 8 9 pub const panic = common.panic; 10 11 comptime { 12 @export(&__fmodh, .{ .name = "__fmodh", .linkage = common.linkage, .visibility = common.visibility }); 13 @export(&fmodf, .{ .name = "fmodf", .linkage = common.linkage, .visibility = common.visibility }); 14 @export(&fmod, .{ .name = "fmod", .linkage = common.linkage, .visibility = common.visibility }); 15 @export(&__fmodx, .{ .name = "__fmodx", .linkage = common.linkage, .visibility = common.visibility }); 16 if (common.want_ppc_abi) { 17 @export(&fmodq, .{ .name = "fmodf128", .linkage = common.linkage, .visibility = common.visibility }); 18 } 19 @export(&fmodq, .{ .name = "fmodq", .linkage = common.linkage, .visibility = common.visibility }); 20 @export(&fmodl, .{ .name = "fmodl", .linkage = common.linkage, .visibility = common.visibility }); 21 } 22 23 pub fn __fmodh(x: f16, y: f16) callconv(.c) f16 { 24 // TODO: more efficient implementation 25 return @floatCast(fmodf(x, y)); 26 } 27 28 pub fn fmodf(x: f32, y: f32) callconv(.c) f32 { 29 return generic_fmod(f32, x, y); 30 } 31 32 pub fn fmod(x: f64, y: f64) callconv(.c) f64 { 33 return generic_fmod(f64, x, y); 34 } 35 36 /// fmodx - floating modulo large, returns the remainder of division for f80 types 37 /// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits 38 pub fn __fmodx(a: f80, b: f80) callconv(.c) f80 { 39 const T = f80; 40 const Z = std.meta.Int(.unsigned, @bitSizeOf(T)); 41 42 const significandBits = math.floatMantissaBits(T); 43 const fractionalBits = math.floatFractionalBits(T); 44 const exponentBits = math.floatExponentBits(T); 45 46 const signBit = (@as(Z, 1) << (significandBits + exponentBits)); 47 const maxExponent = ((1 << exponentBits) - 1); 48 49 var aRep: Z = @bitCast(a); 50 var bRep: Z = @bitCast(b); 51 52 const signA = aRep & signBit; 53 var expA: i32 = @intCast((@as(Z, @bitCast(a)) >> significandBits) & maxExponent); 54 var expB: i32 = @intCast((@as(Z, @bitCast(b)) >> significandBits) & maxExponent); 55 56 // There are 3 cases where the answer is undefined, check for: 57 // - fmodx(val, 0) 58 // - fmodx(val, NaN) 59 // - fmodx(inf, val) 60 // The sign on checked values does not matter. 61 // Doing (a * b) / (a * b) produces undefined results 62 // because the three cases always produce undefined calculations: 63 // - 0 / 0 64 // - val * NaN 65 // - inf / inf 66 if (b == 0 or math.isNan(b) or expA == maxExponent) { 67 return (a * b) / (a * b); 68 } 69 70 // Remove the sign from both 71 aRep &= ~signBit; 72 bRep &= ~signBit; 73 if (aRep <= bRep) { 74 if (aRep == bRep) { 75 return 0 * a; 76 } 77 return a; 78 } 79 80 if (expA == 0) expA = normalize(f80, &aRep); 81 if (expB == 0) expB = normalize(f80, &bRep); 82 83 var highA: u64 = 0; 84 const highB: u64 = 0; 85 var lowA: u64 = @truncate(aRep); 86 const lowB: u64 = @truncate(bRep); 87 88 while (expA > expB) : (expA -= 1) { 89 var high = highA -% highB; 90 const low = lowA -% lowB; 91 if (lowA < lowB) { 92 high -%= 1; 93 } 94 if (high >> 63 == 0) { 95 if ((high | low) == 0) { 96 return 0 * a; 97 } 98 highA = 2 *% high + (low >> 63); 99 lowA = 2 *% low; 100 } else { 101 highA = 2 *% highA + (lowA >> 63); 102 lowA = 2 *% lowA; 103 } 104 } 105 106 var high = highA -% highB; 107 const low = lowA -% lowB; 108 if (lowA < lowB) { 109 high -%= 1; 110 } 111 if (high >> 63 == 0) { 112 if ((high | low) == 0) { 113 return 0 * a; 114 } 115 highA = high; 116 lowA = low; 117 } 118 119 while ((lowA >> fractionalBits) == 0) { 120 lowA = 2 *% lowA; 121 expA = expA - 1; 122 } 123 124 // Combine the exponent with the sign and significand, normalize if happened to be denormalized 125 if (expA < -fractionalBits) { 126 return @bitCast(signA); 127 } else if (expA <= 0) { 128 return @bitCast((lowA >> @intCast(1 - expA)) | signA); 129 } else { 130 return @bitCast(lowA | (@as(Z, @as(u16, @intCast(expA))) << significandBits) | signA); 131 } 132 } 133 134 /// fmodq - floating modulo large, returns the remainder of division for f128 types 135 /// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits 136 pub fn fmodq(a: f128, b: f128) callconv(.c) f128 { 137 var amod = a; 138 var bmod = b; 139 const aPtr_u64: [*]u64 = @ptrCast(&amod); 140 const bPtr_u64: [*]u64 = @ptrCast(&bmod); 141 const aPtr_u16: [*]u16 = @ptrCast(&amod); 142 const bPtr_u16: [*]u16 = @ptrCast(&bmod); 143 144 const exp_and_sign_index = comptime switch (builtin.target.cpu.arch.endian()) { 145 .little => 7, 146 .big => 0, 147 }; 148 const low_index = comptime switch (builtin.target.cpu.arch.endian()) { 149 .little => 0, 150 .big => 1, 151 }; 152 const high_index = comptime switch (builtin.target.cpu.arch.endian()) { 153 .little => 1, 154 .big => 0, 155 }; 156 157 const signA = aPtr_u16[exp_and_sign_index] & 0x8000; 158 var expA: i32 = @intCast((aPtr_u16[exp_and_sign_index] & 0x7fff)); 159 var expB: i32 = @intCast((bPtr_u16[exp_and_sign_index] & 0x7fff)); 160 161 // There are 3 cases where the answer is undefined, check for: 162 // - fmodq(val, 0) 163 // - fmodq(val, NaN) 164 // - fmodq(inf, val) 165 // The sign on checked values does not matter. 166 // Doing (a * b) / (a * b) produces undefined results 167 // because the three cases always produce undefined calculations: 168 // - 0 / 0 169 // - val * NaN 170 // - inf / inf 171 if (b == 0 or std.math.isNan(b) or expA == 0x7fff) { 172 return (a * b) / (a * b); 173 } 174 175 // Remove the sign from both 176 aPtr_u16[exp_and_sign_index] = @bitCast(@as(i16, @intCast(expA))); 177 bPtr_u16[exp_and_sign_index] = @bitCast(@as(i16, @intCast(expB))); 178 if (amod <= bmod) { 179 if (amod == bmod) { 180 return 0 * a; 181 } 182 return a; 183 } 184 185 if (expA == 0) { 186 amod *= 0x1p120; 187 expA = @as(i32, aPtr_u16[exp_and_sign_index]) - 120; 188 } 189 190 if (expB == 0) { 191 bmod *= 0x1p120; 192 expB = @as(i32, bPtr_u16[exp_and_sign_index]) - 120; 193 } 194 195 // OR in extra non-stored mantissa digit 196 var highA: u64 = (aPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48; 197 const highB: u64 = (bPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48; 198 var lowA: u64 = aPtr_u64[low_index]; 199 const lowB: u64 = bPtr_u64[low_index]; 200 201 while (expA > expB) : (expA -= 1) { 202 var high = highA -% highB; 203 const low = lowA -% lowB; 204 if (lowA < lowB) { 205 high -%= 1; 206 } 207 if (high >> 63 == 0) { 208 if ((high | low) == 0) { 209 return 0 * a; 210 } 211 highA = 2 *% high + (low >> 63); 212 lowA = 2 *% low; 213 } else { 214 highA = 2 *% highA + (lowA >> 63); 215 lowA = 2 *% lowA; 216 } 217 } 218 219 var high = highA -% highB; 220 const low = lowA -% lowB; 221 if (lowA < lowB) { 222 high -= 1; 223 } 224 if (high >> 63 == 0) { 225 if ((high | low) == 0) { 226 return 0 * a; 227 } 228 highA = high; 229 lowA = low; 230 } 231 232 while (highA >> 48 == 0) { 233 highA = 2 *% highA + (lowA >> 63); 234 lowA = 2 *% lowA; 235 expA = expA - 1; 236 } 237 238 // Overwrite the current amod with the values in highA and lowA 239 aPtr_u64[high_index] = highA; 240 aPtr_u64[low_index] = lowA; 241 242 // Combine the exponent with the sign, normalize if happened to be denormalized 243 if (expA <= 0) { 244 aPtr_u16[exp_and_sign_index] = @as(u16, @truncate(@as(u32, @bitCast((expA +% 120))))) | signA; 245 amod *= 0x1p-120; 246 } else { 247 aPtr_u16[exp_and_sign_index] = @as(u16, @truncate(@as(u32, @bitCast(expA)))) | signA; 248 } 249 250 return amod; 251 } 252 253 pub fn fmodl(a: c_longdouble, b: c_longdouble) callconv(.c) c_longdouble { 254 switch (@typeInfo(c_longdouble).float.bits) { 255 16 => return __fmodh(a, b), 256 32 => return fmodf(a, b), 257 64 => return fmod(a, b), 258 80 => return __fmodx(a, b), 259 128 => return fmodq(a, b), 260 else => @compileError("unreachable"), 261 } 262 } 263 264 inline fn generic_fmod(comptime T: type, x: T, y: T) T { 265 const bits = @typeInfo(T).float.bits; 266 const uint = std.meta.Int(.unsigned, bits); 267 comptime assert(T == f32 or T == f64); 268 const digits = if (T == f32) 23 else 52; 269 const exp_bits = if (T == f32) 9 else 12; 270 const bits_minus_1 = bits - 1; 271 const mask = if (T == f32) 0xff else 0x7ff; 272 var ux: uint = @bitCast(x); 273 var uy: uint = @bitCast(y); 274 var ex: i32 = @intCast((ux >> digits) & mask); 275 var ey: i32 = @intCast((uy >> digits) & mask); 276 const sx = if (T == f32) @as(u32, @intCast(ux & 0x80000000)) else @as(i32, @intCast(ux >> bits_minus_1)); 277 var i: uint = undefined; 278 279 if (uy << 1 == 0 or math.isNan(@as(T, @bitCast(uy))) or ex == mask) 280 return (x * y) / (x * y); 281 282 if (ux << 1 <= uy << 1) { 283 if (ux << 1 == uy << 1) 284 return 0 * x; 285 return x; 286 } 287 288 // normalize x and y 289 if (ex == 0) { 290 i = ux << exp_bits; 291 while (i >> bits_minus_1 == 0) : ({ 292 ex -= 1; 293 i <<= 1; 294 }) {} 295 ux <<= @intCast(@as(u32, @bitCast(-ex + 1))); 296 } else { 297 ux &= math.maxInt(uint) >> exp_bits; 298 ux |= 1 << digits; 299 } 300 if (ey == 0) { 301 i = uy << exp_bits; 302 while (i >> bits_minus_1 == 0) : ({ 303 ey -= 1; 304 i <<= 1; 305 }) {} 306 uy <<= @intCast(@as(u32, @bitCast(-ey + 1))); 307 } else { 308 uy &= math.maxInt(uint) >> exp_bits; 309 uy |= 1 << digits; 310 } 311 312 // x mod y 313 while (ex > ey) : (ex -= 1) { 314 i = ux -% uy; 315 if (i >> bits_minus_1 == 0) { 316 if (i == 0) 317 return 0 * x; 318 ux = i; 319 } 320 ux <<= 1; 321 } 322 i = ux -% uy; 323 if (i >> bits_minus_1 == 0) { 324 if (i == 0) 325 return 0 * x; 326 ux = i; 327 } 328 while (ux >> digits == 0) : ({ 329 ux <<= 1; 330 ex -= 1; 331 }) {} 332 333 // scale result up 334 if (ex > 0) { 335 ux -%= 1 << digits; 336 ux |= @as(uint, @as(u32, @bitCast(ex))) << digits; 337 } else { 338 ux >>= @intCast(@as(u32, @bitCast(-ex + 1))); 339 } 340 if (T == f32) { 341 ux |= sx; 342 } else { 343 ux |= @as(uint, @intCast(sx)) << bits_minus_1; 344 } 345 return @bitCast(ux); 346 } 347 348 test "fmodf" { 349 const nan_val = math.nan(f32); 350 const inf_val = math.inf(f32); 351 352 try std.testing.expect(math.isNan(fmodf(nan_val, 1.0))); 353 try std.testing.expect(math.isNan(fmodf(1.0, nan_val))); 354 try std.testing.expect(math.isNan(fmodf(inf_val, 1.0))); 355 try std.testing.expect(math.isNan(fmodf(0.0, 0.0))); 356 try std.testing.expect(math.isNan(fmodf(1.0, 0.0))); 357 358 try std.testing.expectEqual(@as(f32, 0.0), fmodf(0.0, 2.0)); 359 try std.testing.expectEqual(@as(f32, -0.0), fmodf(-0.0, 2.0)); 360 361 try std.testing.expectEqual(@as(f32, -2.0), fmodf(-32.0, 10.0)); 362 try std.testing.expectEqual(@as(f32, -2.0), fmodf(-32.0, -10.0)); 363 try std.testing.expectEqual(@as(f32, 2.0), fmodf(32.0, 10.0)); 364 try std.testing.expectEqual(@as(f32, 2.0), fmodf(32.0, -10.0)); 365 } 366 367 test "fmod" { 368 const nan_val = math.nan(f64); 369 const inf_val = math.inf(f64); 370 371 try std.testing.expect(math.isNan(fmod(nan_val, 1.0))); 372 try std.testing.expect(math.isNan(fmod(1.0, nan_val))); 373 try std.testing.expect(math.isNan(fmod(inf_val, 1.0))); 374 try std.testing.expect(math.isNan(fmod(0.0, 0.0))); 375 try std.testing.expect(math.isNan(fmod(1.0, 0.0))); 376 377 try std.testing.expectEqual(@as(f64, 0.0), fmod(0.0, 2.0)); 378 try std.testing.expectEqual(@as(f64, -0.0), fmod(-0.0, 2.0)); 379 380 try std.testing.expectEqual(@as(f64, -2.0), fmod(-32.0, 10.0)); 381 try std.testing.expectEqual(@as(f64, -2.0), fmod(-32.0, -10.0)); 382 try std.testing.expectEqual(@as(f64, 2.0), fmod(32.0, 10.0)); 383 try std.testing.expectEqual(@as(f64, 2.0), fmod(32.0, -10.0)); 384 } 385 386 test { 387 _ = @import("fmodq_test.zig"); 388 _ = @import("fmodx_test.zig"); 389 }