gomp-20050608-branch/libgcc-math/dbl-64/mpexp.c

   1
   2 /*
   3  * IBM Accurate Mathematical Library
   4  * written by International Business Machines Corp.
   5  * Copyright (C) 2001 Free Software Foundation
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU  Lesser General Public License as published by
   9  * the Free Software Foundation; either version 2.1 of the License, or
  10  * (at your option) any later version.
  11  *
  12  * This program is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public License
  18  * along with this program; if not, write to the Free Software
  19  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  20  */
  21 /*************************************************************************/
  22 /*   MODULE_NAME:mpexp.c                                                 */
  23 /*                                                                       */
  24 /*   FUNCTIONS: mpexp                                                    */
  25 /*                                                                       */
  26 /*   FILES NEEDED: mpa.h endian.h mpexp.h                                */
  27 /*                 mpa.c                                                 */
  28 /*                                                                       */
  29 /* Multi-Precision exponential function subroutine                       */
  30 /*   (  for p >= 4, 2**(-55) <= abs(x) <= 1024     ).                    */
  31 /*************************************************************************/
  32
  33 #include "endian.h"
  34 #include "mpa.h"
  35 #include "mpa2.h"
  36 #include "mpexp.h"
  37
  38 /* Multi-Precision exponential function subroutine (for p >= 4,          */
  39 /* 2**(-55) <= abs(x) <= 1024).                                          */
  40 void __mpexp(mp_no *x, mp_no *y, int p) {
  41
  42   int i,j,k,m,m1,m2,n;
  43   double a,b;
  44   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,
  45                              6,6,6,6,7,7,7,7,8,8,8,8,8};
  46   static const int m1p[33]= {0,0,0,0,17,23,23,28,27,38,42,39,43,47,43,47,50,54,
  47                                57,60,64,67,71,74,68,71,74,77,70,73,76,78,81};
  48   static const int m1np[7][18] = {
  49                  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  50                  { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  51                  { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0},
  52                  { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0},
  53                  { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0},
  54                  { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63},
  55                  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}};
  56   mp_no mpone = {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,
  57                     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,
  58                     0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
  59   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,
  60                     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}};
  62   mp_no mps,mpak,mpt1,mpt2;
  63
  64   /* Choose m,n and compute a=2**(-m) */
  65   n = np[p];    m1 = m1p[p];    a = twomm1[p].d;
  66   for (i=0; i<EX; i++)  a *= RADIXI;
  67   for (   ; i>EX; i--)  a *= RADIX;
  68   b = X[1]*RADIXI;   m2 = 24*EX;
  69   for (; b<HALF; m2--)  { a *= TWO;   b *= TWO; }
  70   if (b == HALF) {
  71     for (i=2; i<=p; i++) { if (X[i]!=ZERO)  break; }
  72     if (i==p+1)  { m2--;  a *= TWO; }
  73   }
  74   if ((m=m1+m2) <= 0) {
  75     m=0;  a=ONE;
  76     for (i=n-1; i>0; i--,n--) { if (m1np[i][p]+m2>0)  break; }
  77   }
  78
  79   /* Compute s=x*2**(-m). Put result in mps */
  80   __dbl_mp(a,&mpt1,p);
  81   __mul(x,&mpt1,&mps,p);
  82
  83   /* Evaluate the polynomial. Put result in mpt2 */
  84   mpone.e=1;   mpone.d[0]=ONE;   mpone.d[1]=ONE;
  85   mpk.e = 1;   mpk.d[0] = ONE;   mpk.d[1]=nn[n].d;
  86   __dvd(&mps,&mpk,&mpt1,p);
  87   __add(&mpone,&mpt1,&mpak,p);
  88   for (k=n-1; k>1; k--) {
  89     __mul(&mps,&mpak,&mpt1,p);
  90     mpk.d[1]=nn[k].d;
  91     __dvd(&mpt1,&mpk,&mpt2,p);
  92     __add(&mpone,&mpt2,&mpak,p);
  93   }
  94   __mul(&mps,&mpak,&mpt1,p);
  95   __add(&mpone,&mpt1,&mpt2,p);
  96
  97   /* Raise polynomial value to the power of 2**m. Put result in y */
  98   for (k=0,j=0; k<m; ) {
  99     __mul(&mpt2,&mpt2,&mpt1,p);  k++;
 100     if (k==m)  { j=1;  break; }
 101     __mul(&mpt1,&mpt1,&mpt2,p);  k++;
 102   }
 103   if (j)  __cpy(&mpt1,y,p);
 104   else    __cpy(&mpt2,y,p);
 105   return;
 106 }