1//===-- Shared helpers for gamma family math functions ----------*- 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_GAMMA_UTIL_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_GAMMA_UTIL_H
11
12#include "src/__support/CPP/bit.h"
13#include "src/__support/FPUtil/FPBits.h"
14#include "src/__support/FPUtil/multiply_add.h"
15#include "src/__support/macros/attributes.h"
16#include "src/__support/macros/config.h"
17
18namespace LIBC_NAMESPACE_DECL {
19
20namespace math {
21
22namespace gamma_internal {
23
24template <typename T> LIBC_INLINE constexpr bool is_integer(T x) {
25 using FPBits = fputil::FPBits<T>;
26 using StorageType = typename FPBits::StorageType;
27 FPBits xbits(x);
28 StorageType x_u = xbits.uintval();
29 unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
30 unsigned lsb = static_cast<unsigned>(
31 cpp::countr_zero(static_cast<StorageType>(x_u | FPBits::EXP_MASK)));
32 constexpr unsigned UNIT_EXPONENT =
33 static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
34 return x_e + lsb >= UNIT_EXPONENT;
35}
36
37// sin(pi*x) for x in (0, 1), written as (0.25 - u^2) * P(u^2) where u = x-0.5.
38// Coefficients for P are a degree-7 polynomial in u^2
39// approximating sin(pi*(u+0.5)) / (pi*(0.25-u^2))
40// with max error ~2^{-55}, generated by Sollya with:
41// > P = fpminimax(sin(pi*(x+0.5))/(pi*(0.25-x^2)),
42// [|0,2,4,6,8,10,12,14|], [|D...|], [-0.5, 0.5]);
43LIBC_INLINE double lg_sinpi(double x) {
44 constexpr double COEFFS[8] = {0x000000000000001p+2, -0x1.de9e64df22ea4p+1,
45 0x1.472be122401f8p+0, -0x1.d4fcd82df91bp-3,
46 0x1.9f05c97e0aab2p-6, -0x1.f3091c427b611p-10,
47 0x1.b22c9bfdca547p-14, -0x1.15484325ef569p-18};
48 double u = x - 0.5;
49 double u2 = u * u, u4 = u2 * u2, u8 = u4 * u4;
50 double p01 = fputil::multiply_add(x: u2, y: COEFFS[1], z: COEFFS[0]);
51 double p23 = fputil::multiply_add(x: u2, y: COEFFS[3], z: COEFFS[2]);
52 double p45 = fputil::multiply_add(x: u2, y: COEFFS[5], z: COEFFS[4]);
53 double p67 = fputil::multiply_add(x: u2, y: COEFFS[7], z: COEFFS[6]);
54 double p03 = fputil::multiply_add(x: u4, y: p23, z: p01);
55 double p47 = fputil::multiply_add(x: u4, y: p67, z: p45);
56 // Compute (0.25 - u^2) = (0.5 - u) * (0.5 + u) to avoid ~10-digit
57 // catastrophic cancellation when |u| ~ 0.5 (i.e. x near 0 or 1).
58 return (0.5 - u) * (0.5 + u) * fputil::multiply_add(x: u8, y: p47, z: p03);
59}
60
61// Natural logarithm of x (x > 0), using 16-entry table + degree-7 polynomial.
62// Range reduction: x = 2^e * m where m in [1, 2), decomposed as
63// m = (1+i/16)*(1+z) so log(x) = e*log(2) + log(1+i/16) + log(1+z)
64// = e*log(2) + IL[i] + z*P(z).
65// P approximates log(1+z)/z on z in [-1/16, 1/16] with max
66// error ~2^{-54}, generated by Sollya with:
67// > P = fpminimax(log(1+x)/x, [|0,1,2,3,4,5,6,7|], [|D...|], [-1/16, 1/16]);
68LIBC_INLINE double lg_ln(double x) {
69 using FPBits = fputil::FPBits<double>;
70 uint64_t u = FPBits(x).uintval();
71 int e = static_cast<int>(FPBits(x).get_biased_exponent()) - 0x3ff;
72
73 // Coefficients for log(1 + z)/z on z in [-1/16, 1/16]
74 constexpr double COEFFS[8] = {0x1.fffffffffff24p-1, -0x1.ffffffffd1d67p-2,
75 0x1.55555537802dep-2, -0x1.ffffeca81b866p-3,
76 0x1.999611761d772p-3, -0x1.54f3e581b61bfp-3,
77 0x1.1e642b4cb5143p-3, -0x1.9115a5af1e1edp-4};
78 // IL[i] = log(1 + i/16) for i = 0..15
79 constexpr double IL[16] = {
80 0x1.59caeec280116p-57, 0x1.f0a30c01162aap-5, 0x1.e27076e2af2ebp-4,
81 0x1.5ff3070a793d6p-3, 0x1.c8ff7c79a9a2p-3, 0x1.1675cababa60fp-2,
82 0x1.4618bc21c5ec2p-2, 0x1.739d7f6bbd007p-2, 0x1.9f323ecbf984dp-2,
83 0x1.c8ff7c79a9a21p-2, 0x1.f128f5faf06ecp-2, 0x1.0be72e4252a83p-1,
84 0x1.1e85f5e7040d1p-1, 0x1.307d7334f10bep-1, 0x1.41d8fe84672afp-1,
85 0x1.52a2d265bc5abp-1};
86 // IX[i] = 1 / (1 + i/16) for i = 0..15
87 constexpr double IX[16] = {
88 0x000000000000001p+0, 0x1.e1e1e1e1e1e1ep-1, 0x1.c71c71c71c71cp-1,
89 0x1.af286bca1af28p-1, 0x1.999999999999ap-1, 0x1.8618618618618p-1,
90 0x1.745d1745d1746p-1, 0x1.642c8590b2164p-1, 0x1.5555555555555p-1,
91 0x1.47ae147ae147bp-1, 0x1.3b13b13b13b14p-1, 0x1.2f684bda12f68p-1,
92 0x1.2492492492492p-1, 0x1.1a7b9611a7b96p-1, 0x1.1111111111111p-1,
93 0x1.0842108421084p-1};
94
95 int i = static_cast<int>((u >> 48) & 0xf);
96 // Reduce to mantissa in [1, 2)
97 uint64_t mant_u = (u & (~uint64_t(0) >> 12)) | (uint64_t(0x3ff) << 52);
98 double mant = FPBits(mant_u).get_val();
99 double z = IX[i] * mant - 1.0, z2 = z * z, z4 = z2 * z2;
100 double q01 = fputil::multiply_add(x: z, y: COEFFS[1], z: COEFFS[0]);
101 double q23 = fputil::multiply_add(x: z, y: COEFFS[3], z: COEFFS[2]);
102 double q45 = fputil::multiply_add(x: z, y: COEFFS[5], z: COEFFS[4]);
103 double q67 = fputil::multiply_add(x: z, y: COEFFS[7], z: COEFFS[6]);
104 double q03 = fputil::multiply_add(x: z2, y: q23, z: q01);
105 double q47 = fputil::multiply_add(x: z2, y: q67, z: q45);
106 return e * 0x1.62e42fefa39efp-1 + IL[i] +
107 z * fputil::multiply_add(x: z4, y: q47, z: q03);
108}
109
110} // namespace gamma_internal
111
112} // namespace math
113
114} // namespace LIBC_NAMESPACE_DECL
115
116#endif // LLVM_LIBC_SRC___SUPPORT_MATH_GAMMA_UTIL_H
117