sysdeps/ieee754/dbl-64/mpexp.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2013 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 "mpexp.h"
  34 #include <assert.h>
  35
  36 #ifndef SECTION
  37 # define SECTION
  38 #endif
  39
  40 /* Multi-Precision exponential function subroutine (for p >= 4,          */
  41 /* 2**(-55) <= abs(x) <= 1024).                                          */
  42 void
  43 SECTION
  44 __mpexp(mp_no *x, mp_no *y, int p) {
  45
  46   int i,j,k,m,m1,m2,n;
  47   double a,b;
  48   static const int np[33] = {0,0,0,0,3,3,4,4,5,4,4,5,5,5,6,6,6,6,6,6,
  49                              6,6,6,6,7,7,7,7,8,8,8,8,8};
  50   static const int m1p[33]= {0,0,0,0,17,23,23,28,27,38,42,39,43,47,43,47,50,54,
  51                                57,60,64,67,71,74,68,71,74,77,70,73,76,78,81};
  52   static const int m1np[7][18] = {
  53                  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  54                  { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  55                  { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0},
  56                  { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0},
  57                  { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0},
  58                  { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63},
  59                  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}};
  60   mp_no mpk   = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
  61                     0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
  62                     0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
  63   mp_no mps,mpak,mpt1,mpt2;
  64
  65   /* Choose m,n and compute a=2**(-m) */
  66   n = np[p];    m1 = m1p[p];    a = __mpexp_twomm1[p].d;
  67   for (i=0; i<EX; i++)  a *= RADIXI;
  68   for (   ; i>EX; i--)  a *= RADIX;
  69   b = X[1]*RADIXI;   m2 = 24*EX;
  70   for (; b<HALF; m2--)  { a *= TWO;   b *= TWO; }
  71   if (b == HALF) {
  72     for (i=2; i<=p; i++) { if (X[i]!=ZERO)  break; }
  73     if (i==p+1)  { m2--;  a *= TWO; }
  74   }
  75
  76   m = m1 + m2;
  77   if (__glibc_unlikely (m <= 0))
  78     {
  79       /* The m1np array which is used to determine if we can reduce the
  80          polynomial expansion iterations, has only 18 elements.  Besides,
  81          numbers smaller than those required by p >= 18 should not come here
  82          at all since the fast phase of exp returns 1.0 for anything less
  83          than 2^-55.  */
  84       assert (p < 18);
  85       m = 0;
  86       a = ONE;
  87       for (i = n - 1; i > 0; i--, n--)
  88         if (m1np[i][p] + m2 > 0)
  89           break;
  90     }
  91
  92   /* Compute s=x*2**(-m). Put result in mps */
  93   __dbl_mp(a,&mpt1,p);
  94   __mul(x,&mpt1,&mps,p);
  95
  96   /* Evaluate the polynomial. Put result in mpt2 */
  97   mpk.e = 1;   mpk.d[0] = ONE;   mpk.d[1]=n;
  98   __dvd(&mps,&mpk,&mpt1,p);
  99   __add(&mpone,&mpt1,&mpak,p);
 100   for (k=n-1; k>1; k--) {
 101     __mul(&mps,&mpak,&mpt1,p);
 102     mpk.d[1] = k;
 103     __dvd(&mpt1,&mpk,&mpt2,p);
 104     __add(&mpone,&mpt2,&mpak,p);
 105   }
 106   __mul(&mps,&mpak,&mpt1,p);
 107   __add(&mpone,&mpt1,&mpt2,p);
 108
 109   /* Raise polynomial value to the power of 2**m. Put result in y */
 110   for (k=0,j=0; k<m; ) {
 111     __mul(&mpt2,&mpt2,&mpt1,p);  k++;
 112     if (k==m)  { j=1;  break; }
 113     __mul(&mpt1,&mpt1,&mpt2,p);  k++;
 114   }
 115   if (j)  __cpy(&mpt1,y,p);
 116   else    __cpy(&mpt2,y,p);
 117   return;
 118 }