summaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/dbl-64/e_exp.c
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/dbl-64/e_exp.c')
-rw-r--r--sysdeps/ieee754/dbl-64/e_exp.c150
1 files changed, 71 insertions, 79 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c
index ad1bc84625..ddd2bcb1c2 100644
--- a/sysdeps/ieee754/dbl-64/e_exp.c
+++ b/sysdeps/ieee754/dbl-64/e_exp.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001-2016 Free Software Foundation, Inc.
+ * Copyright (C) 2001-2018 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
@@ -23,10 +23,9 @@
/* 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 */
+/* it computes an almost correctly rounded (to nearest) value of e^x */
/* Assumption: Machine arithmetic operations are performed in */
/* round to nearest mode of IEEE 754 standard. */
/* */
@@ -38,25 +37,25 @@
#include "mydefs.h"
#include "MathLib.h"
#include "uexp.tbl"
+#include <math-barriers.h>
#include <math_private.h>
#include <fenv.h>
#include <float.h>
+#include "eexp.tbl"
#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;
+ double z;
mynumber junk1, junk2, binexp = {{0, 0}};
int4 i, j, m, n, ex;
+ int4 k;
double retval;
{
@@ -66,7 +65,42 @@ __ieee754_exp (double x)
m = junk1.i[HIGH_HALF];
n = m & hugeint;
- if (n > smallint && n < bigint)
+ if (n < 0x3ff0a2b2) /* |x| < 1.03972053527832 */
+ {
+ if (n < 0x3f862e42) /* |x| < 3/2 ln 2 */
+ {
+ if (n < 0x3ed00000) /* |x| < 1/64 ln 2 */
+ {
+ if (n < 0x3e300000) /* |x| < 2^18 */
+ {
+ retval = one + junk1.x;
+ goto ret;
+ }
+ retval = one + junk1.x * (one + half * junk1.x);
+ goto ret;
+ }
+ t = junk1.x * junk1.x;
+ retval = junk1.x + (t * (half + junk1.x * t2) +
+ (t * t) * (t3 + junk1.x * t4 + t * t5));
+ retval = one + retval;
+ goto ret;
+ }
+
+ /* Find the multiple of 2^-6 nearest x. */
+ k = n >> 20;
+ j = (0x00100000 | (n & 0x000fffff)) >> (0x40c - k);
+ j = (j - 1) & ~1;
+ if (m < 0)
+ j += 134;
+ z = junk1.x - TBL2[j];
+ t = z * z;
+ retval = z + (t * (half + (z * t2))
+ + (t * t) * (t3 + z * t4 + t * t5));
+ retval = TBL2[j + 1] + TBL2[j + 1] * retval;
+ goto ret;
+ }
+
+ if (n < bigint) /* && |x| >= 1.03972053527832 */
{
y = x * log2e.x + three51.x;
bexp = y - three51.x; /* multiply the result by 2**bexp */
@@ -93,22 +127,9 @@ __ieee754_exp (double x)
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;
+ /* Maximum relative error is 7.8e-22 (70.1 bits).
+ Maximum ULP error is 0.500007. */
+ retval = res * binexp.x;
goto ret;
}
@@ -166,38 +187,22 @@ __ieee754_exp (double x)
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 check_uflow_ret;
- } /*if error is over bound */
+ /* Does not underflow: res >= 1.0, binexp >= 0x1p-1022
+ Maximum relative error is 7.8e-22 (70.1 bits).
+ Maximum ULP error is 0.500007. */
+ retval = res * binexp.x;
+ goto ret;
}
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 check_uflow_ret;
- }
- else
- {
- retval = __slowexp (x);
- goto check_uflow_ret;
- } /* if error is over bound */
- check_uflow_ret:
+ /* Maximum ULP error is 0.5000035. */
+ binexp.i[HIGH_HALF] = 0x00100000;
+ retval = (res - 1.0) * binexp.x;
if (retval < DBL_MIN)
{
double force_underflow = tiny * tiny;
@@ -210,10 +215,9 @@ __ieee754_exp (double x)
else
{
binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
- if (res == (res + cor * err_0))
- retval = res * binexp.x * t256.x;
- else
- retval = __slowexp (x);
+ /* Maximum relative error is 7.8e-22 (70.1 bits).
+ Maximum ULP error is 0.500007. */
+ retval = res * binexp.x * t256.x;
if (isinf (retval))
goto ret_huge;
else
@@ -233,13 +237,10 @@ ret:
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. */
+/* Compute e^(x+xx). */
double
SECTION
-__exp1 (double x, double xx, double error)
+__exp1 (double x, double xx)
{
double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
mynumber junk1, junk2, binexp = {{0, 0}};
@@ -249,6 +250,7 @@ __exp1 (double x, double xx, double error)
m = junk1.i[HIGH_HALF];
n = m & hugeint; /* no sign */
+ /* fabs (x) > 5.551112e-17 and fabs (x) < 7.080010e+02. */
if (n > smallint && n < bigint)
{
y = x * log2e.x + three51.x;
@@ -276,11 +278,9 @@ __exp1 (double x, double xx, double error)
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;
+ /* Maximum relative error before rounding is 8.8e-22 (69.9 bits).
+ Maximum ULP error is 0.500008. */
+ return res * binexp.x;
}
if (n <= smallint)
@@ -318,6 +318,7 @@ __exp1 (double x, double xx, double error)
cor = (al - res) + rem;
if (m >> 31)
{
+ /* x < 0. */
ex = junk1.i[LOW_HALF];
if (res < 1.0)
{
@@ -328,34 +329,25 @@ __exp1 (double x, double xx, double error)
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;
+ /* Maximum ULP error is 0.500008. */
+ return res * binexp.x;
}
+ /* Denormal case - ex < -1022. */
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;
+ binexp.i[HIGH_HALF] = 0x00100000;
+ /* Maximum ULP error is 0.500004. */
+ return (res - 1.0) * binexp.x;
}
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;
+ /* Maximum ULP error is 0.500008. */
+ return res * binexp.x * t256.x;
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-04-30 17:20:11 +0000