1/* Single-precision SVE log1p
2
3 Copyright (C) 2023-2024 Free Software Foundation, Inc.
4 This file is part of the GNU C Library.
5
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
10
11 The GNU C Library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with the GNU C Library; if not, see
18 <https://www.gnu.org/licenses/>. */
19
20#include "sv_math.h"
21#include "poly_sve_f32.h"
22
23static const struct data
24{
25 float poly[8];
26 float ln2, exp_bias;
27 uint32_t four, three_quarters;
28} data = {.poly = {/* Do not store first term of polynomial, which is -0.5, as
29 this can be fmov-ed directly instead of including it in
30 the main load-and-mla polynomial schedule. */
31 0x1.5555aap-2f, -0x1.000038p-2f, 0x1.99675cp-3f,
32 -0x1.54ef78p-3f, 0x1.28a1f4p-3f, -0x1.0da91p-3f,
33 0x1.abcb6p-4f, -0x1.6f0d5ep-5f},
34 .ln2 = 0x1.62e43p-1f,
35 .exp_bias = 0x1p-23f,
36 .four = 0x40800000,
37 .three_quarters = 0x3f400000};
38
39#define SignExponentMask 0xff800000
40
41static svfloat32_t NOINLINE
42special_case (svfloat32_t x, svfloat32_t y, svbool_t special)
43{
44 return sv_call_f32 (f: log1pf, x, y, cmp: special);
45}
46
47/* Vector log1pf approximation using polynomial on reduced interval. Worst-case
48 error is 1.27 ULP very close to 0.5.
49 _ZGVsMxv_log1pf(0x1.fffffep-2) got 0x1.9f324p-2
50 want 0x1.9f323ep-2. */
51svfloat32_t SV_NAME_F1 (log1p) (svfloat32_t x, svbool_t pg)
52{
53 const struct data *d = ptr_barrier (&data);
54 /* x < -1, Inf/Nan. */
55 svbool_t special = svcmpeq (pg, svreinterpret_u32 (x), 0x7f800000);
56 special = svorn_z (pg, special, svcmpge (pg, x, -1));
57
58 /* With x + 1 = t * 2^k (where t = m + 1 and k is chosen such that m
59 is in [-0.25, 0.5]):
60 log1p(x) = log(t) + log(2^k) = log1p(m) + k*log(2).
61
62 We approximate log1p(m) with a polynomial, then scale by
63 k*log(2). Instead of doing this directly, we use an intermediate
64 scale factor s = 4*k*log(2) to ensure the scale is representable
65 as a normalised fp32 number. */
66 svfloat32_t m = svadd_x (pg, x, 1);
67
68 /* Choose k to scale x to the range [-1/4, 1/2]. */
69 svint32_t k
70 = svand_x (pg, svsub_x (pg, svreinterpret_s32 (m), d->three_quarters),
71 sv_s32 (SignExponentMask));
72
73 /* Scale x by exponent manipulation. */
74 svfloat32_t m_scale = svreinterpret_f32 (
75 svsub_x (pg, svreinterpret_u32 (x), svreinterpret_u32 (k)));
76
77 /* Scale up to ensure that the scale factor is representable as normalised
78 fp32 number, and scale m down accordingly. */
79 svfloat32_t s = svreinterpret_f32 (svsubr_x (pg, k, d->four));
80 m_scale = svadd_x (pg, m_scale, svmla_x (pg, sv_f32 (x: -1), s, 0.25));
81
82 /* Evaluate polynomial on reduced interval. */
83 svfloat32_t ms2 = svmul_x (pg, m_scale, m_scale),
84 ms4 = svmul_x (pg, ms2, ms2);
85 svfloat32_t p = sv_estrin_7_f32_x (pg, x: m_scale, x2: ms2, x4: ms4, poly: d->poly);
86 p = svmad_x (pg, m_scale, p, -0.5);
87 p = svmla_x (pg, m_scale, m_scale, svmul_x (pg, m_scale, p));
88
89 /* The scale factor to be applied back at the end - by multiplying float(k)
90 by 2^-23 we get the unbiased exponent of k. */
91 svfloat32_t scale_back = svmul_x (pg, svcvt_f32_x (pg, k), d->exp_bias);
92
93 /* Apply the scaling back. */
94 svfloat32_t y = svmla_x (pg, p, scale_back, d->ln2);
95
96 if (__glibc_unlikely (svptest_any (pg, special)))
97 return special_case (x, y, special);
98
99 return y;
100}
101

source code of glibc/sysdeps/aarch64/fpu/log1pf_sve.c