newlib/libm/common/exp.c

   1 /* Double-precision e^x function.
   2    Copyright (c) 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 #define N (1 << EXP_TABLE_BITS)
  37 #define InvLn2N __exp_data.invln2N
  38 #define NegLn2hiN __exp_data.negln2hiN
  39 #define NegLn2loN __exp_data.negln2loN
  40 #define Shift __exp_data.shift
  41 #define T __exp_data.tab
  42 #define C2 __exp_data.poly[5 - EXP_POLY_ORDER]
  43 #define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
  44 #define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
  45 #define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
  46 #define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
  47
  48 /* Handle cases that may overflow or underflow when computing the result that
  49    is scale*(1+TMP) without intermediate rounding.  The bit representation of
  50    scale is in SBITS, however it has a computed exponent that may have
  51    overflown into the sign bit so that needs to be adjusted before using it as
  52    a double.  (int32_t)KI is the k used in the argument reduction and exponent
  53    adjustment of scale, positive k here means the result may overflow and
  54    negative k means the result may underflow.  */
  55 static inline double
  56 specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
  57 {
  58   double_t scale, y;
  59
  60   if ((ki & 0x80000000) == 0)
  61     {
  62       /* k > 0, the exponent of scale might have overflowed by <= 460.  */
  63       sbits -= 1009ull << 52;
  64       scale = asdouble (sbits);
  65       y = 0x1p1009 * (scale + scale * tmp);
  66       return check_oflow (y);
  67     }
  68   /* k < 0, need special care in the subnormal range.  */
  69   sbits += 1022ull << 52;
  70   scale = asdouble (sbits);
  71   y = scale + scale * tmp;
  72   if (y < 1.0)
  73     {
  74       /* Round y to the right precision before scaling it into the subnormal
  75          range to avoid double rounding that can cause 0.5+E/2 ulp error where
  76          E is the worst-case ulp error outside the subnormal range.  So this
  77          is only useful if the goal is better than 1 ulp worst-case error.  */
  78       double_t hi, lo;
  79       lo = scale - y + scale * tmp;
  80       hi = 1.0 + y;
  81       lo = 1.0 - hi + y + lo;
  82       y = eval_as_double (hi + lo) - 1.0;
  83       /* Avoid -0.0 with downward rounding.  */
  84       if (WANT_ROUNDING && y == 0.0)
  85         y = 0.0;
  86       /* The underflow exception needs to be signaled explicitly.  */
  87       force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);
  88     }
  89   y = 0x1p-1022 * y;
  90   return check_uflow (y);
  91 }
  92
  93 /* Top 12 bits of a double (sign and exponent bits).  */
  94 static inline uint32_t
  95 top12 (double x)
  96 {
  97   return asuint64 (x) >> 52;
  98 }
  99
 100 double
 101 exp (double x)
 102 {
 103   uint32_t abstop;
 104   uint64_t ki, idx, top, sbits;
 105   /* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
 106   double_t kd, z, r, r2, scale, tail, tmp;
 107
 108   abstop = top12 (x) & 0x7ff;
 109   if (unlikely (abstop - top12 (0x1p-54) >= top12 (512.0) - top12 (0x1p-54)))
 110     {
 111       if (abstop - top12 (0x1p-54) >= 0x80000000)
 112         /* Avoid spurious underflow for tiny x.  */
 113         /* Note: 0 is common input.  */
 114         return WANT_ROUNDING ? 1.0 + x : 1.0;
 115       if (abstop >= top12 (1024.0))
 116         {
 117           if (asuint64 (x) == asuint64 (-INFINITY))
 118             return 0.0;
 119           if (abstop >= top12 (INFINITY))
 120             return 1.0 + x;
 121           if (asuint64 (x) >> 63)
 122             return __math_uflow (0);
 123           else
 124             return __math_oflow (0);
 125         }
 126       /* Large x is special cased below.  */
 127       abstop = 0;
 128     }
 129
 130   /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)].  */
 131   /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N].  */
 132   z = InvLn2N * x;
 133 #if TOINT_INTRINSICS
 134   kd = roundtoint (z);
 135   ki = converttoint (z);
 136 #elif EXP_USE_TOINT_NARROW
 137   /* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes.  */
 138   kd = eval_as_double (z + Shift);
 139   ki = asuint64 (kd) >> 16;
 140   kd = (double_t) (int32_t) ki;
 141 #else
 142   /* z - kd is in [-1, 1] in non-nearest rounding modes.  */
 143   kd = eval_as_double (z + Shift);
 144   ki = asuint64 (kd);
 145   kd -= Shift;
 146 #endif
 147   r = x + kd * NegLn2hiN + kd * NegLn2loN;
 148   /* 2^(k/N) ~= scale * (1 + tail).  */
 149   idx = 2 * (ki % N);
 150   top = ki << (52 - EXP_TABLE_BITS);
 151   tail = asdouble (T[idx]);
 152   /* This is only a valid scale when -1023*N < k < 1024*N.  */
 153   sbits = T[idx + 1] + top;
 154   /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1).  */
 155   /* Evaluation is optimized assuming superscalar pipelined execution.  */
 156   r2 = r * r;
 157   /* Without fma the worst case error is 0.25/N ulp larger.  */
 158   /* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp.  */
 159 #if EXP_POLY_ORDER == 4
 160   tmp = tail + r + r2 * C2 + r * r2 * (C3 + r * C4);
 161 #elif EXP_POLY_ORDER == 5
 162   tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
 163 #elif EXP_POLY_ORDER == 6
 164   tmp = tail + r + r2 * (0.5 + r * C3) + r2 * r2 * (C4 + r * C5 + r2 * C6);
 165 #endif
 166   if (unlikely (abstop == 0))
 167     return specialcase (tmp, sbits, ki);
 168   scale = asdouble (sbits);
 169   /* Note: tmp == 0 or |tmp| > 2^-65 and scale > 2^-739, so there
 170      is no spurious underflow here even without fma.  */
 171   return scale + scale * tmp;
 172 }
 173 #endif