mirror/glibc
2
0
Fork 0
mirror of git://sourceware.org/git/glibc.git synced 2024-11-27 03:41:23 +08:00
Code Issues Packages Projects Releases Wiki Activity
1d70fd6c75
glibc/sysdeps/ieee754/dbl-64/e_exp.c
Siddhesh Poyarekar e7b2d1dd62 Format e_exp.c
2013-10-08 16:22:28 +05:30

343 lines
9.3 KiB
C
Raw Blame History

/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
* Copyright (C) 2001-2013 Free Software Foundation, Inc.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; either version 2.1 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program; if not, see <http://www.gnu.org/licenses/>.
*/
/***************************************************************************/
/* MODULE_NAME:uexp.c */
/* */
/* FUNCTION:uexp */
/* exp1 */
/* */
/* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h */
/* mpa.c mpexp.x slowexp.c */
/* */
/* An ultimate exp routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of e^x */
/* Assumption: Machine arithmetic operations are performed in */
/* round to nearest mode of IEEE 754 standard. */
/* */
/***************************************************************************/
#include "endian.h"
#include "uexp.h"
#include "mydefs.h"
#include "MathLib.h"
#include "uexp.tbl"
#include <math_private.h>
#include <fenv.h>
#ifndef SECTION
# define SECTION
#endif
double __slowexp (double);
/* An ultimate exp routine. Given an IEEE double machine number x it computes
the correctly rounded (to nearest) value of e^x. */
double
SECTION
__ieee754_exp (double x)
{
double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
mynumber junk1, junk2, binexp = {{0, 0}};
int4 i, j, m, n, ex;
double retval;
SET_RESTORE_ROUND (FE_TONEAREST);
junk1.x = x;
m = junk1.i[HIGH_HALF];
n = m & hugeint;
if (n > smallint && n < bigint)
{
y = x * log2e.x + three51.x;
bexp = y - three51.x; /* multiply the result by 2**bexp */
junk1.x = y;
eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */
t = x - bexp * ln_two1.x;
y = t + three33.x;
base = y - three33.x; /* t rounded to a multiple of 2**-18 */
junk2.x = y;
del = (t - base) - eps; /* x = bexp*ln(2) + base + del */
eps = del + del * del * (p3.x * del + p2.x);
binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
j = (junk2.i[LOW_HALF] & 511) << 1;
al = coar.x[i] * fine.x[j];
bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
+ coar.x[i + 1] * fine.x[j + 1]);
rem = (bet + bet * eps) + al * eps;
res = al + rem;
cor = (al - res) + rem;
if (res == (res + cor * err_0))
{
retval = res * binexp.x;
goto ret;
}
else
{
retval = __slowexp (x);
goto ret;
} /*if error is over bound */
}
if (n <= smallint)
{
retval = 1.0;
goto ret;
}
if (n >= badint)
{
if (n > infint)
{
retval = x + x;
goto ret;
} /* x is NaN */
if (n < infint)
{
retval = (x > 0) ? (hhuge * hhuge) : (tiny * tiny);
goto ret;
}
/* x is finite, cause either overflow or underflow */
if (junk1.i[LOW_HALF] != 0)
{
retval = x + x;
goto ret;
} /* x is NaN */
retval = (x > 0) ? inf.x : zero; /* |x| = inf; return either inf or 0 */
goto ret;
}
y = x * log2e.x + three51.x;
bexp = y - three51.x;
junk1.x = y;
eps = bexp * ln_two2.x;
t = x - bexp * ln_two1.x;
y = t + three33.x;
base = y - three33.x;
junk2.x = y;
del = (t - base) - eps;
eps = del + del * del * (p3.x * del + p2.x);
i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
j = (junk2.i[LOW_HALF] & 511) << 1;
al = coar.x[i] * fine.x[j];
bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
+ coar.x[i + 1] * fine.x[j + 1]);
rem = (bet + bet * eps) + al * eps;
res = al + rem;
cor = (al - res) + rem;
if (m >> 31)
{
ex = junk1.i[LOW_HALF];
if (res < 1.0)
{
res += res;
cor += cor;
ex -= 1;
}
if (ex >= -1022)
{
binexp.i[HIGH_HALF] = (1023 + ex) << 20;
if (res == (res + cor * err_0))
{
retval = res * binexp.x;
goto ret;
}
else
{
retval = __slowexp (x);
goto ret;
} /*if error is over bound */
}
ex = -(1022 + ex);
binexp.i[HIGH_HALF] = (1023 - ex) << 20;
res *= binexp.x;
cor *= binexp.x;
eps = 1.0000000001 + err_0 * binexp.x;
t = 1.0 + res;
y = ((1.0 - t) + res) + cor;
res = t + y;
cor = (t - res) + y;
if (res == (res + eps * cor))
{
binexp.i[HIGH_HALF] = 0x00100000;
retval = (res - 1.0) * binexp.x;
goto ret;
}
else
{
retval = __slowexp (x);
goto ret;
} /* if error is over bound */
}
else
{
binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
if (res == (res + cor * err_0))
{
retval = res * binexp.x * t256.x;
goto ret;
}
else
{
retval = __slowexp (x);
goto ret;
}
}
ret:
return retval;
}
#ifndef __ieee754_exp
strong_alias (__ieee754_exp, __exp_finite)
#endif
/* Compute e^(x+xx). The routine also receives bound of error of previous
calculation. If after computing exp the error exceeds the allowed bounds,
the routine returns a non-positive number. Otherwise it returns the
computed result, which is always positive. */
double
SECTION
__exp1 (double x, double xx, double error)
{
double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
mynumber junk1, junk2, binexp = {{0, 0}};
int4 i, j, m, n, ex;
junk1.x = x;
m = junk1.i[HIGH_HALF];
n = m & hugeint; /* no sign */
if (n > smallint && n < bigint)
{
y = x * log2e.x + three51.x;
bexp = y - three51.x; /* multiply the result by 2**bexp */
junk1.x = y;
eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */
t = x - bexp * ln_two1.x;
y = t + three33.x;
base = y - three33.x; /* t rounded to a multiple of 2**-18 */
junk2.x = y;
del = (t - base) + (xx - eps); /* x = bexp*ln(2) + base + del */
eps = del + del * del * (p3.x * del + p2.x);
binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
j = (junk2.i[LOW_HALF] & 511) << 1;
al = coar.x[i] * fine.x[j];
bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
+ coar.x[i + 1] * fine.x[j + 1]);
rem = (bet + bet * eps) + al * eps;
res = al + rem;
cor = (al - res) + rem;
if (res == (res + cor * (1.0 + error + err_1)))
return res * binexp.x;
else
return -10.0;
}
if (n <= smallint)
return 1.0; /* if x->0 e^x=1 */
if (n >= badint)
{
if (n > infint)
return (zero / zero); /* x is NaN, return invalid */
if (n < infint)
return ((x > 0) ? (hhuge * hhuge) : (tiny * tiny));
/* x is finite, cause either overflow or underflow */
if (junk1.i[LOW_HALF] != 0)
return (zero / zero); /* x is NaN */
return ((x > 0) ? inf.x : zero); /* |x| = inf; return either inf or 0 */
}
y = x * log2e.x + three51.x;
bexp = y - three51.x;
junk1.x = y;
eps = bexp * ln_two2.x;
t = x - bexp * ln_two1.x;
y = t + three33.x;
base = y - three33.x;
junk2.x = y;
del = (t - base) + (xx - eps);
eps = del + del * del * (p3.x * del + p2.x);
i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
j = (junk2.i[LOW_HALF] & 511) << 1;
al = coar.x[i] * fine.x[j];
bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
+ coar.x[i + 1] * fine.x[j + 1]);
rem = (bet + bet * eps) + al * eps;
res = al + rem;
cor = (al - res) + rem;
if (m >> 31)
{
ex = junk1.i[LOW_HALF];
if (res < 1.0)
{
res += res;
cor += cor;
ex -= 1;
}
if (ex >= -1022)
{
binexp.i[HIGH_HALF] = (1023 + ex) << 20;
if (res == (res + cor * (1.0 + error + err_1)))
return res * binexp.x;
else
return -10.0;
}
ex = -(1022 + ex);
binexp.i[HIGH_HALF] = (1023 - ex) << 20;
res *= binexp.x;
cor *= binexp.x;
eps = 1.00000000001 + (error + err_1) * binexp.x;
t = 1.0 + res;
y = ((1.0 - t) + res) + cor;
res = t + y;
cor = (t - res) + y;
if (res == (res + eps * cor))
{
binexp.i[HIGH_HALF] = 0x00100000;
return (res - 1.0) * binexp.x;
}
else
return -10.0;
}
else
{
binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
if (res == (res + cor * (1.0 + error + err_1)))
return res * binexp.x * t256.x;
else
return -10.0;
}
}
Reference in New Issue View Git Blame Copy Permalink