1 //===-- lib/divsf3.c - Single-precision division ------------------*- C -*-===//
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
7 //===----------------------------------------------------------------------===//
9 // This file implements single-precision soft-float division
10 // with the IEEE-754 default rounding (to nearest, ties to even).
12 // For simplicity, this implementation currently flushes denormals to zero.
13 // It should be a fairly straightforward exercise to implement gradual
14 // underflow with correct rounding.
16 //===----------------------------------------------------------------------===//
18 #define SINGLE_PRECISION
21 COMPILER_RT_ABI fp_t __divsf3(fp_t a, fp_t b) {
23 const unsigned int aExponent = toRep(a) >> significandBits & maxExponent;
24 const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
25 const rep_t quotientSign = (toRep(a) ^ toRep(b)) & signBit;
27 rep_t aSignificand = toRep(a) & significandMask;
28 rep_t bSignificand = toRep(b) & significandMask;
31 // Detect if a or b is zero, denormal, infinity, or NaN.
32 if (aExponent - 1U >= maxExponent - 1U ||
33 bExponent - 1U >= maxExponent - 1U) {
35 const rep_t aAbs = toRep(a) & absMask;
36 const rep_t bAbs = toRep(b) & absMask;
38 // NaN / anything = qNaN
40 return fromRep(toRep(a) | quietBit);
41 // anything / NaN = qNaN
43 return fromRep(toRep(b) | quietBit);
46 // infinity / infinity = NaN
48 return fromRep(qnanRep);
49 // infinity / anything else = +/- infinity
51 return fromRep(aAbs | quotientSign);
54 // anything else / infinity = +/- 0
56 return fromRep(quotientSign);
61 return fromRep(qnanRep);
62 // zero / anything else = +/- zero
64 return fromRep(quotientSign);
66 // anything else / zero = +/- infinity
68 return fromRep(infRep | quotientSign);
70 // One or both of a or b is denormal. The other (if applicable) is a
71 // normal number. Renormalize one or both of a and b, and set scale to
72 // include the necessary exponent adjustment.
73 if (aAbs < implicitBit)
74 scale += normalize(&aSignificand);
75 if (bAbs < implicitBit)
76 scale -= normalize(&bSignificand);
79 // Set the implicit significand bit. If we fell through from the
80 // denormal path it was already set by normalize( ), but setting it twice
81 // won't hurt anything.
82 aSignificand |= implicitBit;
83 bSignificand |= implicitBit;
84 int quotientExponent = aExponent - bExponent + scale;
85 // 0x7504F333 / 2^32 + 1 = 3/4 + 1/sqrt(2)
87 // Align the significand of b as a Q31 fixed-point number in the range
88 // [1, 2.0) and get a Q32 approximate reciprocal using a small minimax
89 // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This
90 // is accurate to about 3.5 binary digits.
91 uint32_t q31b = bSignificand << 8;
92 uint32_t reciprocal = UINT32_C(0x7504f333) - q31b;
94 // Now refine the reciprocal estimate using a Newton-Raphson iteration:
96 // x1 = x0 * (2 - x0 * b)
98 // This doubles the number of correct binary digits in the approximation
99 // with each iteration.
101 correction = -((uint64_t)reciprocal * q31b >> 32);
102 reciprocal = (uint64_t)reciprocal * correction >> 31;
103 correction = -((uint64_t)reciprocal * q31b >> 32);
104 reciprocal = (uint64_t)reciprocal * correction >> 31;
105 correction = -((uint64_t)reciprocal * q31b >> 32);
106 reciprocal = (uint64_t)reciprocal * correction >> 31;
108 // Adust the final 32-bit reciprocal estimate downward to ensure that it is
109 // strictly smaller than the infinitely precise exact reciprocal. Because
110 // the computation of the Newton-Raphson step is truncating at every step,
111 // this adjustment is small; most of the work is already done.
114 // The numerical reciprocal is accurate to within 2^-28, lies in the
115 // interval [0x1.000000eep-1, 0x1.fffffffcp-1], and is strictly smaller
116 // than the true reciprocal of b. Multiplying a by this reciprocal thus
117 // gives a numerical q = a/b in Q24 with the following properties:
120 // 2. q is in the interval [0x1.000000eep-1, 0x1.fffffffcp0)
121 // 3. The error in q is at most 2^-24 + 2^-27 -- the 2^24 term comes
122 // from the fact that we truncate the product, and the 2^27 term
123 // is the error in the reciprocal of b scaled by the maximum
124 // possible value of a. As a consequence of this error bound,
125 // either q or nextafter(q) is the correctly rounded.
126 rep_t quotient = (uint64_t)reciprocal * (aSignificand << 1) >> 32;
128 // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).
129 // In either case, we are going to compute a residual of the form
133 // We know from the construction of q that r satisfies:
137 // If r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we
138 // already have the correct result. The exact halfway case cannot occur.
139 // We also take this time to right shift quotient if it falls in the [1,2)
140 // range and adjust the exponent accordingly.
142 if (quotient < (implicitBit << 1)) {
143 residual = (aSignificand << 24) - quotient * bSignificand;
147 residual = (aSignificand << 23) - quotient * bSignificand;
150 const int writtenExponent = quotientExponent + exponentBias;
152 if (writtenExponent >= maxExponent) {
153 // If we have overflowed the exponent, return infinity.
154 return fromRep(infRep | quotientSign);
157 else if (writtenExponent < 1) {
158 if (writtenExponent == 0) {
159 // Check whether the rounded result is normal.
160 const bool round = (residual << 1) > bSignificand;
161 // Clear the implicit bit.
162 rep_t absResult = quotient & significandMask;
165 if (absResult & ~significandMask) {
166 // The rounded result is normal; return it.
167 return fromRep(absResult | quotientSign);
170 // Flush denormals to zero. In the future, it would be nice to add
171 // code to round them correctly.
172 return fromRep(quotientSign);
176 const bool round = (residual << 1) > bSignificand;
177 // Clear the implicit bit.
178 rep_t absResult = quotient & significandMask;
179 // Insert the exponent.
180 absResult |= (rep_t)writtenExponent << significandBits;
183 // Insert the sign and return.
184 return fromRep(absResult | quotientSign);
188 #if defined(__ARM_EABI__)
189 #if defined(COMPILER_RT_ARMHF_TARGET)
190 AEABI_RTABI fp_t __aeabi_fdiv(fp_t a, fp_t b) { return __divsf3(a, b); }
192 COMPILER_RT_ALIAS(__divsf3, __aeabi_fdiv)