1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001 Free Software Foundation
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19 */
20/*********************************************************************/
21/*                                                                   */
22/*      MODULE_NAME:ulog.c                                           */
23/*                                                                   */
24/*      FUNCTION:ulog                                                */
25/*                                                                   */
26/*      FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h           */
27/*                    mpexp.c mplog.c mpa.c                          */
28/*                    ulog.tbl                                       */
29/*                                                                   */
30/* An ultimate log routine. Given an IEEE double machine number x    */
31/* it computes the correctly rounded (to nearest) value of log(x).   */
32/* Assumption: Machine arithmetic operations are performed in        */
33/* round to nearest mode of IEEE 754 standard.                       */
34/*                                                                   */
35/*********************************************************************/
36
37
38#include "endian.h"
39#include "dla.h"
40#include "mpa.h"
41#include "MathLib.h"
42#include "math_private.h"
43
44void __mplog(mp_no *, mp_no *, int);
45
46/*********************************************************************/
47/* An ultimate log routine. Given an IEEE double machine number x     */
48/* it computes the correctly rounded (to nearest) value of log(x).   */
49/*********************************************************************/
50double __ieee754_log(double x) {
51#define M 4
52  static const int pr[M]={8,10,18,32};
53  int i,j,n,ux,dx,p;
54#if 0
55  int k;
56#endif
57  double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj,
58         sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb,
59         t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww,
60         a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c;
61  number num;
62  mp_no mpx,mpy,mpy1,mpy2,mperr;
63
64#include "ulog.tbl"
65#include "ulog.h"
66
67  /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */
68
69  num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
70  n=0;
71  if (ux < 0x00100000) {
72    if (((ux & 0x7fffffff) | dx) == 0)  return MHALF/ZERO; /* return -INF */
73    if (ux < 0) return (x-x)/ZERO;                         /* return NaN  */
74    n -= 54;    x *= two54.d;                              /* scale x     */
75    num.d = x;
76  }
77  if (ux >= 0x7ff00000) return x+x;                        /* INF or NaN  */
78
79  /* Regular values of x */
80
81  w = x-ONE;
82  if (ABS(w) > U03) { goto case_03; }
83
84
85  /*--- Stage I, the case abs(x-1) < 0.03 */
86
87  t8 = MHALF*w;
88  EMULV(t8,w,a,aa,t1,t2,t3,t4,t5)
89  EADD(w,a,b,bb)
90
91  /* Evaluate polynomial II */
92  polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+
93          w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w;
94  c = (aa+bb)+polII;
95
96  /* End stage I, case abs(x-1) < 0.03 */
97  if ((y=b+(c+b*E2)) == b+(c-b*E2))  return y;
98
99  /*--- Stage II, the case abs(x-1) < 0.03 */
100
101  a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+
102            w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d))))))));
103  EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5)
104  ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2)
105  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
106  ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2)
107  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
108  ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2)
109  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
110  ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2)
111  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
112  ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2)
113  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
114  ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2)
115  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
116  ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2)
117  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
118  ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2)
119  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
120  ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2)
121  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
122  MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
123  ADD2(w,ZERO,    s3,ss3, b, bb,t1,t2)
124
125  /* End stage II, case abs(x-1) < 0.03 */
126  if ((y=b+(bb+b*E4)) == b+(bb-b*E4))  return y;
127  goto stage_n;
128
129  /*--- Stage I, the case abs(x-1) > 0.03 */
130  case_03:
131
132  /* Find n,u such that x = u*2**n,   1/sqrt(2) < u < sqrt(2)  */
133  n += (num.i[HIGH_HALF] >> 20) - 1023;
134  num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000;
135  if (num.d > SQRT_2) { num.d *= HALF;  n++; }
136  u = num.d;  dbl_n = (double) n;
137
138  /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */
139  num.d += h1.d;
140  i = (num.i[HIGH_HALF] & 0x000fffff) >> 12;
141
142  /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */
143  num.d = u*Iu[i].d + h2.d;
144  j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;
145
146  /* Compute w=(u-ui*vj)/(ui*vj) */
147  p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V);
148  q=u-p0;   r0=Iu[i].d*Iv[j].d;   w=q*r0;
149
150  /* Evaluate polynomial I */
151  polI = w+(a2.d+a3.d*w)*w*w;
152
153  /* Add up everything */
154  nln2a = dbl_n*LN2A;
155  luai  = Lu[i][0].d;   lubi  = Lu[i][1].d;
156  lvaj  = Lv[j][0].d;   lvbj  = Lv[j][1].d;
157  EADD(luai,lvaj,sij,ssij)
158  EADD(nln2a,sij,A  ,ttij)
159  B0 = (((lubi+lvbj)+ssij)+ttij)+dbl_n*LN2B;
160  B  = polI+B0;
161
162  /* End stage I, case abs(x-1) >= 0.03 */
163  if ((y=A+(B+E1)) == A+(B-E1))  return y;
164
165
166  /*--- Stage II, the case abs(x-1) > 0.03 */
167
168  /* Improve the accuracy of r0 */
169  EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5)
170  t=r0*((ONE-sa)-sb);
171  EADD(r0,t,ra,rb)
172
173  /* Compute w */
174  MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8)
175
176  EADD(A,B0,a0,aa0)
177
178  /* Evaluate polynomial III */
179  s1 = (c3.d+(c4.d+c5.d*w)*w)*w;
180  EADD(c2.d,s1,s2,ss2)
181  MUL2(s2,ss2,w,ww,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
182  MUL2(s3,ss3,w,ww,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
183  ADD2(s2,ss2,w,ww,s3,ss3,t1,t2)
184  ADD2(s3,ss3,a0,aa0,a1,aa1,t1,t2)
185
186  /* End stage II, case abs(x-1) >= 0.03 */
187  if ((y=a1+(aa1+E3)) == a1+(aa1-E3)) return y;
188
189
190  /* Final stages. Use multi-precision arithmetic. */
191  stage_n:
192
193  for (i=0; i<M; i++) {
194    p = pr[i];
195    __dbl_mp(x,&mpx,p);  __dbl_mp(y,&mpy,p);
196    __mplog(&mpx,&mpy,p);
197    __dbl_mp(e[i].d,&mperr,p);
198    __add(&mpy,&mperr,&mpy1,p);  __sub(&mpy,&mperr,&mpy2,p);
199    __mp_dbl(&mpy1,&y1,p);       __mp_dbl(&mpy2,&y2,p);
200    if (y1==y2)   return y1;
201  }
202  return y1;
203}
204