1996-03-05 22:41:30 +01:00
|
|
|
/*
|
2001-03-12 01:04:52 +01:00
|
|
|
* IBM Accurate Mathematical Library
|
2002-07-06 08:36:39 +02:00
|
|
|
* written by International Business Machines Corp.
|
2018-01-01 01:32:25 +01:00
|
|
|
* Copyright (C) 2001-2018 Free Software Foundation, Inc.
|
1996-03-05 22:41:30 +01:00
|
|
|
*
|
2001-03-12 01:04:52 +01:00
|
|
|
* 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
|
2002-08-27 00:40:48 +02:00
|
|
|
* the Free Software Foundation; either version 2.1 of the License, or
|
2001-03-12 01:04:52 +01:00
|
|
|
* (at your option) any later version.
|
1996-03-05 22:41:30 +01:00
|
|
|
*
|
2001-03-12 01:04:52 +01:00
|
|
|
* 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
|
2002-08-20 23:51:55 +02:00
|
|
|
* GNU Lesser General Public License for more details.
|
1996-03-05 22:41:30 +01:00
|
|
|
*
|
2001-03-12 01:04:52 +01:00
|
|
|
* You should have received a copy of the GNU Lesser General Public License
|
2012-02-10 00:18:22 +01:00
|
|
|
* along with this program; if not, see <http://www.gnu.org/licenses/>.
|
1996-03-05 22:41:30 +01:00
|
|
|
*/
|
2001-03-12 01:04:52 +01:00
|
|
|
/***************************************************************************/
|
|
|
|
/* MODULE_NAME: upow.c */
|
|
|
|
/* */
|
|
|
|
/* FUNCTIONS: upow */
|
|
|
|
/* log1 */
|
|
|
|
/* checkint */
|
|
|
|
/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h */
|
|
|
|
/* root.tbl uexp.tbl upow.tbl */
|
|
|
|
/* An ultimate power routine. Given two IEEE double machine numbers y,x */
|
|
|
|
/* it computes the correctly rounded (to nearest) value of x^y. */
|
|
|
|
/* Assumption: Machine arithmetic operations are performed in */
|
|
|
|
/* round to nearest mode of IEEE 754 standard. */
|
|
|
|
/* */
|
|
|
|
/***************************************************************************/
|
2014-06-23 22:12:33 +02:00
|
|
|
#include <math.h>
|
2001-03-12 01:04:52 +01:00
|
|
|
#include "endian.h"
|
|
|
|
#include "upow.h"
|
2011-10-23 18:50:28 +02:00
|
|
|
#include <dla.h>
|
2001-03-12 01:04:52 +01:00
|
|
|
#include "mydefs.h"
|
|
|
|
#include "MathLib.h"
|
|
|
|
#include "upow.tbl"
|
2012-03-09 20:29:16 +01:00
|
|
|
#include <math_private.h>
|
Do not include fenv_private.h in math_private.h.
Continuing the clean-up related to the catch-all math_private.h
header, this patch stops math_private.h from including fenv_private.h.
Instead, fenv_private.h is included directly from those users of
math_private.h that also used interfaces from fenv_private.h. No
attempt is made to remove unused includes of math_private.h, but that
is a natural followup.
(However, since math_private.h sometimes defines optimized versions of
math.h interfaces or __* variants thereof, as well as defining its own
interfaces, I think it might make sense to get all those optimized
versions included from include/math.h, not requiring a separate header
at all, before eliminating unused math_private.h includes - that
avoids a file quietly becoming less-optimized if someone adds a call
to one of those interfaces without restoring a math_private.h include
to that file.)
There is still a pitfall that if code uses plain fe* and __fe*
interfaces, but only includes fenv.h and not fenv_private.h or (before
this patch) math_private.h, it will compile on platforms with
exceptions and rounding modes but not get the optimized versions (and
possibly not compile) on platforms without exception and rounding mode
support, so making it easy to break the build for such platforms
accidentally.
I think it would be most natural to move the inlines / macros for fe*
and __fe* in the case of no exceptions and rounding modes into
include/fenv.h, so that all code including fenv.h with _ISOMAC not
defined automatically gets them. Then fenv_private.h would be purely
the header for the libc_fe*, SET_RESTORE_ROUND etc. internal
interfaces and the risk of breaking the build on other platforms than
the one you tested on because of a missing fenv_private.h include
would be much reduced (and there would be some unused fenv_private.h
includes to remove along with unused math_private.h includes).
Tested for x86_64 and x86, and tested with build-many-glibcs.py that
installed stripped shared libraries are unchanged by this patch.
* sysdeps/generic/math_private.h: Do not include <fenv_private.h>.
* math/fromfp.h: Include <fenv_private.h>.
* math/math-narrow.h: Likewise.
* math/s_cexp_template.c: Likewise.
* math/s_csin_template.c: Likewise.
* math/s_csinh_template.c: Likewise.
* math/s_ctan_template.c: Likewise.
* math/s_ctanh_template.c: Likewise.
* math/s_iseqsig_template.c: Likewise.
* math/w_acos_compat.c: Likewise.
* math/w_acosf_compat.c: Likewise.
* math/w_acosl_compat.c: Likewise.
* math/w_asin_compat.c: Likewise.
* math/w_asinf_compat.c: Likewise.
* math/w_asinl_compat.c: Likewise.
* math/w_ilogb_template.c: Likewise.
* math/w_j0_compat.c: Likewise.
* math/w_j0f_compat.c: Likewise.
* math/w_j0l_compat.c: Likewise.
* math/w_j1_compat.c: Likewise.
* math/w_j1f_compat.c: Likewise.
* math/w_j1l_compat.c: Likewise.
* math/w_jn_compat.c: Likewise.
* math/w_jnf_compat.c: Likewise.
* math/w_llogb_template.c: Likewise.
* math/w_log10_compat.c: Likewise.
* math/w_log10f_compat.c: Likewise.
* math/w_log10l_compat.c: Likewise.
* math/w_log2_compat.c: Likewise.
* math/w_log2f_compat.c: Likewise.
* math/w_log2l_compat.c: Likewise.
* math/w_log_compat.c: Likewise.
* math/w_logf_compat.c: Likewise.
* math/w_logl_compat.c: Likewise.
* sysdeps/aarch64/fpu/feholdexcpt.c: Likewise.
* sysdeps/aarch64/fpu/fesetround.c: Likewise.
* sysdeps/aarch64/fpu/fgetexcptflg.c: Likewise.
* sysdeps/aarch64/fpu/ftestexcept.c: Likewise.
* sysdeps/ieee754/dbl-64/e_atan2.c: Likewise.
* sysdeps/ieee754/dbl-64/e_exp.c: Likewise.
* sysdeps/ieee754/dbl-64/e_exp2.c: Likewise.
* sysdeps/ieee754/dbl-64/e_gamma_r.c: Likewise.
* sysdeps/ieee754/dbl-64/e_jn.c: Likewise.
* sysdeps/ieee754/dbl-64/e_pow.c: Likewise.
* sysdeps/ieee754/dbl-64/e_remainder.c: Likewise.
* sysdeps/ieee754/dbl-64/e_sqrt.c: Likewise.
* sysdeps/ieee754/dbl-64/gamma_product.c: Likewise.
* sysdeps/ieee754/dbl-64/lgamma_neg.c: Likewise.
* sysdeps/ieee754/dbl-64/s_atan.c: Likewise.
* sysdeps/ieee754/dbl-64/s_fma.c: Likewise.
* sysdeps/ieee754/dbl-64/s_fmaf.c: Likewise.
* sysdeps/ieee754/dbl-64/s_llrint.c: Likewise.
* sysdeps/ieee754/dbl-64/s_llround.c: Likewise.
* sysdeps/ieee754/dbl-64/s_lrint.c: Likewise.
* sysdeps/ieee754/dbl-64/s_lround.c: Likewise.
* sysdeps/ieee754/dbl-64/s_nearbyint.c: Likewise.
* sysdeps/ieee754/dbl-64/s_sin.c: Likewise.
* sysdeps/ieee754/dbl-64/s_sincos.c: Likewise.
* sysdeps/ieee754/dbl-64/s_tan.c: Likewise.
* sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c: Likewise.
* sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c: Likewise.
* sysdeps/ieee754/dbl-64/x2y2m1.c: Likewise.
* sysdeps/ieee754/float128/float128_private.h: Likewise.
* sysdeps/ieee754/flt-32/e_gammaf_r.c: Likewise.
* sysdeps/ieee754/flt-32/e_j1f.c: Likewise.
* sysdeps/ieee754/flt-32/e_jnf.c: Likewise.
* sysdeps/ieee754/flt-32/lgamma_negf.c: Likewise.
* sysdeps/ieee754/flt-32/s_llrintf.c: Likewise.
* sysdeps/ieee754/flt-32/s_llroundf.c: Likewise.
* sysdeps/ieee754/flt-32/s_lrintf.c: Likewise.
* sysdeps/ieee754/flt-32/s_lroundf.c: Likewise.
* sysdeps/ieee754/flt-32/s_nearbyintf.c: Likewise.
* sysdeps/ieee754/k_standardl.c: Likewise.
* sysdeps/ieee754/ldbl-128/e_expl.c: Likewise.
* sysdeps/ieee754/ldbl-128/e_gammal_r.c: Likewise.
* sysdeps/ieee754/ldbl-128/e_j1l.c: Likewise.
* sysdeps/ieee754/ldbl-128/e_jnl.c: Likewise.
* sysdeps/ieee754/ldbl-128/gamma_productl.c: Likewise.
* sysdeps/ieee754/ldbl-128/lgamma_negl.c: Likewise.
* sysdeps/ieee754/ldbl-128/s_fmal.c: Likewise.
* sysdeps/ieee754/ldbl-128/s_llrintl.c: Likewise.
* sysdeps/ieee754/ldbl-128/s_llroundl.c: Likewise.
* sysdeps/ieee754/ldbl-128/s_lrintl.c: Likewise.
* sysdeps/ieee754/ldbl-128/s_lroundl.c: Likewise.
* sysdeps/ieee754/ldbl-128/s_nearbyintl.c: Likewise.
* sysdeps/ieee754/ldbl-128/x2y2m1l.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_expl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_j1l.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_jnl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_fmal.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_llrintl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_llroundl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_lrintl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_lroundl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_rintl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c: Likewise.
* sysdeps/ieee754/ldbl-96/e_gammal_r.c: Likewise.
* sysdeps/ieee754/ldbl-96/e_jnl.c: Likewise.
* sysdeps/ieee754/ldbl-96/gamma_productl.c: Likewise.
* sysdeps/ieee754/ldbl-96/lgamma_negl.c: Likewise.
* sysdeps/ieee754/ldbl-96/s_fma.c: Likewise.
* sysdeps/ieee754/ldbl-96/s_fmal.c: Likewise.
* sysdeps/ieee754/ldbl-96/s_llrintl.c: Likewise.
* sysdeps/ieee754/ldbl-96/s_llroundl.c: Likewise.
* sysdeps/ieee754/ldbl-96/s_lrintl.c: Likewise.
* sysdeps/ieee754/ldbl-96/s_lroundl.c: Likewise.
* sysdeps/ieee754/ldbl-96/x2y2m1l.c: Likewise.
* sysdeps/powerpc/fpu/e_sqrt.c: Likewise.
* sysdeps/powerpc/fpu/e_sqrtf.c: Likewise.
* sysdeps/riscv/rv64/rvd/s_ceil.c: Likewise.
* sysdeps/riscv/rv64/rvd/s_floor.c: Likewise.
* sysdeps/riscv/rv64/rvd/s_nearbyint.c: Likewise.
* sysdeps/riscv/rv64/rvd/s_round.c: Likewise.
* sysdeps/riscv/rv64/rvd/s_roundeven.c: Likewise.
* sysdeps/riscv/rv64/rvd/s_trunc.c: Likewise.
* sysdeps/riscv/rvd/s_finite.c: Likewise.
* sysdeps/riscv/rvd/s_fmax.c: Likewise.
* sysdeps/riscv/rvd/s_fmin.c: Likewise.
* sysdeps/riscv/rvd/s_fpclassify.c: Likewise.
* sysdeps/riscv/rvd/s_isinf.c: Likewise.
* sysdeps/riscv/rvd/s_isnan.c: Likewise.
* sysdeps/riscv/rvd/s_issignaling.c: Likewise.
* sysdeps/riscv/rvf/fegetround.c: Likewise.
* sysdeps/riscv/rvf/feholdexcpt.c: Likewise.
* sysdeps/riscv/rvf/fesetenv.c: Likewise.
* sysdeps/riscv/rvf/fesetround.c: Likewise.
* sysdeps/riscv/rvf/feupdateenv.c: Likewise.
* sysdeps/riscv/rvf/fgetexcptflg.c: Likewise.
* sysdeps/riscv/rvf/ftestexcept.c: Likewise.
* sysdeps/riscv/rvf/s_ceilf.c: Likewise.
* sysdeps/riscv/rvf/s_finitef.c: Likewise.
* sysdeps/riscv/rvf/s_floorf.c: Likewise.
* sysdeps/riscv/rvf/s_fmaxf.c: Likewise.
* sysdeps/riscv/rvf/s_fminf.c: Likewise.
* sysdeps/riscv/rvf/s_fpclassifyf.c: Likewise.
* sysdeps/riscv/rvf/s_isinff.c: Likewise.
* sysdeps/riscv/rvf/s_isnanf.c: Likewise.
* sysdeps/riscv/rvf/s_issignalingf.c: Likewise.
* sysdeps/riscv/rvf/s_nearbyintf.c: Likewise.
* sysdeps/riscv/rvf/s_roundevenf.c: Likewise.
* sysdeps/riscv/rvf/s_roundf.c: Likewise.
* sysdeps/riscv/rvf/s_truncf.c: Likewise.
2018-09-03 23:09:04 +02:00
|
|
|
#include <fenv_private.h>
|
2018-05-10 02:53:04 +02:00
|
|
|
#include <math-underflow.h>
|
2012-03-05 13:22:46 +01:00
|
|
|
#include <fenv.h>
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2011-10-25 06:56:33 +02:00
|
|
|
#ifndef SECTION
|
|
|
|
# define SECTION
|
|
|
|
#endif
|
|
|
|
|
2012-04-09 11:43:18 +02:00
|
|
|
static const double huge = 1.0e300, tiny = 1.0e-300;
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2018-02-12 11:42:42 +01:00
|
|
|
double __exp1 (double x, double xx);
|
|
|
|
static double log1 (double x, double *delta);
|
2013-10-08 12:53:16 +02:00
|
|
|
static int checkint (double x);
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
/* An ultimate power routine. Given two IEEE double machine numbers y, x it
|
|
|
|
computes the correctly rounded (to nearest) value of X^y. */
|
2011-10-25 06:56:33 +02:00
|
|
|
double
|
|
|
|
SECTION
|
2013-10-08 12:53:16 +02:00
|
|
|
__ieee754_pow (double x, double y)
|
|
|
|
{
|
2018-02-12 11:42:42 +01:00
|
|
|
double z, a, aa, t, a1, a2, y1, y2;
|
2013-10-08 12:53:16 +02:00
|
|
|
mynumber u, v;
|
2001-03-12 01:04:52 +01:00
|
|
|
int k;
|
2013-10-08 12:53:16 +02:00
|
|
|
int4 qx, qy;
|
|
|
|
v.x = y;
|
|
|
|
u.x = x;
|
|
|
|
if (v.i[LOW_HALF] == 0)
|
|
|
|
{ /* of y */
|
|
|
|
qx = u.i[HIGH_HALF] & 0x7fffffff;
|
|
|
|
/* Is x a NaN? */
|
Fix pow (qNaN, 0) result with -lieee (bug 20919), remove dead parts of wrappers.
The dbl-64 implementation of __ieee754_pow returns a NaN for pow
(qNaN, 0) when it should return 1. Normally this is covered up by the
wrappers ending up calling __kernel_standard which fixes up the result
for this case, but for -lieee the wrappers are bypassed and the bad
result gets through as a return value.
Now, the wrappers fixing this are dealing with variant error handling
that wants a result of NaN for pow (qNaN, 0), and only ever call
__kernel_standard for this case if NaN resulted from __ieee754_pow.
This leads to a question of whether the dbl-64 code might be
deliberately returning NaN in order to use those code paths. However,
I can find no sign that this is deliberate. If it were deliberate one
would expect other implementations to do the same, and would expect
the return of NaN to be very old, but it appears it came in by
accident when the present e_pow.c implementation replaced an fdlibm
implementation in 2001. So it appears to be unintended that this path
through the pow wrapper could be used at all.
So this patch fixes the implementation to return 1 in this case as
expected. This is consistent with all the other implementations. The
relevant path through the wrappers is now unreachable, so is removed
(which is the main motivation of this patch: to avoid that path
becoming accidentally reachable when implementing TS 18661-1 semantics
that pow (sNaN, 0) should return qNaN with "invalid" raised). Another
path that would require __ieee754_pow (0, 0) to return 0 is also
unreachable (as all implementations return 1, in accordance with C99
semantics), so is removed as well.
Note: we don't have anything set up to test -lieee, which in any case
is obsolescent (at some point we should remove the ability for new
programs to access _LIB_VERSION or define matherr and have it called
by glibc). So testing will be implicit through sNaN tests added when
making sNaN inputs work correctly for pow functions.
Tested for x86_64 and x86.
[BZ #20919]
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Do not return
NaN first argument when raised to power 0.
* math/w_pow.c (__pow): Do not check for NaN or zero results from
raising to power zero.
* math/w_powf.c (__powf): Likewise.
* math/w_powl.c (__powl): Likewise.
* sysdeps/ieee754/k_standard.c (__kernel_standard): Do not handle
pow (0, 0) or pow (NaN, 0).
2016-12-02 23:50:46 +01:00
|
|
|
if ((((qx == 0x7ff00000) && (u.i[LOW_HALF] != 0)) || (qx > 0x7ff00000))
|
Fix sysdeps/ieee754 pow handling of sNaN arguments (bug 20916).
Various pow function implementations mishandle sNaN arguments in
various ways. This includes returning sNaN instead of qNaN for sNaN
arguments. For arguments (1, sNaN) and (sNaN, 0), TS 18661-1
semantics are also that the result should be qNaN, whereas with a qNaN
argument there the result should be 1, but for the dbl-64
implementation of pow there are issues with sNaN arguments beyond not
implementing the TS 18661-1 semantics in those special cases.
This patch makes the implementations in sysdeps/ieee754 follow the TS
18661-1 semantics consistently. Because x86 / x86_64 implementations
still need fixing, testcases are not included with this patch; they
will be included with the fix for the x86 / x86_64 versions.
Tested for x86_64, x86, mips64 and powerpc (with such testcases, which
pass in the mips64 and powerpc cases).
[BZ #20916]
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Do not return 1
for arguments (sNaN, 0) or (1, sNaN). Do arithmetic on NaN
arguments to compute result.
* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Do not return
1 for arguments (sNaN, 0) or (1, sNaN).
* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Likewise.
2016-12-03 00:21:15 +01:00
|
|
|
&& (y != 0 || issignaling (x)))
|
|
|
|
return x + x;
|
2013-10-08 12:53:16 +02:00
|
|
|
if (y == 1.0)
|
|
|
|
return x;
|
|
|
|
if (y == 2.0)
|
|
|
|
return x * x;
|
|
|
|
if (y == -1.0)
|
|
|
|
return 1.0 / x;
|
|
|
|
if (y == 0)
|
|
|
|
return 1.0;
|
|
|
|
}
|
2001-03-12 01:04:52 +01:00
|
|
|
/* else */
|
2013-10-08 12:53:16 +02:00
|
|
|
if (((u.i[HIGH_HALF] > 0 && u.i[HIGH_HALF] < 0x7ff00000) || /* x>0 and not x->0 */
|
|
|
|
(u.i[HIGH_HALF] == 0 && u.i[LOW_HALF] != 0)) &&
|
|
|
|
/* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
|
|
|
|
(v.i[HIGH_HALF] & 0x7fffffff) < 0x4ff00000)
|
|
|
|
{ /* if y<-1 or y>1 */
|
|
|
|
double retval;
|
2012-03-05 13:22:46 +01:00
|
|
|
|
2014-06-23 22:12:33 +02:00
|
|
|
{
|
|
|
|
SET_RESTORE_ROUND (FE_TONEAREST);
|
2012-03-05 13:22:46 +01:00
|
|
|
|
2014-06-23 22:12:33 +02:00
|
|
|
/* Avoid internal underflow for tiny y. The exact value of y does
|
|
|
|
not matter if |y| <= 2**-64. */
|
2015-05-15 12:53:55 +02:00
|
|
|
if (fabs (y) < 0x1p-64)
|
2014-06-23 22:12:33 +02:00
|
|
|
y = y < 0 ? -0x1p-64 : 0x1p-64;
|
2018-02-12 11:42:42 +01:00
|
|
|
z = log1 (x, &aa); /* x^y =e^(y log (X)) */
|
2014-06-23 22:12:33 +02:00
|
|
|
t = y * CN;
|
|
|
|
y1 = t - (t - y);
|
|
|
|
y2 = y - y1;
|
|
|
|
t = z * CN;
|
|
|
|
a1 = t - (t - z);
|
|
|
|
a2 = (z - a1) + aa;
|
|
|
|
a = y1 * a1;
|
|
|
|
aa = y2 * a1 + y * a2;
|
|
|
|
a1 = a + aa;
|
|
|
|
a2 = (a - a1) + aa;
|
2018-02-12 11:42:42 +01:00
|
|
|
|
|
|
|
/* Maximum relative error RElog of log1 is 1.0e-21 (69.7 bits).
|
|
|
|
Maximum relative error REexp of __exp1 is 8.8e-22 (69.9 bits).
|
|
|
|
We actually compute exp ((1 + RElog) * log (x) * y) * (1 + REexp).
|
|
|
|
Since RElog/REexp are tiny and log (x) * y is at most log (DBL_MAX),
|
|
|
|
this is equivalent to pow (x, y) * (1 + 710 * RElog + REexp).
|
|
|
|
So the relative error is 710 * 1.0e-21 + 8.8e-22 = 7.1e-19
|
|
|
|
(60.2 bits). The worst-case ULP error is 0.5064. */
|
|
|
|
|
|
|
|
retval = __exp1 (a1, a2);
|
2014-06-23 22:12:33 +02:00
|
|
|
}
|
2012-03-05 13:22:46 +01:00
|
|
|
|
2015-06-03 16:36:34 +02:00
|
|
|
if (isinf (retval))
|
2014-06-23 22:12:33 +02:00
|
|
|
retval = huge * huge;
|
|
|
|
else if (retval == 0)
|
|
|
|
retval = tiny * tiny;
|
2015-09-26 00:29:10 +02:00
|
|
|
else
|
|
|
|
math_check_force_underflow_nonneg (retval);
|
2013-10-08 12:53:16 +02:00
|
|
|
return retval;
|
|
|
|
}
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
if (x == 0)
|
|
|
|
{
|
|
|
|
if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
|
|
|
|
|| (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000) /* NaN */
|
Fix sysdeps/ieee754 pow handling of sNaN arguments (bug 20916).
Various pow function implementations mishandle sNaN arguments in
various ways. This includes returning sNaN instead of qNaN for sNaN
arguments. For arguments (1, sNaN) and (sNaN, 0), TS 18661-1
semantics are also that the result should be qNaN, whereas with a qNaN
argument there the result should be 1, but for the dbl-64
implementation of pow there are issues with sNaN arguments beyond not
implementing the TS 18661-1 semantics in those special cases.
This patch makes the implementations in sysdeps/ieee754 follow the TS
18661-1 semantics consistently. Because x86 / x86_64 implementations
still need fixing, testcases are not included with this patch; they
will be included with the fix for the x86 / x86_64 versions.
Tested for x86_64, x86, mips64 and powerpc (with such testcases, which
pass in the mips64 and powerpc cases).
[BZ #20916]
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Do not return 1
for arguments (sNaN, 0) or (1, sNaN). Do arithmetic on NaN
arguments to compute result.
* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Do not return
1 for arguments (sNaN, 0) or (1, sNaN).
* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Likewise.
2016-12-03 00:21:15 +01:00
|
|
|
return y + y;
|
2015-05-15 12:53:55 +02:00
|
|
|
if (fabs (y) > 1.0e20)
|
2013-10-08 12:53:16 +02:00
|
|
|
return (y > 0) ? 0 : 1.0 / 0.0;
|
|
|
|
k = checkint (y);
|
|
|
|
if (k == -1)
|
|
|
|
return y < 0 ? 1.0 / x : x;
|
|
|
|
else
|
|
|
|
return y < 0 ? 1.0 / 0.0 : 0.0; /* return 0 */
|
|
|
|
}
|
2007-03-05 20:38:56 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
qx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign */
|
|
|
|
qy = v.i[HIGH_HALF] & 0x7fffffff; /* no sign */
|
2007-03-05 20:38:56 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0)) /* NaN */
|
Fix sysdeps/ieee754 pow handling of sNaN arguments (bug 20916).
Various pow function implementations mishandle sNaN arguments in
various ways. This includes returning sNaN instead of qNaN for sNaN
arguments. For arguments (1, sNaN) and (sNaN, 0), TS 18661-1
semantics are also that the result should be qNaN, whereas with a qNaN
argument there the result should be 1, but for the dbl-64
implementation of pow there are issues with sNaN arguments beyond not
implementing the TS 18661-1 semantics in those special cases.
This patch makes the implementations in sysdeps/ieee754 follow the TS
18661-1 semantics consistently. Because x86 / x86_64 implementations
still need fixing, testcases are not included with this patch; they
will be included with the fix for the x86 / x86_64 versions.
Tested for x86_64, x86, mips64 and powerpc (with such testcases, which
pass in the mips64 and powerpc cases).
[BZ #20916]
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Do not return 1
for arguments (sNaN, 0) or (1, sNaN). Do arithmetic on NaN
arguments to compute result.
* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Do not return
1 for arguments (sNaN, 0) or (1, sNaN).
* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Likewise.
2016-12-03 00:21:15 +01:00
|
|
|
return x + y;
|
2013-10-08 12:53:16 +02:00
|
|
|
if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0)) /* NaN */
|
Fix sysdeps/ieee754 pow handling of sNaN arguments (bug 20916).
Various pow function implementations mishandle sNaN arguments in
various ways. This includes returning sNaN instead of qNaN for sNaN
arguments. For arguments (1, sNaN) and (sNaN, 0), TS 18661-1
semantics are also that the result should be qNaN, whereas with a qNaN
argument there the result should be 1, but for the dbl-64
implementation of pow there are issues with sNaN arguments beyond not
implementing the TS 18661-1 semantics in those special cases.
This patch makes the implementations in sysdeps/ieee754 follow the TS
18661-1 semantics consistently. Because x86 / x86_64 implementations
still need fixing, testcases are not included with this patch; they
will be included with the fix for the x86 / x86_64 versions.
Tested for x86_64, x86, mips64 and powerpc (with such testcases, which
pass in the mips64 and powerpc cases).
[BZ #20916]
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Do not return 1
for arguments (sNaN, 0) or (1, sNaN). Do arithmetic on NaN
arguments to compute result.
* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Do not return
1 for arguments (sNaN, 0) or (1, sNaN).
* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Likewise.
2016-12-03 00:21:15 +01:00
|
|
|
return x == 1.0 && !issignaling (y) ? 1.0 : y + y;
|
2007-03-05 20:38:56 +01:00
|
|
|
|
2001-03-12 01:04:52 +01:00
|
|
|
/* if x<0 */
|
2013-10-08 12:53:16 +02:00
|
|
|
if (u.i[HIGH_HALF] < 0)
|
|
|
|
{
|
|
|
|
k = checkint (y);
|
|
|
|
if (k == 0)
|
|
|
|
{
|
|
|
|
if (qy == 0x7ff00000)
|
|
|
|
{
|
|
|
|
if (x == -1.0)
|
|
|
|
return 1.0;
|
|
|
|
else if (x > -1.0)
|
|
|
|
return v.i[HIGH_HALF] < 0 ? INF.x : 0.0;
|
|
|
|
else
|
|
|
|
return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x;
|
|
|
|
}
|
|
|
|
else if (qx == 0x7ff00000)
|
|
|
|
return y < 0 ? 0.0 : INF.x;
|
|
|
|
return (x - x) / (x - x); /* y not integer and x<0 */
|
|
|
|
}
|
2007-03-05 20:38:56 +01:00
|
|
|
else if (qx == 0x7ff00000)
|
2013-10-08 12:53:16 +02:00
|
|
|
{
|
|
|
|
if (k < 0)
|
|
|
|
return y < 0 ? nZERO.x : nINF.x;
|
|
|
|
else
|
|
|
|
return y < 0 ? 0.0 : INF.x;
|
|
|
|
}
|
|
|
|
/* if y even or odd */
|
2014-06-23 22:12:33 +02:00
|
|
|
if (k == 1)
|
|
|
|
return __ieee754_pow (-x, y);
|
|
|
|
else
|
|
|
|
{
|
|
|
|
double retval;
|
|
|
|
{
|
|
|
|
SET_RESTORE_ROUND (FE_TONEAREST);
|
|
|
|
retval = -__ieee754_pow (-x, y);
|
|
|
|
}
|
2015-06-03 16:36:34 +02:00
|
|
|
if (isinf (retval))
|
2014-06-23 22:12:33 +02:00
|
|
|
retval = -huge * huge;
|
|
|
|
else if (retval == 0)
|
|
|
|
retval = -tiny * tiny;
|
|
|
|
return retval;
|
|
|
|
}
|
2001-03-13 03:01:34 +01:00
|
|
|
}
|
2001-03-12 01:04:52 +01:00
|
|
|
/* x>0 */
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
if (qx == 0x7ff00000) /* x= 2^-0x3ff */
|
2013-04-03 21:10:02 +02:00
|
|
|
return y > 0 ? x : 0;
|
2001-02-20 00:16:05 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
if (qy > 0x45f00000 && qy < 0x7ff00000)
|
|
|
|
{
|
|
|
|
if (x == 1.0)
|
|
|
|
return 1.0;
|
|
|
|
if (y > 0)
|
|
|
|
return (x > 1.0) ? huge * huge : tiny * tiny;
|
|
|
|
if (y < 0)
|
|
|
|
return (x < 1.0) ? huge * huge : tiny * tiny;
|
|
|
|
}
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
if (x == 1.0)
|
|
|
|
return 1.0;
|
|
|
|
if (y > 0)
|
|
|
|
return (x > 1.0) ? INF.x : 0;
|
|
|
|
if (y < 0)
|
|
|
|
return (x < 1.0) ? INF.x : 0;
|
|
|
|
return 0; /* unreachable, to make the compiler happy */
|
2001-03-12 01:04:52 +01:00
|
|
|
}
|
2013-10-08 12:53:16 +02:00
|
|
|
|
2011-10-25 02:19:17 +02:00
|
|
|
#ifndef __ieee754_pow
|
2011-10-12 17:27:51 +02:00
|
|
|
strong_alias (__ieee754_pow, __pow_finite)
|
2011-10-25 02:19:17 +02:00
|
|
|
#endif
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
/* Compute log(x) (x is left argument). The result is the returned double + the
|
2018-02-12 11:42:42 +01:00
|
|
|
parameter DELTA. */
|
2011-10-25 06:56:33 +02:00
|
|
|
static double
|
|
|
|
SECTION
|
2018-02-12 11:42:42 +01:00
|
|
|
log1 (double x, double *delta)
|
2013-10-08 12:53:16 +02:00
|
|
|
{
|
2016-01-04 16:28:52 +01:00
|
|
|
unsigned int i, j;
|
|
|
|
int m;
|
2013-10-08 12:53:16 +02:00
|
|
|
double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0;
|
|
|
|
mynumber u, v;
|
2002-08-24 03:36:09 +02:00
|
|
|
#ifdef BIG_ENDI
|
2013-10-08 12:53:16 +02:00
|
|
|
mynumber /**/ two52 = {{0x43300000, 0x00000000}}; /* 2**52 */
|
2002-08-24 03:36:09 +02:00
|
|
|
#else
|
2013-10-08 12:53:16 +02:00
|
|
|
# ifdef LITTLE_ENDI
|
|
|
|
mynumber /**/ two52 = {{0x00000000, 0x43300000}}; /* 2**52 */
|
|
|
|
# endif
|
2002-08-24 03:36:09 +02:00
|
|
|
#endif
|
update from main archive 961030
Thu Oct 31 00:01:39 1996 Ulrich Drepper <drepper@cygnus.com>
* signal/Makefile (routines): Add sigwait.
* signal/signal.h: Add prototype for sigwait.
* sysdeps/posix/sigwait.c: New file. Implementation of sigwait
function from POSIX.1c.
* sysdeps/stub/sigwait.c: New file. Stub version of sigwait.
Wed Oct 30 02:01:17 1996 Richard Henderson <rth@tamu.edu>
* sunrpc/xdr_float.c (xdr_float): Handle sizeof(float)!=sizeof(long),
but don't bother going farther than sizeof(float)==sizeof(int).
(xdr_double): Handle little-endian machines! Handle sizeof(double)
!= 2*sizeof(long), though again don't bother with more than int.
Thu Oct 29 16:09:42 1996 Craig Metz <cmetz@inner.net>
* sysdeps/posix/getaddrinfo.c: Use buffer limits for inet_ntop
function.
Tue Oct 29 12:37:22 1996 Ulrich Drepper <drepper@cygnus.com>
* Makerules: Create symbolic links for linking in $(libdir).
(make-link): Use absolute path for destination if this is not in
the same directory.
* elf/rtld.c (dl_main): When verifying don't check the name of
the dynamic linker.
* shlib-versions: Change entries for Hurd specific libs from
*-*-gnu* to *-*-gnu?* so that i586-pc-linux-gnu does not match
these entries.
* assert/assert.h: Reformat copyright.
Change reference to ANSI into reference to ISO C.
* ctype/ctype.h: Likewise.
* errno.h: Likewise.
* limits.h: Likewise.
* math/math.h: Likewise.
* setjmp/setjmp.h: Likewise.
* stdio/stdio.h: Likewise.
* libio/stdio.h: Likewise.
* stdlib/stdlib.h: Likewise.
* string/string.h: Likewise.
* time/time.h: Likewise.
* string/argz.h: Use __const is definitions.
* elf/dlfcn.h: Use __const and __P. Reformat copyright.
* misc/err.h: Likewise.
* wctype/wctype.h (wctrans_t): Use __const instead of const.
* Makeconfig ($(common-objpfx)soversions.mk): Generate list of
sonames for versioned libraries.
* Makefile: Remove code to generate libc-version.h.
Generate gnu/lib-names.h with info from soversions.mk.
* features.h: Define __GLIBC__ and __GLIBC_MINOR__.
* dirent/tst-seekdir.c: Initialize save3.
* grp/testgrp.c: Initialize my_group.
* grp/fgetgrent_r.c: Change interface to follow POSIX.1c.
* grp/grp.h: Likewise.
* nss/getXXbyYY.c: Likewise.
* nss/getXXbyYY_r.c: Likewise.
* nss/getXXent.c: Likewise.
* nss/getXXent_r.c: Likewise.
* pwd/fgetpwent_r.c: Likewise.
* pwd/pwd.h: Likewise.
* shadow/fgetspent_r.c: Likewise.
* shadow/sgetspent.c: Likewise.
* shadow/sgetspent_r.c: Likewise.
* grp/fgetgrent.c: Adapt for change in interface of fgetgrent_r.
* pwd/fgetpwent.c: Likewise, for fgetpwent_r.c.
* shadow/fgetspent.c: Likewise, for fgetpwent_r.c.
* resolv/netdb.h: Adapt prototypes for reentrant functions to
follow POSIX.1c.
* sunrpc/rpc/netdb.h: Likewise,
* shadow/shadow.h: Likewise.
* inet/getnetgrent_r.c: Follow change in pwd/grp function interface.
* sysdeps/unix/getlogin_r.c: Return ERANGE when buffer is too small.
* inet/herrno.c: Don't define __h_errno. Only h_errno otherwise the
ELF aliasing creates strange situations.
* sysdeps/unix/sysv/linux/errnos.H: Define __set_errno as inline
function.
* sysdeps/unix/sysv/linux/i386/sysdep.S: Don't define __errno.
* sysdeps/unix/sysv/linux/m68k/sysdep.S: Likewise.
* libio/libio.h: Don't declare _IO_flockfile and _IO_funlockfile
weak.
* locale/programs/charmap.c: Add casts to prevent warnings.
* locale/programs/linereader.h: Likewise.
* locale/programs/ld-collate.c: Likewise.
* locale/programs/stringtrans.c: Likewise.
Change types for various variables to prevent warnings.
* locale/programs/ld-ctype.c: Likewise.
* locale/programs/linereader.h (lr_ungetc): Likewise.
* locale/programs/charset.h (struct charset): Use `unsigned int'
as type for width_default.
* posix/regex.c: Change type of `this_reg' variables.
* stdio-common/Makefile: Use -Wno-format for tstdiomisc.c.
* stdio-common/bug5.c: De-ANSI-fy. Use correct types for
variables.
* stdio-common/printf_fp.c: Initialize to_shift.
* stdio-common/test_rdwr.c: Add cast.
* stdio-common/vfprintf.c: Add casts and use correct types to
prevent warnings.
* stdio-common/vfscanf.c: Initialize str and strptr.
* sysdeps/libm-ieee754/e_jnf.c: Use correct types to prevent warnings.
* sysdeps/libm-ieee754/e_pow.c: Likewise.
* sysdeps/libm-ieee754/e_powf.c: Likewise.
* sysdeps/libm-ieee754/e_rem_pio2f.c: Likewise.
* time/test-tz.c: Likewise.
* manual/creature.texi: Document _REENTRANT and _THREAD_SAFE.
* manual/libc.texinfo: Prevent makeinfo failure by avoiding
libc.cp index. This must be fixed.
* manual/nss.texi: Adapt for correct POSIX.1c interface of
reentrant functions.
* manual/users.texi: Document netgroup functions.
* po/es.po: Updated.
* po/fr.po: Updated.
* posix/fnmatch.c: Change to match libit version.
* posix/unistd.h: Change prototype for ttyname_r to match POSIX.1c.
* sysdep/posix/ttyname_r.c: Likewise.
* stdlib/atexit.h (__new_exitfn): Add internal locking.
* stdlib/exit.c: De-ANSI-fy. Handle new ef_us value for flavor.
* stdlib/exit.h: De-ANSI-fy. Define new ef_us value for flavor.
* stdlib/random.c (__srandom): Add internal locking.
(__initstate): Likewise.
(__setstate): Likewise.
(__random): Likewise.
Mon Oct 28 22:28:37 1996 NIIBE Yutaka <gniibe@mri.co.jp>
* sysdeps/generic/crypt-entry.c (crypt_r): Use __set_errno.
(crypt): Likewise.
* resolv/gethnamaddr.c (gethostbyname2): Likewise.
* sysdeps/generic/uname.c: Likewise.
* sysdeps/posix/rename.c: Likewise.
* sysdeps/stub/setrlimit.c: Likewise.
* nss/nss_db/db-netgrp.c (_nss_db_setnetgrent): Fix typo.
Sun Oct 27 11:12:50 1996 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* locale/programs/ld-collate.c (collate_order_elem): Fix format
string.
(collate_element_to): Cast field width argument to `int' for
format string.
(collate_symbol): Likewise.
(collate_order_elem): Likewise.
(collate_weight_bsymbol): Likewise.
(collate_simple_weight): Likewise.
* locale/programs/ld-time.c (STRARR_ELEM): Fix format string.
* locale/programs/ld-ctype.c (ctype_class_newP): Add missing
argument for format string.
(ctype_map_newP): Likewise.
(set_class_defaults): Fix format string.
* locale/programs/localedef.c (construct_output_path): Putting an
explicit \0 into the format string does not work, use %c.
Sat Oct 26 20:38:36 1996 Richard Henderson <rth@tamu.edu>
* Makerules: Install all shared libraries in $(slibdir).
* login/Makefile: Build libutil.so in others pass after
libc.so is created.
* misc/mntent.h: Include <paths.h> for _PATH_MNTTAB & _PATH_MOUNTED.
* string/stratcliff.c: Allocate 3 pages instead of one, then use
mprotect so that we know that the adjacent pages are inaccessible.
* resource/sys/resource.h: Move all structures and enums to ...
* sysdeps/generic/resourcebits.h: ... here ...
* sysdeps/unix/bsd/sun/sunos4/resourcebits.h: ... and here.
* sysdeps/unix/sysv/linux/alpha/resourcebits.h: Remove.
* sysdeps/unix/sysv/linux/i386/resourcebits.h: Remove.
* sysdeps/unix/sysv/linux/m68k/resourcebits.h: Remove.
* sysdeps/unix/sysv/linux/mips/resourcebits.h: Remove.
* sysdeps/unix/sysv/linux/resourcebits.h: New file. Use kernel
header for RLIMIT_* definitions. The members of struct rlimit
are longs.
Thu Oct 24 17:43:34 1996 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* MakeTAGS (sysdep-dirs): Fix typo.
Wed Oct 23 03:45:22 1996 Ulrich Drepper <drepper@cygnus.com>
* Makefile (headers): Don't mention libc-version.h.
(install-others): ...but here.
* time/strptime.c: Recognize %s, %u, %g, and %G format.
nothing is found. This guarantees all subsequent calls behave
* sysdeps/unix/sysv/linux/syscalls.list: Change function name for
* io/getwd.c (getwd) [! PATH_MAX]: Don't assume that the user's
buffer is any longer than the amount necessary to hold the
filename; the Hurd getcwd uses the *entire* contents of the
buffer, however long it is specified to be.
* posix/getconf.c: De-ANSI-fy. Recognize POSIX.2 constant names.
since these do not depend on the platform.
1996-10-31 03:57:12 +01:00
|
|
|
|
2001-03-12 01:04:52 +01:00
|
|
|
u.x = x;
|
|
|
|
m = u.i[HIGH_HALF];
|
2018-02-12 11:42:42 +01:00
|
|
|
if (m < 0x00100000) /* Handle denormal x. */
|
2013-10-08 12:53:16 +02:00
|
|
|
{
|
|
|
|
x = x * t52.x;
|
|
|
|
add = -52.0;
|
|
|
|
u.x = x;
|
|
|
|
m = u.i[HIGH_HALF];
|
|
|
|
}
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
if ((m & 0x000fffff) < 0x0006a09e)
|
|
|
|
{
|
|
|
|
u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3ff00000;
|
|
|
|
two52.i[LOW_HALF] = (m >> 20);
|
|
|
|
}
|
2001-03-12 01:04:52 +01:00
|
|
|
else
|
2013-10-08 12:53:16 +02:00
|
|
|
{
|
|
|
|
u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3fe00000;
|
|
|
|
two52.i[LOW_HALF] = (m >> 20) + 1;
|
|
|
|
}
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2001-03-12 01:04:52 +01:00
|
|
|
v.x = u.x + bigu.x;
|
|
|
|
uu = v.x - bigu.x;
|
2013-10-08 12:53:16 +02:00
|
|
|
i = (v.i[LOW_HALF] & 0x000003ff) << 2;
|
2018-02-12 11:42:42 +01:00
|
|
|
if (two52.i[LOW_HALF] == 1023) /* Exponent of x is 0. */
|
2013-10-08 12:53:16 +02:00
|
|
|
{
|
|
|
|
if (i > 1192 && i < 1208) /* |x-1| < 1.5*2**-10 */
|
|
|
|
{
|
2001-03-12 01:04:52 +01:00
|
|
|
t = x - 1.0;
|
2013-10-08 12:53:16 +02:00
|
|
|
t1 = (t + 5.0e6) - 5.0e6;
|
|
|
|
t2 = t - t1;
|
|
|
|
e1 = t - 0.5 * t1 * t1;
|
|
|
|
e2 = (t * t * t * (r3 + t * (r4 + t * (r5 + t * (r6 + t
|
|
|
|
* (r7 + t * r8)))))
|
|
|
|
- 0.5 * t2 * (t + t1));
|
|
|
|
res = e1 + e2;
|
|
|
|
*delta = (e1 - res) + e2;
|
2018-02-12 11:42:42 +01:00
|
|
|
/* Max relative error is 1.464844e-24, so accurate to 79.1 bits. */
|
2001-03-12 01:04:52 +01:00
|
|
|
return res;
|
2013-10-08 12:53:16 +02:00
|
|
|
} /* |x-1| < 1.5*2**-10 */
|
2001-03-12 01:04:52 +01:00
|
|
|
else
|
2013-10-08 12:53:16 +02:00
|
|
|
{
|
|
|
|
v.x = u.x * (ui.x[i] + ui.x[i + 1]) + bigv.x;
|
|
|
|
vv = v.x - bigv.x;
|
|
|
|
j = v.i[LOW_HALF] & 0x0007ffff;
|
|
|
|
j = j + j + j;
|
|
|
|
eps = u.x - uu * vv;
|
|
|
|
e1 = eps * ui.x[i];
|
|
|
|
e2 = eps * (ui.x[i + 1] + vj.x[j] * (ui.x[i] + ui.x[i + 1]));
|
|
|
|
e = e1 + e2;
|
|
|
|
e2 = ((e1 - e) + e2);
|
|
|
|
t = ui.x[i + 2] + vj.x[j + 1];
|
|
|
|
t1 = t + e;
|
|
|
|
t2 = ((((t - t1) + e) + (ui.x[i + 3] + vj.x[j + 2])) + e2 + e * e
|
|
|
|
* (p2 + e * (p3 + e * p4)));
|
|
|
|
res = t1 + t2;
|
|
|
|
*delta = (t1 - res) + t2;
|
2018-02-12 11:42:42 +01:00
|
|
|
/* Max relative error is 1.0e-24, so accurate to 79.7 bits. */
|
2001-03-12 01:04:52 +01:00
|
|
|
return res;
|
2013-10-08 12:53:16 +02:00
|
|
|
}
|
2018-02-12 11:42:42 +01:00
|
|
|
}
|
|
|
|
else /* Exponent of x != 0. */
|
2013-10-08 12:53:16 +02:00
|
|
|
{
|
2001-03-12 01:04:52 +01:00
|
|
|
eps = u.x - uu;
|
2013-10-08 12:53:16 +02:00
|
|
|
nx = (two52.x - two52e.x) + add;
|
|
|
|
e1 = eps * ui.x[i];
|
|
|
|
e2 = eps * ui.x[i + 1];
|
|
|
|
e = e1 + e2;
|
|
|
|
e2 = (e1 - e) + e2;
|
|
|
|
t = nx * ln2a.x + ui.x[i + 2];
|
|
|
|
t1 = t + e;
|
|
|
|
t2 = ((((t - t1) + e) + nx * ln2b.x + ui.x[i + 3] + e2) + e * e
|
|
|
|
* (q2 + e * (q3 + e * (q4 + e * (q5 + e * q6)))));
|
|
|
|
res = t1 + t2;
|
|
|
|
*delta = (t1 - res) + t2;
|
2018-02-12 11:42:42 +01:00
|
|
|
/* Max relative error is 1.0e-21, so accurate to 69.7 bits. */
|
2001-03-12 01:04:52 +01:00
|
|
|
return res;
|
2013-10-08 12:53:16 +02:00
|
|
|
}
|
2001-03-12 01:04:52 +01:00
|
|
|
}
|
1996-03-05 22:41:30 +01:00
|
|
|
|
2018-02-12 11:42:42 +01:00
|
|
|
|
2013-10-08 12:53:16 +02:00
|
|
|
/* This function receives a double x and checks if it is an integer. If not,
|
|
|
|
it returns 0, else it returns 1 if even or -1 if odd. */
|
2011-10-25 06:56:33 +02:00
|
|
|
static int
|
|
|
|
SECTION
|
2013-10-08 12:53:16 +02:00
|
|
|
checkint (double x)
|
|
|
|
{
|
|
|
|
union
|
|
|
|
{
|
|
|
|
int4 i[2];
|
|
|
|
double x;
|
|
|
|
} u;
|
2017-12-19 19:41:01 +01:00
|
|
|
int k;
|
|
|
|
unsigned int m, n;
|
2001-03-12 01:04:52 +01:00
|
|
|
u.x = x;
|
2013-10-08 12:53:16 +02:00
|
|
|
m = u.i[HIGH_HALF] & 0x7fffffff; /* no sign */
|
|
|
|
if (m >= 0x7ff00000)
|
|
|
|
return 0; /* x is +/-inf or NaN */
|
|
|
|
if (m >= 0x43400000)
|
|
|
|
return 1; /* |x| >= 2**53 */
|
|
|
|
if (m < 0x40000000)
|
|
|
|
return 0; /* |x| < 2, can not be 0 or 1 */
|
2001-03-12 01:04:52 +01:00
|
|
|
n = u.i[LOW_HALF];
|
2013-10-08 12:53:16 +02:00
|
|
|
k = (m >> 20) - 1023; /* 1 <= k <= 52 */
|
|
|
|
if (k == 52)
|
|
|
|
return (n & 1) ? -1 : 1; /* odd or even */
|
|
|
|
if (k > 20)
|
|
|
|
{
|
2016-10-14 21:53:27 +02:00
|
|
|
if (n << (k - 20) != 0)
|
2013-10-08 12:53:16 +02:00
|
|
|
return 0; /* if not integer */
|
2016-10-14 21:53:27 +02:00
|
|
|
return (n << (k - 21) != 0) ? -1 : 1;
|
2013-10-08 12:53:16 +02:00
|
|
|
}
|
|
|
|
if (n)
|
|
|
|
return 0; /*if not integer */
|
|
|
|
if (k == 20)
|
|
|
|
return (m & 1) ? -1 : 1;
|
2016-10-14 21:53:27 +02:00
|
|
|
if (m << (k + 12) != 0)
|
2013-10-08 12:53:16 +02:00
|
|
|
return 0;
|
2016-10-14 21:53:27 +02:00
|
|
|
return (m << (k + 11) != 0) ? -1 : 1;
|
1996-03-05 22:41:30 +01:00
|
|
|
}
|