sysdeps/aarch64/fpu/expm1f_advsimd.c

   1 /* Single-precision AdvSIMD expm1
   2
   3    Copyright (C) 2023 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 "v_math.h"
  21 #include "poly_advsimd_f32.h"
  22
  23 static const struct data
  24 {
  25   float32x4_t poly[5];
  26   float32x4_t invln2, ln2_lo, ln2_hi, shift;
  27   int32x4_t exponent_bias;
  28 #if WANT_SIMD_EXCEPT
  29   uint32x4_t thresh;
  30 #else
  31   float32x4_t oflow_bound;
  32 #endif
  33 } data = {
  34   /* Generated using fpminimax with degree=5 in [-log(2)/2, log(2)/2].  */
  35   .poly = { V4 (0x1.fffffep-2), V4 (0x1.5554aep-3), V4 (0x1.555736p-5),
  36             V4 (0x1.12287cp-7), V4 (0x1.6b55a2p-10) },
  37   .invln2 = V4 (0x1.715476p+0f),
  38   .ln2_hi = V4 (0x1.62e4p-1f),
  39   .ln2_lo = V4 (0x1.7f7d1cp-20f),
  40   .shift = V4 (0x1.8p23f),
  41   .exponent_bias = V4 (0x3f800000),
  42 #if !WANT_SIMD_EXCEPT
  43   /* Value above which expm1f(x) should overflow. Absolute value of the
  44      underflow bound is greater than this, so it catches both cases - there is
  45      a small window where fallbacks are triggered unnecessarily.  */
  46   .oflow_bound = V4 (0x1.5ebc4p+6),
  47 #else
  48   /* asuint(oflow_bound) - asuint(0x1p-23), shifted left by 1 for absolute
  49      compare.  */
  50   .thresh = V4 (0x1d5ebc40),
  51 #endif
  52 };
  53
  54 /* asuint(0x1p-23), shifted by 1 for abs compare.  */
  55 #define TinyBound v_u32 (0x34000000 << 1)
  56
  57 static float32x4_t VPCS_ATTR NOINLINE
  58 special_case (float32x4_t x, float32x4_t y, uint32x4_t special)
  59 {
  60   return v_call_f32 (expm1f, x, y, special);
  61 }
  62
  63 /* Single-precision vector exp(x) - 1 function.
  64    The maximum error is 1.51 ULP:
  65    _ZGVnN4v_expm1f (0x1.8baa96p-2) got 0x1.e2fb9p-2
  66                                   want 0x1.e2fb94p-2.  */
  67 float32x4_t VPCS_ATTR V_NAME_F1 (expm1) (float32x4_t x)
  68 {
  69   const struct data *d = ptr_barrier (&data);
  70   uint32x4_t ix = vreinterpretq_u32_f32 (x);
  71
  72 #if WANT_SIMD_EXCEPT
  73   /* If fp exceptions are to be triggered correctly, fall back to scalar for
  74      |x| < 2^-23, |x| > oflow_bound, Inf & NaN. Add ix to itself for
  75      shift-left by 1, and compare with thresh which was left-shifted offline -
  76      this is effectively an absolute compare.  */
  77   uint32x4_t special
  78       = vcgeq_u32 (vsubq_u32 (vaddq_u32 (ix, ix), TinyBound), d->thresh);
  79   if (__glibc_unlikely (v_any_u32 (special)))
  80     x = v_zerofy_f32 (x, special);
  81 #else
  82   /* Handles very large values (+ve and -ve), +/-NaN, +/-Inf.  */
  83   uint32x4_t special = vceqzq_u32 (vcaltq_f32 (x, d->oflow_bound));
  84 #endif
  85
  86   /* Reduce argument to smaller range:
  87      Let i = round(x / ln2)
  88      and f = x - i * ln2, then f is in [-ln2/2, ln2/2].
  89      exp(x) - 1 = 2^i * (expm1(f) + 1) - 1
  90      where 2^i is exact because i is an integer.  */
  91   float32x4_t j = vsubq_f32 (vfmaq_f32 (d->shift, d->invln2, x), d->shift);
  92   int32x4_t i = vcvtq_s32_f32 (j);
  93   float32x4_t f = vfmsq_f32 (x, j, d->ln2_hi);
  94   f = vfmsq_f32 (f, j, d->ln2_lo);
  95
  96   /* Approximate expm1(f) using polynomial.
  97      Taylor expansion for expm1(x) has the form:
  98          x + ax^2 + bx^3 + cx^4 ....
  99      So we calculate the polynomial P(f) = a + bf + cf^2 + ...
 100      and assemble the approximation expm1(f) ~= f + f^2 * P(f).  */
 101   float32x4_t p = v_horner_4_f32 (f, d->poly);
 102   p = vfmaq_f32 (f, vmulq_f32 (f, f), p);
 103
 104   /* Assemble the result.
 105      expm1(x) ~= 2^i * (p + 1) - 1
 106      Let t = 2^i.  */
 107   int32x4_t u = vaddq_s32 (vshlq_n_s32 (i, 23), d->exponent_bias);
 108   float32x4_t t = vreinterpretq_f32_s32 (u);
 109
 110   if (__glibc_unlikely (v_any_u32 (special)))
 111     return special_case (vreinterpretq_f32_u32 (ix),
 112                          vfmaq_f32 (vsubq_f32 (t, v_f32 (1.0f)), p, t),
 113                          special);
 114
 115   /* expm1(x) ~= p * t + (t - 1).  */
 116   return vfmaq_f32 (vsubq_f32 (t, v_f32 (1.0f)), p, t);
 117 }