1 //===-- lib/divsf3.c - Single-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 single-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 SINGLE_PRECISION
19 #include "fp_lib.h"
20 
21 COMPILER_RT_ABI fp_t __divsf3(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   // 0x7504F333 / 2^32 + 1 = 3/4 + 1/sqrt(2)
86 
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;
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 correction;
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;
107 
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.
112   reciprocal -= 2;
113 
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:
118   //
119   //    1. q < a/b
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;
127 
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
130   //
131   //     r = a - q*b
132   //
133   // We know from the construction of q that r satisfies:
134   //
135   //     0 <= r < ulp(q)*b
136   //
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.
141   rep_t residual;
142   if (quotient < (implicitBit << 1)) {
143     residual = (aSignificand << 24) - quotient * bSignificand;
144     quotientExponent--;
145   } else {
146     quotient >>= 1;
147     residual = (aSignificand << 23) - quotient * bSignificand;
148   }
149 
150   const int writtenExponent = quotientExponent + exponentBias;
151 
152   if (writtenExponent >= maxExponent) {
153     // If we have overflowed the exponent, return infinity.
154     return fromRep(infRep | quotientSign);
155   }
156 
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;
163       // Round.
164       absResult += round;
165       if (absResult & ~significandMask) {
166         // The rounded result is normal; return it.
167         return fromRep(absResult | quotientSign);
168       }
169     }
170     // Flush denormals to zero.  In the future, it would be nice to add
171     // code to round them correctly.
172     return fromRep(quotientSign);
173   }
174 
175   else {
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;
181     // Round.
182     absResult += round;
183     // Insert the sign and return.
184     return fromRep(absResult | quotientSign);
185   }
186 }
187 
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); }
191 #else
192 COMPILER_RT_ALIAS(__divsf3, __aeabi_fdiv)
193 #endif
194 #endif
195