1//! Ported from:
2//!
3//! https://github.com/llvm/llvm-project/commit/d674d96bc56c0f377879d01c9d8dfdaaa7859cdb/compiler-rt/lib/builtins/divsf3.c
4
5const std = @import("std");
6
7const compiler_rt = @import("../compiler_rt.zig");
8const symbol = compiler_rt.symbol;
9const normalize = compiler_rt.normalize;
10
11comptime {
12 symbol(&__divhf3, "__divhf3");
13 if (compiler_rt.want_aeabi) {
14 symbol(&__aeabi_fdiv, "__aeabi_fdiv");
15 } else {
16 symbol(&__divsf3, "__divsf3");
17 }
18}
19
20fn __divhf3(a: compiler_rt.f16.Abi, b: compiler_rt.f16.Abi) callconv(.c) compiler_rt.f16.Abi {
21 return compiler_rt.f16.toAbi(div_f16(compiler_rt.f16.fromAbi(a), compiler_rt.f16.fromAbi(b)));
22}
23pub fn div_f16(a: f16, b: f16) f16 {
24 // TODO: more efficient implementation
25 return @floatCast(div_f32(a, b));
26}
27
28fn __divsf3(a: compiler_rt.f32.Abi, b: compiler_rt.f32.Abi) callconv(.c) compiler_rt.f32.Abi {
29 return compiler_rt.f32.toAbi(div_f32(compiler_rt.f32.fromAbi(a), compiler_rt.f32.fromAbi(b)));
30}
31
32fn __aeabi_fdiv(a: f32, b: f32) callconv(.{ .arm_aapcs = .{} }) f32 {
33 return div_f32(a, b);
34}
35
36pub fn div_f32(a: f32, b: f32) f32 {
37 const Z = @Int(.unsigned, 32);
38
39 const significandBits = std.math.floatMantissaBits(f32);
40 const exponentBits = std.math.floatExponentBits(f32);
41
42 const signBit = (@as(Z, 1) << (significandBits + exponentBits));
43 const maxExponent = ((1 << exponentBits) - 1);
44 const exponentBias = (maxExponent >> 1);
45
46 const implicitBit = (@as(Z, 1) << significandBits);
47 const quietBit = implicitBit >> 1;
48 const significandMask = implicitBit - 1;
49
50 const absMask = signBit - 1;
51 const exponentMask = absMask ^ significandMask;
52 const qnanRep = exponentMask | quietBit;
53 const infRep: Z = @bitCast(std.math.inf(f32));
54
55 const aExponent: u32 = @truncate((@as(Z, @bitCast(a)) >> significandBits) & maxExponent);
56 const bExponent: u32 = @truncate((@as(Z, @bitCast(b)) >> significandBits) & maxExponent);
57 const quotientSign: Z = (@as(Z, @bitCast(a)) ^ @as(Z, @bitCast(b))) & signBit;
58
59 var aSignificand: Z = @as(Z, @bitCast(a)) & significandMask;
60 var bSignificand: Z = @as(Z, @bitCast(b)) & significandMask;
61 var scale: i32 = 0;
62
63 // Detect if a or b is zero, denormal, infinity, or NaN.
64 if (aExponent -% 1 >= maxExponent - 1 or bExponent -% 1 >= maxExponent - 1) {
65 const aAbs: Z = @as(Z, @bitCast(a)) & absMask;
66 const bAbs: Z = @as(Z, @bitCast(b)) & absMask;
67
68 // NaN / anything = qNaN
69 if (aAbs > infRep) return @bitCast(@as(Z, @bitCast(a)) | quietBit);
70 // anything / NaN = qNaN
71 if (bAbs > infRep) return @bitCast(@as(Z, @bitCast(b)) | quietBit);
72
73 if (aAbs == infRep) {
74 // infinity / infinity = NaN
75 if (bAbs == infRep) {
76 return @bitCast(qnanRep);
77 }
78 // infinity / anything else = +/- infinity
79 else {
80 return @bitCast(aAbs | quotientSign);
81 }
82 }
83
84 // anything else / infinity = +/- 0
85 if (bAbs == infRep) return @bitCast(quotientSign);
86
87 if (aAbs == 0) {
88 // zero / zero = NaN
89 if (bAbs == 0) {
90 return @bitCast(qnanRep);
91 }
92 // zero / anything else = +/- zero
93 else {
94 return @bitCast(quotientSign);
95 }
96 }
97 // anything else / zero = +/- infinity
98 if (bAbs == 0) return @bitCast(infRep | quotientSign);
99
100 // one or both of a or b is denormal, the other (if applicable) is a
101 // normal number. Renormalize one or both of a and b, and set scale to
102 // include the necessary exponent adjustment.
103 if (aAbs < implicitBit) scale +%= normalize(f32, &aSignificand);
104 if (bAbs < implicitBit) scale -%= normalize(f32, &bSignificand);
105 }
106
107 // Or in the implicit significand bit. (If we fell through from the
108 // denormal path it was already set by normalize( ), but setting it twice
109 // won't hurt anything.)
110 aSignificand |= implicitBit;
111 bSignificand |= implicitBit;
112 var quotientExponent: i32 = @as(i32, @bitCast(aExponent -% bExponent)) +% scale;
113
114 // Align the significand of b as a Q31 fixed-point number in the range
115 // [1, 2.0) and get a Q32 approximate reciprocal using a small minimax
116 // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This
117 // is accurate to about 3.5 binary digits.
118 const q31b = bSignificand << 8;
119 var reciprocal = @as(u32, 0x7504f333) -% q31b;
120
121 // Now refine the reciprocal estimate using a Newton-Raphson iteration:
122 //
123 // x1 = x0 * (2 - x0 * b)
124 //
125 // This doubles the number of correct binary digits in the approximation
126 // with each iteration, so after three iterations, we have about 28 binary
127 // digits of accuracy.
128 var correction: u32 = undefined;
129 correction = @truncate(~(@as(u64, reciprocal) *% q31b >> 32) +% 1);
130 reciprocal = @truncate(@as(u64, reciprocal) *% correction >> 31);
131 correction = @truncate(~(@as(u64, reciprocal) *% q31b >> 32) +% 1);
132 reciprocal = @truncate(@as(u64, reciprocal) *% correction >> 31);
133 correction = @truncate(~(@as(u64, reciprocal) *% q31b >> 32) +% 1);
134 reciprocal = @truncate(@as(u64, reciprocal) *% correction >> 31);
135
136 // Exhaustive testing shows that the error in reciprocal after three steps
137 // is in the interval [-0x1.f58108p-31, 0x1.d0e48cp-29], in line with our
138 // expectations. We bump the reciprocal by a tiny value to force the error
139 // to be strictly positive (in the range [0x1.4fdfp-37,0x1.287246p-29], to
140 // be specific). This also causes 1/1 to give a sensible approximation
141 // instead of zero (due to overflow).
142 reciprocal -%= 2;
143
144 // The numerical reciprocal is accurate to within 2^-28, lies in the
145 // interval [0x1.000000eep-1, 0x1.fffffffcp-1], and is strictly smaller
146 // than the true reciprocal of b. Multiplying a by this reciprocal thus
147 // gives a numerical q = a/b in Q24 with the following properties:
148 //
149 // 1. q < a/b
150 // 2. q is in the interval [0x1.000000eep-1, 0x1.fffffffcp0)
151 // 3. the error in q is at most 2^-24 + 2^-27 -- the 2^24 term comes
152 // from the fact that we truncate the product, and the 2^27 term
153 // is the error in the reciprocal of b scaled by the maximum
154 // possible value of a. As a consequence of this error bound,
155 // either q or nextafter(q) is the correctly rounded
156 var quotient: Z = @truncate(@as(u64, reciprocal) *% (aSignificand << 1) >> 32);
157
158 // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).
159 // In either case, we are going to compute a residual of the form
160 //
161 // r = a - q*b
162 //
163 // We know from the construction of q that r satisfies:
164 //
165 // 0 <= r < ulp(q)*b
166 //
167 // if r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we
168 // already have the correct result. The exact halfway case cannot occur.
169 // We also take this time to right shift quotient if it falls in the [1,2)
170 // range and adjust the exponent accordingly.
171 var residual: Z = undefined;
172 if (quotient < (implicitBit << 1)) {
173 residual = (aSignificand << 24) -% quotient *% bSignificand;
174 quotientExponent -%= 1;
175 } else {
176 quotient >>= 1;
177 residual = (aSignificand << 23) -% quotient *% bSignificand;
178 }
179
180 const writtenExponent = quotientExponent +% exponentBias;
181
182 const round = @intFromBool((residual << 1) >= bSignificand);
183
184 if (writtenExponent >= maxExponent) {
185 // If we have overflowed the exponent, return infinity.
186 return @bitCast(infRep | quotientSign);
187 } else if (writtenExponent < 1) {
188 if (writtenExponent == 0) {
189 // Check whether the rounded result is normal.
190 var absResult = quotient & significandMask;
191 // Round.
192 absResult += round;
193 if ((absResult & ~significandMask) > 0) {
194 // The rounded result is normal; return it.
195 return @bitCast(absResult | quotientSign);
196 }
197 }
198
199 const roundedQuotient = quotient +% round;
200 const shiftAmount: u32 = @intCast(1 - writtenExponent);
201 if (shiftAmount > significandBits + 1) {
202 return @bitCast(quotientSign);
203 }
204
205 const denormQuotient = roundedQuotient >> @as(std.math.Log2Int(Z), @intCast(shiftAmount));
206 return @bitCast((denormQuotient & significandMask) | quotientSign);
207 } else {
208 // Clear the implicit bit
209 var absResult = quotient & significandMask;
210 // Insert the exponent
211 absResult |= @as(Z, @bitCast(writtenExponent)) << significandBits;
212 // Round
213 absResult +%= round;
214 // Insert the sign and return
215 return @bitCast(absResult | quotientSign);
216 }
217}
218
219test {
220 _ = @import("divsf3_test.zig");
221}