sysdeps/aarch64/fpu/log1pf_sve.c

   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
  23 static 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
  41 static svfloat32_t NOINLINE
  42 special_case (svfloat32_t x, svfloat32_t y, svbool_t special)
  43 {
  44   return sv_call_f32 (log1pf, x, y, 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.  */
  51 svfloat32_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 (-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, m_scale, ms2, ms4, 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 }