1  /* Single-precision 2^x function.
       2     Copyright (C) 2017-2023 Free Software Foundation, Inc.
       3     This file is part of the GNU C Library.
       4  
       5     The GNU C Library is free software; you can redistribute it and/or
       6     modify it under the terms of the GNU Lesser General Public
       7     License as published by the Free Software Foundation; either
       8     version 2.1 of the License, or (at your option) any later version.
       9  
      10     The GNU C Library is distributed in the hope that it will be useful,
      11     but WITHOUT ANY WARRANTY; without even the implied warranty of
      12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      13     Lesser General Public License for more details.
      14  
      15     You should have received a copy of the GNU Lesser General Public
      16     License along with the GNU C Library; if not, see
      17     <https://www.gnu.org/licenses/>.  */
      18  
      19  #include <math.h>
      20  #include <math-narrow-eval.h>
      21  #include <stdint.h>
      22  #include <libm-alias-finite.h>
      23  #include <libm-alias-float.h>
      24  #include "math_config.h"
      25  
      26  /*
      27  EXP2F_TABLE_BITS = 5
      28  EXP2F_POLY_ORDER = 3
      29  
      30  ULP error: 0.502 (nearest rounding.)
      31  Relative error: 1.69 * 2^-34 in [-1/64, 1/64] (before rounding.)
      32  Wrong count: 168353 (all nearest rounding wrong results with fma.)
      33  Non-nearest ULP error: 1 (rounded ULP error)
      34  */
      35  
      36  #define N (1 << EXP2F_TABLE_BITS)
      37  #define T __exp2f_data.tab
      38  #define C __exp2f_data.poly
      39  #define SHIFT __exp2f_data.shift_scaled
      40  
      41  static inline uint32_t
      42  top12 (float x)
      43  {
      44    return asuint (x) >> 20;
      45  }
      46  
      47  float
      48  __exp2f (float x)
      49  {
      50    uint32_t abstop;
      51    uint64_t ki, t;
      52    /* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
      53    double_t kd, xd, z, r, r2, y, s;
      54  
      55    xd = (double_t) x;
      56    abstop = top12 (x) & 0x7ff;
      57    if (__glibc_unlikely (abstop >= top12 (128.0f)))
      58      {
      59        /* |x| >= 128 or x is nan.  */
      60        if (asuint (x) == asuint (-INFINITY))
      61  	return 0.0f;
      62        if (abstop >= top12 (INFINITY))
      63  	return x + x;
      64        if (x > 0.0f)
      65  	return __math_oflowf (0);
      66        if (x <= -150.0f)
      67  	return __math_uflowf (0);
      68  #if WANT_ERRNO_UFLOW
      69        if (x < -149.0f)
      70  	return __math_may_uflowf (0);
      71  #endif
      72      }
      73  
      74    /* x = k/N + r with r in [-1/(2N), 1/(2N)] and int k.  */
      75    kd = math_narrow_eval ((double) (xd + SHIFT)); /* Needs to be double.  */
      76    ki = asuint64 (kd);
      77    kd -= SHIFT; /* k/N for int k.  */
      78    r = xd - kd;
      79  
      80    /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */
      81    t = T[ki % N];
      82    t += ki << (52 - EXP2F_TABLE_BITS);
      83    s = asdouble (t);
      84    z = C[0] * r + C[1];
      85    r2 = r * r;
      86    y = C[2] * r + 1;
      87    y = z * r2 + y;
      88    y = y * s;
      89    return (float) y;
      90  }
      91  #ifndef __exp2f
      92  strong_alias (__exp2f, __ieee754_exp2f)
      93  libm_alias_finite (__ieee754_exp2f, __exp2f)
      94  versioned_symbol (libm, __exp2f, exp2f, GLIBC_2_27);
      95  libm_alias_float_other (__exp2, exp2)
      96  #endif