Back to index

glibc  2.9
mplog.c
Go to the documentation of this file.
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 02111-1307, 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 /* Multi-Precision logarithm function subroutine (for precision p >= 4, */
00031 /* 2**(-1024) < x < 2**1024) and x is outside of the interval           */
00032 /* [1-2**(-54),1+2**(-54)]. Upon entry, x should be set to the          */
00033 /* multi-precision 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 }