1 //===-- lib/divdf3.c - Double-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 double-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 DOUBLE_PRECISION
21 COMPILER_RT_ABI fp_t __divdf3(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;
86 // Align the significand of b as a Q31 fixed-point number in the range
87 // [1, 2.0) and get a Q32 approximate reciprocal using a small minimax
88 // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This
89 // is accurate to about 3.5 binary digits.
90 const uint32_t q31b = bSignificand >> 21;
91 uint32_t recip32 = UINT32_C(0x7504f333) - q31b;
92 // 0x7504F333 / 2^32 + 1 = 3/4 + 1/sqrt(2)
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.
100 uint32_t correction32;
101 correction32 = -((uint64_t)recip32 * q31b >> 32);
102 recip32 = (uint64_t)recip32 * correction32 >> 31;
103 correction32 = -((uint64_t)recip32 * q31b >> 32);
104 recip32 = (uint64_t)recip32 * correction32 >> 31;
105 correction32 = -((uint64_t)recip32 * q31b >> 32);
106 recip32 = (uint64_t)recip32 * correction32 >> 31;
108 // The reciprocal may have overflowed to zero if the upper half of b is
109 // exactly 1.0. This would sabatoge the full-width final stage of the
110 // computation that follows, so we adjust the reciprocal down by one bit.
113 // We need to perform one more iteration to get us to 56 binary digits.
114 // The last iteration needs to happen with extra precision.
115 const uint32_t q63blo = bSignificand << 11;
116 uint64_t correction, reciprocal;
117 correction = -((uint64_t)recip32 * q31b + ((uint64_t)recip32 * q63blo >> 32));
118 uint32_t cHi = correction >> 32;
119 uint32_t cLo = correction;
120 reciprocal = (uint64_t)recip32 * cHi + ((uint64_t)recip32 * cLo >> 32);
122 // Adjust the final 64-bit reciprocal estimate downward to ensure that it is
123 // strictly smaller than the infinitely precise exact reciprocal. Because
124 // the computation of the Newton-Raphson step is truncating at every step,
125 // this adjustment is small; most of the work is already done.
128 // The numerical reciprocal is accurate to within 2^-56, lies in the
129 // interval [0.5, 1.0), and is strictly smaller than the true reciprocal
130 // of b. Multiplying a by this reciprocal thus gives a numerical q = a/b
131 // in Q53 with the following properties:
134 // 2. q is in the interval [0.5, 2.0)
135 // 3. The error in q is bounded away from 2^-53 (actually, we have a
136 // couple of bits to spare, but this is all we need).
138 // We need a 64 x 64 multiply high to compute q, which isn't a basic
139 // operation in C, so we need to be a little bit fussy.
140 rep_t quotient, quotientLo;
141 wideMultiply(aSignificand << 2, reciprocal, "ient, "ientLo);
143 // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).
144 // In either case, we are going to compute a residual of the form
148 // We know from the construction of q that r satisfies:
152 // If r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we
153 // already have the correct result. The exact halfway case cannot occur.
154 // We also take this time to right shift quotient if it falls in the [1,2)
155 // range and adjust the exponent accordingly.
157 if (quotient < (implicitBit << 1)) {
158 residual = (aSignificand << 53) - quotient * bSignificand;
162 residual = (aSignificand << 52) - quotient * bSignificand;
165 const int writtenExponent = quotientExponent + exponentBias;
167 if (writtenExponent >= maxExponent) {
168 // If we have overflowed the exponent, return infinity.
169 return fromRep(infRep | quotientSign);
172 else if (writtenExponent < 1) {
173 if (writtenExponent == 0) {
174 // Check whether the rounded result is normal.
175 const bool round = (residual << 1) > bSignificand;
176 // Clear the implicit bit.
177 rep_t absResult = quotient & significandMask;
180 if (absResult & ~significandMask) {
181 // The rounded result is normal; return it.
182 return fromRep(absResult | quotientSign);
185 // Flush denormals to zero. In the future, it would be nice to add
186 // code to round them correctly.
187 return fromRep(quotientSign);
191 const bool round = (residual << 1) > bSignificand;
192 // Clear the implicit bit.
193 rep_t absResult = quotient & significandMask;
194 // Insert the exponent.
195 absResult |= (rep_t)writtenExponent << significandBits;
198 // Insert the sign and return.
199 const double result = fromRep(absResult | quotientSign);
204 #if defined(__ARM_EABI__)
205 #if defined(COMPILER_RT_ARMHF_TARGET)
206 AEABI_RTABI fp_t __aeabi_ddiv(fp_t a, fp_t b) { return __divdf3(a, b); }
208 COMPILER_RT_ALIAS(__divdf3, __aeabi_ddiv)