glibc
2.9

00001 00002 /* 00003 * IBM Accurate Mathematical Library 00004 * written by International Business Machines Corp. 00005 * Copyright (C) 2001 Free Software Foundation 00006 * 00007 * This program is free software; you can redistribute it and/or modify 00008 * it under the terms of the GNU Lesser General Public License as published by 00009 * the Free Software Foundation; either version 2.1 of the License, or 00010 * (at your option) any later version. 00011 * 00012 * This program is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU Lesser General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU Lesser General Public License 00018 * along with this program; if not, write to the Free Software 00019 * Foundation, Inc., 59 Temple Place  Suite 330, Boston, MA 021111307, USA. 00020 */ 00021 /*************************************************************************/ 00022 /* MODULE_NAME:mpexp.c */ 00023 /* */ 00024 /* FUNCTIONS: mpexp */ 00025 /* */ 00026 /* FILES NEEDED: mpa.h endian.h mpexp.h */ 00027 /* mpa.c */ 00028 /* */ 00029 /* MultiPrecision exponential function subroutine */ 00030 /* ( for p >= 4, 2**(55) <= abs(x) <= 1024 ). */ 00031 /*************************************************************************/ 00032 00033 #include "endian.h" 00034 #include "mpa.h" 00035 #include "mpexp.h" 00036 00037 /* MultiPrecision exponential function subroutine (for p >= 4, */ 00038 /* 2**(55) <= abs(x) <= 1024). */ 00039 void __mpexp(mp_no *x, mp_no *y, int p) { 00040 00041 int i,j,k,m,m1,m2,n; 00042 double a,b; 00043 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, 00044 6,6,6,6,7,7,7,7,8,8,8,8,8}; 00045 static const int m1p[33]= {0,0,0,0,17,23,23,28,27,38,42,39,43,47,43,47,50,54, 00046 57,60,64,67,71,74,68,71,74,77,70,73,76,78,81}; 00047 static const int m1np[7][18] = { 00048 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 00049 { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 00050 { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0}, 00051 { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0}, 00052 { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0}, 00053 { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63}, 00054 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}}; 00055 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, 00056 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, 00057 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; 00058 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, 00059 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, 00060 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; 00061 mp_no mps,mpak,mpt1,mpt2; 00062 00063 /* Choose m,n and compute a=2**(m) */ 00064 n = np[p]; m1 = m1p[p]; a = twomm1[p].d; 00065 for (i=0; i<EX; i++) a *= RADIXI; 00066 for ( ; i>EX; i) a *= RADIX; 00067 b = X[1]*RADIXI; m2 = 24*EX; 00068 for (; b<HALF; m2) { a *= TWO; b *= TWO; } 00069 if (b == HALF) { 00070 for (i=2; i<=p; i++) { if (X[i]!=ZERO) break; } 00071 if (i==p+1) { m2; a *= TWO; } 00072 } 00073 if ((m=m1+m2) <= 0) { 00074 m=0; a=ONE; 00075 for (i=n1; i>0; i,n) { if (m1np[i][p]+m2>0) break; } 00076 } 00077 00078 /* Compute s=x*2**(m). Put result in mps */ 00079 __dbl_mp(a,&mpt1,p); 00080 __mul(x,&mpt1,&mps,p); 00081 00082 /* Evaluate the polynomial. Put result in mpt2 */ 00083 mpone.e=1; mpone.d[0]=ONE; mpone.d[1]=ONE; 00084 mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=nn[n].d; 00085 __dvd(&mps,&mpk,&mpt1,p); 00086 __add(&mpone,&mpt1,&mpak,p); 00087 for (k=n1; k>1; k) { 00088 __mul(&mps,&mpak,&mpt1,p); 00089 mpk.d[1]=nn[k].d; 00090 __dvd(&mpt1,&mpk,&mpt2,p); 00091 __add(&mpone,&mpt2,&mpak,p); 00092 } 00093 __mul(&mps,&mpak,&mpt1,p); 00094 __add(&mpone,&mpt1,&mpt2,p); 00095 00096 /* Raise polynomial value to the power of 2**m. Put result in y */ 00097 for (k=0,j=0; k<m; ) { 00098 __mul(&mpt2,&mpt2,&mpt1,p); k++; 00099 if (k==m) { j=1; break; } 00100 __mul(&mpt1,&mpt1,&mpt2,p); k++; 00101 } 00102 if (j) __cpy(&mpt1,y,p); 00103 else __cpy(&mpt2,y,p); 00104 return; 00105 }