blob: 656ade4f7735354cbc045ccc167d7fa91c12424e (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
|
//===-- Square root of x86 long double numbers ------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_SQRT_80_BIT_LONG_DOUBLE_H
#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_SQRT_80_BIT_LONG_DOUBLE_H
#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/rounding_mode.h"
#include "src/__support/UInt128.h"
#include "src/__support/common.h"
namespace LIBC_NAMESPACE {
namespace fputil {
namespace x86 {
LIBC_INLINE void normalize(int &exponent, UInt128 &mantissa) {
const unsigned int shift = static_cast<unsigned int>(
cpp::countl_zero(static_cast<uint64_t>(mantissa)) -
(8 * sizeof(uint64_t) - 1 - FPBits<long double>::FRACTION_LEN));
exponent -= shift;
mantissa <<= shift;
}
// if constexpr statement in sqrt.h still requires x86::sqrt to be declared
// even when it's not used.
LIBC_INLINE long double sqrt(long double x);
// Correctly rounded SQRT for all rounding modes.
// Shift-and-add algorithm.
#if defined(LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80)
LIBC_INLINE long double sqrt(long double x) {
using LDBits = FPBits<long double>;
using StorageType = typename LDBits::StorageType;
constexpr StorageType ONE = StorageType(1) << int(LDBits::FRACTION_LEN);
constexpr auto LDNAN = LDBits::quiet_nan().get_val();
LDBits bits(x);
if (bits == LDBits::inf(Sign::POS) || bits.is_zero() || bits.is_nan()) {
// sqrt(+Inf) = +Inf
// sqrt(+0) = +0
// sqrt(-0) = -0
// sqrt(NaN) = NaN
// sqrt(-NaN) = -NaN
return x;
} else if (bits.is_neg()) {
// sqrt(-Inf) = NaN
// sqrt(-x) = NaN
return LDNAN;
} else {
int x_exp = bits.get_explicit_exponent();
StorageType x_mant = bits.get_mantissa();
// Step 1a: Normalize denormal input
if (bits.get_implicit_bit()) {
x_mant |= ONE;
} else if (bits.is_subnormal()) {
normalize(x_exp, x_mant);
}
// Step 1b: Make sure the exponent is even.
if (x_exp & 1) {
--x_exp;
x_mant <<= 1;
}
// After step 1b, x = 2^(x_exp) * x_mant, where x_exp is even, and
// 1 <= x_mant < 4. So sqrt(x) = 2^(x_exp / 2) * y, with 1 <= y < 2.
// Notice that the output of sqrt is always in the normal range.
// To perform shift-and-add algorithm to find y, let denote:
// y(n) = 1.y_1 y_2 ... y_n, we can define the nth residue to be:
// r(n) = 2^n ( x_mant - y(n)^2 ).
// That leads to the following recurrence formula:
// r(n) = 2*r(n-1) - y_n*[ 2*y(n-1) + 2^(-n-1) ]
// with the initial conditions: y(0) = 1, and r(0) = x - 1.
// So the nth digit y_n of the mantissa of sqrt(x) can be found by:
// y_n = 1 if 2*r(n-1) >= 2*y(n - 1) + 2^(-n-1)
// 0 otherwise.
StorageType y = ONE;
StorageType r = x_mant - ONE;
for (StorageType current_bit = ONE >> 1; current_bit; current_bit >>= 1) {
r <<= 1;
StorageType tmp = (y << 1) + current_bit; // 2*y(n - 1) + 2^(-n-1)
if (r >= tmp) {
r -= tmp;
y += current_bit;
}
}
// We compute one more iteration in order to round correctly.
bool lsb = static_cast<bool>(y & 1); // Least significant bit
bool rb = false; // Round bit
r <<= 2;
StorageType tmp = (y << 2) + 1;
if (r >= tmp) {
r -= tmp;
rb = true;
}
// Append the exponent field.
x_exp = ((x_exp >> 1) + LDBits::EXP_BIAS);
y |= (static_cast<StorageType>(x_exp) << (LDBits::FRACTION_LEN + 1));
switch (quick_get_round()) {
case FE_TONEAREST:
// Round to nearest, ties to even
if (rb && (lsb || (r != 0)))
++y;
break;
case FE_UPWARD:
if (rb || (r != 0))
++y;
break;
}
// Extract output
FPBits<long double> out(0.0L);
out.set_biased_exponent(x_exp);
out.set_implicit_bit(1);
out.set_mantissa((y & (ONE - 1)));
return out.get_val();
}
}
#endif // LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80
} // namespace x86
} // namespace fputil
} // namespace LIBC_NAMESPACE
#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_SQRT_80_BIT_LONG_DOUBLE_H
|
