sysdeps/ieee754/dbl-64/mpexp.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2016 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:mpexp.c                                                 */
  21 /*                                                                       */
  22 /*   FUNCTIONS: mpexp                                                    */
  23 /*                                                                       */
  24 /*   FILES NEEDED: mpa.h endian.h mpexp.h                                */
  25 /*                 mpa.c                                                 */
  26 /*                                                                       */
  27 /* Multi-Precision exponential function subroutine                       */
  28 /*   (  for p >= 4, 2**(-55) <= abs(x) <= 1024     ).                    */
  29 /*************************************************************************/
  30
  31 #include "endian.h"
  32 #include "mpa.h"
  33 #include <assert.h>
  34
  35 #ifndef SECTION
  36 # define SECTION
  37 #endif
  38
  39 /* Multi-Precision exponential function subroutine (for p >= 4,
  40    2**(-55) <= abs(x) <= 1024).  */
  41 void
  42 SECTION
  43 __mpexp (mp_no *x, mp_no *y, int p)
  44 {
  45   int i, j, k, m, m1, m2, n;
  46   mantissa_t b;
  47   static const int np[33] =
  48     {
  49       0, 0, 0, 0, 3, 3, 4, 4, 5, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6,
  50       6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8
  51     };
  52
  53   static const int m1p[33] =
  54     {
  55       0, 0, 0, 0,
  56       17, 23, 23, 28,
  57       27, 38, 42, 39,
  58       43, 47, 43, 47,
  59       50, 54, 57, 60,
  60       64, 67, 71, 74,
  61       68, 71, 74, 77,
  62       70, 73, 76, 78,
  63       81
  64     };
  65   static const int m1np[7][18] =
  66     {
  67       {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  68       {0, 0, 0, 0, 36, 48, 60, 72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  69       {0, 0, 0, 0, 24, 32, 40, 48, 56, 64, 72, 0, 0, 0, 0, 0, 0, 0},
  70       {0, 0, 0, 0, 17, 23, 29, 35, 41, 47, 53, 59, 65, 0, 0, 0, 0, 0},
  71       {0, 0, 0, 0, 0, 0, 23, 28, 33, 38, 42, 47, 52, 57, 62, 66, 0, 0},
  72       {0, 0, 0, 0, 0, 0, 0, 0, 27, 0, 0, 39, 43, 47, 51, 55, 59, 63},
  73       {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 43, 47, 50, 54}
  74     };
  75   mp_no mps, mpk, mpt1, mpt2;
  76
  77   /* Choose m,n and compute a=2**(-m).  */
  78   n = np[p];
  79   m1 = m1p[p];
  80   b = X[1];
  81   m2 = 24 * EX;
  82   for (; b < HALFRAD; m2--)
  83     b *= 2;
  84   if (b == HALFRAD)
  85     {
  86       for (i = 2; i <= p; i++)
  87         {
  88           if (X[i] != 0)
  89             break;
  90         }
  91       if (i == p + 1)
  92         m2--;
  93     }
  94
  95   m = m1 + m2;
  96   if (__glibc_unlikely (m <= 0))
  97     {
  98       /* The m1np array which is used to determine if we can reduce the
  99          polynomial expansion iterations, has only 18 elements.  Besides,
 100          numbers smaller than those required by p >= 18 should not come here
 101          at all since the fast phase of exp returns 1.0 for anything less
 102          than 2^-55.  */
 103       assert (p < 18);
 104       m = 0;
 105       for (i = n - 1; i > 0; i--, n--)
 106         if (m1np[i][p] + m2 > 0)
 107           break;
 108     }
 109
 110   /* Compute s=x*2**(-m). Put result in mps.  This is the range-reduced input
 111      that we will use to compute e^s.  For the final result, simply raise it
 112      to 2^m.  */
 113   __pow_mp (-m, &mpt1, p);
 114   __mul (x, &mpt1, &mps, p);
 115
 116   /* Compute the Taylor series for e^s:
 117
 118          1 + x/1! + x^2/2! + x^3/3! ...
 119
 120      for N iterations.  We compute this as:
 121
 122          e^x = 1 + (x * n!/1! + x^2 * n!/2! + x^3 * n!/3!) / n!
 123              = 1 + (x * (n!/1! + x * (n!/2! + x * (n!/3! + x ...)))) / n!
 124
 125      k! is computed on the fly as KF and at the end of the polynomial loop, KF
 126      is n!, which can be used directly.  */
 127   __cpy (&mps, &mpt2, p);
 128
 129   double kf = 1.0;
 130
 131   /* Evaluate the rest.  The result will be in mpt2.  */
 132   for (k = n - 1; k > 0; k--)
 133     {
 134       /* n! / k! = n * (n - 1) ... * (n - k + 1) */
 135       kf *= k + 1;
 136
 137       __dbl_mp (kf, &mpk, p);
 138       __add (&mpt2, &mpk, &mpt1, p);
 139       __mul (&mps, &mpt1, &mpt2, p);
 140     }
 141   __dbl_mp (kf, &mpk, p);
 142   __dvd (&mpt2, &mpk, &mpt1, p);
 143   __add (&__mpone, &mpt1, &mpt2, p);
 144
 145   /* Raise polynomial value to the power of 2**m. Put result in y.  */
 146   for (k = 0, j = 0; k < m;)
 147     {
 148       __sqr (&mpt2, &mpt1, p);
 149       k++;
 150       if (k == m)
 151         {
 152           j = 1;
 153           break;
 154         }
 155       __sqr (&mpt1, &mpt2, p);
 156       k++;
 157     }
 158   if (j)
 159     __cpy (&mpt1, y, p);
 160   else
 161     __cpy (&mpt2, y, p);
 162   return;
 163 }