master
1const builtin = @import("builtin");
2const std = @import("std");
3const math = std.math;
4const assert = std.debug.assert;
5const arch = builtin.cpu.arch;
6const common = @import("common.zig");
7const normalize = common.normalize;
8
9pub const panic = common.panic;
10
11comptime {
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
23pub fn __fmodh(x: f16, y: f16) callconv(.c) f16 {
24 // TODO: more efficient implementation
25 return @floatCast(fmodf(x, y));
26}
27
28pub fn fmodf(x: f32, y: f32) callconv(.c) f32 {
29 return generic_fmod(f32, x, y);
30}
31
32pub 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
38pub 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
136pub 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
253pub 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
264inline 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
348test "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
367test "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
386test {
387 _ = @import("fmodq_test.zig");
388 _ = @import("fmodx_test.zig");
389}