1const std = @import("std");
2const builtin = @import("builtin");
3const arch = builtin.cpu.arch;
4
5const compiler_rt = @import("../compiler_rt.zig");
6const symbol = compiler_rt.symbol;
7const normalize = compiler_rt.normalize;
8const wideMultiply = compiler_rt.wideMultiply;
9
10comptime {
11 symbol(&__divxf3, "__divxf3");
12}
13
14fn __divxf3(a: compiler_rt.f80.Abi, b: compiler_rt.f80.Abi) callconv(.c) compiler_rt.f80.Abi {
15 return compiler_rt.f80.toAbi(div_f80(compiler_rt.f80.fromAbi(a), compiler_rt.f80.fromAbi(b)));
16}
17pub fn div_f80(a: f80, b: f80) f80 {
18 const T = f80;
19 const Z = @Int(.unsigned, @bitSizeOf(T));
20
21 const significandBits = std.math.floatMantissaBits(T);
22 const fractionalBits = std.math.floatFractionalBits(T);
23 const exponentBits = std.math.floatExponentBits(T);
24
25 const signBit = (@as(Z, 1) << (significandBits + exponentBits));
26 const maxExponent = ((1 << exponentBits) - 1);
27 const exponentBias = (maxExponent >> 1);
28
29 const integerBit = (@as(Z, 1) << fractionalBits);
30 const quietBit = integerBit >> 1;
31 const significandMask = (@as(Z, 1) << significandBits) - 1;
32
33 const absMask = signBit - 1;
34 const qnanRep = @as(Z, @bitCast(std.math.nan(T))) | quietBit;
35 const infRep: Z = @bitCast(std.math.inf(T));
36
37 const aExponent: u32 = @truncate((@as(Z, @bitCast(a)) >> significandBits) & maxExponent);
38 const bExponent: u32 = @truncate((@as(Z, @bitCast(b)) >> significandBits) & maxExponent);
39 const quotientSign: Z = (@as(Z, @bitCast(a)) ^ @as(Z, @bitCast(b))) & signBit;
40
41 var aSignificand: Z = @as(Z, @bitCast(a)) & significandMask;
42 var bSignificand: Z = @as(Z, @bitCast(b)) & significandMask;
43 var scale: i32 = 0;
44
45 // Detect if a or b is zero, denormal, infinity, or NaN.
46 if (aExponent -% 1 >= maxExponent - 1 or bExponent -% 1 >= maxExponent - 1) {
47 const aAbs: Z = @as(Z, @bitCast(a)) & absMask;
48 const bAbs: Z = @as(Z, @bitCast(b)) & absMask;
49
50 // NaN / anything = qNaN
51 if (aAbs > infRep) return @bitCast(@as(Z, @bitCast(a)) | quietBit);
52 // anything / NaN = qNaN
53 if (bAbs > infRep) return @bitCast(@as(Z, @bitCast(b)) | quietBit);
54
55 if (aAbs == infRep) {
56 // infinity / infinity = NaN
57 if (bAbs == infRep) {
58 return @bitCast(qnanRep);
59 }
60 // infinity / anything else = +/- infinity
61 else {
62 return @bitCast(aAbs | quotientSign);
63 }
64 }
65
66 // anything else / infinity = +/- 0
67 if (bAbs == infRep) return @bitCast(quotientSign);
68
69 if (aAbs == 0) {
70 // zero / zero = NaN
71 if (bAbs == 0) {
72 return @bitCast(qnanRep);
73 }
74 // zero / anything else = +/- zero
75 else {
76 return @bitCast(quotientSign);
77 }
78 }
79 // anything else / zero = +/- infinity
80 if (bAbs == 0) return @bitCast(infRep | quotientSign);
81
82 // one or both of a or b is denormal, the other (if applicable) is a
83 // normal number. Renormalize one or both of a and b, and set scale to
84 // include the necessary exponent adjustment.
85 if (aAbs < integerBit) scale +%= normalize(T, &aSignificand);
86 if (bAbs < integerBit) scale -%= normalize(T, &bSignificand);
87 }
88 var quotientExponent: i32 = @as(i32, @bitCast(aExponent -% bExponent)) +% scale;
89
90 // Align the significand of b as a Q63 fixed-point number in the range
91 // [1, 2.0) and get a Q64 approximate reciprocal using a small minimax
92 // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This
93 // is accurate to about 3.5 binary digits.
94 const q63b: u64 = @intCast(bSignificand);
95 var recip64 = @as(u64, 0x7504f333F9DE6484) -% q63b;
96 // 0x7504f333F9DE6484 / 2^64 + 1 = 3/4 + 1/sqrt(2)
97
98 // Now refine the reciprocal estimate using a Newton-Raphson iteration:
99 //
100 // x1 = x0 * (2 - x0 * b)
101 //
102 // This doubles the number of correct binary digits in the approximation
103 // with each iteration.
104 var correction64: u64 = undefined;
105 correction64 = @truncate(~(@as(u128, recip64) *% q63b >> 64) +% 1);
106 recip64 = @truncate(@as(u128, recip64) *% correction64 >> 63);
107 correction64 = @truncate(~(@as(u128, recip64) *% q63b >> 64) +% 1);
108 recip64 = @truncate(@as(u128, recip64) *% correction64 >> 63);
109 correction64 = @truncate(~(@as(u128, recip64) *% q63b >> 64) +% 1);
110 recip64 = @truncate(@as(u128, recip64) *% correction64 >> 63);
111 correction64 = @truncate(~(@as(u128, recip64) *% q63b >> 64) +% 1);
112 recip64 = @truncate(@as(u128, recip64) *% correction64 >> 63);
113 correction64 = @truncate(~(@as(u128, recip64) *% q63b >> 64) +% 1);
114 recip64 = @truncate(@as(u128, recip64) *% correction64 >> 63);
115
116 // The reciprocal may have overflowed to zero if the upper half of b is
117 // exactly 1.0. This would sabatoge the full-width final stage of the
118 // computation that follows, so we adjust the reciprocal down by one bit.
119 recip64 -%= 1;
120
121 // We need to perform one more iteration to get us to 112 binary digits;
122 // The last iteration needs to happen with extra precision.
123
124 // NOTE: This operation is equivalent to __multi3, which is not implemented
125 // in some architechures
126 var reciprocal: u128 = undefined;
127 var correction: u128 = undefined;
128 var dummy: u128 = undefined;
129 wideMultiply(u128, recip64, q63b, &dummy, &correction);
130
131 correction = -%correction;
132
133 const cHi: u64 = @truncate(correction >> 64);
134 const cLo: u64 = @truncate(correction);
135
136 var r64cH: u128 = undefined;
137 var r64cL: u128 = undefined;
138 wideMultiply(u128, recip64, cHi, &dummy, &r64cH);
139 wideMultiply(u128, recip64, cLo, &dummy, &r64cL);
140
141 reciprocal = r64cH + (r64cL >> 64);
142
143 // Adjust the final 128-bit reciprocal estimate downward to ensure that it
144 // is strictly smaller than the infinitely precise exact reciprocal. Because
145 // the computation of the Newton-Raphson step is truncating at every step,
146 // this adjustment is small; most of the work is already done.
147 reciprocal -%= 2;
148
149 // The numerical reciprocal is accurate to within 2^-112, lies in the
150 // interval [0.5, 1.0), and is strictly smaller than the true reciprocal
151 // of b. Multiplying a by this reciprocal thus gives a numerical q = a/b
152 // in Q127 with the following properties:
153 //
154 // 1. q < a/b
155 // 2. q is in the interval [0.5, 2.0)
156 // 3. The error in q is bounded away from 2^-63 (actually, we have
157 // many bits to spare, but this is all we need).
158
159 // We need a 128 x 128 multiply high to compute q.
160 var quotient128: u128 = undefined;
161 var quotientLo: u128 = undefined;
162 wideMultiply(u128, aSignificand << 2, reciprocal, &quotient128, &quotientLo);
163
164 // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).
165 // Right shift the quotient if it falls in the [1,2) range and adjust the
166 // exponent accordingly.
167 const quotient: u64 = if (quotient128 < (integerBit << 1)) b: {
168 quotientExponent -= 1;
169 break :b @intCast(quotient128);
170 } else @intCast(quotient128 >> 1);
171
172 // We are going to compute a residual of the form
173 //
174 // r = a - q*b
175 //
176 // We know from the construction of q that r satisfies:
177 //
178 // 0 <= r < ulp(q)*b
179 //
180 // If r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we
181 // already have the correct result. The exact halfway case cannot occur.
182 const residual: u64 = -%(quotient *% q63b);
183
184 const writtenExponent = quotientExponent + exponentBias;
185 if (writtenExponent >= maxExponent) {
186 // If we have overflowed the exponent, return infinity.
187 return @bitCast(infRep | quotientSign);
188 } else if (writtenExponent < 1) {
189 if (writtenExponent == 0) {
190 // Check whether the rounded result is normal.
191 if (residual > (bSignificand >> 1)) { // round
192 if (quotient == (integerBit - 1)) // If the rounded result is normal, return it
193 return @bitCast(@as(Z, @bitCast(std.math.floatMin(T))) | quotientSign);
194 }
195 }
196 // Flush denormals to zero. In the future, it would be nice to add
197 // code to round them correctly.
198 return @bitCast(quotientSign);
199 } else {
200 const round = @intFromBool(residual > (bSignificand >> 1));
201 // Insert the exponent
202 var absResult = quotient | (@as(Z, @intCast(writtenExponent)) << significandBits);
203 // Round
204 absResult +%= round;
205 // Insert the sign and return
206 return @bitCast(absResult | quotientSign | integerBit);
207 }
208}
209
210test {
211 _ = @import("divxf3_test.zig");
212}