1*e8d8bef9SDimitry Andric//===-- fp_div_impl.inc - Floating point division -----------------*- C -*-===// 2*e8d8bef9SDimitry Andric// 3*e8d8bef9SDimitry Andric// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. 4*e8d8bef9SDimitry Andric// See https://llvm.org/LICENSE.txt for license information. 5*e8d8bef9SDimitry Andric// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception 6*e8d8bef9SDimitry Andric// 7*e8d8bef9SDimitry Andric//===----------------------------------------------------------------------===// 8*e8d8bef9SDimitry Andric// 9*e8d8bef9SDimitry Andric// This file implements soft-float division with the IEEE-754 default 10*e8d8bef9SDimitry Andric// rounding (to nearest, ties to even). 11*e8d8bef9SDimitry Andric// 12*e8d8bef9SDimitry Andric//===----------------------------------------------------------------------===// 13*e8d8bef9SDimitry Andric 14*e8d8bef9SDimitry Andric#include "fp_lib.h" 15*e8d8bef9SDimitry Andric 16*e8d8bef9SDimitry Andric// The __divXf3__ function implements Newton-Raphson floating point division. 17*e8d8bef9SDimitry Andric// It uses 3 iterations for float32, 4 for float64 and 5 for float128, 18*e8d8bef9SDimitry Andric// respectively. Due to number of significant bits being roughly doubled 19*e8d8bef9SDimitry Andric// every iteration, the two modes are supported: N full-width iterations (as 20*e8d8bef9SDimitry Andric// it is done for float32 by default) and (N-1) half-width iteration plus one 21*e8d8bef9SDimitry Andric// final full-width iteration. It is expected that half-width integer 22*e8d8bef9SDimitry Andric// operations (w.r.t rep_t size) can be performed faster for some hardware but 23*e8d8bef9SDimitry Andric// they require error estimations to be computed separately due to larger 24*e8d8bef9SDimitry Andric// computational errors caused by truncating intermediate results. 25*e8d8bef9SDimitry Andric 26*e8d8bef9SDimitry Andric// Half the bit-size of rep_t 27*e8d8bef9SDimitry Andric#define HW (typeWidth / 2) 28*e8d8bef9SDimitry Andric// rep_t-sized bitmask with lower half of bits set to ones 29*e8d8bef9SDimitry Andric#define loMask (REP_C(-1) >> HW) 30*e8d8bef9SDimitry Andric 31*e8d8bef9SDimitry Andric#if NUMBER_OF_FULL_ITERATIONS < 1 32*e8d8bef9SDimitry Andric#error At least one full iteration is required 33*e8d8bef9SDimitry Andric#endif 34*e8d8bef9SDimitry Andric 35*e8d8bef9SDimitry Andricstatic __inline fp_t __divXf3__(fp_t a, fp_t b) { 36*e8d8bef9SDimitry Andric 37*e8d8bef9SDimitry Andric const unsigned int aExponent = toRep(a) >> significandBits & maxExponent; 38*e8d8bef9SDimitry Andric const unsigned int bExponent = toRep(b) >> significandBits & maxExponent; 39*e8d8bef9SDimitry Andric const rep_t quotientSign = (toRep(a) ^ toRep(b)) & signBit; 40*e8d8bef9SDimitry Andric 41*e8d8bef9SDimitry Andric rep_t aSignificand = toRep(a) & significandMask; 42*e8d8bef9SDimitry Andric rep_t bSignificand = toRep(b) & significandMask; 43*e8d8bef9SDimitry Andric int scale = 0; 44*e8d8bef9SDimitry Andric 45*e8d8bef9SDimitry Andric // Detect if a or b is zero, denormal, infinity, or NaN. 46*e8d8bef9SDimitry Andric if (aExponent - 1U >= maxExponent - 1U || 47*e8d8bef9SDimitry Andric bExponent - 1U >= maxExponent - 1U) { 48*e8d8bef9SDimitry Andric 49*e8d8bef9SDimitry Andric const rep_t aAbs = toRep(a) & absMask; 50*e8d8bef9SDimitry Andric const rep_t bAbs = toRep(b) & absMask; 51*e8d8bef9SDimitry Andric 52*e8d8bef9SDimitry Andric // NaN / anything = qNaN 53*e8d8bef9SDimitry Andric if (aAbs > infRep) 54*e8d8bef9SDimitry Andric return fromRep(toRep(a) | quietBit); 55*e8d8bef9SDimitry Andric // anything / NaN = qNaN 56*e8d8bef9SDimitry Andric if (bAbs > infRep) 57*e8d8bef9SDimitry Andric return fromRep(toRep(b) | quietBit); 58*e8d8bef9SDimitry Andric 59*e8d8bef9SDimitry Andric if (aAbs == infRep) { 60*e8d8bef9SDimitry Andric // infinity / infinity = NaN 61*e8d8bef9SDimitry Andric if (bAbs == infRep) 62*e8d8bef9SDimitry Andric return fromRep(qnanRep); 63*e8d8bef9SDimitry Andric // infinity / anything else = +/- infinity 64*e8d8bef9SDimitry Andric else 65*e8d8bef9SDimitry Andric return fromRep(aAbs | quotientSign); 66*e8d8bef9SDimitry Andric } 67*e8d8bef9SDimitry Andric 68*e8d8bef9SDimitry Andric // anything else / infinity = +/- 0 69*e8d8bef9SDimitry Andric if (bAbs == infRep) 70*e8d8bef9SDimitry Andric return fromRep(quotientSign); 71*e8d8bef9SDimitry Andric 72*e8d8bef9SDimitry Andric if (!aAbs) { 73*e8d8bef9SDimitry Andric // zero / zero = NaN 74*e8d8bef9SDimitry Andric if (!bAbs) 75*e8d8bef9SDimitry Andric return fromRep(qnanRep); 76*e8d8bef9SDimitry Andric // zero / anything else = +/- zero 77*e8d8bef9SDimitry Andric else 78*e8d8bef9SDimitry Andric return fromRep(quotientSign); 79*e8d8bef9SDimitry Andric } 80*e8d8bef9SDimitry Andric // anything else / zero = +/- infinity 81*e8d8bef9SDimitry Andric if (!bAbs) 82*e8d8bef9SDimitry Andric return fromRep(infRep | quotientSign); 83*e8d8bef9SDimitry Andric 84*e8d8bef9SDimitry Andric // One or both of a or b is denormal. The other (if applicable) is a 85*e8d8bef9SDimitry Andric // normal number. Renormalize one or both of a and b, and set scale to 86*e8d8bef9SDimitry Andric // include the necessary exponent adjustment. 87*e8d8bef9SDimitry Andric if (aAbs < implicitBit) 88*e8d8bef9SDimitry Andric scale += normalize(&aSignificand); 89*e8d8bef9SDimitry Andric if (bAbs < implicitBit) 90*e8d8bef9SDimitry Andric scale -= normalize(&bSignificand); 91*e8d8bef9SDimitry Andric } 92*e8d8bef9SDimitry Andric 93*e8d8bef9SDimitry Andric // Set the implicit significand bit. If we fell through from the 94*e8d8bef9SDimitry Andric // denormal path it was already set by normalize( ), but setting it twice 95*e8d8bef9SDimitry Andric // won't hurt anything. 96*e8d8bef9SDimitry Andric aSignificand |= implicitBit; 97*e8d8bef9SDimitry Andric bSignificand |= implicitBit; 98*e8d8bef9SDimitry Andric 99*e8d8bef9SDimitry Andric int writtenExponent = (aExponent - bExponent + scale) + exponentBias; 100*e8d8bef9SDimitry Andric 101*e8d8bef9SDimitry Andric const rep_t b_UQ1 = bSignificand << (typeWidth - significandBits - 1); 102*e8d8bef9SDimitry Andric 103*e8d8bef9SDimitry Andric // Align the significand of b as a UQ1.(n-1) fixed-point number in the range 104*e8d8bef9SDimitry Andric // [1.0, 2.0) and get a UQ0.n approximate reciprocal using a small minimax 105*e8d8bef9SDimitry Andric // polynomial approximation: x0 = 3/4 + 1/sqrt(2) - b/2. 106*e8d8bef9SDimitry Andric // The max error for this approximation is achieved at endpoints, so 107*e8d8bef9SDimitry Andric // abs(x0(b) - 1/b) <= abs(x0(1) - 1/1) = 3/4 - 1/sqrt(2) = 0.04289..., 108*e8d8bef9SDimitry Andric // which is about 4.5 bits. 109*e8d8bef9SDimitry Andric // The initial approximation is between x0(1.0) = 0.9571... and x0(2.0) = 0.4571... 110*e8d8bef9SDimitry Andric 111*e8d8bef9SDimitry Andric // Then, refine the reciprocal estimate using a quadratically converging 112*e8d8bef9SDimitry Andric // Newton-Raphson iteration: 113*e8d8bef9SDimitry Andric // x_{n+1} = x_n * (2 - x_n * b) 114*e8d8bef9SDimitry Andric // 115*e8d8bef9SDimitry Andric // Let b be the original divisor considered "in infinite precision" and 116*e8d8bef9SDimitry Andric // obtained from IEEE754 representation of function argument (with the 117*e8d8bef9SDimitry Andric // implicit bit set). Corresponds to rep_t-sized b_UQ1 represented in 118*e8d8bef9SDimitry Andric // UQ1.(W-1). 119*e8d8bef9SDimitry Andric // 120*e8d8bef9SDimitry Andric // Let b_hw be an infinitely precise number obtained from the highest (HW-1) 121*e8d8bef9SDimitry Andric // bits of divisor significand (with the implicit bit set). Corresponds to 122*e8d8bef9SDimitry Andric // half_rep_t-sized b_UQ1_hw represented in UQ1.(HW-1) that is a **truncated** 123*e8d8bef9SDimitry Andric // version of b_UQ1. 124*e8d8bef9SDimitry Andric // 125*e8d8bef9SDimitry Andric // Let e_n := x_n - 1/b_hw 126*e8d8bef9SDimitry Andric // E_n := x_n - 1/b 127*e8d8bef9SDimitry Andric // abs(E_n) <= abs(e_n) + (1/b_hw - 1/b) 128*e8d8bef9SDimitry Andric // = abs(e_n) + (b - b_hw) / (b*b_hw) 129*e8d8bef9SDimitry Andric // <= abs(e_n) + 2 * 2^-HW 130*e8d8bef9SDimitry Andric 131*e8d8bef9SDimitry Andric // rep_t-sized iterations may be slower than the corresponding half-width 132*e8d8bef9SDimitry Andric // variant depending on the handware and whether single/double/quad precision 133*e8d8bef9SDimitry Andric // is selected. 134*e8d8bef9SDimitry Andric // NB: Using half-width iterations increases computation errors due to 135*e8d8bef9SDimitry Andric // rounding, so error estimations have to be computed taking the selected 136*e8d8bef9SDimitry Andric // mode into account! 137*e8d8bef9SDimitry Andric#if NUMBER_OF_HALF_ITERATIONS > 0 138*e8d8bef9SDimitry Andric // Starting with (n-1) half-width iterations 139*e8d8bef9SDimitry Andric const half_rep_t b_UQ1_hw = bSignificand >> (significandBits + 1 - HW); 140*e8d8bef9SDimitry Andric 141*e8d8bef9SDimitry Andric // C is (3/4 + 1/sqrt(2)) - 1 truncated to W0 fractional bits as UQ0.HW 142*e8d8bef9SDimitry Andric // with W0 being either 16 or 32 and W0 <= HW. 143*e8d8bef9SDimitry Andric // That is, C is the aforementioned 3/4 + 1/sqrt(2) constant (from which 144*e8d8bef9SDimitry Andric // b/2 is subtracted to obtain x0) wrapped to [0, 1) range. 145*e8d8bef9SDimitry Andric#if defined(SINGLE_PRECISION) 146*e8d8bef9SDimitry Andric // Use 16-bit initial estimation in case we are using half-width iterations 147*e8d8bef9SDimitry Andric // for float32 division. This is expected to be useful for some 16-bit 148*e8d8bef9SDimitry Andric // targets. Not used by default as it requires performing more work during 149*e8d8bef9SDimitry Andric // rounding and would hardly help on regular 32- or 64-bit targets. 150*e8d8bef9SDimitry Andric const half_rep_t C_hw = HALF_REP_C(0x7504); 151*e8d8bef9SDimitry Andric#else 152*e8d8bef9SDimitry Andric // HW is at least 32. Shifting into the highest bits if needed. 153*e8d8bef9SDimitry Andric const half_rep_t C_hw = HALF_REP_C(0x7504F333) << (HW - 32); 154*e8d8bef9SDimitry Andric#endif 155*e8d8bef9SDimitry Andric 156*e8d8bef9SDimitry Andric // b >= 1, thus an upper bound for 3/4 + 1/sqrt(2) - b/2 is about 0.9572, 157*e8d8bef9SDimitry Andric // so x0 fits to UQ0.HW without wrapping. 158*e8d8bef9SDimitry Andric half_rep_t x_UQ0_hw = C_hw - (b_UQ1_hw /* exact b_hw/2 as UQ0.HW */); 159*e8d8bef9SDimitry Andric // An e_0 error is comprised of errors due to 160*e8d8bef9SDimitry Andric // * x0 being an inherently imprecise first approximation of 1/b_hw 161*e8d8bef9SDimitry Andric // * C_hw being some (irrational) number **truncated** to W0 bits 162*e8d8bef9SDimitry Andric // Please note that e_0 is calculated against the infinitely precise 163*e8d8bef9SDimitry Andric // reciprocal of b_hw (that is, **truncated** version of b). 164*e8d8bef9SDimitry Andric // 165*e8d8bef9SDimitry Andric // e_0 <= 3/4 - 1/sqrt(2) + 2^-W0 166*e8d8bef9SDimitry Andric 167*e8d8bef9SDimitry Andric // By construction, 1 <= b < 2 168*e8d8bef9SDimitry Andric // f(x) = x * (2 - b*x) = 2*x - b*x^2 169*e8d8bef9SDimitry Andric // f'(x) = 2 * (1 - b*x) 170*e8d8bef9SDimitry Andric // 171*e8d8bef9SDimitry Andric // On the [0, 1] interval, f(0) = 0, 172*e8d8bef9SDimitry Andric // then it increses until f(1/b) = 1 / b, maximum on (0, 1), 173*e8d8bef9SDimitry Andric // then it decreses to f(1) = 2 - b 174*e8d8bef9SDimitry Andric // 175*e8d8bef9SDimitry Andric // Let g(x) = x - f(x) = b*x^2 - x. 176*e8d8bef9SDimitry Andric // On (0, 1/b), g(x) < 0 <=> f(x) > x 177*e8d8bef9SDimitry Andric // On (1/b, 1], g(x) > 0 <=> f(x) < x 178*e8d8bef9SDimitry Andric // 179*e8d8bef9SDimitry Andric // For half-width iterations, b_hw is used instead of b. 180*e8d8bef9SDimitry Andric REPEAT_N_TIMES(NUMBER_OF_HALF_ITERATIONS, { 181*e8d8bef9SDimitry Andric // corr_UQ1_hw can be **larger** than 2 - b_hw*x by at most 1*Ulp 182*e8d8bef9SDimitry Andric // of corr_UQ1_hw. 183*e8d8bef9SDimitry Andric // "0.0 - (...)" is equivalent to "2.0 - (...)" in UQ1.(HW-1). 184*e8d8bef9SDimitry Andric // On the other hand, corr_UQ1_hw should not overflow from 2.0 to 0.0 provided 185*e8d8bef9SDimitry Andric // no overflow occurred earlier: ((rep_t)x_UQ0_hw * b_UQ1_hw >> HW) is 186*e8d8bef9SDimitry Andric // expected to be strictly positive because b_UQ1_hw has its highest bit set 187*e8d8bef9SDimitry Andric // and x_UQ0_hw should be rather large (it converges to 1/2 < 1/b_hw <= 1). 188*e8d8bef9SDimitry Andric half_rep_t corr_UQ1_hw = 0 - ((rep_t)x_UQ0_hw * b_UQ1_hw >> HW); 189*e8d8bef9SDimitry Andric 190*e8d8bef9SDimitry Andric // Now, we should multiply UQ0.HW and UQ1.(HW-1) numbers, naturally 191*e8d8bef9SDimitry Andric // obtaining an UQ1.(HW-1) number and proving its highest bit could be 192*e8d8bef9SDimitry Andric // considered to be 0 to be able to represent it in UQ0.HW. 193*e8d8bef9SDimitry Andric // From the above analysis of f(x), if corr_UQ1_hw would be represented 194*e8d8bef9SDimitry Andric // without any intermediate loss of precision (that is, in twice_rep_t) 195*e8d8bef9SDimitry Andric // x_UQ0_hw could be at most [1.]000... if b_hw is exactly 1.0 and strictly 196*e8d8bef9SDimitry Andric // less otherwise. On the other hand, to obtain [1.]000..., one have to pass 197*e8d8bef9SDimitry Andric // 1/b_hw == 1.0 to f(x), so this cannot occur at all without overflow (due 198*e8d8bef9SDimitry Andric // to 1.0 being not representable as UQ0.HW). 199*e8d8bef9SDimitry Andric // The fact corr_UQ1_hw was virtually round up (due to result of 200*e8d8bef9SDimitry Andric // multiplication being **first** truncated, then negated - to improve 201*e8d8bef9SDimitry Andric // error estimations) can increase x_UQ0_hw by up to 2*Ulp of x_UQ0_hw. 202*e8d8bef9SDimitry Andric x_UQ0_hw = (rep_t)x_UQ0_hw * corr_UQ1_hw >> (HW - 1); 203*e8d8bef9SDimitry Andric // Now, either no overflow occurred or x_UQ0_hw is 0 or 1 in its half_rep_t 204*e8d8bef9SDimitry Andric // representation. In the latter case, x_UQ0_hw will be either 0 or 1 after 205*e8d8bef9SDimitry Andric // any number of iterations, so just subtract 2 from the reciprocal 206*e8d8bef9SDimitry Andric // approximation after last iteration. 207*e8d8bef9SDimitry Andric 208*e8d8bef9SDimitry Andric // In infinite precision, with 0 <= eps1, eps2 <= U = 2^-HW: 209*e8d8bef9SDimitry Andric // corr_UQ1_hw = 2 - (1/b_hw + e_n) * b_hw + 2*eps1 210*e8d8bef9SDimitry Andric // = 1 - e_n * b_hw + 2*eps1 211*e8d8bef9SDimitry Andric // x_UQ0_hw = (1/b_hw + e_n) * (1 - e_n*b_hw + 2*eps1) - eps2 212*e8d8bef9SDimitry Andric // = 1/b_hw - e_n + 2*eps1/b_hw + e_n - e_n^2*b_hw + 2*e_n*eps1 - eps2 213*e8d8bef9SDimitry Andric // = 1/b_hw + 2*eps1/b_hw - e_n^2*b_hw + 2*e_n*eps1 - eps2 214*e8d8bef9SDimitry Andric // e_{n+1} = -e_n^2*b_hw + 2*eps1/b_hw + 2*e_n*eps1 - eps2 215*e8d8bef9SDimitry Andric // = 2*e_n*eps1 - (e_n^2*b_hw + eps2) + 2*eps1/b_hw 216*e8d8bef9SDimitry Andric // \------ >0 -------/ \-- >0 ---/ 217*e8d8bef9SDimitry Andric // abs(e_{n+1}) <= 2*abs(e_n)*U + max(2*e_n^2 + U, 2 * U) 218*e8d8bef9SDimitry Andric }) 219*e8d8bef9SDimitry Andric // For initial half-width iterations, U = 2^-HW 220*e8d8bef9SDimitry Andric // Let abs(e_n) <= u_n * U, 221*e8d8bef9SDimitry Andric // then abs(e_{n+1}) <= 2 * u_n * U^2 + max(2 * u_n^2 * U^2 + U, 2 * U) 222*e8d8bef9SDimitry Andric // u_{n+1} <= 2 * u_n * U + max(2 * u_n^2 * U + 1, 2) 223*e8d8bef9SDimitry Andric 224*e8d8bef9SDimitry Andric // Account for possible overflow (see above). For an overflow to occur for the 225*e8d8bef9SDimitry Andric // first time, for "ideal" corr_UQ1_hw (that is, without intermediate 226*e8d8bef9SDimitry Andric // truncation), the result of x_UQ0_hw * corr_UQ1_hw should be either maximum 227*e8d8bef9SDimitry Andric // value representable in UQ0.HW or less by 1. This means that 1/b_hw have to 228*e8d8bef9SDimitry Andric // be not below that value (see g(x) above), so it is safe to decrement just 229*e8d8bef9SDimitry Andric // once after the final iteration. On the other hand, an effective value of 230*e8d8bef9SDimitry Andric // divisor changes after this point (from b_hw to b), so adjust here. 231*e8d8bef9SDimitry Andric x_UQ0_hw -= 1U; 232*e8d8bef9SDimitry Andric rep_t x_UQ0 = (rep_t)x_UQ0_hw << HW; 233*e8d8bef9SDimitry Andric x_UQ0 -= 1U; 234*e8d8bef9SDimitry Andric 235*e8d8bef9SDimitry Andric#else 236*e8d8bef9SDimitry Andric // C is (3/4 + 1/sqrt(2)) - 1 truncated to 32 fractional bits as UQ0.n 237*e8d8bef9SDimitry Andric const rep_t C = REP_C(0x7504F333) << (typeWidth - 32); 238*e8d8bef9SDimitry Andric rep_t x_UQ0 = C - b_UQ1; 239*e8d8bef9SDimitry Andric // E_0 <= 3/4 - 1/sqrt(2) + 2 * 2^-32 240*e8d8bef9SDimitry Andric#endif 241*e8d8bef9SDimitry Andric 242*e8d8bef9SDimitry Andric // Error estimations for full-precision iterations are calculated just 243*e8d8bef9SDimitry Andric // as above, but with U := 2^-W and taking extra decrementing into account. 244*e8d8bef9SDimitry Andric // We need at least one such iteration. 245*e8d8bef9SDimitry Andric 246*e8d8bef9SDimitry Andric#ifdef USE_NATIVE_FULL_ITERATIONS 247*e8d8bef9SDimitry Andric REPEAT_N_TIMES(NUMBER_OF_FULL_ITERATIONS, { 248*e8d8bef9SDimitry Andric rep_t corr_UQ1 = 0 - ((twice_rep_t)x_UQ0 * b_UQ1 >> typeWidth); 249*e8d8bef9SDimitry Andric x_UQ0 = (twice_rep_t)x_UQ0 * corr_UQ1 >> (typeWidth - 1); 250*e8d8bef9SDimitry Andric }) 251*e8d8bef9SDimitry Andric#else 252*e8d8bef9SDimitry Andric#if NUMBER_OF_FULL_ITERATIONS != 1 253*e8d8bef9SDimitry Andric#error Only a single emulated full iteration is supported 254*e8d8bef9SDimitry Andric#endif 255*e8d8bef9SDimitry Andric#if !(NUMBER_OF_HALF_ITERATIONS > 0) 256*e8d8bef9SDimitry Andric // Cannot normally reach here: only one full-width iteration is requested and 257*e8d8bef9SDimitry Andric // the total number of iterations should be at least 3 even for float32. 258*e8d8bef9SDimitry Andric#error Check NUMBER_OF_HALF_ITERATIONS, NUMBER_OF_FULL_ITERATIONS and USE_NATIVE_FULL_ITERATIONS. 259*e8d8bef9SDimitry Andric#endif 260*e8d8bef9SDimitry Andric // Simulating operations on a twice_rep_t to perform a single final full-width 261*e8d8bef9SDimitry Andric // iteration. Using ad-hoc multiplication implementations to take advantage 262*e8d8bef9SDimitry Andric // of particular structure of operands. 263*e8d8bef9SDimitry Andric rep_t blo = b_UQ1 & loMask; 264*e8d8bef9SDimitry Andric // x_UQ0 = x_UQ0_hw * 2^HW - 1 265*e8d8bef9SDimitry Andric // x_UQ0 * b_UQ1 = (x_UQ0_hw * 2^HW) * (b_UQ1_hw * 2^HW + blo) - b_UQ1 266*e8d8bef9SDimitry Andric // 267*e8d8bef9SDimitry Andric // <--- higher half ---><--- lower half ---> 268*e8d8bef9SDimitry Andric // [x_UQ0_hw * b_UQ1_hw] 269*e8d8bef9SDimitry Andric // + [ x_UQ0_hw * blo ] 270*e8d8bef9SDimitry Andric // - [ b_UQ1 ] 271*e8d8bef9SDimitry Andric // = [ result ][.... discarded ...] 272*e8d8bef9SDimitry Andric rep_t corr_UQ1 = 0U - ( (rep_t)x_UQ0_hw * b_UQ1_hw 273*e8d8bef9SDimitry Andric + ((rep_t)x_UQ0_hw * blo >> HW) 274*e8d8bef9SDimitry Andric - REP_C(1)); // account for *possible* carry 275*e8d8bef9SDimitry Andric rep_t lo_corr = corr_UQ1 & loMask; 276*e8d8bef9SDimitry Andric rep_t hi_corr = corr_UQ1 >> HW; 277*e8d8bef9SDimitry Andric // x_UQ0 * corr_UQ1 = (x_UQ0_hw * 2^HW) * (hi_corr * 2^HW + lo_corr) - corr_UQ1 278*e8d8bef9SDimitry Andric x_UQ0 = ((rep_t)x_UQ0_hw * hi_corr << 1) 279*e8d8bef9SDimitry Andric + ((rep_t)x_UQ0_hw * lo_corr >> (HW - 1)) 280*e8d8bef9SDimitry Andric - REP_C(2); // 1 to account for the highest bit of corr_UQ1 can be 1 281*e8d8bef9SDimitry Andric // 1 to account for possible carry 282*e8d8bef9SDimitry Andric // Just like the case of half-width iterations but with possibility 283*e8d8bef9SDimitry Andric // of overflowing by one extra Ulp of x_UQ0. 284*e8d8bef9SDimitry Andric x_UQ0 -= 1U; 285*e8d8bef9SDimitry Andric // ... and then traditional fixup by 2 should work 286*e8d8bef9SDimitry Andric 287*e8d8bef9SDimitry Andric // On error estimation: 288*e8d8bef9SDimitry Andric // abs(E_{N-1}) <= (u_{N-1} + 2 /* due to conversion e_n -> E_n */) * 2^-HW 289*e8d8bef9SDimitry Andric // + (2^-HW + 2^-W)) 290*e8d8bef9SDimitry Andric // abs(E_{N-1}) <= (u_{N-1} + 3.01) * 2^-HW 291*e8d8bef9SDimitry Andric 292*e8d8bef9SDimitry Andric // Then like for the half-width iterations: 293*e8d8bef9SDimitry Andric // With 0 <= eps1, eps2 < 2^-W 294*e8d8bef9SDimitry Andric // E_N = 4 * E_{N-1} * eps1 - (E_{N-1}^2 * b + 4 * eps2) + 4 * eps1 / b 295*e8d8bef9SDimitry Andric // abs(E_N) <= 2^-W * [ 4 * abs(E_{N-1}) + max(2 * abs(E_{N-1})^2 * 2^W + 4, 8)) ] 296*e8d8bef9SDimitry Andric // abs(E_N) <= 2^-W * [ 4 * (u_{N-1} + 3.01) * 2^-HW + max(4 + 2 * (u_{N-1} + 3.01)^2, 8) ] 297*e8d8bef9SDimitry Andric#endif 298*e8d8bef9SDimitry Andric 299*e8d8bef9SDimitry Andric // Finally, account for possible overflow, as explained above. 300*e8d8bef9SDimitry Andric x_UQ0 -= 2U; 301*e8d8bef9SDimitry Andric 302*e8d8bef9SDimitry Andric // u_n for different precisions (with N-1 half-width iterations): 303*e8d8bef9SDimitry Andric // W0 is the precision of C 304*e8d8bef9SDimitry Andric // u_0 = (3/4 - 1/sqrt(2) + 2^-W0) * 2^HW 305*e8d8bef9SDimitry Andric 306*e8d8bef9SDimitry Andric // Estimated with bc: 307*e8d8bef9SDimitry Andric // define half1(un) { return 2.0 * (un + un^2) / 2.0^hw + 1.0; } 308*e8d8bef9SDimitry Andric // define half2(un) { return 2.0 * un / 2.0^hw + 2.0; } 309*e8d8bef9SDimitry Andric // define full1(un) { return 4.0 * (un + 3.01) / 2.0^hw + 2.0 * (un + 3.01)^2 + 4.0; } 310*e8d8bef9SDimitry Andric // define full2(un) { return 4.0 * (un + 3.01) / 2.0^hw + 8.0; } 311*e8d8bef9SDimitry Andric 312*e8d8bef9SDimitry Andric // | f32 (0 + 3) | f32 (2 + 1) | f64 (3 + 1) | f128 (4 + 1) 313*e8d8bef9SDimitry Andric // u_0 | < 184224974 | < 2812.1 | < 184224974 | < 791240234244348797 314*e8d8bef9SDimitry Andric // u_1 | < 15804007 | < 242.7 | < 15804007 | < 67877681371350440 315*e8d8bef9SDimitry Andric // u_2 | < 116308 | < 2.81 | < 116308 | < 499533100252317 316*e8d8bef9SDimitry Andric // u_3 | < 7.31 | | < 7.31 | < 27054456580 317*e8d8bef9SDimitry Andric // u_4 | | | | < 80.4 318*e8d8bef9SDimitry Andric // Final (U_N) | same as u_3 | < 72 | < 218 | < 13920 319*e8d8bef9SDimitry Andric 320*e8d8bef9SDimitry Andric // Add 2 to U_N due to final decrement. 321*e8d8bef9SDimitry Andric 322*e8d8bef9SDimitry Andric#if defined(SINGLE_PRECISION) && NUMBER_OF_HALF_ITERATIONS == 2 && NUMBER_OF_FULL_ITERATIONS == 1 323*e8d8bef9SDimitry Andric#define RECIPROCAL_PRECISION REP_C(74) 324*e8d8bef9SDimitry Andric#elif defined(SINGLE_PRECISION) && NUMBER_OF_HALF_ITERATIONS == 0 && NUMBER_OF_FULL_ITERATIONS == 3 325*e8d8bef9SDimitry Andric#define RECIPROCAL_PRECISION REP_C(10) 326*e8d8bef9SDimitry Andric#elif defined(DOUBLE_PRECISION) && NUMBER_OF_HALF_ITERATIONS == 3 && NUMBER_OF_FULL_ITERATIONS == 1 327*e8d8bef9SDimitry Andric#define RECIPROCAL_PRECISION REP_C(220) 328*e8d8bef9SDimitry Andric#elif defined(QUAD_PRECISION) && NUMBER_OF_HALF_ITERATIONS == 4 && NUMBER_OF_FULL_ITERATIONS == 1 329*e8d8bef9SDimitry Andric#define RECIPROCAL_PRECISION REP_C(13922) 330*e8d8bef9SDimitry Andric#else 331*e8d8bef9SDimitry Andric#error Invalid number of iterations 332*e8d8bef9SDimitry Andric#endif 333*e8d8bef9SDimitry Andric 334*e8d8bef9SDimitry Andric // Suppose 1/b - P * 2^-W < x < 1/b + P * 2^-W 335*e8d8bef9SDimitry Andric x_UQ0 -= RECIPROCAL_PRECISION; 336*e8d8bef9SDimitry Andric // Now 1/b - (2*P) * 2^-W < x < 1/b 337*e8d8bef9SDimitry Andric // FIXME Is x_UQ0 still >= 0.5? 338*e8d8bef9SDimitry Andric 339*e8d8bef9SDimitry Andric rep_t quotient_UQ1, dummy; 340*e8d8bef9SDimitry Andric wideMultiply(x_UQ0, aSignificand << 1, "ient_UQ1, &dummy); 341*e8d8bef9SDimitry Andric // Now, a/b - 4*P * 2^-W < q < a/b for q=<quotient_UQ1:dummy> in UQ1.(SB+1+W). 342*e8d8bef9SDimitry Andric 343*e8d8bef9SDimitry Andric // quotient_UQ1 is in [0.5, 2.0) as UQ1.(SB+1), 344*e8d8bef9SDimitry Andric // adjust it to be in [1.0, 2.0) as UQ1.SB. 345*e8d8bef9SDimitry Andric rep_t residualLo; 346*e8d8bef9SDimitry Andric if (quotient_UQ1 < (implicitBit << 1)) { 347*e8d8bef9SDimitry Andric // Highest bit is 0, so just reinterpret quotient_UQ1 as UQ1.SB, 348*e8d8bef9SDimitry Andric // effectively doubling its value as well as its error estimation. 349*e8d8bef9SDimitry Andric residualLo = (aSignificand << (significandBits + 1)) - quotient_UQ1 * bSignificand; 350*e8d8bef9SDimitry Andric writtenExponent -= 1; 351*e8d8bef9SDimitry Andric aSignificand <<= 1; 352*e8d8bef9SDimitry Andric } else { 353*e8d8bef9SDimitry Andric // Highest bit is 1 (the UQ1.(SB+1) value is in [1, 2)), convert it 354*e8d8bef9SDimitry Andric // to UQ1.SB by right shifting by 1. Least significant bit is omitted. 355*e8d8bef9SDimitry Andric quotient_UQ1 >>= 1; 356*e8d8bef9SDimitry Andric residualLo = (aSignificand << significandBits) - quotient_UQ1 * bSignificand; 357*e8d8bef9SDimitry Andric } 358*e8d8bef9SDimitry Andric // NB: residualLo is calculated above for the normal result case. 359*e8d8bef9SDimitry Andric // It is re-computed on denormal path that is expected to be not so 360*e8d8bef9SDimitry Andric // performance-sensitive. 361*e8d8bef9SDimitry Andric 362*e8d8bef9SDimitry Andric // Now, q cannot be greater than a/b and can differ by at most 8*P * 2^-W + 2^-SB 363*e8d8bef9SDimitry Andric // Each NextAfter() increments the floating point value by at least 2^-SB 364*e8d8bef9SDimitry Andric // (more, if exponent was incremented). 365*e8d8bef9SDimitry Andric // Different cases (<---> is of 2^-SB length, * = a/b that is shown as a midpoint): 366*e8d8bef9SDimitry Andric // q 367*e8d8bef9SDimitry Andric // | | * | | | | | 368*e8d8bef9SDimitry Andric // <---> 2^t 369*e8d8bef9SDimitry Andric // | | | | | * | | 370*e8d8bef9SDimitry Andric // q 371*e8d8bef9SDimitry Andric // To require at most one NextAfter(), an error should be less than 1.5 * 2^-SB. 372*e8d8bef9SDimitry Andric // (8*P) * 2^-W + 2^-SB < 1.5 * 2^-SB 373*e8d8bef9SDimitry Andric // (8*P) * 2^-W < 0.5 * 2^-SB 374*e8d8bef9SDimitry Andric // P < 2^(W-4-SB) 375*e8d8bef9SDimitry Andric // Generally, for at most R NextAfter() to be enough, 376*e8d8bef9SDimitry Andric // P < (2*R - 1) * 2^(W-4-SB) 377*e8d8bef9SDimitry Andric // For f32 (0+3): 10 < 32 (OK) 378*e8d8bef9SDimitry Andric // For f32 (2+1): 32 < 74 < 32 * 3, so two NextAfter() are required 379*e8d8bef9SDimitry Andric // For f64: 220 < 256 (OK) 380*e8d8bef9SDimitry Andric // For f128: 4096 * 3 < 13922 < 4096 * 5 (three NextAfter() are required) 381*e8d8bef9SDimitry Andric 382*e8d8bef9SDimitry Andric // If we have overflowed the exponent, return infinity 383*e8d8bef9SDimitry Andric if (writtenExponent >= maxExponent) 384*e8d8bef9SDimitry Andric return fromRep(infRep | quotientSign); 385*e8d8bef9SDimitry Andric 386*e8d8bef9SDimitry Andric // Now, quotient_UQ1_SB <= the correctly-rounded result 387*e8d8bef9SDimitry Andric // and may need taking NextAfter() up to 3 times (see error estimates above) 388*e8d8bef9SDimitry Andric // r = a - b * q 389*e8d8bef9SDimitry Andric rep_t absResult; 390*e8d8bef9SDimitry Andric if (writtenExponent > 0) { 391*e8d8bef9SDimitry Andric // Clear the implicit bit 392*e8d8bef9SDimitry Andric absResult = quotient_UQ1 & significandMask; 393*e8d8bef9SDimitry Andric // Insert the exponent 394*e8d8bef9SDimitry Andric absResult |= (rep_t)writtenExponent << significandBits; 395*e8d8bef9SDimitry Andric residualLo <<= 1; 396*e8d8bef9SDimitry Andric } else { 397*e8d8bef9SDimitry Andric // Prevent shift amount from being negative 398*e8d8bef9SDimitry Andric if (significandBits + writtenExponent < 0) 399*e8d8bef9SDimitry Andric return fromRep(quotientSign); 400*e8d8bef9SDimitry Andric 401*e8d8bef9SDimitry Andric absResult = quotient_UQ1 >> (-writtenExponent + 1); 402*e8d8bef9SDimitry Andric 403*e8d8bef9SDimitry Andric // multiplied by two to prevent shift amount to be negative 404*e8d8bef9SDimitry Andric residualLo = (aSignificand << (significandBits + writtenExponent)) - (absResult * bSignificand << 1); 405*e8d8bef9SDimitry Andric } 406*e8d8bef9SDimitry Andric 407*e8d8bef9SDimitry Andric // Round 408*e8d8bef9SDimitry Andric residualLo += absResult & 1; // tie to even 409*e8d8bef9SDimitry Andric // The above line conditionally turns the below LT comparison into LTE 410*e8d8bef9SDimitry Andric absResult += residualLo > bSignificand; 411*e8d8bef9SDimitry Andric#if defined(QUAD_PRECISION) || (defined(SINGLE_PRECISION) && NUMBER_OF_HALF_ITERATIONS > 0) 412*e8d8bef9SDimitry Andric // Do not round Infinity to NaN 413*e8d8bef9SDimitry Andric absResult += absResult < infRep && residualLo > (2 + 1) * bSignificand; 414*e8d8bef9SDimitry Andric#endif 415*e8d8bef9SDimitry Andric#if defined(QUAD_PRECISION) 416*e8d8bef9SDimitry Andric absResult += absResult < infRep && residualLo > (4 + 1) * bSignificand; 417*e8d8bef9SDimitry Andric#endif 418*e8d8bef9SDimitry Andric return fromRep(absResult | quotientSign); 419*e8d8bef9SDimitry Andric} 420