|
| 1 | +//===-- Implementation of hypotf function ---------------------------------===// |
| 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 | +#include "src/__support/common.h" |
| 9 | +#include "utils/FPUtil/BasicOperations.h" |
| 10 | +#include "utils/FPUtil/FPBits.h" |
| 11 | + |
| 12 | +namespace __llvm_libc { |
| 13 | + |
| 14 | +using namespace fputil; |
| 15 | + |
| 16 | +uint32_t findLeadingOne(uint32_t mant, int &shift_length) { |
| 17 | + shift_length = 0; |
| 18 | + constexpr int nsteps = 5; |
| 19 | + constexpr uint32_t bounds[nsteps] = {1 << 16, 1 << 8, 1 << 4, 1 << 2, 1 << 1}; |
| 20 | + constexpr int shifts[nsteps] = {16, 8, 4, 2, 1}; |
| 21 | + for (int i = 0; i < nsteps; ++i) { |
| 22 | + if (mant >= bounds[i]) { |
| 23 | + shift_length += shifts[i]; |
| 24 | + mant >>= shifts[i]; |
| 25 | + } |
| 26 | + } |
| 27 | + return 1U << shift_length; |
| 28 | +} |
| 29 | + |
| 30 | +// Correctly rounded IEEE 754 HYPOT(x, y) with round to nearest, ties to even. |
| 31 | +// |
| 32 | +// Algorithm: |
| 33 | +// - Let a = max(|x|, |y|), b = min(|x|, |y|), then we have that: |
| 34 | +// a <= sqrt(a^2 + b^2) <= min(a + b, a*sqrt(2)) |
| 35 | +// 1. So if b < eps(a)/2, then HYPOT(x, y) = a. |
| 36 | +// |
| 37 | +// - Moreover, the exponent part of HYPOT(x, y) is either the same or 1 more |
| 38 | +// than the exponent part of a. |
| 39 | +// |
| 40 | +// 2. For the remaining cases, we will use the digit-by-digit (shift-and-add) |
| 41 | +// algorithm to compute SQRT(Z): |
| 42 | +// |
| 43 | +// - For Y = y0.y1...yn... = SQRT(Z), |
| 44 | +// let Y(n) = y0.y1...yn be the first n fractional digits of Y. |
| 45 | +// |
| 46 | +// - The nth scaled residual R(n) is defined to be: |
| 47 | +// R(n) = 2^n * (Z - Y(n)^2) |
| 48 | +// |
| 49 | +// - Since Y(n) = Y(n - 1) + yn * 2^(-n), the scaled residual |
| 50 | +// satisfies the following recurrence formula: |
| 51 | +// R(n) = 2*R(n - 1) - yn*(2*Y(n - 1) + 2^(-n)), |
| 52 | +// with the initial conditions: |
| 53 | +// Y(0) = y0, and R(0) = Z - y0. |
| 54 | +// |
| 55 | +// - So the nth fractional digit of Y = SQRT(Z) can be decided by: |
| 56 | +// yn = 1 if 2*R(n - 1) >= 2*Y(n - 1) + 2^(-n), |
| 57 | +// 0 otherwise. |
| 58 | +// |
| 59 | +// 3. Precision analysis: |
| 60 | +// |
| 61 | +// - Notice that in the decision function: |
| 62 | +// 2*R(n - 1) >= 2*Y(n - 1) + 2^(-n), |
| 63 | +// the right hand side only uses up to the 2^(-n)-bit, and both sides are |
| 64 | +// non-negative, so R(n - 1) can be truncated at the 2^(-(n + 1))-bit, so |
| 65 | +// that 2*R(n - 1) is corrected up to the 2^(-n)-bit. |
| 66 | +// |
| 67 | +// - Thus, in order to round SQRT(a^2 + b^2) correctly up to n-fractional |
| 68 | +// bits, we need to perform the summation (a^2 + b^2) correctly up to (2n + |
| 69 | +// 2)-fractional bits, and the remaining bits are sticky bits (i.e. we only |
| 70 | +// care if they are 0 or > 0), and the comparisons, additions/subtractions |
| 71 | +// can be done in n-fractional bits precision. |
| 72 | +// |
| 73 | +// - For single precision (float), we can use uint64_t to store the sum a^2 + |
| 74 | +// b^2 exact up to (2n + 2)-fractional bits. |
| 75 | +// |
| 76 | +// - Then we can feed this sum into the digit-by-digit algorithm for SQRT(Z) |
| 77 | +// described above. |
| 78 | +// |
| 79 | +// |
| 80 | +// Special cases: |
| 81 | +// - HYPOT(x, y) is +Inf if x or y is +Inf or -Inf; else |
| 82 | +// - HYPOT(x, y) is NaN if x or y is NaN. |
| 83 | +// |
| 84 | +float LLVM_LIBC_ENTRYPOINT(hypotf)(float x, float y) { |
| 85 | + FPBits<float> x_bits(x), y_bits(y); |
| 86 | + |
| 87 | + if (x_bits.isInf() || y_bits.isInf()) { |
| 88 | + return FPBits<float>::inf(); |
| 89 | + } |
| 90 | + if (x_bits.isNaN()) { |
| 91 | + return x; |
| 92 | + } |
| 93 | + if (y_bits.isNaN()) { |
| 94 | + return y; |
| 95 | + } |
| 96 | + |
| 97 | + uint16_t a_exp, b_exp, out_exp; |
| 98 | + uint32_t a_mant, b_mant; |
| 99 | + uint64_t a_mant_sq, b_mant_sq; |
| 100 | + bool sticky_bits; |
| 101 | + |
| 102 | + if ((x_bits.exponent >= y_bits.exponent + MantissaWidth<float>::value + 2) || |
| 103 | + (y == 0)) { |
| 104 | + return abs(x); |
| 105 | + } else if ((y_bits.exponent >= |
| 106 | + x_bits.exponent + MantissaWidth<float>::value + 2) || |
| 107 | + (x == 0)) { |
| 108 | + y_bits.sign = 0; |
| 109 | + return abs(y); |
| 110 | + } |
| 111 | + |
| 112 | + if (x >= y) { |
| 113 | + a_exp = x_bits.exponent; |
| 114 | + a_mant = x_bits.mantissa; |
| 115 | + b_exp = y_bits.exponent; |
| 116 | + b_mant = y_bits.mantissa; |
| 117 | + } else { |
| 118 | + a_exp = y_bits.exponent; |
| 119 | + a_mant = y_bits.mantissa; |
| 120 | + b_exp = x_bits.exponent; |
| 121 | + b_mant = x_bits.mantissa; |
| 122 | + } |
| 123 | + |
| 124 | + out_exp = a_exp; |
| 125 | + |
| 126 | + // Add an extra bit to simplify the final rounding bit computation. |
| 127 | + constexpr uint32_t one = 1U << (MantissaWidth<float>::value + 1); |
| 128 | + |
| 129 | + a_mant <<= 1; |
| 130 | + b_mant <<= 1; |
| 131 | + |
| 132 | + uint32_t leading_one; |
| 133 | + int y_mant_width; |
| 134 | + if (a_exp != 0) { |
| 135 | + leading_one = one; |
| 136 | + a_mant |= one; |
| 137 | + y_mant_width = MantissaWidth<float>::value + 1; |
| 138 | + } else { |
| 139 | + leading_one = findLeadingOne(a_mant, y_mant_width); |
| 140 | + } |
| 141 | + |
| 142 | + if (b_exp != 0) { |
| 143 | + b_mant |= one; |
| 144 | + } |
| 145 | + |
| 146 | + a_mant_sq = static_cast<uint64_t>(a_mant) * a_mant; |
| 147 | + b_mant_sq = static_cast<uint64_t>(b_mant) * b_mant; |
| 148 | + |
| 149 | + // At this point, a_exp >= b_exp > a_exp - 25, so in order to line up aSqMant |
| 150 | + // and bSqMant, we need to shift bSqMant to the right by (a_exp - b_exp) bits. |
| 151 | + // But before that, remember to store the losing bits to sticky. |
| 152 | + // The shift length is for a^2 and b^2, so it's double of the exponent |
| 153 | + // difference between a and b. |
| 154 | + uint16_t shift_length = 2 * (a_exp - b_exp); |
| 155 | + sticky_bits = ((b_mant_sq & ((1ULL << shift_length) - 1)) != 0); |
| 156 | + b_mant_sq >>= shift_length; |
| 157 | + |
| 158 | + uint64_t sum = a_mant_sq + b_mant_sq; |
| 159 | + if (sum >= (1ULL << (2 * y_mant_width + 2))) { |
| 160 | + // a^2 + b^2 >= 4* leading_one^2, so we will need an extra bit to the left. |
| 161 | + if (leading_one == one) { |
| 162 | + // For normal result, we discard the last 2 bits of the sum and increase |
| 163 | + // the exponent. |
| 164 | + sticky_bits = sticky_bits || ((sum & 0x3U) != 0); |
| 165 | + sum >>= 2; |
| 166 | + ++out_exp; |
| 167 | + if (out_exp >= FPBits<float>::maxExponent) { |
| 168 | + return FPBits<float>::inf(); |
| 169 | + } |
| 170 | + } else { |
| 171 | + // For denormal result, we simply move the leading bit of the result to |
| 172 | + // the left by 1. |
| 173 | + leading_one <<= 1; |
| 174 | + ++y_mant_width; |
| 175 | + } |
| 176 | + } |
| 177 | + |
| 178 | + uint32_t Y = leading_one; |
| 179 | + uint32_t R = static_cast<uint32_t>(sum >> y_mant_width) - leading_one; |
| 180 | + uint32_t tailBits = static_cast<uint32_t>(sum) & (leading_one - 1); |
| 181 | + |
| 182 | + for (uint32_t current_bit = leading_one >> 1; current_bit; |
| 183 | + current_bit >>= 1) { |
| 184 | + R = (R << 1) + ((tailBits & current_bit) ? 1 : 0); |
| 185 | + uint32_t tmp = (Y << 1) + current_bit; // 2*y(n - 1) + 2^(-n) |
| 186 | + if (R >= tmp) { |
| 187 | + R -= tmp; |
| 188 | + Y += current_bit; |
| 189 | + } |
| 190 | + } |
| 191 | + |
| 192 | + bool round_bit = Y & 1U; |
| 193 | + bool lsb = Y & 2U; |
| 194 | + |
| 195 | + if (Y >= one) { |
| 196 | + Y -= one; |
| 197 | + |
| 198 | + if (out_exp == 0) { |
| 199 | + out_exp = 1; |
| 200 | + } |
| 201 | + } |
| 202 | + |
| 203 | + Y >>= 1; |
| 204 | + |
| 205 | + // Round to the nearest, tie to even. |
| 206 | + if (round_bit && (lsb || sticky_bits || (R != 0))) { |
| 207 | + ++Y; |
| 208 | + } |
| 209 | + |
| 210 | + if (Y >= (one >> 1)) { |
| 211 | + Y -= one >> 1; |
| 212 | + ++out_exp; |
| 213 | + if (out_exp >= FPBits<float>::maxExponent) { |
| 214 | + return FPBits<float>::inf(); |
| 215 | + } |
| 216 | + } |
| 217 | + |
| 218 | + Y |= static_cast<uint32_t>(out_exp) << MantissaWidth<float>::value; |
| 219 | + return *reinterpret_cast<float *>(&Y); |
| 220 | +} |
| 221 | + |
| 222 | +} // namespace __llvm_libc |
0 commit comments