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