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