|
Packit |
6c4009 |
/*
|
|
Packit |
6c4009 |
* IBM Accurate Mathematical Library
|
|
Packit |
6c4009 |
* written by International Business Machines Corp.
|
|
Packit |
6c4009 |
* Copyright (C) 2001-2018 Free Software Foundation, Inc.
|
|
Packit |
6c4009 |
*
|
|
Packit |
6c4009 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
6c4009 |
* it under the terms of the GNU Lesser General Public License as published by
|
|
Packit |
6c4009 |
* the Free Software Foundation; either version 2.1 of the License, or
|
|
Packit |
6c4009 |
* (at your option) any later version.
|
|
Packit |
6c4009 |
*
|
|
Packit |
6c4009 |
* This program is distributed in the hope that it will be useful,
|
|
Packit |
6c4009 |
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
6c4009 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
Packit |
6c4009 |
* GNU Lesser General Public License for more details.
|
|
Packit |
6c4009 |
*
|
|
Packit |
6c4009 |
* You should have received a copy of the GNU Lesser General Public License
|
|
Packit |
6c4009 |
* along with this program; if not, see <http://www.gnu.org/licenses/>.
|
|
Packit |
6c4009 |
*/
|
|
Packit |
6c4009 |
/*********************************************************************/
|
|
Packit |
6c4009 |
/* */
|
|
Packit |
6c4009 |
/* MODULE_NAME:ulog.c */
|
|
Packit |
6c4009 |
/* */
|
|
Packit |
6c4009 |
/* FUNCTION:ulog */
|
|
Packit |
6c4009 |
/* */
|
|
Packit |
6c4009 |
/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h */
|
|
Packit |
6c4009 |
/* ulog.tbl */
|
|
Packit |
6c4009 |
/* */
|
|
Packit |
6c4009 |
/* An ultimate log routine. Given an IEEE double machine number x */
|
|
Packit |
6c4009 |
/* it computes the rounded (to nearest) value of log(x). */
|
|
Packit |
6c4009 |
/* Assumption: Machine arithmetic operations are performed in */
|
|
Packit |
6c4009 |
/* round to nearest mode of IEEE 754 standard. */
|
|
Packit |
6c4009 |
/* */
|
|
Packit |
6c4009 |
/*********************************************************************/
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
#include "endian.h"
|
|
Packit |
6c4009 |
#include <dla.h>
|
|
Packit |
6c4009 |
#include "mpa.h"
|
|
Packit |
6c4009 |
#include "MathLib.h"
|
|
Packit |
6c4009 |
#include <math.h>
|
|
Packit |
6c4009 |
#include <math_private.h>
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
#ifndef SECTION
|
|
Packit |
6c4009 |
# define SECTION
|
|
Packit |
6c4009 |
#endif
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/*********************************************************************/
|
|
Packit |
6c4009 |
/* An ultimate log routine. Given an IEEE double machine number x */
|
|
Packit |
6c4009 |
/* it computes the rounded (to nearest) value of log(x). */
|
|
Packit |
6c4009 |
/*********************************************************************/
|
|
Packit |
6c4009 |
double
|
|
Packit |
6c4009 |
SECTION
|
|
Packit |
6c4009 |
__ieee754_log (double x)
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
int i, j, n, ux, dx;
|
|
Packit |
6c4009 |
double dbl_n, u, p0, q, r0, w, nln2a, luai, lubi, lvaj, lvbj,
|
|
Packit |
6c4009 |
sij, ssij, ttij, A, B, B0, polI, polII, t8, a, aa, b, bb, c;
|
|
Packit |
6c4009 |
#ifndef DLA_FMS
|
|
Packit |
6c4009 |
double t1, t2, t3, t4, t5;
|
|
Packit |
6c4009 |
#endif
|
|
Packit |
6c4009 |
number num;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
#include "ulog.tbl"
|
|
Packit |
6c4009 |
#include "ulog.h"
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
num.d = x;
|
|
Packit |
6c4009 |
ux = num.i[HIGH_HALF];
|
|
Packit |
6c4009 |
dx = num.i[LOW_HALF];
|
|
Packit |
6c4009 |
n = 0;
|
|
Packit |
6c4009 |
if (__glibc_unlikely (ux < 0x00100000))
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
if (__glibc_unlikely (((ux & 0x7fffffff) | dx) == 0))
|
|
Packit |
6c4009 |
return MHALF / 0.0; /* return -INF */
|
|
Packit |
6c4009 |
if (__glibc_unlikely (ux < 0))
|
|
Packit |
6c4009 |
return (x - x) / 0.0; /* return NaN */
|
|
Packit |
6c4009 |
n -= 54;
|
|
Packit |
6c4009 |
x *= two54.d; /* scale x */
|
|
Packit |
6c4009 |
num.d = x;
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
if (__glibc_unlikely (ux >= 0x7ff00000))
|
|
Packit |
6c4009 |
return x + x; /* INF or NaN */
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Regular values of x */
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
w = x - 1;
|
|
Packit |
6c4009 |
if (__glibc_likely (fabs (w) > U03))
|
|
Packit |
6c4009 |
goto case_03;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* log (1) is +0 in all rounding modes. */
|
|
Packit |
6c4009 |
if (w == 0.0)
|
|
Packit |
6c4009 |
return 0.0;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/*--- The case abs(x-1) < 0.03 */
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
t8 = MHALF * w;
|
|
Packit |
6c4009 |
EMULV (t8, w, a, aa, t1, t2, t3, t4, t5);
|
|
Packit |
6c4009 |
EADD (w, a, b, bb);
|
|
Packit |
6c4009 |
/* Evaluate polynomial II */
|
|
Packit |
6c4009 |
polII = b7.d + w * b8.d;
|
|
Packit |
6c4009 |
polII = b6.d + w * polII;
|
|
Packit |
6c4009 |
polII = b5.d + w * polII;
|
|
Packit |
6c4009 |
polII = b4.d + w * polII;
|
|
Packit |
6c4009 |
polII = b3.d + w * polII;
|
|
Packit |
6c4009 |
polII = b2.d + w * polII;
|
|
Packit |
6c4009 |
polII = b1.d + w * polII;
|
|
Packit |
6c4009 |
polII = b0.d + w * polII;
|
|
Packit |
6c4009 |
polII *= w * w * w;
|
|
Packit |
6c4009 |
c = (aa + bb) + polII;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Here b contains the high part of the result, and c the low part.
|
|
Packit |
6c4009 |
Maximum error is b * 2.334e-19, so accuracy is >61 bits.
|
|
Packit |
6c4009 |
Therefore max ULP error of b + c is ~0.502. */
|
|
Packit |
6c4009 |
return b + c;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/*--- The case abs(x-1) > 0.03 */
|
|
Packit |
6c4009 |
case_03:
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Find n,u such that x = u*2**n, 1/sqrt(2) < u < sqrt(2) */
|
|
Packit |
6c4009 |
n += (num.i[HIGH_HALF] >> 20) - 1023;
|
|
Packit |
6c4009 |
num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000;
|
|
Packit |
6c4009 |
if (num.d > SQRT_2)
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
num.d *= HALF;
|
|
Packit |
6c4009 |
n++;
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
u = num.d;
|
|
Packit |
6c4009 |
dbl_n = (double) n;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */
|
|
Packit |
6c4009 |
num.d += h1.d;
|
|
Packit |
6c4009 |
i = (num.i[HIGH_HALF] & 0x000fffff) >> 12;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */
|
|
Packit |
6c4009 |
num.d = u * Iu[i].d + h2.d;
|
|
Packit |
6c4009 |
j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Compute w=(u-ui*vj)/(ui*vj) */
|
|
Packit |
6c4009 |
p0 = (1 + (i - 75) * DEL_U) * (1 + (j - 180) * DEL_V);
|
|
Packit |
6c4009 |
q = u - p0;
|
|
Packit |
6c4009 |
r0 = Iu[i].d * Iv[j].d;
|
|
Packit |
6c4009 |
w = q * r0;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Evaluate polynomial I */
|
|
Packit |
6c4009 |
polI = w + (a2.d + a3.d * w) * w * w;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Add up everything */
|
|
Packit |
6c4009 |
nln2a = dbl_n * LN2A;
|
|
Packit |
6c4009 |
luai = Lu[i][0].d;
|
|
Packit |
6c4009 |
lubi = Lu[i][1].d;
|
|
Packit |
6c4009 |
lvaj = Lv[j][0].d;
|
|
Packit |
6c4009 |
lvbj = Lv[j][1].d;
|
|
Packit |
6c4009 |
EADD (luai, lvaj, sij, ssij);
|
|
Packit |
6c4009 |
EADD (nln2a, sij, A, ttij);
|
|
Packit |
6c4009 |
B0 = (((lubi + lvbj) + ssij) + ttij) + dbl_n * LN2B;
|
|
Packit |
6c4009 |
B = polI + B0;
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
/* Here A contains the high part of the result, and B the low part.
|
|
Packit |
6c4009 |
Maximum abs error is 6.095e-21 and min log (x) is 0.0295 since x > 1.03.
|
|
Packit |
6c4009 |
Therefore max ULP error of A + B is ~0.502. */
|
|
Packit |
6c4009 |
return A + B;
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
#ifndef __ieee754_log
|
|
Packit |
6c4009 |
strong_alias (__ieee754_log, __log_finite)
|
|
Packit |
6c4009 |
#endif
|