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 /* */ 00023 /* MODULE_NAME:mplog.c */ 00024 /* */ 00025 /* FUNCTIONS: mplog */ 00026 /* */ 00027 /* FILES NEEDED: endian.h mpa.h mplog.h */ 00028 /* mpexp.c */ 00029 /* */ 00030 /* MultiPrecision logarithm function subroutine (for precision p >= 4, */ 00031 /* 2**(1024) < x < 2**1024) and x is outside of the interval */ 00032 /* [12**(54),1+2**(54)]. Upon entry, x should be set to the */ 00033 /* multiprecision value of the input and y should be set into a multi */ 00034 /* precision value of an approximation of log(x) with relative error */ 00035 /* bound of at most 2**(52). The routine improves the accuracy of y. */ 00036 /* */ 00037 /************************************************************************/ 00038 #include "endian.h" 00039 #include "mpa.h" 00040 00041 void __mpexp(mp_no *, mp_no *, int); 00042 00043 void __mplog(mp_no *x, mp_no *y, int p) { 00044 #include "mplog.h" 00045 int i,m; 00046 #if 0 00047 int j,k,m1,m2,n; 00048 double a,b; 00049 #endif 00050 static const int mp[33] = {0,0,0,0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3, 00051 4,4,4,4,4,4,4,4,4,4,4,4,4,4}; 00052 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, 00053 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, 00054 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; 00055 mp_no mpt1,mpt2; 00056 00057 /* Choose m and initiate mpone */ 00058 m = mp[p]; mpone.e = 1; mpone.d[0]=mpone.d[1]=ONE; 00059 00060 /* Perform m newton iterations to solve for y: exp(y)x=0. */ 00061 /* The iterations formula is: y(n+1)=y(n)+(x*exp(y(n))1). */ 00062 __cpy(y,&mpt1,p); 00063 for (i=0; i<m; i++) { 00064 mpt1.d[0]=mpt1.d[0]; 00065 __mpexp(&mpt1,&mpt2,p); 00066 __mul(x,&mpt2,&mpt1,p); 00067 __sub(&mpt1,&mpone,&mpt2,p); 00068 __add(y,&mpt2,&mpt1,p); 00069 __cpy(&mpt1,y,p); 00070 } 00071 return; 00072 }