1/*
2 * Single-precision SVE log10 function.
3 *
4 * Copyright (c) 2022-2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8#include "sv_math.h"
9#include "pl_sig.h"
10#include "pl_test.h"
11
12static const struct data
13{
14  float poly_0246[4];
15  float poly_1357[4];
16  float ln2, inv_ln10;
17} data = {
18  .poly_1357 = {
19    /* Coefficients copied from the AdvSIMD routine, then rearranged so that coeffs
20       1, 3, 5 and 7 can be loaded as a single quad-word, hence used with _lane
21       variant of MLA intrinsic.  */
22    0x1.2879c8p-3f, 0x1.6408f8p-4f, 0x1.f0e514p-5f, 0x1.f5f76ap-5f
23  },
24  .poly_0246 = { -0x1.bcb79cp-3f, -0x1.bcd472p-4f, -0x1.246f8p-4f,
25		 -0x1.0fc92cp-4f },
26  .ln2 = 0x1.62e43p-1f,
27  .inv_ln10 = 0x1.bcb7b2p-2f,
28};
29
30#define Min 0x00800000
31#define Max 0x7f800000
32#define Thres 0x7f000000  /* Max - Min.  */
33#define Offset 0x3f2aaaab /* 0.666667.  */
34#define MantissaMask 0x007fffff
35
36static svfloat32_t NOINLINE
37special_case (svfloat32_t x, svfloat32_t y, svbool_t special)
38{
39  return sv_call_f32 (log10f, x, y, special);
40}
41
42/* Optimised implementation of SVE log10f using the same algorithm and
43   polynomial as AdvSIMD log10f.
44   Maximum error is 3.31ulps:
45   SV_NAME_F1 (log10)(0x1.555c16p+0) got 0x1.ffe2fap-4
46				    want 0x1.ffe2f4p-4.  */
47svfloat32_t SV_NAME_F1 (log10) (svfloat32_t x, const svbool_t pg)
48{
49  const struct data *d = ptr_barrier (&data);
50  svuint32_t ix = svreinterpret_u32 (x);
51  svbool_t special = svcmpge (pg, svsub_x (pg, ix, Min), Thres);
52
53  /* x = 2^n * (1+r), where 2/3 < 1+r < 4/3.  */
54  ix = svsub_x (pg, ix, Offset);
55  svfloat32_t n = svcvt_f32_x (
56      pg, svasr_x (pg, svreinterpret_s32 (ix), 23)); /* signextend.  */
57  ix = svand_x (pg, ix, MantissaMask);
58  ix = svadd_x (pg, ix, Offset);
59  svfloat32_t r = svsub_x (pg, svreinterpret_f32 (ix), 1.0f);
60
61  /* y = log10(1+r) + n*log10(2)
62     log10(1+r) ~ r * InvLn(10) + P(r)
63     where P(r) is a polynomial. Use order 9 for log10(1+x), i.e. order 8 for
64     log10(1+x)/x, with x in [-1/3, 1/3] (offset=2/3).  */
65  svfloat32_t r2 = svmul_x (pg, r, r);
66  svfloat32_t r4 = svmul_x (pg, r2, r2);
67  svfloat32_t p_1357 = svld1rq (svptrue_b32 (), &d->poly_1357[0]);
68  svfloat32_t q_01 = svmla_lane (sv_f32 (d->poly_0246[0]), r, p_1357, 0);
69  svfloat32_t q_23 = svmla_lane (sv_f32 (d->poly_0246[1]), r, p_1357, 1);
70  svfloat32_t q_45 = svmla_lane (sv_f32 (d->poly_0246[2]), r, p_1357, 2);
71  svfloat32_t q_67 = svmla_lane (sv_f32 (d->poly_0246[3]), r, p_1357, 3);
72  svfloat32_t q_47 = svmla_x (pg, q_45, r2, q_67);
73  svfloat32_t q_03 = svmla_x (pg, q_01, r2, q_23);
74  svfloat32_t y = svmla_x (pg, q_03, r4, q_47);
75
76  /* Using hi = Log10(2)*n + r*InvLn(10) is faster but less accurate.  */
77  svfloat32_t hi = svmla_x (pg, r, n, d->ln2);
78  hi = svmul_x (pg, hi, d->inv_ln10);
79
80  if (unlikely (svptest_any (pg, special)))
81    return special_case (x, svmla_x (svnot_z (pg, special), hi, r2, y),
82			 special);
83  return svmla_x (pg, hi, r2, y);
84}
85
86PL_SIG (SV, F, 1, log10, 0.01, 11.1)
87PL_TEST_ULP (SV_NAME_F1 (log10), 2.82)
88PL_TEST_INTERVAL (SV_NAME_F1 (log10), -0.0, -0x1p126, 100)
89PL_TEST_INTERVAL (SV_NAME_F1 (log10), 0x1p-149, 0x1p-126, 4000)
90PL_TEST_INTERVAL (SV_NAME_F1 (log10), 0x1p-126, 0x1p-23, 50000)
91PL_TEST_INTERVAL (SV_NAME_F1 (log10), 0x1p-23, 1.0, 50000)
92PL_TEST_INTERVAL (SV_NAME_F1 (log10), 1.0, 100, 50000)
93PL_TEST_INTERVAL (SV_NAME_F1 (log10), 100, inf, 50000)
94