Home | History | Annotate | Download | only in single 
      1 /*
      2  * e_expf.c - single-precision exp function
      3  *
      4  * Copyright (c) 2009-2018, Arm Limited.
      5  * SPDX-License-Identifier: MIT
      6  */
      7 
      8 /*
      9  * Algorithm was once taken from Cody & Waite, but has been munged
     10  * out of all recognition by SGT.
     11  */
     12 
     13 #include <math.h>
     14 #include <errno.h>
     15 #include "math_private.h"
     16 
     17 float
     18 expf(float X)
     19 {
     20   int N; float XN, g, Rg, Result;
     21   unsigned ix = fai(X), edgecaseflag = 0;
     22 
     23   /*
     24    * Handle infinities, NaNs and big numbers.
     25    */
     26   if (__builtin_expect((ix << 1) - 0x67000000 > 0x85500000 - 0x67000000, 0)) {
     27     if (!(0x7f800000 & ~ix)) {
     28       if (ix == 0xff800000)
     29         return 0.0f;
     30       else
     31         return FLOAT_INFNAN(X);/* do the right thing with both kinds of NaN and with +inf */
     32     } else if ((ix << 1) < 0x67000000) {
     33       return 1.0f;               /* magnitude so small the answer can't be distinguished from 1 */
     34     } else if ((ix << 1) > 0x85a00000) {
     35       __set_errno(ERANGE);
     36       if (ix & 0x80000000) {
     37         return FLOAT_UNDERFLOW;
     38       } else {
     39         return FLOAT_OVERFLOW;
     40       }
     41     } else {
     42       edgecaseflag = 1;
     43     }
     44   }
     45 
     46   /*
     47    * Split the input into an integer multiple of log(2)/4, and a
     48    * fractional part.
     49    */
     50   XN = X * (4.0f*1.4426950408889634074f);
     51 #ifdef __TARGET_FPU_SOFTVFP
     52   XN = _frnd(XN);
     53   N = (int)XN;
     54 #else
     55   N = (int)(XN + (ix & 0x80000000 ? -0.5f : 0.5f));
     56   XN = N;
     57 #endif
     58   g = (X - XN * 0x1.62ep-3F) - XN * 0x1.0bfbe8p-17F;  /* prec-and-a-half representation of log(2)/4 */
     59 
     60   /*
     61    * Now we compute exp(X) in, conceptually, three parts:
     62    *  - a pure power of two which we get from N>>2
     63    *  - exp(g) for g in [-log(2)/8,+log(2)/8], which we compute
     64    *    using a Remez-generated polynomial approximation
     65    *  - exp(k*log(2)/4) (aka 2^(k/4)) for k in [0..3], which we
     66    *    get from a lookup table in precision-and-a-half and
     67    *    multiply by g.
     68    *
     69    * We gain a bit of extra precision by the fact that actually
     70    * our polynomial approximation gives us exp(g)-1, and we add
     71    * the 1 back on by tweaking the prec-and-a-half multiplication
     72    * step.
     73    *
     74    * Coefficients generated by the command
     75 
     76 ./auxiliary/remez.jl --variable=g --suffix=f -- '-log(BigFloat(2))/8' '+log(BigFloat(2))/8' 3 0 '(expm1(x))/x'
     77 
     78   */
     79   Rg = g * (
     80             9.999999412829185331953781321128516523408059996430919985217971370689774264850229e-01f+g*(4.999999608551332693833317084753864837160947932961832943901913087652889900683833e-01f+g*(1.667292360203016574303631953046104769969440903672618034272397630620346717392378e-01f+g*(4.168230895653321517750133783431970715648192153539929404872173693978116154823859e-02f)))
     81             );
     82 
     83   /*
     84    * Do the table lookup and combine it with Rg, to get our final
     85    * answer apart from the exponent.
     86    */
     87   {
     88     static const float twotokover4top[4] = { 0x1p+0F, 0x1.306p+0F, 0x1.6ap+0F, 0x1.ae8p+0F };
     89     static const float twotokover4bot[4] = { 0x0p+0F, 0x1.fc1464p-13F, 0x1.3cccfep-13F, 0x1.3f32b6p-13F };
     90     static const float twotokover4all[4] = { 0x1p+0F, 0x1.306fep+0F, 0x1.6a09e6p+0F, 0x1.ae89fap+0F };
     91     int index = (N & 3);
     92     Rg = twotokover4top[index] + (twotokover4bot[index] + twotokover4all[index]*Rg);
     93     N >>= 2;
     94   }
     95 
     96   /*
     97    * Combine the output exponent and mantissa, and return.
     98    */
     99   if (__builtin_expect(edgecaseflag, 0)) {
    100     Result = fhex(((N/2) << 23) + 0x3f800000);
    101     Result *= Rg;
    102     Result *= fhex(((N-N/2) << 23) + 0x3f800000);
    103     /*
    104      * Step not mentioned in C&W: set errno reliably.
    105      */
    106     if (fai(Result) == 0)
    107       return MATHERR_EXPF_UFL(Result);
    108     if (fai(Result) == 0x7f800000)
    109       return MATHERR_EXPF_OFL(Result);
    110     return FLOAT_CHECKDENORM(Result);
    111   } else {
    112     Result = fhex(N * 8388608.0f + (float)0x3f800000);
    113     Result *= Rg;
    114   }
    115 
    116   return Result;
    117 }
    118