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
21COMPILER_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)
190AEABI_RTABI fp_t __aeabi_fdiv(fp_t a, fp_t b) { return __divsf3(a, b); }
191#else
192COMPILER_RT_ALIAS(__divsf3, __aeabi_fdiv)
193#endif
194#endif
195