sysdeps/ieee754/dbl-64/mplog.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 /*                                                                      */
  23 /*     MODULE_NAME:mplog.c                                              */
  24 /*                                                                      */
  25 /*     FUNCTIONS:  mplog                                                */
  26 /*                                                                      */
  27 /*     FILES NEEDED: endian.h mpa.h  mplog.h                            */
  28 /*                    mpexp.c                                           */
  29 /*                                                                      */
  30 /* Multi-Precision logarithm function subroutine (for precision p >= 4, */
  31 /* 2**(-1024) < x < 2**1024) and x is outside of the interval           */
  32 /* [1-2**(-54),1+2**(-54)]. Upon entry, x should be set to the          */
  33 /* multi-precision value of the input and y should be set into a multi- */
  34 /* precision value of an approximation of log(x) with relative error    */
  35 /* bound of at most 2**(-52). The routine improves the accuracy of y.   */
  36 /*                                                                      */
  37 /************************************************************************/
  38 #include "endian.h"
  39 #include "mpa.h"
  40
  41 void __mpexp(mp_no *, mp_no *, int);
  42
  43 void __mplog(mp_no *x, mp_no *y, int p) {
  44 #include "mplog.h"
  45   int i,m;
  46 #if 0
  47   int j,k,m1,m2,n;
  48   double a,b;
  49 #endif
  50   static const int mp[33] = {0,0,0,0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,
  51                              4,4,4,4,4,4,4,4,4,4,4,4,4,4};
  52   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,
  53                     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,
  54                     0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
  55   mp_no mpt1,mpt2;
  56
  57   /* Choose m and initiate mpone */
  58   m = mp[p];  mpone.e = 1;  mpone.d[0]=mpone.d[1]=ONE;
  59
  60   /* Perform m newton iterations to solve for y: exp(y)-x=0.     */
  61   /* The iterations formula is:  y(n+1)=y(n)+(x*exp(-y(n))-1).   */
  62   __cpy(y,&mpt1,p);
  63   for (i=0; i<m; i++) {
  64     mpt1.d[0]=-mpt1.d[0];
  65     __mpexp(&mpt1,&mpt2,p);
  66     __mul(x,&mpt2,&mpt1,p);
  67     __sub(&mpt1,&mpone,&mpt2,p);
  68     __add(y,&mpt2,&mpt1,p);
  69     __cpy(&mpt1,y,p);
  70   }
  71   return;
  72 }