blob: 614aa399fcc21b8227c7bae0a9c2a67ca461e3c9 (
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
|
//===-- Implementation of hypotf function ---------------------------------===//
//
// 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
//
//===----------------------------------------------------------------------===//
#include "src/math/hypotf.h"
#include "src/__support/FPUtil/BasicOperations.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/sqrt.h"
#include "src/__support/common.h"
namespace LIBC_NAMESPACE {
LLVM_LIBC_FUNCTION(float, hypotf, (float x, float y)) {
using DoubleBits = fputil::FPBits<double>;
using FPBits = fputil::FPBits<float>;
FPBits x_bits(x), y_bits(y);
uint16_t x_exp = x_bits.get_biased_exponent();
uint16_t y_exp = y_bits.get_biased_exponent();
uint16_t exp_diff = (x_exp > y_exp) ? (x_exp - y_exp) : (y_exp - x_exp);
if (exp_diff >= FPBits::FRACTION_LEN + 2) {
return fputil::abs(x) + fputil::abs(y);
}
double xd = static_cast<double>(x);
double yd = static_cast<double>(y);
// These squares are exact.
double x_sq = xd * xd;
double y_sq = yd * yd;
// Compute the sum of squares.
double sum_sq = x_sq + y_sq;
// Compute the rounding error with Fast2Sum algorithm:
// x_sq + y_sq = sum_sq - err
double err = (x_sq >= y_sq) ? (sum_sq - x_sq) - y_sq : (sum_sq - y_sq) - x_sq;
// Take sqrt in double precision.
DoubleBits result(fputil::sqrt(sum_sq));
if (!DoubleBits(sum_sq).is_inf_or_nan()) {
// Correct rounding.
double r_sq = result.get_val() * result.get_val();
double diff = sum_sq - r_sq;
constexpr uint64_t mask = 0x0000'0000'3FFF'FFFFULL;
uint64_t lrs = result.uintval() & mask;
if (lrs == 0x0000'0000'1000'0000ULL && err < diff) {
result.set_uintval(result.uintval() | 1ULL);
} else if (lrs == 0x0000'0000'3000'0000ULL && err > diff) {
result.set_uintval(result.uintval() - 1ULL);
}
} else {
FPBits bits_x(x), bits_y(y);
if (bits_x.is_inf_or_nan() || bits_y.is_inf_or_nan()) {
if (bits_x.is_inf() || bits_y.is_inf())
return FPBits::inf().get_val();
if (bits_x.is_nan())
return x;
return y;
}
}
return static_cast<float>(result.get_val());
}
} // namespace LIBC_NAMESPACE
|
