sysdeps/ieee754/dbl-64/slowexp.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2015 Free Software Foundation, Inc.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program 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
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  18  */
  19 /**************************************************************************/
  20 /*  MODULE_NAME:slowexp.c                                                 */
  21 /*                                                                        */
  22 /*  FUNCTION:slowexp                                                      */
  23 /*                                                                        */
  24 /*  FILES NEEDED:mpa.h                                                    */
  25 /*               mpa.c mpexp.c                                            */
  26 /*                                                                        */
  27 /*Converting from double precision to Multi-precision and calculating     */
  28 /* e^x                                                                    */
  29 /**************************************************************************/
  30 #include <math_private.h>
  31
  32 #include <stap-probe.h>
  33
  34 #ifndef USE_LONG_DOUBLE_FOR_MP
  35 # include "mpa.h"
  36 void __mpexp (mp_no *x, mp_no *y, int p);
  37 #endif
  38
  39 #ifndef SECTION
  40 # define SECTION
  41 #endif
  42
  43 /*Converting from double precision to Multi-precision and calculating  e^x */
  44 double
  45 SECTION
  46 __slowexp (double x)
  47 {
  48 #ifndef USE_LONG_DOUBLE_FOR_MP
  49   double w, z, res, eps = 3.0e-26;
  50   int p;
  51   mp_no mpx, mpy, mpz, mpw, mpeps, mpcor;
  52
  53   /* Use the multiple precision __MPEXP function to compute the exponential
  54      First at 144 bits and if it is not accurate enough, at 768 bits.  */
  55   p = 6;
  56   __dbl_mp (x, &mpx, p);
  57   __mpexp (&mpx, &mpy, p);
  58   __dbl_mp (eps, &mpeps, p);
  59   __mul (&mpeps, &mpy, &mpcor, p);
  60   __add (&mpy, &mpcor, &mpw, p);
  61   __sub (&mpy, &mpcor, &mpz, p);
  62   __mp_dbl (&mpw, &w, p);
  63   __mp_dbl (&mpz, &z, p);
  64   if (w == z)
  65     {
  66       /* Track how often we get to the slow exp code plus
  67          its input/output values.  */
  68       LIBC_PROBE (slowexp_p6, 2, &x, &w);
  69       return w;
  70     }
  71   else
  72     {
  73       p = 32;
  74       __dbl_mp (x, &mpx, p);
  75       __mpexp (&mpx, &mpy, p);
  76       __mp_dbl (&mpy, &res, p);
  77
  78       /* Track how often we get to the uber-slow exp code plus
  79          its input/output values.  */
  80       LIBC_PROBE (slowexp_p32, 2, &x, &res);
  81       return res;
  82     }
  83 #else
  84   return (double) __ieee754_expl((long double)x);
  85 #endif
  86 }