authorgravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2022-04-26 13:05:30-04:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2022-04-26 13:05:30-04:00
log0324975d1d6111b2631aa4c49831afb105f8e33d
tree63d82c1a738b99674a22a7de7a6b5cf741cecd94
parent5f2d0d414dc44af7bda0e8d805d038e8f1a6f9d3
parent28c05b0a53e451cbe06d321f4e2f92cf3b1a7be4
signaturebadge-question-mark Signed by PGP key 4AEE18F83AFDEB23

Merge pull request #11526 from topolarity/f80-div

compiler_rt: Implement `__divxf3` and `fmodx`

17 files changed, 668 insertions(+), 209 deletions(-)

lib/std/fmt.zig+4-3
...@@ -1124,6 +1124,7 @@ pub fn formatFloatHexadecimal(...@@ -1124,6 +1124,7 @@ pub fn formatFloatHexadecimal(
1124 const TU = std.meta.Int(.unsigned, std.meta.bitCount(T));1124 const TU = std.meta.Int(.unsigned, std.meta.bitCount(T));
11251125
1126 const mantissa_bits = math.floatMantissaBits(T);1126 const mantissa_bits = math.floatMantissaBits(T);
1127 const fractional_bits = math.floatFractionalBits(T);
1127 const exponent_bits = math.floatExponentBits(T);1128 const exponent_bits = math.floatExponentBits(T);
1128 const mantissa_mask = (1 << mantissa_bits) - 1;1129 const mantissa_mask = (1 << mantissa_bits) - 1;
1129 const exponent_mask = (1 << exponent_bits) - 1;1130 const exponent_mask = (1 << exponent_bits) - 1;
...@@ -1155,14 +1156,14 @@ pub fn formatFloatHexadecimal(...@@ -1155,14 +1156,14 @@ pub fn formatFloatHexadecimal(
1155 // Adjust the exponent for printing.1156 // Adjust the exponent for printing.
1156 exponent += 1;1157 exponent += 1;
1157 } else {1158 } else {
1158 // Add the implicit 1.1159 if (fractional_bits == mantissa_bits)
1159 mantissa |= 1 << mantissa_bits;1160 mantissa |= 1 << fractional_bits; // Add the implicit integer bit.
1160 }1161 }
11611162
1162 // Fill in zeroes to round the mantissa width to a multiple of 4.1163 // Fill in zeroes to round the mantissa width to a multiple of 4.
1163 if (T == f16) mantissa <<= 2 else if (T == f32) mantissa <<= 1;1164 if (T == f16) mantissa <<= 2 else if (T == f32) mantissa <<= 1;
11641165
1165 const mantissa_digits = (mantissa_bits + 3) / 4;1166 const mantissa_digits = (fractional_bits + 3) / 4;
11661167
1167 if (options.precision) |precision| {1168 if (options.precision) |precision| {
1168 // Round if needed.1169 // Round if needed.
lib/std/math/signbit.zig+6
...@@ -9,6 +9,7 @@ pub fn signbit(x: anytype) bool {...@@ -9,6 +9,7 @@ pub fn signbit(x: anytype) bool {
9 f16 => signbit16(x),9 f16 => signbit16(x),
10 f32 => signbit32(x),10 f32 => signbit32(x),
11 f64 => signbit64(x),11 f64 => signbit64(x),
12 f80 => signbit80(x),
12 f128 => signbit128(x),13 f128 => signbit128(x),
13 else => @compileError("signbit not implemented for " ++ @typeName(T)),14 else => @compileError("signbit not implemented for " ++ @typeName(T)),
14 };15 };
...@@ -29,6 +30,11 @@ fn signbit64(x: f64) bool {...@@ -29,6 +30,11 @@ fn signbit64(x: f64) bool {
29 return bits >> 63 != 0;30 return bits >> 63 != 0;
30}31}
3132
33fn signbit80(x: f80) bool {
34 const bits = @bitCast(u80, x);
35 return bits >> 79 != 0;
36}
37
32fn signbit128(x: f128) bool {38fn signbit128(x: f128) bool {
33 const bits = @bitCast(u128, x);39 const bits = @bitCast(u128, x);
34 return bits >> 127 != 0;40 return bits >> 127 != 0;
lib/std/special/compiler_rt.zig+14-5
...@@ -253,6 +253,8 @@ comptime {...@@ -253,6 +253,8 @@ comptime {
253 @export(__divsf3, .{ .name = "__divsf3", .linkage = linkage });253 @export(__divsf3, .{ .name = "__divsf3", .linkage = linkage });
254 const __divdf3 = @import("compiler_rt/divdf3.zig").__divdf3;254 const __divdf3 = @import("compiler_rt/divdf3.zig").__divdf3;
255 @export(__divdf3, .{ .name = "__divdf3", .linkage = linkage });255 @export(__divdf3, .{ .name = "__divdf3", .linkage = linkage });
256 const __divxf3 = @import("compiler_rt/divxf3.zig").__divxf3;
257 @export(__divxf3, .{ .name = "__divxf3", .linkage = linkage });
256 const __divtf3 = @import("compiler_rt/divtf3.zig").__divtf3;258 const __divtf3 = @import("compiler_rt/divtf3.zig").__divtf3;
257 @export(__divtf3, .{ .name = "__divtf3", .linkage = linkage });259 @export(__divtf3, .{ .name = "__divtf3", .linkage = linkage });
258260
...@@ -725,12 +727,18 @@ comptime {...@@ -725,12 +727,18 @@ comptime {
725 }727 }
726728
727 if (!is_test) {729 if (!is_test) {
728 @export(fmodl, .{ .name = "fmodl", .linkage = linkage });730 if (long_double_is_f80) {
729 if (long_double_is_f128) {731 @export(fmodx, .{ .name = "fmodl", .linkage = linkage });
730 @export(fmodl, .{ .name = "fmodq", .linkage = linkage });732 } else if (long_double_is_f128) {
733 @export(fmodq, .{ .name = "fmodl", .linkage = linkage });
731 } else {734 } else {
732 @export(fmodq, .{ .name = "fmodq", .linkage = linkage });735 @export(fmodl, .{ .name = "fmodl", .linkage = linkage });
736 }
737 if (long_double_is_f80 or builtin.zig_backend == .stage1) {
738 // TODO: https://github.com/ziglang/zig/issues/11161
739 @export(fmodx, .{ .name = "fmodx", .linkage = linkage });
733 }740 }
741 @export(fmodq, .{ .name = "fmodq", .linkage = linkage });
734742
735 @export(floorf, .{ .name = "floorf", .linkage = linkage });743 @export(floorf, .{ .name = "floorf", .linkage = linkage });
736 @export(floor, .{ .name = "floor", .linkage = linkage });744 @export(floor, .{ .name = "floor", .linkage = linkage });
...@@ -884,7 +892,8 @@ fn ceill(x: c_longdouble) callconv(.C) c_longdouble {...@@ -884,7 +892,8 @@ fn ceill(x: c_longdouble) callconv(.C) c_longdouble {
884 return math.ceil(x);892 return math.ceil(x);
885}893}
886894
887const fmodq = @import("compiler_rt/floatfmodq.zig").fmodq;895const fmodq = @import("compiler_rt/fmodq.zig").fmodq;
896const fmodx = @import("compiler_rt/fmodx.zig").fmodx;
888fn fmodl(x: c_longdouble, y: c_longdouble) callconv(.C) c_longdouble {897fn fmodl(x: c_longdouble, y: c_longdouble) callconv(.C) c_longdouble {
889 if (!long_double_is_f128) {898 if (!long_double_is_f128) {
890 @panic("TODO implement this");899 @panic("TODO implement this");
lib/std/special/compiler_rt/README.md+17-17
...@@ -152,53 +152,53 @@ Bugs should be solved by trying to duplicate the bug upstream, if possible....@@ -152,53 +152,53 @@ Bugs should be solved by trying to duplicate the bug upstream, if possible.
152- todo todo __fixsfsi // convert a to i32, rounding towards zero152- todo todo __fixsfsi // convert a to i32, rounding towards zero
153- todo todo __fixdfsi //153- todo todo __fixdfsi //
154- todo todo __fixtfsi //154- todo todo __fixtfsi //
155- none none __fixxfsi // missing155- todo todo __fixxfsi //
156- todo todo __fixsfdi // convert a to i64, rounding towards zero156- todo todo __fixsfdi // convert a to i64, rounding towards zero
157- todo todo __fixdfdi //157- todo todo __fixdfdi //
158- todo todo __fixtfdi //158- todo todo __fixtfdi //
159- none none __fixxfdi // missing159- todo todo __fixxfdi //
160- todo todo __fixsfti // convert a to i128, rounding towards zero160- todo todo __fixsfti // convert a to i128, rounding towards zero
161- todo todo __fixdfti //161- todo todo __fixdfti //
162- todo todo __fixtfdi //162- todo todo __fixtfdi //
163- none none __fixxfti // missing163- todo todo __fixxfti //
164164
165- __fixunssfsi // convert to u32, rounding towards zero. negative values become 0.165- __fixunssfsi // convert to u32, rounding towards zero. negative values become 0.
166- __fixunsdfsi //166- __fixunsdfsi //
167- __fixunstfsi //167- __fixunstfsi //
168- __fixunsxfsi // missing168- __fixunsxfsi //
169- __fixunssfdi // convert to u64, rounding towards zero. negative values become 0.169- __fixunssfdi // convert to u64, rounding towards zero. negative values become 0.
170- __fixunsdfdi //170- __fixunsdfdi //
171- __fixunstfdi //171- __fixunstfdi //
172- __fixunsxfdi // missing172- __fixunsxfdi //
173- __fixunssfti // convert to u128, rounding towards zero. negative values become 0.173- __fixunssfti // convert to u128, rounding towards zero. negative values become 0.
174- __fixunsdfti //174- __fixunsdfti //
175- __fixunstfdi //175- __fixunstfdi //
176- __fixunsxfti // missing176- __fixunsxfti //
177177
178- __floatsisf // convert i32 to floating point178- __floatsisf // convert i32 to floating point
179- __floatsidf //179- __floatsidf //
180- __floatsitf //180- __floatsitf //
181- __floatsixf // missing181- __floatsixf //
182- __floatdisf // convert i64 to floating point182- __floatdisf // convert i64 to floating point
183- __floatdidf //183- __floatdidf //
184- __floatditf //184- __floatditf //
185- __floatdixf // missing185- __floatdixf //
186- __floattisf // convert i128 to floating point186- __floattisf // convert i128 to floating point
187- __floattidf //187- __floattidf //
188- __floattixf // missing188- __floattixf //
189189
190- __floatunsisf // convert i32 to floating point190- __floatunsisf // convert u32 to floating point
191- __floatunsidf //191- __floatunsidf //
192- __floatunsitf //192- __floatunsitf //
193- __floatunsixf // missing193- __floatunsixf //
194- __floatundisf // convert i64 to floating point194- __floatundisf // convert u64 to floating point
195- __floatundidf //195- __floatundidf //
196- __floatunditf //196- __floatunditf //
197- __floatundixf // missing197- __floatundixf //
198- __floatuntisf // convert i128 to floating point198- __floatuntisf // convert u128 to floating point
199- __floatuntidf //199- __floatuntidf //
200- __floatuntitf //200- __floatuntitf //
201- __floatuntixf // missing201- __floatuntixf //
202202
203#### Float Comparison203#### Float Comparison
204- __cmpsf2 // return (a<b)=>-1,(a==b)=>0,(a>b)=>1,Nan=>1 dont rely on this204- __cmpsf2 // return (a<b)=>-1,(a==b)=>0,(a>b)=>1,Nan=>1 dont rely on this
...@@ -242,11 +242,11 @@ Bugs should be solved by trying to duplicate the bug upstream, if possible....@@ -242,11 +242,11 @@ Bugs should be solved by trying to duplicate the bug upstream, if possible.
242- __mulsf3 // a * b242- __mulsf3 // a * b
243- __muldf3 // a * b243- __muldf3 // a * b
244- __multf3 // a * b244- __multf3 // a * b
245- __mulxf3 // a * b missing245- __mulxf3 // a * b
246- __divsf3 // a / b246- __divsf3 // a / b
247- __divdf3 // a / b247- __divdf3 // a / b
248- __divtf3 // a / b248- __divtf3 // a / b
249- __divxf3 // a / b missing249- __divxf3 // a / b
250- __negsf2 // -a symbol-level compatibility: libgcc uses this for the rl78250- __negsf2 // -a symbol-level compatibility: libgcc uses this for the rl78
251- __negdf2 // -a unnecessary: can be lowered directly to a xor251- __negdf2 // -a unnecessary: can be lowered directly to a xor
252- __negtf2 // -a252- __negtf2 // -a
lib/std/special/compiler_rt/addXf3.zig+5-7
...@@ -74,7 +74,7 @@ fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T...@@ -74,7 +74,7 @@ fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T
7474
75 const shift = @clz(std.meta.Int(.unsigned, bits), significand.*) - @clz(Z, integerBit);75 const shift = @clz(std.meta.Int(.unsigned, bits), significand.*) - @clz(Z, integerBit);
76 significand.* <<= @intCast(S, shift);76 significand.* <<= @intCast(S, shift);
77 return 1 - shift;77 return @as(i32, 1) - shift;
78}78}
7979
80// TODO: restore inline keyword, see: https://github.com/ziglang/zig/issues/215480// TODO: restore inline keyword, see: https://github.com/ziglang/zig/issues/2154
...@@ -210,12 +210,10 @@ fn addXf3(comptime T: type, a: T, b: T) T {...@@ -210,12 +210,10 @@ fn addXf3(comptime T: type, a: T, b: T) T {
210 if (aExponent >= maxExponent) return @bitCast(T, infRep | resultSign);210 if (aExponent >= maxExponent) return @bitCast(T, infRep | resultSign);
211211
212 if (aExponent <= 0) {212 if (aExponent <= 0) {
213 // Result is denormal before rounding; the exponent is zero and we213 // Result is denormal; the exponent and round/sticky bits are zero.
214 // need to shift the significand.214 // All we need to do is shift the significand and apply the correct sign.
215 const shift = @intCast(Z, 1 - aExponent);215 aSignificand >>= @intCast(S, 4 - aExponent);
216 const sticky = if (aSignificand << @intCast(S, typeWidth - shift) != 0) @as(Z, 1) else 0;216 return @bitCast(T, resultSign | aSignificand);
217 aSignificand = aSignificand >> @intCast(S, shift | sticky);
218 aExponent = 0;
219 }217 }
220218
221 // Low three bits are round, guard, and sticky.219 // Low three bits are round, guard, and sticky.
lib/std/special/compiler_rt/addXf3_test.zig+2
...@@ -152,4 +152,6 @@ test "addxf3" {...@@ -152,4 +152,6 @@ test "addxf3" {
152 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.8p-63, 0x3FFF_8800000000000000); // round down to even152 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.8p-63, 0x3FFF_8800000000000000); // round down to even
153 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.9p-63, 0x3FFF_8800000000000001); // round up153 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.9p-63, 0x3FFF_8800000000000001); // round up
154 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x2.0p-63, 0x3FFF_8800000000000001); // exact154 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x2.0p-63, 0x3FFF_8800000000000001); // exact
155 try test__addxf3(0x0.ffff_ffff_ffff_fffcp-16382, 0x0.0000_0000_0000_0002p-16382, 0x0000_7FFFFFFFFFFFFFFF); // exact
156 try test__addxf3(0x0.1fff_ffff_ffff_fffcp-16382, 0x0.0000_0000_0000_0002p-16382, 0x0000_0FFFFFFFFFFFFFFF); // exact
155}157}
lib/std/special/compiler_rt/divdf3.zig+3-4
...@@ -314,12 +314,11 @@ pub fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void {...@@ -314,12 +314,11 @@ pub fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void {
314pub fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T).Float.bits)) i32 {314pub fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T).Float.bits)) i32 {
315 @setRuntimeSafety(builtin.is_test);315 @setRuntimeSafety(builtin.is_test);
316 const Z = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);316 const Z = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
317 const significandBits = std.math.floatMantissaBits(T);317 const integerBit = @as(Z, 1) << std.math.floatFractionalBits(T);
318 const implicitBit = @as(Z, 1) << significandBits;
319318
320 const shift = @clz(Z, significand.*) - @clz(Z, implicitBit);319 const shift = @clz(Z, significand.*) - @clz(Z, integerBit);
321 significand.* <<= @intCast(std.math.Log2Int(Z), shift);320 significand.* <<= @intCast(std.math.Log2Int(Z), shift);
322 return 1 - shift;321 return @as(i32, 1) - shift;
323}322}
324323
325pub fn __aeabi_ddiv(a: f64, b: f64) callconv(.AAPCS) f64 {324pub fn __aeabi_ddiv(a: f64, b: f64) callconv(.AAPCS) f64 {
lib/std/special/compiler_rt/divxf3.zig created+202
...@@ -0,0 +1,202 @@
1const std = @import("std");
2const builtin = @import("builtin");
3const normalize = @import("divdf3.zig").normalize;
4const wideMultiply = @import("divdf3.zig").wideMultiply;
5
6pub fn __divxf3(a: f80, b: f80) callconv(.C) f80 {
7 @setRuntimeSafety(builtin.is_test);
8 const T = f80;
9 const Z = std.meta.Int(.unsigned, @bitSizeOf(T));
10
11 const significandBits = std.math.floatMantissaBits(T);
12 const fractionalBits = std.math.floatFractionalBits(T);
13 const exponentBits = std.math.floatExponentBits(T);
14
15 const signBit = (@as(Z, 1) << (significandBits + exponentBits));
16 const maxExponent = ((1 << exponentBits) - 1);
17 const exponentBias = (maxExponent >> 1);
18
19 const integerBit = (@as(Z, 1) << fractionalBits);
20 const quietBit = integerBit >> 1;
21 const significandMask = (@as(Z, 1) << significandBits) - 1;
22
23 const absMask = signBit - 1;
24 const qnanRep = @bitCast(Z, std.math.nan(T)) | quietBit;
25 const infRep = @bitCast(Z, std.math.inf(T));
26
27 const aExponent = @truncate(u32, (@bitCast(Z, a) >> significandBits) & maxExponent);
28 const bExponent = @truncate(u32, (@bitCast(Z, b) >> significandBits) & maxExponent);
29 const quotientSign: Z = (@bitCast(Z, a) ^ @bitCast(Z, b)) & signBit;
30
31 var aSignificand: Z = @bitCast(Z, a) & significandMask;
32 var bSignificand: Z = @bitCast(Z, b) & significandMask;
33 var scale: i32 = 0;
34
35 // Detect if a or b is zero, denormal, infinity, or NaN.
36 if (aExponent -% 1 >= maxExponent - 1 or bExponent -% 1 >= maxExponent - 1) {
37 const aAbs: Z = @bitCast(Z, a) & absMask;
38 const bAbs: Z = @bitCast(Z, b) & absMask;
39
40 // NaN / anything = qNaN
41 if (aAbs > infRep) return @bitCast(T, @bitCast(Z, a) | quietBit);
42 // anything / NaN = qNaN
43 if (bAbs > infRep) return @bitCast(T, @bitCast(Z, b) | quietBit);
44
45 if (aAbs == infRep) {
46 // infinity / infinity = NaN
47 if (bAbs == infRep) {
48 return @bitCast(T, qnanRep);
49 }
50 // infinity / anything else = +/- infinity
51 else {
52 return @bitCast(T, aAbs | quotientSign);
53 }
54 }
55
56 // anything else / infinity = +/- 0
57 if (bAbs == infRep) return @bitCast(T, quotientSign);
58
59 if (aAbs == 0) {
60 // zero / zero = NaN
61 if (bAbs == 0) {
62 return @bitCast(T, qnanRep);
63 }
64 // zero / anything else = +/- zero
65 else {
66 return @bitCast(T, quotientSign);
67 }
68 }
69 // anything else / zero = +/- infinity
70 if (bAbs == 0) return @bitCast(T, infRep | quotientSign);
71
72 // one or both of a or b is denormal, the other (if applicable) is a
73 // normal number. Renormalize one or both of a and b, and set scale to
74 // include the necessary exponent adjustment.
75 if (aAbs < integerBit) scale +%= normalize(T, &aSignificand);
76 if (bAbs < integerBit) scale -%= normalize(T, &bSignificand);
77 }
78 var quotientExponent: i32 = @bitCast(i32, aExponent -% bExponent) +% scale;
79
80 // Align the significand of b as a Q63 fixed-point number in the range
81 // [1, 2.0) and get a Q64 approximate reciprocal using a small minimax
82 // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This
83 // is accurate to about 3.5 binary digits.
84 const q63b = @intCast(u64, bSignificand);
85 var recip64 = @as(u64, 0x7504f333F9DE6484) -% q63b;
86 // 0x7504f333F9DE6484 / 2^64 + 1 = 3/4 + 1/sqrt(2)
87
88 // Now refine the reciprocal estimate using a Newton-Raphson iteration:
89 //
90 // x1 = x0 * (2 - x0 * b)
91 //
92 // This doubles the number of correct binary digits in the approximation
93 // with each iteration.
94 var correction64: u64 = undefined;
95 correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
96 recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
97 correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
98 recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
99 correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
100 recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
101 correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
102 recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
103 correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
104 recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
105
106 // The reciprocal may have overflowed to zero if the upper half of b is
107 // exactly 1.0. This would sabatoge the full-width final stage of the
108 // computation that follows, so we adjust the reciprocal down by one bit.
109 recip64 -%= 1;
110
111 // We need to perform one more iteration to get us to 112 binary digits;
112 // The last iteration needs to happen with extra precision.
113
114 // NOTE: This operation is equivalent to __multi3, which is not implemented
115 // in some architechures
116 var reciprocal: u128 = undefined;
117 var correction: u128 = undefined;
118 var dummy: u128 = undefined;
119 wideMultiply(u128, recip64, q63b, &dummy, &correction);
120
121 correction = -%correction;
122
123 const cHi = @truncate(u64, correction >> 64);
124 const cLo = @truncate(u64, correction);
125
126 var r64cH: u128 = undefined;
127 var r64cL: u128 = undefined;
128 wideMultiply(u128, recip64, cHi, &dummy, &r64cH);
129 wideMultiply(u128, recip64, cLo, &dummy, &r64cL);
130
131 reciprocal = r64cH + (r64cL >> 64);
132
133 // Adjust the final 128-bit reciprocal estimate downward to ensure that it
134 // is strictly smaller than the infinitely precise exact reciprocal. Because
135 // the computation of the Newton-Raphson step is truncating at every step,
136 // this adjustment is small; most of the work is already done.
137 reciprocal -%= 2;
138
139 // The numerical reciprocal is accurate to within 2^-112, lies in the
140 // interval [0.5, 1.0), and is strictly smaller than the true reciprocal
141 // of b. Multiplying a by this reciprocal thus gives a numerical q = a/b
142 // in Q127 with the following properties:
143 //
144 // 1. q < a/b
145 // 2. q is in the interval [0.5, 2.0)
146 // 3. The error in q is bounded away from 2^-63 (actually, we have
147 // many bits to spare, but this is all we need).
148
149 // We need a 128 x 128 multiply high to compute q.
150 var quotient128: u128 = undefined;
151 var quotientLo: u128 = undefined;
152 wideMultiply(u128, aSignificand << 2, reciprocal, &quotient128, &quotientLo);
153
154 // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).
155 // Right shift the quotient if it falls in the [1,2) range and adjust the
156 // exponent accordingly.
157 var quotient: u64 = if (quotient128 < (integerBit << 1)) b: {
158 quotientExponent -= 1;
159 break :b @intCast(u64, quotient128);
160 } else @intCast(u64, quotient128 >> 1);
161
162 // We are going to compute a residual of the form
163 //
164 // r = a - q*b
165 //
166 // We know from the construction of q that r satisfies:
167 //
168 // 0 <= r < ulp(q)*b
169 //
170 // If r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we
171 // already have the correct result. The exact halfway case cannot occur.
172 var residual: u64 = -%(quotient *% q63b);
173
174 const writtenExponent = quotientExponent + exponentBias;
175 if (writtenExponent >= maxExponent) {
176 // If we have overflowed the exponent, return infinity.
177 return @bitCast(T, infRep | quotientSign);
178 } else if (writtenExponent < 1) {
179 if (writtenExponent == 0) {
180 // Check whether the rounded result is normal.
181 if (residual > (bSignificand >> 1)) { // round
182 if (quotient == (integerBit - 1)) // If the rounded result is normal, return it
183 return @bitCast(T, @bitCast(Z, std.math.floatMin(T)) | quotientSign);
184 }
185 }
186 // Flush denormals to zero. In the future, it would be nice to add
187 // code to round them correctly.
188 return @bitCast(T, quotientSign);
189 } else {
190 const round = @boolToInt(residual > (bSignificand >> 1));
191 // Insert the exponent
192 var absResult = quotient | (@intCast(Z, writtenExponent) << significandBits);
193 // Round
194 absResult +%= round;
195 // Insert the sign and return
196 return @bitCast(T, absResult | quotientSign | integerBit);
197 }
198}
199
200test {
201 _ = @import("divxf3_test.zig");
202}
lib/std/special/compiler_rt/divxf3_test.zig created+65
...@@ -0,0 +1,65 @@
1const std = @import("std");
2const math = std.math;
3const testing = std.testing;
4
5const __divxf3 = @import("divxf3.zig").__divxf3;
6
7fn compareResult(result: f80, expected: u80) bool {
8 const rep = @bitCast(u80, result);
9
10 if (rep == expected) return true;
11 // test other possible NaN representations (signal NaN)
12 if (math.isNan(result) and math.isNan(@bitCast(f80, expected))) return true;
13
14 return false;
15}
16
17fn expect__divxf3_result(a: f80, b: f80, expected: u80) !void {
18 const x = __divxf3(a, b);
19 const ret = compareResult(x, expected);
20 try testing.expect(ret == true);
21}
22
23fn test__divxf3(a: f80, b: f80) !void {
24 const integerBit = 1 << math.floatFractionalBits(f80);
25 const x = __divxf3(a, b);
26
27 // Next float (assuming normal, non-zero result)
28 const x_plus_eps = @bitCast(f80, (@bitCast(u80, x) + 1) | integerBit);
29 // Prev float (assuming normal, non-zero result)
30 const x_minus_eps = @bitCast(f80, (@bitCast(u80, x) - 1) | integerBit);
31
32 // Make sure result is more accurate than the adjacent floats
33 const err_x = std.math.fabs(@mulAdd(f80, x, b, -a));
34 const err_x_plus_eps = std.math.fabs(@mulAdd(f80, x_plus_eps, b, -a));
35 const err_x_minus_eps = std.math.fabs(@mulAdd(f80, x_minus_eps, b, -a));
36
37 try testing.expect(err_x_minus_eps > err_x);
38 try testing.expect(err_x_plus_eps > err_x);
39}
40
41test "divxf3" {
42 // qNaN / any = qNaN
43 try expect__divxf3_result(math.qnan_f80, 0x1.23456789abcdefp+5, 0x7fffC000000000000000);
44 // NaN / any = NaN
45 try expect__divxf3_result(math.nan_f80, 0x1.23456789abcdefp+5, 0x7fffC000000000000000);
46 // inf / any(except inf and nan) = inf
47 try expect__divxf3_result(math.inf(f80), 0x1.23456789abcdefp+5, 0x7fff8000000000000000);
48 // inf / inf = nan
49 try expect__divxf3_result(math.inf(f80), math.inf(f80), 0x7fffC000000000000000);
50 // inf / nan = nan
51 try expect__divxf3_result(math.inf(f80), math.nan(f80), 0x7fffC000000000000000);
52
53 try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.eedcbaba3a94546558237654321fp-1);
54 try test__divxf3(0x1.a2b34c56d745382f9abf2c3dfeffp-50, 0x1.ed2c3ba15935332532287654321fp-9);
55 try test__divxf3(0x1.2345f6aaaa786555f42432abcdefp+456, 0x1.edacbba9874f765463544dd3621fp+6400);
56 try test__divxf3(0x1.2d3456f789ba6322bc665544edefp-234, 0x1.eddcdba39f3c8b7a36564354321fp-4455);
57 try test__divxf3(0x1.2345f6b77b7a8953365433abcdefp+234, 0x1.edcba987d6bb3aa467754354321fp-4055);
58 try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.a2b34c56d745382f9abf2c3dfeffp-50);
59 try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.1234567890abcdef987654321123p0);
60 try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.12394205810257120adae8929f23p+16);
61 try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.febdcefa1231245f9abf2c3dfeffp-50);
62
63 // Result rounds down to zero
64 try expect__divxf3_result(6.72420628622418701252535563464350521E-4932, 2.0, 0x0);
65}
lib/std/special/compiler_rt/floatfmodq.zig deleted-126
...@@ -1,126 +0,0 @@
1const builtin = @import("builtin");
2const std = @import("std");
3
4// fmodq - floating modulo large, returns the remainder of division for f128 types
5// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits
6pub fn fmodq(a: f128, b: f128) callconv(.C) f128 {
7 @setRuntimeSafety(builtin.is_test);
8 var amod = a;
9 var bmod = b;
10 const aPtr_u64 = @ptrCast([*]u64, &amod);
11 const bPtr_u64 = @ptrCast([*]u64, &bmod);
12 const aPtr_u16 = @ptrCast([*]u16, &amod);
13 const bPtr_u16 = @ptrCast([*]u16, &bmod);
14
15 const exp_and_sign_index = comptime switch (builtin.target.cpu.arch.endian()) {
16 .Little => 7,
17 .Big => 0,
18 };
19 const low_index = comptime switch (builtin.target.cpu.arch.endian()) {
20 .Little => 0,
21 .Big => 1,
22 };
23 const high_index = comptime switch (builtin.target.cpu.arch.endian()) {
24 .Little => 1,
25 .Big => 0,
26 };
27
28 const signA = aPtr_u16[exp_and_sign_index] & 0x8000;
29 var expA = @intCast(i32, (aPtr_u16[exp_and_sign_index] & 0x7fff));
30 var expB = bPtr_u16[exp_and_sign_index] & 0x7fff;
31
32 // There are 3 cases where the answer is undefined, check for:
33 // - fmodq(val, 0)
34 // - fmodq(val, NaN)
35 // - fmodq(inf, val)
36 // The sign on checked values does not matter.
37 // Doing (a * b) / (a * b) procudes undefined results
38 // because the three cases always produce undefined calculations:
39 // - 0 / 0
40 // - val * NaN
41 // - inf / inf
42 if (b == 0 or std.math.isNan(b) or expA == 0x7fff) {
43 return (a * b) / (a * b);
44 }
45
46 // Remove the sign from both
47 aPtr_u16[exp_and_sign_index] = @bitCast(u16, @intCast(i16, expA));
48 bPtr_u16[exp_and_sign_index] = @bitCast(u16, @intCast(i16, expB));
49 if (amod <= bmod) {
50 if (amod == bmod) {
51 return 0 * a;
52 }
53 return a;
54 }
55
56 if (expA == 0) {
57 amod *= 0x1p120;
58 expA = aPtr_u16[exp_and_sign_index] -% 120;
59 }
60
61 if (expB == 0) {
62 bmod *= 0x1p120;
63 expB = bPtr_u16[exp_and_sign_index] -% 120;
64 }
65
66 // OR in extra non-stored mantissa digit
67 var highA: u64 = (aPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48;
68 var highB: u64 = (bPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48;
69 var lowA: u64 = aPtr_u64[low_index];
70 var lowB: u64 = bPtr_u64[low_index];
71
72 while (expA > expB) : (expA -= 1) {
73 var high = highA -% highB;
74 var low = lowA -% lowB;
75 if (lowA < lowB) {
76 high = highA -% 1;
77 }
78 if (high >> 63 == 0) {
79 if ((high | low) == 0) {
80 return 0 * a;
81 }
82 highA = 2 *% high + (low >> 63);
83 lowA = 2 *% low;
84 } else {
85 highA = 2 *% highA + (lowA >> 63);
86 lowA = 2 *% lowA;
87 }
88 }
89
90 var high = highA -% highB;
91 var low = lowA -% lowB;
92 if (lowA < lowB) {
93 high -= 1;
94 }
95 if (high >> 63 == 0) {
96 if ((high | low) == 0) {
97 return 0 * a;
98 }
99 highA = high;
100 lowA = low;
101 }
102
103 while (highA >> 48 == 0) {
104 highA = 2 *% highA + (lowA >> 63);
105 lowA = 2 *% lowA;
106 expA = expA - 1;
107 }
108
109 // Overwrite the current amod with the values in highA and lowA
110 aPtr_u64[high_index] = highA;
111 aPtr_u64[low_index] = lowA;
112
113 // Combine the exponent with the sign, normalize if happend to be denormalized
114 if (expA <= 0) {
115 aPtr_u16[exp_and_sign_index] = @truncate(u16, @bitCast(u32, (expA +% 120))) | signA;
116 amod *= 0x1p-120;
117 } else {
118 aPtr_u16[exp_and_sign_index] = @truncate(u16, @bitCast(u32, expA)) | signA;
119 }
120
121 return amod;
122}
123
124test {
125 _ = @import("floatfmodq_test.zig");
126}
lib/std/special/compiler_rt/floatfmodq_test.zig deleted-46
...@@ -1,46 +0,0 @@
1const std = @import("std");
2const fmodq = @import("floatfmodq.zig");
3const testing = std.testing;
4
5fn test_fmodq(a: f128, b: f128, exp: f128) !void {
6 const res = fmodq.fmodq(a, b);
7 try testing.expect(exp == res);
8}
9
10fn test_fmodq_nans() !void {
11 try testing.expect(std.math.isNan(fmodq.fmodq(1.0, std.math.nan(f128))));
12 try testing.expect(std.math.isNan(fmodq.fmodq(1.0, -std.math.nan(f128))));
13 try testing.expect(std.math.isNan(fmodq.fmodq(std.math.nan(f128), 1.0)));
14 try testing.expect(std.math.isNan(fmodq.fmodq(-std.math.nan(f128), 1.0)));
15}
16
17fn test_fmodq_infs() !void {
18 try testing.expect(fmodq.fmodq(1.0, std.math.inf(f128)) == 1.0);
19 try testing.expect(fmodq.fmodq(1.0, -std.math.inf(f128)) == 1.0);
20 try testing.expect(std.math.isNan(fmodq.fmodq(std.math.inf(f128), 1.0)));
21 try testing.expect(std.math.isNan(fmodq.fmodq(-std.math.inf(f128), 1.0)));
22}
23
24test "fmodq" {
25 try test_fmodq(6.8, 4.0, 2.8);
26 try test_fmodq(6.8, -4.0, 2.8);
27 try test_fmodq(-6.8, 4.0, -2.8);
28 try test_fmodq(-6.8, -4.0, -2.8);
29 try test_fmodq(3.0, 2.0, 1.0);
30 try test_fmodq(-5.0, 3.0, -2.0);
31 try test_fmodq(3.0, 2.0, 1.0);
32 try test_fmodq(1.0, 2.0, 1.0);
33 try test_fmodq(0.0, 1.0, 0.0);
34 try test_fmodq(-0.0, 1.0, -0.0);
35 try test_fmodq(7046119.0, 5558362.0, 1487757.0);
36 try test_fmodq(9010357.0, 1957236.0, 1181413.0);
37
38 // Denormals
39 const a: f128 = 0xedcb34a235253948765432134674p-16494;
40 const b: f128 = 0x5d2e38791cfbc0737402da5a9518p-16494;
41 const exp: f128 = 0x336ec3affb2db8618e4e7d5e1c44p-16494;
42 try test_fmodq(a, b, exp);
43
44 try test_fmodq_nans();
45 try test_fmodq_infs();
46}
lib/std/special/compiler_rt/fmodq.zig created+126
...@@ -0,0 +1,126 @@
1const builtin = @import("builtin");
2const std = @import("std");
3
4// fmodq - floating modulo large, returns the remainder of division for f128 types
5// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits
6pub fn fmodq(a: f128, b: f128) callconv(.C) f128 {
7 @setRuntimeSafety(builtin.is_test);
8 var amod = a;
9 var bmod = b;
10 const aPtr_u64 = @ptrCast([*]u64, &amod);
11 const bPtr_u64 = @ptrCast([*]u64, &bmod);
12 const aPtr_u16 = @ptrCast([*]u16, &amod);
13 const bPtr_u16 = @ptrCast([*]u16, &bmod);
14
15 const exp_and_sign_index = comptime switch (builtin.target.cpu.arch.endian()) {
16 .Little => 7,
17 .Big => 0,
18 };
19 const low_index = comptime switch (builtin.target.cpu.arch.endian()) {
20 .Little => 0,
21 .Big => 1,
22 };
23 const high_index = comptime switch (builtin.target.cpu.arch.endian()) {
24 .Little => 1,
25 .Big => 0,
26 };
27
28 const signA = aPtr_u16[exp_and_sign_index] & 0x8000;
29 var expA = @intCast(i32, (aPtr_u16[exp_and_sign_index] & 0x7fff));
30 var expB = @intCast(i32, (bPtr_u16[exp_and_sign_index] & 0x7fff));
31
32 // There are 3 cases where the answer is undefined, check for:
33 // - fmodq(val, 0)
34 // - fmodq(val, NaN)
35 // - fmodq(inf, val)
36 // The sign on checked values does not matter.
37 // Doing (a * b) / (a * b) procudes undefined results
38 // because the three cases always produce undefined calculations:
39 // - 0 / 0
40 // - val * NaN
41 // - inf / inf
42 if (b == 0 or std.math.isNan(b) or expA == 0x7fff) {
43 return (a * b) / (a * b);
44 }
45
46 // Remove the sign from both
47 aPtr_u16[exp_and_sign_index] = @bitCast(u16, @intCast(i16, expA));
48 bPtr_u16[exp_and_sign_index] = @bitCast(u16, @intCast(i16, expB));
49 if (amod <= bmod) {
50 if (amod == bmod) {
51 return 0 * a;
52 }
53 return a;
54 }
55
56 if (expA == 0) {
57 amod *= 0x1p120;
58 expA = @as(i32, aPtr_u16[exp_and_sign_index]) - 120;
59 }
60
61 if (expB == 0) {
62 bmod *= 0x1p120;
63 expB = @as(i32, bPtr_u16[exp_and_sign_index]) - 120;
64 }
65
66 // OR in extra non-stored mantissa digit
67 var highA: u64 = (aPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48;
68 var highB: u64 = (bPtr_u64[high_index] & (std.math.maxInt(u64) >> 16)) | 1 << 48;
69 var lowA: u64 = aPtr_u64[low_index];
70 var lowB: u64 = bPtr_u64[low_index];
71
72 while (expA > expB) : (expA -= 1) {
73 var high = highA -% highB;
74 var low = lowA -% lowB;
75 if (lowA < lowB) {
76 high -%= 1;
77 }
78 if (high >> 63 == 0) {
79 if ((high | low) == 0) {
80 return 0 * a;
81 }
82 highA = 2 *% high + (low >> 63);
83 lowA = 2 *% low;
84 } else {
85 highA = 2 *% highA + (lowA >> 63);
86 lowA = 2 *% lowA;
87 }
88 }
89
90 var high = highA -% highB;
91 var low = lowA -% lowB;
92 if (lowA < lowB) {
93 high -= 1;
94 }
95 if (high >> 63 == 0) {
96 if ((high | low) == 0) {
97 return 0 * a;
98 }
99 highA = high;
100 lowA = low;
101 }
102
103 while (highA >> 48 == 0) {
104 highA = 2 *% highA + (lowA >> 63);
105 lowA = 2 *% lowA;
106 expA = expA - 1;
107 }
108
109 // Overwrite the current amod with the values in highA and lowA
110 aPtr_u64[high_index] = highA;
111 aPtr_u64[low_index] = lowA;
112
113 // Combine the exponent with the sign, normalize if happend to be denormalized
114 if (expA <= 0) {
115 aPtr_u16[exp_and_sign_index] = @truncate(u16, @bitCast(u32, (expA +% 120))) | signA;
116 amod *= 0x1p-120;
117 } else {
118 aPtr_u16[exp_and_sign_index] = @truncate(u16, @bitCast(u32, expA)) | signA;
119 }
120
121 return amod;
122}
123
124test {
125 _ = @import("fmodq_test.zig");
126}
lib/std/special/compiler_rt/fmodq_test.zig created+52
...@@ -0,0 +1,52 @@
1const std = @import("std");
2const fmodq = @import("fmodq.zig");
3const testing = std.testing;
4
5fn test_fmodq(a: f128, b: f128, exp: f128) !void {
6 const res = fmodq.fmodq(a, b);
7 try testing.expect(exp == res);
8}
9
10fn test_fmodq_nans() !void {
11 try testing.expect(std.math.isNan(fmodq.fmodq(1.0, std.math.nan(f128))));
12 try testing.expect(std.math.isNan(fmodq.fmodq(1.0, -std.math.nan(f128))));
13 try testing.expect(std.math.isNan(fmodq.fmodq(std.math.nan(f128), 1.0)));
14 try testing.expect(std.math.isNan(fmodq.fmodq(-std.math.nan(f128), 1.0)));
15}
16
17fn test_fmodq_infs() !void {
18 try testing.expect(fmodq.fmodq(1.0, std.math.inf(f128)) == 1.0);
19 try testing.expect(fmodq.fmodq(1.0, -std.math.inf(f128)) == 1.0);
20 try testing.expect(std.math.isNan(fmodq.fmodq(std.math.inf(f128), 1.0)));
21 try testing.expect(std.math.isNan(fmodq.fmodq(-std.math.inf(f128), 1.0)));
22}
23
24test "fmodq" {
25 try test_fmodq(6.8, 4.0, 2.8);
26 try test_fmodq(6.8, -4.0, 2.8);
27 try test_fmodq(-6.8, 4.0, -2.8);
28 try test_fmodq(-6.8, -4.0, -2.8);
29 try test_fmodq(3.0, 2.0, 1.0);
30 try test_fmodq(-5.0, 3.0, -2.0);
31 try test_fmodq(3.0, 2.0, 1.0);
32 try test_fmodq(1.0, 2.0, 1.0);
33 try test_fmodq(0.0, 1.0, 0.0);
34 try test_fmodq(-0.0, 1.0, -0.0);
35 try test_fmodq(7046119.0, 5558362.0, 1487757.0);
36 try test_fmodq(9010357.0, 1957236.0, 1181413.0);
37 try test_fmodq(5192296858534827628530496329220095, 10.0, 5.0);
38 try test_fmodq(5192296858534827628530496329220095, 922337203681230954775807, 220474884073715748246157);
39
40 // Denormals
41 const a1: f128 = 0xedcb34a235253948765432134674p-16494;
42 const b1: f128 = 0x5d2e38791cfbc0737402da5a9518p-16494;
43 const exp1: f128 = 0x336ec3affb2db8618e4e7d5e1c44p-16494;
44 try test_fmodq(a1, b1, exp1);
45 const a2: f128 = 0x0.7654_3210_fdec_ba98_7654_3210_fdecp-16382;
46 const b2: f128 = 0x0.0012_fdac_bdef_1234_fdec_3222_1111p-16382;
47 const exp2: f128 = 0x0.0001_aecd_9d66_4a6e_67b7_d7d0_a901p-16382;
48 try test_fmodq(a2, b2, exp2);
49
50 try test_fmodq_nans();
51 try test_fmodq_infs();
52}
lib/std/special/compiler_rt/fmodx.zig created+108
...@@ -0,0 +1,108 @@
1const builtin = @import("builtin");
2const std = @import("std");
3const math = std.math;
4const normalize = @import("divdf3.zig").normalize;
5
6// fmodx - floating modulo large, returns the remainder of division for f80 types
7// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits
8pub fn fmodx(a: f80, b: f80) callconv(.C) f80 {
9 @setRuntimeSafety(builtin.is_test);
10
11 const T = f80;
12 const Z = std.meta.Int(.unsigned, @bitSizeOf(T));
13
14 const significandBits = math.floatMantissaBits(T);
15 const fractionalBits = math.floatFractionalBits(T);
16 const exponentBits = math.floatExponentBits(T);
17
18 const signBit = (@as(Z, 1) << (significandBits + exponentBits));
19 const maxExponent = ((1 << exponentBits) - 1);
20
21 var aRep = @bitCast(Z, a);
22 var bRep = @bitCast(Z, b);
23
24 const signA = aRep & signBit;
25 var expA = @intCast(i32, (@bitCast(Z, a) >> significandBits) & maxExponent);
26 var expB = @intCast(i32, (@bitCast(Z, b) >> significandBits) & maxExponent);
27
28 // There are 3 cases where the answer is undefined, check for:
29 // - fmodx(val, 0)
30 // - fmodx(val, NaN)
31 // - fmodx(inf, val)
32 // The sign on checked values does not matter.
33 // Doing (a * b) / (a * b) procudes undefined results
34 // because the three cases always produce undefined calculations:
35 // - 0 / 0
36 // - val * NaN
37 // - inf / inf
38 if (b == 0 or math.isNan(b) or expA == maxExponent) {
39 return (a * b) / (a * b);
40 }
41
42 // Remove the sign from both
43 aRep &= ~signBit;
44 bRep &= ~signBit;
45 if (aRep <= bRep) {
46 if (aRep == bRep) {
47 return 0 * a;
48 }
49 return a;
50 }
51
52 if (expA == 0) expA = normalize(f80, &aRep);
53 if (expB == 0) expB = normalize(f80, &bRep);
54
55 var highA: u64 = 0;
56 var highB: u64 = 0;
57 var lowA: u64 = @truncate(u64, aRep);
58 var lowB: u64 = @truncate(u64, bRep);
59
60 while (expA > expB) : (expA -= 1) {
61 var high = highA -% highB;
62 var low = lowA -% lowB;
63 if (lowA < lowB) {
64 high -%= 1;
65 }
66 if (high >> 63 == 0) {
67 if ((high | low) == 0) {
68 return 0 * a;
69 }
70 highA = 2 *% high + (low >> 63);
71 lowA = 2 *% low;
72 } else {
73 highA = 2 *% highA + (lowA >> 63);
74 lowA = 2 *% lowA;
75 }
76 }
77
78 var high = highA -% highB;
79 var low = lowA -% lowB;
80 if (lowA < lowB) {
81 high -%= 1;
82 }
83 if (high >> 63 == 0) {
84 if ((high | low) == 0) {
85 return 0 * a;
86 }
87 highA = high;
88 lowA = low;
89 }
90
91 while ((lowA >> fractionalBits) == 0) {
92 lowA = 2 *% lowA;
93 expA = expA - 1;
94 }
95
96 // Combine the exponent with the sign and significand, normalize if happened to be denormalized
97 if (expA < -fractionalBits) {
98 return @bitCast(T, signA);
99 } else if (expA <= 0) {
100 return @bitCast(T, (lowA >> @intCast(math.Log2Int(u64), 1 - expA)) | signA);
101 } else {
102 return @bitCast(T, lowA | (@as(Z, @intCast(u16, expA)) << significandBits) | signA);
103 }
104}
105
106test {
107 _ = @import("fmodx_test.zig");
108}
lib/std/special/compiler_rt/fmodx_test.zig created+51
...@@ -0,0 +1,51 @@
1const std = @import("std");
2const fmodx = @import("fmodx.zig");
3const testing = std.testing;
4
5fn test_fmodx(a: f80, b: f80, exp: f80) !void {
6 const res = fmodx.fmodx(a, b);
7 try testing.expect(exp == res);
8}
9
10fn test_fmodx_nans() !void {
11 try testing.expect(std.math.isNan(fmodx.fmodx(1.0, std.math.nan(f80))));
12 try testing.expect(std.math.isNan(fmodx.fmodx(1.0, -std.math.nan(f80))));
13 try testing.expect(std.math.isNan(fmodx.fmodx(std.math.nan(f80), 1.0)));
14 try testing.expect(std.math.isNan(fmodx.fmodx(-std.math.nan(f80), 1.0)));
15}
16
17fn test_fmodx_infs() !void {
18 try testing.expect(fmodx.fmodx(1.0, std.math.inf(f80)) == 1.0);
19 try testing.expect(fmodx.fmodx(1.0, -std.math.inf(f80)) == 1.0);
20 try testing.expect(std.math.isNan(fmodx.fmodx(std.math.inf(f80), 1.0)));
21 try testing.expect(std.math.isNan(fmodx.fmodx(-std.math.inf(f80), 1.0)));
22}
23
24test "fmodx" {
25 try test_fmodx(6.4, 4.0, 2.4);
26 try test_fmodx(6.4, -4.0, 2.4);
27 try test_fmodx(-6.4, 4.0, -2.4);
28 try test_fmodx(-6.4, -4.0, -2.4);
29 try test_fmodx(3.0, 2.0, 1.0);
30 try test_fmodx(-5.0, 3.0, -2.0);
31 try test_fmodx(3.0, 2.0, 1.0);
32 try test_fmodx(1.0, 2.0, 1.0);
33 try test_fmodx(0.0, 1.0, 0.0);
34 try test_fmodx(-0.0, 1.0, -0.0);
35 try test_fmodx(7046119.0, 5558362.0, 1487757.0);
36 try test_fmodx(9010357.0, 1957236.0, 1181413.0);
37 try test_fmodx(9223372036854775807, 10.0, 7.0);
38
39 // Denormals
40 const a1: f80 = 0x0.76e5_9a51_1a92_9ca4p-16381;
41 const b1: f80 = 0x0.2e97_1c3c_8e7d_e03ap-16381;
42 const exp1: f80 = 0x0.19b7_61d7_fd96_dc30p-16381;
43 try test_fmodx(a1, b1, exp1);
44 const a2: f80 = 0x0.76e5_9a51_1a92_9ca4p-16381;
45 const b2: f80 = 0x0.0e97_1c3c_8e7d_e03ap-16381;
46 const exp2: f80 = 0x0.022c_b86c_a6a3_9ad4p-16381;
47 try test_fmodx(a2, b2, exp2);
48
49 try test_fmodx_nans();
50 try test_fmodx_infs();
51}
lib/std/special/compiler_rt/mulXf3.zig+9-1
...@@ -152,6 +152,10 @@ fn mulXf3(comptime T: type, a: T, b: T) T {...@@ -152,6 +152,10 @@ fn mulXf3(comptime T: type, a: T, b: T) T {
152 const sticky = wideShrWithTruncation(ZSignificand, &productHi, &productLo, shift);152 const sticky = wideShrWithTruncation(ZSignificand, &productHi, &productLo, shift);
153 productLo |= @boolToInt(sticky);153 productLo |= @boolToInt(sticky);
154 result = productHi;154 result = productHi;
155
156 // We include the integer bit so that rounding will carry to the exponent,
157 // but it will be removed later if the result is still denormal
158 if (significandBits != fractionalBits) result |= integerBit;
155 } else {159 } else {
156 // Result is normal before rounding; insert the exponent.160 // Result is normal before rounding; insert the exponent.
157 result = productHi & significandMask;161 result = productHi & significandMask;
...@@ -166,7 +170,11 @@ fn mulXf3(comptime T: type, a: T, b: T) T {...@@ -166,7 +170,11 @@ fn mulXf3(comptime T: type, a: T, b: T) T {
166170
167 // Restore any explicit integer bit, if it was rounded off171 // Restore any explicit integer bit, if it was rounded off
168 if (significandBits != fractionalBits) {172 if (significandBits != fractionalBits) {
169 if ((result >> significandBits) != 0) result |= integerBit;173 if ((result >> significandBits) != 0) {
174 result |= integerBit;
175 } else {
176 result &= ~integerBit;
177 }
170 }178 }
171179
172 // Insert the sign of the result:180 // Insert the sign of the result:
lib/std/special/compiler_rt/mulXf3_test.zig+4
...@@ -105,6 +105,7 @@ test "multf3" {...@@ -105,6 +105,7 @@ test "multf3" {
105105
106 try test__multf3(0x1.0000_0000_0000_0000_0000_0000_0001p+0, 0x1.8p+5, 0x4004_8000_0000_0000, 0x0000_0000_0000_0002);106 try test__multf3(0x1.0000_0000_0000_0000_0000_0000_0001p+0, 0x1.8p+5, 0x4004_8000_0000_0000, 0x0000_0000_0000_0002);
107 try test__multf3(0x1.0000_0000_0000_0000_0000_0000_0002p+0, 0x1.8p+5, 0x4004_8000_0000_0000, 0x0000_0000_0000_0003);107 try test__multf3(0x1.0000_0000_0000_0000_0000_0000_0002p+0, 0x1.8p+5, 0x4004_8000_0000_0000, 0x0000_0000_0000_0003);
108 try test__multf3(2.0, math.floatTrueMin(f128), 0x0000_0000_0000_0000, 0x0000_0000_0000_0002);
108}109}
109110
110const qnan80 = @bitCast(f80, @bitCast(u80, math.nan(f80)) | (1 << (math.floatFractionalBits(f80) - 1)));111const qnan80 = @bitCast(f80, @bitCast(u80, math.nan(f80)) | (1 << (math.floatFractionalBits(f80) - 1)));
...@@ -164,4 +165,7 @@ test "mulxf3" {...@@ -164,4 +165,7 @@ test "mulxf3" {
164165
165 try test__mulxf3(0x1.0000_0001p+0, 0x1.0000_0001p+0, 0x3FFF_8000_0001_0000_0000); // round down to even166 try test__mulxf3(0x1.0000_0001p+0, 0x1.0000_0001p+0, 0x3FFF_8000_0001_0000_0000); // round down to even
166 try test__mulxf3(0x1.0000_0001p+0, 0x1.0000_0001_0002p+0, 0x3FFF_8000_0001_0001_0001); // round up167 try test__mulxf3(0x1.0000_0001p+0, 0x1.0000_0001_0002p+0, 0x3FFF_8000_0001_0001_0001); // round up
168 try test__mulxf3(0x0.8000_0000_0000_0000p-16382, 2.0, 0x0001_8000_0000_0000_0000); // denormal -> normal
169 try test__mulxf3(0x0.7fff_ffff_ffff_fffep-16382, 0x2.0000_0000_0000_0008p0, 0x0001_8000_0000_0000_0000); // denormal -> normal
170 try test__mulxf3(0x0.7fff_ffff_ffff_fffep-16382, 0x1.0000_0000_0000_0000p0, 0x0000_3FFF_FFFF_FFFF_FFFF); // denormal -> denormal
167}171}