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
41void __mpexp(mp_no *, mp_no *, int);
42
43void __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}
73