1//===-- Implementation header for log1pf ------------------------*- 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#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_LOG1PF_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_LOG1PF_H
11
12#include "src/__support/FPUtil/FEnvImpl.h"
13#include "src/__support/FPUtil/FMA.h"
14#include "src/__support/FPUtil/FPBits.h"
15#include "src/__support/FPUtil/PolyEval.h"
16#include "src/__support/FPUtil/except_value_utils.h"
17#include "src/__support/FPUtil/multiply_add.h"
18#include "src/__support/common.h"
19#include "src/__support/macros/config.h"
20#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
21#include "src/__support/macros/properties/cpu_features.h"
22#include "src/__support/math/acoshf_utils.h"
23
24// This is an algorithm for log(1+x) in single precision which is
25// correctly rounded for all rounding modes.
26// - An exhaustive test show that when x >= 2^45, log1pf(x) == logf(x)
27// for all rounding modes.
28// - When 2^(-6) <= |x| < 2^45, the sum (double(x) + 1.0) is exact,
29// so we can adapt the correctly rounded algorithm of logf to compute
30// log(double(x) + 1.0) correctly. For more information about the logf
31// algorithm, see `libc/src/math/generic/logf.cpp`.
32// - When |x| < 2^(-6), we use a degree-8 polynomial in double precision
33// generated with Sollya using the following command:
34// fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]);
35
36namespace LIBC_NAMESPACE_DECL {
37
38namespace math {
39
40LIBC_INLINE float log1pf(float x) {
41 using FPBits = typename fputil::FPBits<float>;
42 FPBits xbits(x);
43 uint32_t x_u = xbits.uintval();
44 uint32_t x_a = x_u & 0x7fff'ffffU;
45 double xd = static_cast<double>(x);
46
47 if (x_a <= 0x3c80'0000U) {
48 // |x| <= 2^-6.
49#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
50 // Hard-to round cases.
51 switch (x_u) {
52 case 0x3540'0003U: // x = 0x1.800006p-21f
53 return fputil::round_result_slightly_down(value_rn: 0x1.7ffffep-21f);
54 case 0x3710'001bU: // x = 0x1.200036p-17f
55 return fputil::round_result_slightly_down(value_rn: 0x1.1fffe6p-17f);
56 case 0xb53f'fffdU: // x = -0x1.7ffffap-21
57 return fputil::round_result_slightly_down(value_rn: -0x1.800002p-21f);
58 case 0xb70f'ffe5U: // x = -0x1.1fffcap-17
59 return fputil::round_result_slightly_down(value_rn: -0x1.20001ap-17f);
60 case 0xbb0e'c8c4U: // x = -0x1.1d9188p-9
61 return fputil::round_result_slightly_up(value_rn: -0x1.1de14ap-9f);
62 }
63#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
64
65 // Polynomial generated by Sollya with:
66 // > P = fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]);
67 // > dirtyinfnorm((log(1 + x) - x*P)/log(1 + x), [-2^-6, 2^-6]);
68 // 0x1.1447755e54a327941f7db7316f8dcd7cf33d15ffp-58
69 constexpr double COEFFS[7] = {-0x1.0000000000000p-1, 0x1.5555555556aadp-2,
70 -0x1.000000000181ap-2, 0x1.999998998124ep-3,
71 -0x1.55555452e2a2bp-3, 0x1.24adb8cde4aa7p-3,
72 -0x1.0019db915ef6fp-3};
73
74 double xsq = xd * xd;
75 double c0 = fputil::multiply_add(x: xd, y: COEFFS[1], z: COEFFS[0]);
76 double c1 = fputil::multiply_add(x: xd, y: COEFFS[3], z: COEFFS[2]);
77 double c2 = fputil::multiply_add(x: xd, y: COEFFS[5], z: COEFFS[4]);
78 double x4 = xsq * xsq;
79 double d0 = fputil::multiply_add(x: xsq, y: c1, z: c0);
80 double d1 = fputil::multiply_add(x: xsq, y: COEFFS[6], z: c2);
81 double d2 = fputil::multiply_add(x: x4, y: d1, z: d0);
82 double r = fputil::multiply_add(x: xsq, y: d2, z: xd);
83
84 return static_cast<float>(r);
85 }
86
87 // Use log1p(x) = log(1 + x) for |x| > 2^-6;
88
89 // Check for exceptional cases.
90 if (x_a >= 0x3f80'0000) {
91 // |x| >= 1.
92 if (LIBC_UNLIKELY(x_u >= 0x7f80'0000)) {
93 // x is inf, nan, or x <= -1.
94 if (x == -1.0f) {
95 // x = -1
96 fputil::set_errno_if_required(ERANGE);
97 fputil::raise_except_if_required(FE_DIVBYZERO);
98 return FPBits::inf(sign: Sign::NEG).get_val();
99 }
100 if (xbits.is_signaling_nan() || x < 1.0f) {
101 // x is signaling NaNs or x < -1
102 if (x < 1.0f)
103 fputil::set_errno_if_required(EDOM);
104 fputil::raise_except_if_required(FE_INVALID);
105 return fputil::FPBits<float>::quiet_nan().get_val();
106 }
107 // x is +inf or quiet NaN
108 return x;
109 }
110#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
111 // Filter hard-to-round cases:
112 if (LIBC_UNLIKELY(x >= 0x1.30bf04p+43f)) {
113 switch (x_u) {
114 case 0x5518'5f82U: // x = 0x1.30bf04p+43f
115 return fputil::round_result_slightly_up(value_rn: 0x1.dfac9p+4f);
116 case 0x5cd6'9e88U: // x = 0x1.ad3d1p+58f
117 return fputil::round_result_slightly_up(value_rn: 0x1.45c146p+5f);
118 case 0x5ee8'984eU: // x = 0x1.d1309cp+62f
119 return fputil::round_result_slightly_up(value_rn: 0x1.5c9442p+5f);
120 case 0x65d8'90d3U: // x = 0x1.b121a6p+76f
121 return fputil::round_result_slightly_down(value_rn: 0x1.a9a3f2p+5f);
122 case 0x6f31'a8ecU: // x = 0x1.6351d8p+95f
123 return fputil::round_result_slightly_down(value_rn: 0x1.08b512p+6f);
124 case 0x7a17'f30aU: // x = 0x1.2fe614p+117f
125 return fputil::round_result_slightly_up(value_rn: 0x1.451436p+6f);
126#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
127 case 0x58f1'9e31U: // x = 0x1.e33c62p+50f
128 return fputil::round_result_slightly_down(value_rn: 0x1.1a576cp+5f);
129 case 0x665e'7ca6U: // x = 0x1.bcf94cp+77f
130 return fputil::round_result_slightly_up(value_rn: 0x1.af66cp+5f);
131 case 0x79e7'ec37U: // x = 0x1.cfd86ep+116f
132 return fputil::round_result_slightly_up(value_rn: 0x1.43ff6ep+6f);
133#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
134 }
135 }
136 } else {
137 if (LIBC_UNLIKELY(x_u == 0x3efd'81adU)) // x = 0x1.fb035ap-2f
138 return fputil::round_result_slightly_up(value_rn: 0x1.9bddc2p-2f);
139#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
140 }
141
142 double r = acoshf_internal::log_eval(x: xd + 1.0);
143 return static_cast<float>(r);
144}
145
146} // namespace math
147} // namespace LIBC_NAMESPACE_DECL
148
149#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LOG1PF_H
150