newlib/libm/common/sf_pow.c

   1 /* Single-precision pow function.
   2    Copyright (c) 2017-2018 Arm Ltd.  All rights reserved.
   3
   4    SPDX-License-Identifier: BSD-3-Clause
   5
   6    Redistribution and use in source and binary forms, with or without
   7    modification, are permitted provided that the following conditions
   8    are met:
   9    1. Redistributions of source code must retain the above copyright
  10       notice, this list of conditions and the following disclaimer.
  11    2. Redistributions in binary form must reproduce the above copyright
  12       notice, this list of conditions and the following disclaimer in the
  13       documentation and/or other materials provided with the distribution.
  14    3. The name of the company may not be used to endorse or promote
  15       products derived from this software without specific prior written
  16       permission.
  17
  18    THIS SOFTWARE IS PROVIDED BY ARM LTD ``AS IS'' AND ANY EXPRESS OR IMPLIED
  19    WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  20    MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  21    IN NO EVENT SHALL ARM LTD BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  22    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
  23    TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  24    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  25    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  26    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  27    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
  28
  29 #include "fdlibm.h"
  30 #if !__OBSOLETE_MATH
  31
  32 #include <math.h>
  33 #include <stdint.h>
  34 #include "math_config.h"
  35
  36 /*
  37 POWF_LOG2_POLY_ORDER = 5
  38 EXP2F_TABLE_BITS = 5
  39
  40 ULP error: 0.82 (~ 0.5 + relerr*2^24)
  41 relerr: 1.27 * 2^-26 (Relative error ~= 128*Ln2*relerr_log2 + relerr_exp2)
  42 relerr_log2: 1.83 * 2^-33 (Relative error of logx.)
  43 relerr_exp2: 1.69 * 2^-34 (Relative error of exp2(ylogx).)
  44 */
  45
  46 #define N (1 << POWF_LOG2_TABLE_BITS)
  47 #define T __powf_log2_data.tab
  48 #define A __powf_log2_data.poly
  49 #define OFF 0x3f330000
  50
  51 /* Subnormal input is normalized so ix has negative biased exponent.
  52    Output is multiplied by N (POWF_SCALE) if TOINT_INTRINICS is set.  */
  53 static inline double_t
  54 log2_inline (uint32_t ix)
  55 {
  56   /* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
  57   double_t z, r, r2, r4, p, q, y, y0, invc, logc;
  58   uint32_t iz, top, tmp;
  59   int k, i;
  60
  61   /* x = 2^k z; where z is in range [OFF,2*OFF] and exact.
  62      The range is split into N subintervals.
  63      The ith subinterval contains z and c is near its center.  */
  64   tmp = ix - OFF;
  65   i = (tmp >> (23 - POWF_LOG2_TABLE_BITS)) % N;
  66   top = tmp & 0xff800000;
  67   iz = ix - top;
  68   k = (int32_t) top >> (23 - POWF_SCALE_BITS); /* arithmetic shift */
  69   invc = T[i].invc;
  70   logc = T[i].logc;
  71   z = (double_t) asfloat (iz);
  72
  73   /* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */
  74   r = z * invc - 1;
  75   y0 = logc + (double_t) k;
  76
  77   /* Pipelined polynomial evaluation to approximate log1p(r)/ln2.  */
  78   r2 = r * r;
  79   y = A[0] * r + A[1];
  80   p = A[2] * r + A[3];
  81   r4 = r2 * r2;
  82   q = A[4] * r + y0;
  83   q = p * r2 + q;
  84   y = y * r4 + q;
  85   return y;
  86 }
  87
  88 #undef N
  89 #undef T
  90 #define N (1 << EXP2F_TABLE_BITS)
  91 #define T __exp2f_data.tab
  92 #define SIGN_BIAS (1 << (EXP2F_TABLE_BITS + 11))
  93
  94 /* The output of log2 and thus the input of exp2 is either scaled by N
  95    (in case of fast toint intrinsics) or not.  The unscaled xd must be
  96    in [-1021,1023], sign_bias sets the sign of the result.  */
  97 static inline double_t
  98 exp2_inline (double_t xd, uint32_t sign_bias)
  99 {
 100   uint64_t ki, ski, t;
 101   /* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
 102   double_t kd, z, r, r2, y, s;
 103
 104 #if TOINT_INTRINSICS
 105 # define C __exp2f_data.poly_scaled
 106   /* N*x = k + r with r in [-1/2, 1/2] */
 107   kd = roundtoint (xd); /* k */
 108   ki = converttoint (xd);
 109 #else
 110 # define C __exp2f_data.poly
 111 # define SHIFT __exp2f_data.shift_scaled
 112   /* x = k/N + r with r in [-1/(2N), 1/(2N)] */
 113   kd = (double) (xd + SHIFT); /* Rounding to double precision is required.  */
 114   ki = asuint64 (kd);
 115   kd -= SHIFT; /* k/N */
 116 #endif
 117   r = xd - kd;
 118
 119   /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */
 120   t = T[ki % N];
 121   ski = ki + sign_bias;
 122   t += ski << (52 - EXP2F_TABLE_BITS);
 123   s = asdouble (t);
 124   z = C[0] * r + C[1];
 125   r2 = r * r;
 126   y = C[2] * r + 1;
 127   y = z * r2 + y;
 128   y = y * s;
 129   return y;
 130 }
 131
 132 /* Returns 0 if not int, 1 if odd int, 2 if even int.  The argument is
 133    the bit representation of a non-zero finite floating-point value.  */
 134 static inline int
 135 checkint (uint32_t iy)
 136 {
 137   int e = iy >> 23 & 0xff;
 138   if (e < 0x7f)
 139     return 0;
 140   if (e > 0x7f + 23)
 141     return 2;
 142   if (iy & ((1 << (0x7f + 23 - e)) - 1))
 143     return 0;
 144   if (iy & (1 << (0x7f + 23 - e)))
 145     return 1;
 146   return 2;
 147 }
 148
 149 static inline int
 150 zeroinfnan (uint32_t ix)
 151 {
 152   return 2 * ix - 1 >= 2u * 0x7f800000 - 1;
 153 }
 154
 155 float
 156 powf (float x, float y)
 157 {
 158   uint32_t sign_bias = 0;
 159   uint32_t ix, iy;
 160
 161   ix = asuint (x);
 162   iy = asuint (y);
 163   if (__builtin_expect (ix - 0x00800000 >= 0x7f800000 - 0x00800000
 164                           || zeroinfnan (iy),
 165                         0))
 166     {
 167       /* Either (x < 0x1p-126 or inf or nan) or (y is 0 or inf or nan).  */
 168       if (__builtin_expect (zeroinfnan (iy), 0))
 169         {
 170           if (2 * iy == 0)
 171             return issignalingf_inline (x) ? x + y : 1.0f;
 172           if (ix == 0x3f800000)
 173             return issignalingf_inline (y) ? x + y : 1.0f;
 174           if (2 * ix > 2u * 0x7f800000 || 2 * iy > 2u * 0x7f800000)
 175             return x + y;
 176           if (2 * ix == 2 * 0x3f800000)
 177             return 1.0f;
 178           if ((2 * ix < 2 * 0x3f800000) == !(iy & 0x80000000))
 179             return 0.0f; /* |x|<1 && y==inf or |x|>1 && y==-inf.  */
 180           return y * y;
 181         }
 182       if (__builtin_expect (zeroinfnan (ix), 0))
 183         {
 184           float_t x2 = x * x;
 185           if (ix & 0x80000000 && checkint (iy) == 1)
 186             {
 187               x2 = -x2;
 188               sign_bias = 1;
 189             }
 190 #if WANT_ERRNO
 191           if (2 * ix == 0 && iy & 0x80000000)
 192             return __math_divzerof (sign_bias);
 193 #endif
 194           return iy & 0x80000000 ? 1 / x2 : x2;
 195         }
 196       /* x and y are non-zero finite.  */
 197       if (ix & 0x80000000)
 198         {
 199           /* Finite x < 0.  */
 200           int yint = checkint (iy);
 201           if (yint == 0)
 202             return __math_invalidf (x);
 203           if (yint == 1)
 204             sign_bias = SIGN_BIAS;
 205           ix &= 0x7fffffff;
 206         }
 207       if (ix < 0x00800000)
 208         {
 209           /* Normalize subnormal x so exponent becomes negative.  */
 210           ix = asuint (x * 0x1p23f);
 211           ix &= 0x7fffffff;
 212           ix -= 23 << 23;
 213         }
 214     }
 215   double_t logx = log2_inline (ix);
 216   double_t ylogx = y * logx; /* Note: cannot overflow, y is single prec.  */
 217   if (__builtin_expect ((asuint64 (ylogx) >> 47 & 0xffff)
 218                           >= asuint64 (126.0 * POWF_SCALE) >> 47,
 219                         0))
 220     {
 221       /* |y*log(x)| >= 126.  */
 222       if (ylogx > 0x1.fffffffd1d571p+6 * POWF_SCALE)
 223         /* |x^y| > 0x1.ffffffp127.  */
 224         return __math_oflowf (sign_bias);
 225       if (WANT_ROUNDING && WANT_ERRNO
 226           && ylogx > 0x1.fffffffa3aae2p+6 * POWF_SCALE)
 227         /* |x^y| > 0x1.fffffep127, check if we round away from 0.  */
 228         if ((!sign_bias
 229              && eval_as_float (1.0f + opt_barrier_float (0x1p-25f)) != 1.0f)
 230             || (sign_bias
 231                 && eval_as_float (-1.0f - opt_barrier_float (0x1p-25f))
 232                      != -1.0f))
 233           return __math_oflowf (sign_bias);
 234       if (ylogx <= -150.0 * POWF_SCALE)
 235         return __math_uflowf (sign_bias);
 236 #if WANT_ERRNO_UFLOW
 237       if (ylogx < -149.0 * POWF_SCALE)
 238         return __math_may_uflowf (sign_bias);
 239 #endif
 240     }
 241   return (float) exp2_inline (ylogx, sign_bias);
 242 }
 243 #endif /* !__OBSOLETE_MATH */