xref: /freebsd/contrib/llvm-project/compiler-rt/lib/builtins/divdf3.c (revision c66ec88fed842fbaad62c30d510644ceb7bd2d71)
1 //===-- lib/divdf3.c - Double-precision division ------------------*- C -*-===//
2 //
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
6 //
7 //===----------------------------------------------------------------------===//
8 //
9 // This file implements double-precision soft-float division
10 // with the IEEE-754 default rounding (to nearest, ties to even).
11 //
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.
15 //
16 //===----------------------------------------------------------------------===//
17 
18 #define DOUBLE_PRECISION
19 #include "fp_lib.h"
20 
21 COMPILER_RT_ABI fp_t __divdf3(fp_t a, fp_t b) {
22 
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;
26 
27   rep_t aSignificand = toRep(a) & significandMask;
28   rep_t bSignificand = toRep(b) & significandMask;
29   int scale = 0;
30 
31   // Detect if a or b is zero, denormal, infinity, or NaN.
32   if (aExponent - 1U >= maxExponent - 1U ||
33       bExponent - 1U >= maxExponent - 1U) {
34 
35     const rep_t aAbs = toRep(a) & absMask;
36     const rep_t bAbs = toRep(b) & absMask;
37 
38     // NaN / anything = qNaN
39     if (aAbs > infRep)
40       return fromRep(toRep(a) | quietBit);
41     // anything / NaN = qNaN
42     if (bAbs > infRep)
43       return fromRep(toRep(b) | quietBit);
44 
45     if (aAbs == infRep) {
46       // infinity / infinity = NaN
47       if (bAbs == infRep)
48         return fromRep(qnanRep);
49       // infinity / anything else = +/- infinity
50       else
51         return fromRep(aAbs | quotientSign);
52     }
53 
54     // anything else / infinity = +/- 0
55     if (bAbs == infRep)
56       return fromRep(quotientSign);
57 
58     if (!aAbs) {
59       // zero / zero = NaN
60       if (!bAbs)
61         return fromRep(qnanRep);
62       // zero / anything else = +/- zero
63       else
64         return fromRep(quotientSign);
65     }
66     // anything else / zero = +/- infinity
67     if (!bAbs)
68       return fromRep(infRep | quotientSign);
69 
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);
77   }
78 
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 
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)
93 
94   // Now refine the reciprocal estimate using a Newton-Raphson iteration:
95   //
96   //     x1 = x0 * (2 - x0 * b)
97   //
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;
107 
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.
111   recip32--;
112 
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);
121 
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.
126   reciprocal -= 2;
127 
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:
132   //
133   //    1. q < a/b
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).
137 
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, &quotient, &quotientLo);
142 
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
145   //
146   //     r = a - q*b
147   //
148   // We know from the construction of q that r satisfies:
149   //
150   //     0 <= r < ulp(q)*b
151   //
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.
156   rep_t residual;
157   if (quotient < (implicitBit << 1)) {
158     residual = (aSignificand << 53) - quotient * bSignificand;
159     quotientExponent--;
160   } else {
161     quotient >>= 1;
162     residual = (aSignificand << 52) - quotient * bSignificand;
163   }
164 
165   const int writtenExponent = quotientExponent + exponentBias;
166 
167   if (writtenExponent >= maxExponent) {
168     // If we have overflowed the exponent, return infinity.
169     return fromRep(infRep | quotientSign);
170   }
171 
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;
178       // Round.
179       absResult += round;
180       if (absResult & ~significandMask) {
181         // The rounded result is normal; return it.
182         return fromRep(absResult | quotientSign);
183       }
184     }
185     // Flush denormals to zero.  In the future, it would be nice to add
186     // code to round them correctly.
187     return fromRep(quotientSign);
188   }
189 
190   else {
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;
196     // Round.
197     absResult += round;
198     // Insert the sign and return.
199     const double result = fromRep(absResult | quotientSign);
200     return result;
201   }
202 }
203 
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); }
207 #else
208 COMPILER_RT_ALIAS(__divdf3, __aeabi_ddiv)
209 #endif
210 #endif
211