mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-21 20:40:05 +00:00
Add new pow implementation
The algorithm is exp(y * log(x)), where log(x) is computed with about 1.3*2^-68 relative error (1.5*2^-68 without fma), returning the result in two doubles, and the exp part uses the same algorithm (and lookup tables) as exp, but takes the input as two doubles and a sign (to handle negative bases with odd integer exponent). The __exp1 internal symbol is no longer necessary. There is separate code path when fma is not available but the worst case error is about 0.54 ULP in both cases. The lookup table and consts for log are 4168 bytes. The .rodata+.text is decreased by 37908 bytes on aarch64. The non-nearest rounding error is less than 1 ULP. Improvements on Cortex-A72 compared to current glibc master: pow thruput: 2.40x in [0.01 11.1]x[0.01 11.1] pow latency: 1.84x in [0.01 11.1]x[0.01 11.1] Tested on aarch64-linux-gnu (defined __FP_FAST_FMA, TOINT_INTRINSICS) and arm-linux-gnueabihf (!defined __FP_FAST_FMA, !TOINT_INTRINSICS) and x86_64-linux-gnu (!defined __FP_FAST_FMA, !TOINT_INTRINSICS) and powerpc64le-linux-gnu (defined __FP_FAST_FMA, !TOINT_INTRINSICS) targets. * NEWS: Mention pow improvements. * math/Makefile (type-double-routines): Add e_pow_log_data. * sysdeps/generic/math_private.h (__exp1): Remove. * sysdeps/i386/fpu/e_pow_log_data.c: New file. * sysdeps/ia64/fpu/e_pow_log_data.c: New file. * sysdeps/ieee754/dbl-64/Makefile (CFLAGS-e_pow.c): Allow fma contraction. * sysdeps/ieee754/dbl-64/e_exp.c (__exp1): Remove. (exp_inline): Remove. (__ieee754_exp): Only single double input is handled. * sysdeps/ieee754/dbl-64/e_pow.c: Rewrite. * sysdeps/ieee754/dbl-64/e_pow_log_data.c: New file. * sysdeps/ieee754/dbl-64/math_config.h (issignaling_inline): Define. (__pow_log_data): Define. * sysdeps/ieee754/dbl-64/upow.h: Remove. * sysdeps/ieee754/dbl-64/upow.tbl: Remove. * sysdeps/m68k/m680x0/fpu/e_pow_log_data.c: New file. * sysdeps/x86_64/fpu/multiarch/Makefile (CFLAGS-e_pow-fma.c): Allow fma contraction. (CFLAGS-e_pow-fma4.c): Likewise.
This commit is contained in:
Szabolcs Nagy
parent
dab9c3488e
commit
424c4f60ed
23
ChangeLog
23
ChangeLog
@ -1,3 +1,26 @@
|
||||
2018-09-19 Szabolcs Nagy <szabolcs.nagy@arm.com>
|
||||
|
||||
* NEWS: Mention pow improvements.
|
||||
* math/Makefile (type-double-routines): Add e_pow_log_data.
|
||||
* sysdeps/generic/math_private.h (__exp1): Remove.
|
||||
* sysdeps/i386/fpu/e_pow_log_data.c: New file.
|
||||
* sysdeps/ia64/fpu/e_pow_log_data.c: New file.
|
||||
* sysdeps/ieee754/dbl-64/Makefile (CFLAGS-e_pow.c): Allow fma
|
||||
contraction.
|
||||
* sysdeps/ieee754/dbl-64/e_exp.c (__exp1): Remove.
|
||||
(exp_inline): Remove.
|
||||
(__ieee754_exp): Only single double input is handled.
|
||||
* sysdeps/ieee754/dbl-64/e_pow.c: Rewrite.
|
||||
* sysdeps/ieee754/dbl-64/e_pow_log_data.c: New file.
|
||||
* sysdeps/ieee754/dbl-64/math_config.h (issignaling_inline): Define.
|
||||
(__pow_log_data): Define.
|
||||
* sysdeps/ieee754/dbl-64/upow.h: Remove.
|
||||
* sysdeps/ieee754/dbl-64/upow.tbl: Remove.
|
||||
* sysdeps/m68k/m680x0/fpu/e_pow_log_data.c: New file.
|
||||
* sysdeps/x86_64/fpu/multiarch/Makefile (CFLAGS-e_pow-fma.c): Allow fma
|
||||
contraction.
|
||||
(CFLAGS-e_pow-fma4.c): Likewise.
|
||||
|
||||
2018-09-18 Paul Eggert <eggert@cs.ucla.edu>
|
||||
|
||||
Simplify tzfile fstat failure code
|
||||
|
||||
2
NEWS
2
NEWS
@ -16,7 +16,7 @@ Major new features:
|
||||
to set the install root if you wish to install into a non-default
|
||||
configured location.
|
||||
|
||||
* Optimized generic exp, exp2, log, log2, sinf, cosf, sincosf and tanf.
|
||||
* Optimized generic exp, exp2, log, log2, pow, sinf, cosf, sincosf and tanf.
|
||||
|
||||
* The reallocarray function is now declared under _DEFAULT_SOURCE, not just
|
||||
for _GNU_SOURCE, to match BSD environments.
|
||||
|
||||
@ -128,7 +128,7 @@ type-double-suffix :=
|
||||
type-double-routines := branred doasin dosincos mpa mpatan2 \
|
||||
k_rem_pio2 mpatan mpsqrt mptan sincos32 \
|
||||
sincostab math_err e_exp_data e_log_data \
|
||||
e_log2_data
|
||||
e_log2_data e_pow_log_data
|
||||
|
||||
# float support
|
||||
type-float-suffix := f
|
||||
|
||||
@ -225,7 +225,6 @@ do { \
|
||||
|
||||
|
||||
/* Prototypes for functions of the IBM Accurate Mathematical Library. */
|
||||
extern double __exp1 (double __x, double __xx);
|
||||
extern double __sin (double __x);
|
||||
extern double __cos (double __x);
|
||||
extern int __branred (double __x, double *__a, double *__aa);
|
||||
|
||||
1
sysdeps/i386/fpu/e_pow_log_data.c
Normal file
1
sysdeps/i386/fpu/e_pow_log_data.c
Normal file
@ -0,0 +1 @@
|
||||
/* Not needed. */
|
||||
1
sysdeps/ia64/fpu/e_pow_log_data.c
Normal file
1
sysdeps/ia64/fpu/e_pow_log_data.c
Normal file
@ -0,0 +1 @@
|
||||
/* Not needed. */
|
||||
@ -2,5 +2,4 @@ ifeq ($(subdir),math)
|
||||
# branred depends on precise IEEE double rounding
|
||||
CFLAGS-branred.c += $(config-cflags-nofma)
|
||||
CFLAGS-e_sqrt.c += $(config-cflags-nofma)
|
||||
CFLAGS-e_pow.c += $(config-cflags-nofma)
|
||||
endif
|
||||
|
||||
@ -85,10 +85,13 @@ top12 (double x)
|
||||
return asuint64 (x) >> 52;
|
||||
}
|
||||
|
||||
/* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
|
||||
If hastail is 0 then xtail is assumed to be 0 too. */
|
||||
static inline double
|
||||
exp_inline (double x, double xtail, int hastail)
|
||||
#ifndef SECTION
|
||||
# define SECTION
|
||||
#endif
|
||||
|
||||
double
|
||||
SECTION
|
||||
__ieee754_exp (double x)
|
||||
{
|
||||
uint32_t abstop;
|
||||
uint64_t ki, idx, top, sbits;
|
||||
@ -131,9 +134,6 @@ exp_inline (double x, double xtail, int hastail)
|
||||
kd -= Shift;
|
||||
#endif
|
||||
r = x + kd * NegLn2hiN + kd * NegLn2loN;
|
||||
/* The code assumes 2^-200 < |xtail| < 2^-8/N. */
|
||||
if (hastail)
|
||||
r += xtail;
|
||||
/* 2^(k/N) ~= scale * (1 + tail). */
|
||||
idx = 2 * (ki % N);
|
||||
top = ki << (52 - EXP_TABLE_BITS);
|
||||
@ -149,29 +149,10 @@ exp_inline (double x, double xtail, int hastail)
|
||||
if (__glibc_unlikely (abstop == 0))
|
||||
return specialcase (tmp, sbits, ki);
|
||||
scale = asdouble (sbits);
|
||||
/* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
|
||||
/* Note: tmp == 0 or |tmp| > 2^-65 and scale > 2^-739, so there
|
||||
is no spurious underflow here even without fma. */
|
||||
return scale + scale * tmp;
|
||||
}
|
||||
|
||||
#ifndef SECTION
|
||||
# define SECTION
|
||||
#endif
|
||||
|
||||
double
|
||||
SECTION
|
||||
__ieee754_exp (double x)
|
||||
{
|
||||
return exp_inline (x, 0, 0);
|
||||
}
|
||||
#ifndef __ieee754_exp
|
||||
strong_alias (__ieee754_exp, __exp_finite)
|
||||
#endif
|
||||
|
||||
/* Compute e^(x+xx). */
|
||||
double
|
||||
SECTION
|
||||
__exp1 (double x, double xx)
|
||||
{
|
||||
return exp_inline (x, xx, 1);
|
||||
}
|
||||
|
||||
@ -1,360 +1,380 @@
|
||||
/*
|
||||
* IBM Accurate Mathematical Library
|
||||
* written by International Business Machines Corp.
|
||||
* 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
|
||||
* 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: 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. */
|
||||
/* */
|
||||
/***************************************************************************/
|
||||
/* Double-precision x^y function.
|
||||
Copyright (C) 2018 Free Software Foundation, Inc.
|
||||
This file is part of the GNU C Library.
|
||||
|
||||
The GNU C Library 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.
|
||||
|
||||
The GNU C Library 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 the GNU C Library; if not, see
|
||||
<http://www.gnu.org/licenses/>. */
|
||||
|
||||
#include <math.h>
|
||||
#include "endian.h"
|
||||
#include "upow.h"
|
||||
#include <dla.h>
|
||||
#include "mydefs.h"
|
||||
#include "MathLib.h"
|
||||
#include "upow.tbl"
|
||||
#include <math_private.h>
|
||||
#include <fenv_private.h>
|
||||
#include <math-underflow.h>
|
||||
#include <fenv.h>
|
||||
#include <stdint.h>
|
||||
#include <math-barriers.h>
|
||||
#include <math-narrow-eval.h>
|
||||
#include "math_config.h"
|
||||
|
||||
/*
|
||||
Worst-case error: 0.54 ULP (~= ulperr_exp + 1024*Ln2*relerr_log*2^53)
|
||||
relerr_log: 1.3 * 2^-68 (Relative error of log, 1.5 * 2^-68 without fma)
|
||||
ulperr_exp: 0.509 ULP (ULP error of exp, 0.511 ULP without fma)
|
||||
*/
|
||||
|
||||
#define T __pow_log_data.tab
|
||||
#define A __pow_log_data.poly
|
||||
#define Ln2hi __pow_log_data.ln2hi
|
||||
#define Ln2lo __pow_log_data.ln2lo
|
||||
#define N (1 << POW_LOG_TABLE_BITS)
|
||||
#define OFF 0x3fe6955500000000
|
||||
|
||||
/* Top 12 bits of a double (sign and exponent bits). */
|
||||
static inline uint32_t
|
||||
top12 (double x)
|
||||
{
|
||||
return asuint64 (x) >> 52;
|
||||
}
|
||||
|
||||
/* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about
|
||||
additional 15 bits precision. IX is the bit representation of x, but
|
||||
normalized in the subnormal range using the sign bit for the exponent. */
|
||||
static inline double_t
|
||||
log_inline (uint64_t ix, double_t *tail)
|
||||
{
|
||||
/* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
|
||||
double_t z, r, y, invc, logc, logctail, kd, hi, t1, t2, lo, lo1, lo2, p;
|
||||
uint64_t iz, tmp;
|
||||
int k, i;
|
||||
|
||||
/* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
|
||||
The range is split into N subintervals.
|
||||
The ith subinterval contains z and c is near its center. */
|
||||
tmp = ix - OFF;
|
||||
i = (tmp >> (52 - POW_LOG_TABLE_BITS)) % N;
|
||||
k = (int64_t) tmp >> 52; /* arithmetic shift */
|
||||
iz = ix - (tmp & 0xfffULL << 52);
|
||||
z = asdouble (iz);
|
||||
kd = (double_t) k;
|
||||
|
||||
/* log(x) = k*Ln2 + log(c) + log1p(z/c-1). */
|
||||
invc = T[i].invc;
|
||||
logc = T[i].logc;
|
||||
logctail = T[i].logctail;
|
||||
|
||||
/* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and
|
||||
|z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */
|
||||
#ifdef __FP_FAST_FMA
|
||||
r = __builtin_fma (z, invc, -1.0);
|
||||
#else
|
||||
/* Split z such that rhi, rlo and rhi*rhi are exact and |rlo| <= |r|. */
|
||||
double_t zhi = asdouble ((iz + (1ULL << 31)) & (-1ULL << 32));
|
||||
double_t zlo = z - zhi;
|
||||
double_t rhi = zhi * invc - 1.0;
|
||||
double_t rlo = zlo * invc;
|
||||
r = rhi + rlo;
|
||||
#endif
|
||||
|
||||
/* k*Ln2 + log(c) + r. */
|
||||
t1 = kd * Ln2hi + logc;
|
||||
t2 = t1 + r;
|
||||
lo1 = kd * Ln2lo + logctail;
|
||||
lo2 = t1 - t2 + r;
|
||||
|
||||
/* Evaluation is optimized assuming superscalar pipelined execution. */
|
||||
double_t ar, ar2, ar3, lo3, lo4;
|
||||
ar = A[0] * r; /* A[0] = -0.5. */
|
||||
ar2 = r * ar;
|
||||
ar3 = r * ar2;
|
||||
/* k*Ln2 + log(c) + r + A[0]*r*r. */
|
||||
#ifdef __FP_FAST_FMA
|
||||
hi = t2 + ar2;
|
||||
lo3 = __builtin_fma (ar, r, -ar2);
|
||||
lo4 = t2 - hi + ar2;
|
||||
#else
|
||||
double_t arhi = A[0] * rhi;
|
||||
double_t arhi2 = rhi * arhi;
|
||||
hi = t2 + arhi2;
|
||||
lo3 = rlo * (ar + arhi);
|
||||
lo4 = t2 - hi + arhi2;
|
||||
#endif
|
||||
/* p = log1p(r) - r - A[0]*r*r. */
|
||||
p = (ar3
|
||||
* (A[1] + r * A[2] + ar2 * (A[3] + r * A[4] + ar2 * (A[5] + r * A[6]))));
|
||||
lo = lo1 + lo2 + lo3 + lo4 + p;
|
||||
y = hi + lo;
|
||||
*tail = hi - y + lo;
|
||||
return y;
|
||||
}
|
||||
|
||||
#undef N
|
||||
#undef T
|
||||
#define N (1 << EXP_TABLE_BITS)
|
||||
#define InvLn2N __exp_data.invln2N
|
||||
#define NegLn2hiN __exp_data.negln2hiN
|
||||
#define NegLn2loN __exp_data.negln2loN
|
||||
#define Shift __exp_data.shift
|
||||
#define T __exp_data.tab
|
||||
#define C2 __exp_data.poly[5 - EXP_POLY_ORDER]
|
||||
#define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
|
||||
#define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
|
||||
#define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
|
||||
#define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
|
||||
|
||||
/* Handle cases that may overflow or underflow when computing the result that
|
||||
is scale*(1+TMP) without intermediate rounding. The bit representation of
|
||||
scale is in SBITS, however it has a computed exponent that may have
|
||||
overflown into the sign bit so that needs to be adjusted before using it as
|
||||
a double. (int32_t)KI is the k used in the argument reduction and exponent
|
||||
adjustment of scale, positive k here means the result may overflow and
|
||||
negative k means the result may underflow. */
|
||||
static inline double
|
||||
specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
|
||||
{
|
||||
double_t scale, y;
|
||||
|
||||
if ((ki & 0x80000000) == 0)
|
||||
{
|
||||
/* k > 0, the exponent of scale might have overflowed by <= 460. */
|
||||
sbits -= 1009ull << 52;
|
||||
scale = asdouble (sbits);
|
||||
y = 0x1p1009 * (scale + scale * tmp);
|
||||
return check_oflow (y);
|
||||
}
|
||||
/* k < 0, need special care in the subnormal range. */
|
||||
sbits += 1022ull << 52;
|
||||
/* Note: sbits is signed scale. */
|
||||
scale = asdouble (sbits);
|
||||
y = scale + scale * tmp;
|
||||
if (fabs (y) < 1.0)
|
||||
{
|
||||
/* Round y to the right precision before scaling it into the subnormal
|
||||
range to avoid double rounding that can cause 0.5+E/2 ulp error where
|
||||
E is the worst-case ulp error outside the subnormal range. So this
|
||||
is only useful if the goal is better than 1 ulp worst-case error. */
|
||||
double_t hi, lo, one = 1.0;
|
||||
if (y < 0.0)
|
||||
one = -1.0;
|
||||
lo = scale - y + scale * tmp;
|
||||
hi = one + y;
|
||||
lo = one - hi + y + lo;
|
||||
y = math_narrow_eval (hi + lo) - one;
|
||||
/* Fix the sign of 0. */
|
||||
if (y == 0.0)
|
||||
y = asdouble (sbits & 0x8000000000000000);
|
||||
/* The underflow exception needs to be signaled explicitly. */
|
||||
math_force_eval (math_opt_barrier (0x1p-1022) * 0x1p-1022);
|
||||
}
|
||||
y = 0x1p-1022 * y;
|
||||
return check_uflow (y);
|
||||
}
|
||||
|
||||
#define SIGN_BIAS (0x800 << EXP_TABLE_BITS)
|
||||
|
||||
/* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
|
||||
The sign_bias argument is SIGN_BIAS or 0 and sets the sign to -1 or 1. */
|
||||
static inline double
|
||||
exp_inline (double x, double xtail, uint32_t sign_bias)
|
||||
{
|
||||
uint32_t abstop;
|
||||
uint64_t ki, idx, top, sbits;
|
||||
/* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
|
||||
double_t kd, z, r, r2, scale, tail, tmp;
|
||||
|
||||
abstop = top12 (x) & 0x7ff;
|
||||
if (__glibc_unlikely (abstop - top12 (0x1p-54)
|
||||
>= top12 (512.0) - top12 (0x1p-54)))
|
||||
{
|
||||
if (abstop - top12 (0x1p-54) >= 0x80000000)
|
||||
{
|
||||
/* Avoid spurious underflow for tiny x. */
|
||||
/* Note: 0 is common input. */
|
||||
double_t one = WANT_ROUNDING ? 1.0 + x : 1.0;
|
||||
return sign_bias ? -one : one;
|
||||
}
|
||||
if (abstop >= top12 (1024.0))
|
||||
{
|
||||
/* Note: inf and nan are already handled. */
|
||||
if (asuint64 (x) >> 63)
|
||||
return __math_uflow (sign_bias);
|
||||
else
|
||||
return __math_oflow (sign_bias);
|
||||
}
|
||||
/* Large x is special cased below. */
|
||||
abstop = 0;
|
||||
}
|
||||
|
||||
/* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */
|
||||
/* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */
|
||||
z = InvLn2N * x;
|
||||
#if TOINT_INTRINSICS
|
||||
/* z - kd is in [-0.5, 0.5] in all rounding modes. */
|
||||
kd = roundtoint (z);
|
||||
ki = converttoint (z);
|
||||
#else
|
||||
/* z - kd is in [-1, 1] in non-nearest rounding modes. */
|
||||
kd = math_narrow_eval (z + Shift);
|
||||
ki = asuint64 (kd);
|
||||
kd -= Shift;
|
||||
#endif
|
||||
r = x + kd * NegLn2hiN + kd * NegLn2loN;
|
||||
/* The code assumes 2^-200 < |xtail| < 2^-8/N. */
|
||||
r += xtail;
|
||||
/* 2^(k/N) ~= scale * (1 + tail). */
|
||||
idx = 2 * (ki % N);
|
||||
top = (ki + sign_bias) << (52 - EXP_TABLE_BITS);
|
||||
tail = asdouble (T[idx]);
|
||||
/* This is only a valid scale when -1023*N < k < 1024*N. */
|
||||
sbits = T[idx + 1] + top;
|
||||
/* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1). */
|
||||
/* Evaluation is optimized assuming superscalar pipelined execution. */
|
||||
r2 = r * r;
|
||||
/* Without fma the worst case error is 0.25/N ulp larger. */
|
||||
/* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp. */
|
||||
tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
|
||||
if (__glibc_unlikely (abstop == 0))
|
||||
return specialcase (tmp, sbits, ki);
|
||||
scale = asdouble (sbits);
|
||||
/* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
|
||||
is no spurious underflow here even without fma. */
|
||||
return scale + scale * tmp;
|
||||
}
|
||||
|
||||
/* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is
|
||||
the bit representation of a non-zero finite floating-point value. */
|
||||
static inline int
|
||||
checkint (uint64_t iy)
|
||||
{
|
||||
int e = iy >> 52 & 0x7ff;
|
||||
if (e < 0x3ff)
|
||||
return 0;
|
||||
if (e > 0x3ff + 52)
|
||||
return 2;
|
||||
if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
|
||||
return 0;
|
||||
if (iy & (1ULL << (0x3ff + 52 - e)))
|
||||
return 1;
|
||||
return 2;
|
||||
}
|
||||
|
||||
/* Returns 1 if input is the bit representation of 0, infinity or nan. */
|
||||
static inline int
|
||||
zeroinfnan (uint64_t i)
|
||||
{
|
||||
return 2 * i - 1 >= 2 * asuint64 (INFINITY) - 1;
|
||||
}
|
||||
|
||||
#ifndef SECTION
|
||||
# define SECTION
|
||||
#endif
|
||||
|
||||
static const double huge = 1.0e300, tiny = 1.0e-300;
|
||||
|
||||
double __exp1 (double x, double xx);
|
||||
static double log1 (double x, double *delta);
|
||||
static int checkint (double x);
|
||||
|
||||
/* An ultimate power routine. Given two IEEE double machine numbers y, x it
|
||||
computes the correctly rounded (to nearest) value of X^y. */
|
||||
double
|
||||
SECTION
|
||||
__ieee754_pow (double x, double y)
|
||||
{
|
||||
double z, a, aa, t, a1, a2, y1, y2;
|
||||
mynumber u, v;
|
||||
int k;
|
||||
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? */
|
||||
if ((((qx == 0x7ff00000) && (u.i[LOW_HALF] != 0)) || (qx > 0x7ff00000))
|
||||
&& (y != 0 || issignaling (x)))
|
||||
return x + x;
|
||||
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;
|
||||
}
|
||||
/* else */
|
||||
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;
|
||||
uint32_t sign_bias = 0;
|
||||
uint64_t ix, iy;
|
||||
uint32_t topx, topy;
|
||||
|
||||
{
|
||||
SET_RESTORE_ROUND (FE_TONEAREST);
|
||||
|
||||
/* Avoid internal underflow for tiny y. The exact value of y does
|
||||
not matter if |y| <= 2**-64. */
|
||||
if (fabs (y) < 0x1p-64)
|
||||
y = y < 0 ? -0x1p-64 : 0x1p-64;
|
||||
z = log1 (x, &aa); /* x^y =e^(y log (X)) */
|
||||
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;
|
||||
|
||||
/* Maximum relative error RElog of log1 is 1.0e-21 (69.7 bits).
|
||||
Maximum relative error REexp of __exp1 is 1.0e-18 (59.8 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 + 1.0e-18 = 1.7e-18
|
||||
(59 bits). The worst-case ULP error is 0.515. */
|
||||
|
||||
retval = __exp1 (a1, a2);
|
||||
}
|
||||
|
||||
if (isinf (retval))
|
||||
retval = huge * huge;
|
||||
else if (retval == 0)
|
||||
retval = tiny * tiny;
|
||||
else
|
||||
math_check_force_underflow_nonneg (retval);
|
||||
return retval;
|
||||
}
|
||||
|
||||
if (x == 0)
|
||||
ix = asuint64 (x);
|
||||
iy = asuint64 (y);
|
||||
topx = top12 (x);
|
||||
topy = top12 (y);
|
||||
if (__glibc_unlikely (topx - 0x001 >= 0x7ff - 0x001
|
||||
|| (topy & 0x7ff) - 0x3be >= 0x43e - 0x3be))
|
||||
{
|
||||
if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
|
||||
|| (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000) /* NaN */
|
||||
return y + y;
|
||||
if (fabs (y) > 1.0e20)
|
||||
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 */
|
||||
}
|
||||
|
||||
qx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign */
|
||||
qy = v.i[HIGH_HALF] & 0x7fffffff; /* no sign */
|
||||
|
||||
if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0)) /* NaN */
|
||||
return x + y;
|
||||
if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0)) /* NaN */
|
||||
return x == 1.0 && !issignaling (y) ? 1.0 : y + y;
|
||||
|
||||
/* if x<0 */
|
||||
if (u.i[HIGH_HALF] < 0)
|
||||
{
|
||||
k = checkint (y);
|
||||
if (k == 0)
|
||||
/* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0
|
||||
and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1. */
|
||||
/* Special cases: (x < 0x1p-126 or inf or nan) or
|
||||
(|y| < 0x1p-65 or |y| >= 0x1p63 or nan). */
|
||||
if (__glibc_unlikely (zeroinfnan (iy)))
|
||||
{
|
||||
if (qy == 0x7ff00000)
|
||||
if (2 * iy == 0)
|
||||
return issignaling_inline (x) ? x + y : 1.0;
|
||||
if (ix == asuint64 (1.0))
|
||||
return issignaling_inline (y) ? x + y : 1.0;
|
||||
if (2 * ix > 2 * asuint64 (INFINITY)
|
||||
|| 2 * iy > 2 * asuint64 (INFINITY))
|
||||
return x + y;
|
||||
if (2 * ix == 2 * asuint64 (1.0))
|
||||
return 1.0;
|
||||
if ((2 * ix < 2 * asuint64 (1.0)) == !(iy >> 63))
|
||||
return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf. */
|
||||
return y * y;
|
||||
}
|
||||
if (__glibc_unlikely (zeroinfnan (ix)))
|
||||
{
|
||||
double_t x2 = x * x;
|
||||
if (ix >> 63 && checkint (iy) == 1)
|
||||
{
|
||||
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;
|
||||
x2 = -x2;
|
||||
sign_bias = 1;
|
||||
}
|
||||
else if (qx == 0x7ff00000)
|
||||
return y < 0 ? 0.0 : INF.x;
|
||||
return (x - x) / (x - x); /* y not integer and x<0 */
|
||||
if (WANT_ERRNO && 2 * ix == 0 && iy >> 63)
|
||||
return __math_divzero (sign_bias);
|
||||
/* Without the barrier some versions of clang hoist the 1/x2 and
|
||||
thus division by zero exception can be signaled spuriously. */
|
||||
return iy >> 63 ? math_opt_barrier (1 / x2) : x2;
|
||||
}
|
||||
else if (qx == 0x7ff00000)
|
||||
/* Here x and y are non-zero finite. */
|
||||
if (ix >> 63)
|
||||
{
|
||||
if (k < 0)
|
||||
return y < 0 ? nZERO.x : nINF.x;
|
||||
else
|
||||
return y < 0 ? 0.0 : INF.x;
|
||||
/* Finite x < 0. */
|
||||
int yint = checkint (iy);
|
||||
if (yint == 0)
|
||||
return __math_invalid (x);
|
||||
if (yint == 1)
|
||||
sign_bias = SIGN_BIAS;
|
||||
ix &= 0x7fffffffffffffff;
|
||||
topx &= 0x7ff;
|
||||
}
|
||||
/* if y even or odd */
|
||||
if (k == 1)
|
||||
return __ieee754_pow (-x, y);
|
||||
else
|
||||
if ((topy & 0x7ff) - 0x3be >= 0x43e - 0x3be)
|
||||
{
|
||||
double retval;
|
||||
{
|
||||
SET_RESTORE_ROUND (FE_TONEAREST);
|
||||
retval = -__ieee754_pow (-x, y);
|
||||
}
|
||||
if (isinf (retval))
|
||||
retval = -huge * huge;
|
||||
else if (retval == 0)
|
||||
retval = -tiny * tiny;
|
||||
return retval;
|
||||
/* Note: sign_bias == 0 here because y is not odd. */
|
||||
if (ix == asuint64 (1.0))
|
||||
return 1.0;
|
||||
if ((topy & 0x7ff) < 0x3be)
|
||||
{
|
||||
/* |y| < 2^-65, x^y ~= 1 + y*log(x). */
|
||||
if (WANT_ROUNDING)
|
||||
return ix > asuint64 (1.0) ? 1.0 + y : 1.0 - y;
|
||||
else
|
||||
return 1.0;
|
||||
}
|
||||
return (ix > asuint64 (1.0)) == (topy < 0x800) ? __math_oflow (0)
|
||||
: __math_uflow (0);
|
||||
}
|
||||
if (topx == 0)
|
||||
{
|
||||
/* Normalize subnormal x so exponent becomes negative. */
|
||||
ix = asuint64 (x * 0x1p52);
|
||||
ix &= 0x7fffffffffffffff;
|
||||
ix -= 52ULL << 52;
|
||||
}
|
||||
}
|
||||
/* x>0 */
|
||||
|
||||
if (qx == 0x7ff00000) /* x= 2^-0x3ff */
|
||||
return y > 0 ? x : 0;
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
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 */
|
||||
double_t lo;
|
||||
double_t hi = log_inline (ix, &lo);
|
||||
double_t ehi, elo;
|
||||
#ifdef __FP_FAST_FMA
|
||||
ehi = y * hi;
|
||||
elo = y * lo + __builtin_fma (y, hi, -ehi);
|
||||
#else
|
||||
double_t yhi = asdouble (iy & -1ULL << 27);
|
||||
double_t ylo = y - yhi;
|
||||
double_t lhi = asdouble (asuint64 (hi) & -1ULL << 27);
|
||||
double_t llo = hi - lhi + lo;
|
||||
ehi = yhi * lhi;
|
||||
elo = ylo * lhi + y * llo; /* |elo| < |ehi| * 2^-25. */
|
||||
#endif
|
||||
return exp_inline (ehi, elo, sign_bias);
|
||||
}
|
||||
|
||||
#ifndef __ieee754_pow
|
||||
strong_alias (__ieee754_pow, __pow_finite)
|
||||
#endif
|
||||
|
||||
/* Compute log(x) (x is left argument). The result is the returned double + the
|
||||
parameter DELTA. */
|
||||
static double
|
||||
SECTION
|
||||
log1 (double x, double *delta)
|
||||
{
|
||||
unsigned int i, j;
|
||||
int m;
|
||||
double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0;
|
||||
mynumber u, v;
|
||||
#ifdef BIG_ENDI
|
||||
mynumber /**/ two52 = {{0x43300000, 0x00000000}}; /* 2**52 */
|
||||
#else
|
||||
# ifdef LITTLE_ENDI
|
||||
mynumber /**/ two52 = {{0x00000000, 0x43300000}}; /* 2**52 */
|
||||
# endif
|
||||
#endif
|
||||
|
||||
u.x = x;
|
||||
m = u.i[HIGH_HALF];
|
||||
if (m < 0x00100000) /* Handle denormal x. */
|
||||
{
|
||||
x = x * t52.x;
|
||||
add = -52.0;
|
||||
u.x = x;
|
||||
m = u.i[HIGH_HALF];
|
||||
}
|
||||
|
||||
if ((m & 0x000fffff) < 0x0006a09e)
|
||||
{
|
||||
u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3ff00000;
|
||||
two52.i[LOW_HALF] = (m >> 20);
|
||||
}
|
||||
else
|
||||
{
|
||||
u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3fe00000;
|
||||
two52.i[LOW_HALF] = (m >> 20) + 1;
|
||||
}
|
||||
|
||||
v.x = u.x + bigu.x;
|
||||
uu = v.x - bigu.x;
|
||||
i = (v.i[LOW_HALF] & 0x000003ff) << 2;
|
||||
if (two52.i[LOW_HALF] == 1023) /* Exponent of x is 0. */
|
||||
{
|
||||
if (i > 1192 && i < 1208) /* |x-1| < 1.5*2**-10 */
|
||||
{
|
||||
t = x - 1.0;
|
||||
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;
|
||||
/* Max relative error is 1.464844e-24, so accurate to 79.1 bits. */
|
||||
return res;
|
||||
} /* |x-1| < 1.5*2**-10 */
|
||||
else
|
||||
{
|
||||
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;
|
||||
/* Max relative error is 1.0e-24, so accurate to 79.7 bits. */
|
||||
return res;
|
||||
}
|
||||
}
|
||||
else /* Exponent of x != 0. */
|
||||
{
|
||||
eps = u.x - uu;
|
||||
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;
|
||||
/* Max relative error is 1.0e-21, so accurate to 69.7 bits. */
|
||||
return res;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* 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. */
|
||||
static int
|
||||
SECTION
|
||||
checkint (double x)
|
||||
{
|
||||
union
|
||||
{
|
||||
int4 i[2];
|
||||
double x;
|
||||
} u;
|
||||
int k;
|
||||
unsigned int m, n;
|
||||
u.x = x;
|
||||
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 */
|
||||
n = u.i[LOW_HALF];
|
||||
k = (m >> 20) - 1023; /* 1 <= k <= 52 */
|
||||
if (k == 52)
|
||||
return (n & 1) ? -1 : 1; /* odd or even */
|
||||
if (k > 20)
|
||||
{
|
||||
if (n << (k - 20) != 0)
|
||||
return 0; /* if not integer */
|
||||
return (n << (k - 21) != 0) ? -1 : 1;
|
||||
}
|
||||
if (n)
|
||||
return 0; /*if not integer */
|
||||
if (k == 20)
|
||||
return (m & 1) ? -1 : 1;
|
||||
if (m << (k + 12) != 0)
|
||||
return 0;
|
||||
return (m << (k + 11) != 0) ? -1 : 1;
|
||||
}
|
||||
|
||||
195
sysdeps/ieee754/dbl-64/e_pow_log_data.c
Normal file
195
sysdeps/ieee754/dbl-64/e_pow_log_data.c
Normal file
@ -0,0 +1,195 @@
|
||||
/* Data for the log part of pow.
|
||||
Copyright (C) 2018 Free Software Foundation, Inc.
|
||||
This file is part of the GNU C Library.
|
||||
|
||||
The GNU C Library 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.
|
||||
|
||||
The GNU C Library 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 the GNU C Library; if not, see
|
||||
<http://www.gnu.org/licenses/>. */
|
||||
|
||||
#include "math_config.h"
|
||||
|
||||
#define N (1 << POW_LOG_TABLE_BITS)
|
||||
|
||||
const struct pow_log_data __pow_log_data = {
|
||||
.ln2hi = 0x1.62e42fefa3800p-1,
|
||||
.ln2lo = 0x1.ef35793c76730p-45,
|
||||
.poly = {
|
||||
#if N == 128 && POW_LOG_POLY_ORDER == 8
|
||||
// relative error: 0x1.11922ap-70
|
||||
// in -0x1.6bp-8 0x1.6bp-8
|
||||
// Coefficients are scaled to match the scaling during evaluation.
|
||||
-0x1p-1,
|
||||
0x1.555555555556p-2 * -2,
|
||||
-0x1.0000000000006p-2 * -2,
|
||||
0x1.999999959554ep-3 * 4,
|
||||
-0x1.555555529a47ap-3 * 4,
|
||||
0x1.2495b9b4845e9p-3 * -8,
|
||||
-0x1.0002b8b263fc3p-3 * -8,
|
||||
#endif
|
||||
},
|
||||
/* Algorithm:
|
||||
|
||||
x = 2^k z
|
||||
log(x) = k ln2 + log(c) + log(z/c)
|
||||
log(z/c) = poly(z/c - 1)
|
||||
|
||||
where z is in [0x1.69555p-1; 0x1.69555p0] which is split into N subintervals
|
||||
and z falls into the ith one, then table entries are computed as
|
||||
|
||||
tab[i].invc = 1/c
|
||||
tab[i].logc = round(0x1p43*log(c))/0x1p43
|
||||
tab[i].logctail = (double)(log(c) - logc)
|
||||
|
||||
where c is chosen near the center of the subinterval such that 1/c has only a
|
||||
few precision bits so z/c - 1 is exactly representible as double:
|
||||
|
||||
1/c = center < 1 ? round(N/center)/N : round(2*N/center)/N/2
|
||||
|
||||
Note: |z/c - 1| < 1/N for the chosen c, |log(c) - logc - logctail| < 0x1p-97,
|
||||
the last few bits of logc are rounded away so k*ln2hi + logc has no rounding
|
||||
error and the interval for z is selected such that near x == 1, where log(x)
|
||||
is tiny, large cancellation error is avoided in logc + poly(z/c - 1). */
|
||||
.tab = {
|
||||
#if N == 128
|
||||
#define A(a,b,c) {a,0,b,c},
|
||||
A(0x1.6a00000000000p+0, -0x1.62c82f2b9c800p-2, 0x1.ab42428375680p-48)
|
||||
A(0x1.6800000000000p+0, -0x1.5d1bdbf580800p-2, -0x1.ca508d8e0f720p-46)
|
||||
A(0x1.6600000000000p+0, -0x1.5767717455800p-2, -0x1.362a4d5b6506dp-45)
|
||||
A(0x1.6400000000000p+0, -0x1.51aad872df800p-2, -0x1.684e49eb067d5p-49)
|
||||
A(0x1.6200000000000p+0, -0x1.4be5f95777800p-2, -0x1.41b6993293ee0p-47)
|
||||
A(0x1.6000000000000p+0, -0x1.4618bc21c6000p-2, 0x1.3d82f484c84ccp-46)
|
||||
A(0x1.5e00000000000p+0, -0x1.404308686a800p-2, 0x1.c42f3ed820b3ap-50)
|
||||
A(0x1.5c00000000000p+0, -0x1.3a64c55694800p-2, 0x1.0b1c686519460p-45)
|
||||
A(0x1.5a00000000000p+0, -0x1.347dd9a988000p-2, 0x1.5594dd4c58092p-45)
|
||||
A(0x1.5800000000000p+0, -0x1.2e8e2bae12000p-2, 0x1.67b1e99b72bd8p-45)
|
||||
A(0x1.5600000000000p+0, -0x1.2895a13de8800p-2, 0x1.5ca14b6cfb03fp-46)
|
||||
A(0x1.5600000000000p+0, -0x1.2895a13de8800p-2, 0x1.5ca14b6cfb03fp-46)
|
||||
A(0x1.5400000000000p+0, -0x1.22941fbcf7800p-2, -0x1.65a242853da76p-46)
|
||||
A(0x1.5200000000000p+0, -0x1.1c898c1699800p-2, -0x1.fafbc68e75404p-46)
|
||||
A(0x1.5000000000000p+0, -0x1.1675cababa800p-2, 0x1.f1fc63382a8f0p-46)
|
||||
A(0x1.4e00000000000p+0, -0x1.1058bf9ae4800p-2, -0x1.6a8c4fd055a66p-45)
|
||||
A(0x1.4c00000000000p+0, -0x1.0a324e2739000p-2, -0x1.c6bee7ef4030ep-47)
|
||||
A(0x1.4a00000000000p+0, -0x1.0402594b4d000p-2, -0x1.036b89ef42d7fp-48)
|
||||
A(0x1.4a00000000000p+0, -0x1.0402594b4d000p-2, -0x1.036b89ef42d7fp-48)
|
||||
A(0x1.4800000000000p+0, -0x1.fb9186d5e4000p-3, 0x1.d572aab993c87p-47)
|
||||
A(0x1.4600000000000p+0, -0x1.ef0adcbdc6000p-3, 0x1.b26b79c86af24p-45)
|
||||
A(0x1.4400000000000p+0, -0x1.e27076e2af000p-3, -0x1.72f4f543fff10p-46)
|
||||
A(0x1.4200000000000p+0, -0x1.d5c216b4fc000p-3, 0x1.1ba91bbca681bp-45)
|
||||
A(0x1.4000000000000p+0, -0x1.c8ff7c79aa000p-3, 0x1.7794f689f8434p-45)
|
||||
A(0x1.4000000000000p+0, -0x1.c8ff7c79aa000p-3, 0x1.7794f689f8434p-45)
|
||||
A(0x1.3e00000000000p+0, -0x1.bc286742d9000p-3, 0x1.94eb0318bb78fp-46)
|
||||
A(0x1.3c00000000000p+0, -0x1.af3c94e80c000p-3, 0x1.a4e633fcd9066p-52)
|
||||
A(0x1.3a00000000000p+0, -0x1.a23bc1fe2b000p-3, -0x1.58c64dc46c1eap-45)
|
||||
A(0x1.3a00000000000p+0, -0x1.a23bc1fe2b000p-3, -0x1.58c64dc46c1eap-45)
|
||||
A(0x1.3800000000000p+0, -0x1.9525a9cf45000p-3, -0x1.ad1d904c1d4e3p-45)
|
||||
A(0x1.3600000000000p+0, -0x1.87fa06520d000p-3, 0x1.bbdbf7fdbfa09p-45)
|
||||
A(0x1.3400000000000p+0, -0x1.7ab890210e000p-3, 0x1.bdb9072534a58p-45)
|
||||
A(0x1.3400000000000p+0, -0x1.7ab890210e000p-3, 0x1.bdb9072534a58p-45)
|
||||
A(0x1.3200000000000p+0, -0x1.6d60fe719d000p-3, -0x1.0e46aa3b2e266p-46)
|
||||
A(0x1.3000000000000p+0, -0x1.5ff3070a79000p-3, -0x1.e9e439f105039p-46)
|
||||
A(0x1.3000000000000p+0, -0x1.5ff3070a79000p-3, -0x1.e9e439f105039p-46)
|
||||
A(0x1.2e00000000000p+0, -0x1.526e5e3a1b000p-3, -0x1.0de8b90075b8fp-45)
|
||||
A(0x1.2c00000000000p+0, -0x1.44d2b6ccb8000p-3, 0x1.70cc16135783cp-46)
|
||||
A(0x1.2c00000000000p+0, -0x1.44d2b6ccb8000p-3, 0x1.70cc16135783cp-46)
|
||||
A(0x1.2a00000000000p+0, -0x1.371fc201e9000p-3, 0x1.178864d27543ap-48)
|
||||
A(0x1.2800000000000p+0, -0x1.29552f81ff000p-3, -0x1.48d301771c408p-45)
|
||||
A(0x1.2600000000000p+0, -0x1.1b72ad52f6000p-3, -0x1.e80a41811a396p-45)
|
||||
A(0x1.2600000000000p+0, -0x1.1b72ad52f6000p-3, -0x1.e80a41811a396p-45)
|
||||
A(0x1.2400000000000p+0, -0x1.0d77e7cd09000p-3, 0x1.a699688e85bf4p-47)
|
||||
A(0x1.2400000000000p+0, -0x1.0d77e7cd09000p-3, 0x1.a699688e85bf4p-47)
|
||||
A(0x1.2200000000000p+0, -0x1.fec9131dbe000p-4, -0x1.575545ca333f2p-45)
|
||||
A(0x1.2000000000000p+0, -0x1.e27076e2b0000p-4, 0x1.a342c2af0003cp-45)
|
||||
A(0x1.2000000000000p+0, -0x1.e27076e2b0000p-4, 0x1.a342c2af0003cp-45)
|
||||
A(0x1.1e00000000000p+0, -0x1.c5e548f5bc000p-4, -0x1.d0c57585fbe06p-46)
|
||||
A(0x1.1c00000000000p+0, -0x1.a926d3a4ae000p-4, 0x1.53935e85baac8p-45)
|
||||
A(0x1.1c00000000000p+0, -0x1.a926d3a4ae000p-4, 0x1.53935e85baac8p-45)
|
||||
A(0x1.1a00000000000p+0, -0x1.8c345d631a000p-4, 0x1.37c294d2f5668p-46)
|
||||
A(0x1.1a00000000000p+0, -0x1.8c345d631a000p-4, 0x1.37c294d2f5668p-46)
|
||||
A(0x1.1800000000000p+0, -0x1.6f0d28ae56000p-4, -0x1.69737c93373dap-45)
|
||||
A(0x1.1600000000000p+0, -0x1.51b073f062000p-4, 0x1.f025b61c65e57p-46)
|
||||
A(0x1.1600000000000p+0, -0x1.51b073f062000p-4, 0x1.f025b61c65e57p-46)
|
||||
A(0x1.1400000000000p+0, -0x1.341d7961be000p-4, 0x1.c5edaccf913dfp-45)
|
||||
A(0x1.1400000000000p+0, -0x1.341d7961be000p-4, 0x1.c5edaccf913dfp-45)
|
||||
A(0x1.1200000000000p+0, -0x1.16536eea38000p-4, 0x1.47c5e768fa309p-46)
|
||||
A(0x1.1000000000000p+0, -0x1.f0a30c0118000p-5, 0x1.d599e83368e91p-45)
|
||||
A(0x1.1000000000000p+0, -0x1.f0a30c0118000p-5, 0x1.d599e83368e91p-45)
|
||||
A(0x1.0e00000000000p+0, -0x1.b42dd71198000p-5, 0x1.c827ae5d6704cp-46)
|
||||
A(0x1.0e00000000000p+0, -0x1.b42dd71198000p-5, 0x1.c827ae5d6704cp-46)
|
||||
A(0x1.0c00000000000p+0, -0x1.77458f632c000p-5, -0x1.cfc4634f2a1eep-45)
|
||||
A(0x1.0c00000000000p+0, -0x1.77458f632c000p-5, -0x1.cfc4634f2a1eep-45)
|
||||
A(0x1.0a00000000000p+0, -0x1.39e87b9fec000p-5, 0x1.502b7f526feaap-48)
|
||||
A(0x1.0a00000000000p+0, -0x1.39e87b9fec000p-5, 0x1.502b7f526feaap-48)
|
||||
A(0x1.0800000000000p+0, -0x1.f829b0e780000p-6, -0x1.980267c7e09e4p-45)
|
||||
A(0x1.0800000000000p+0, -0x1.f829b0e780000p-6, -0x1.980267c7e09e4p-45)
|
||||
A(0x1.0600000000000p+0, -0x1.7b91b07d58000p-6, -0x1.88d5493faa639p-45)
|
||||
A(0x1.0400000000000p+0, -0x1.fc0a8b0fc0000p-7, -0x1.f1e7cf6d3a69cp-50)
|
||||
A(0x1.0400000000000p+0, -0x1.fc0a8b0fc0000p-7, -0x1.f1e7cf6d3a69cp-50)
|
||||
A(0x1.0200000000000p+0, -0x1.fe02a6b100000p-8, -0x1.9e23f0dda40e4p-46)
|
||||
A(0x1.0200000000000p+0, -0x1.fe02a6b100000p-8, -0x1.9e23f0dda40e4p-46)
|
||||
A(0x1.0000000000000p+0, 0x0.0000000000000p+0, 0x0.0000000000000p+0)
|
||||
A(0x1.0000000000000p+0, 0x0.0000000000000p+0, 0x0.0000000000000p+0)
|
||||
A(0x1.fc00000000000p-1, 0x1.0101575890000p-7, -0x1.0c76b999d2be8p-46)
|
||||
A(0x1.f800000000000p-1, 0x1.0205658938000p-6, -0x1.3dc5b06e2f7d2p-45)
|
||||
A(0x1.f400000000000p-1, 0x1.8492528c90000p-6, -0x1.aa0ba325a0c34p-45)
|
||||
A(0x1.f000000000000p-1, 0x1.0415d89e74000p-5, 0x1.111c05cf1d753p-47)
|
||||
A(0x1.ec00000000000p-1, 0x1.466aed42e0000p-5, -0x1.c167375bdfd28p-45)
|
||||
A(0x1.e800000000000p-1, 0x1.894aa149fc000p-5, -0x1.97995d05a267dp-46)
|
||||
A(0x1.e400000000000p-1, 0x1.ccb73cdddc000p-5, -0x1.a68f247d82807p-46)
|
||||
A(0x1.e200000000000p-1, 0x1.eea31c006c000p-5, -0x1.e113e4fc93b7bp-47)
|
||||
A(0x1.de00000000000p-1, 0x1.1973bd1466000p-4, -0x1.5325d560d9e9bp-45)
|
||||
A(0x1.da00000000000p-1, 0x1.3bdf5a7d1e000p-4, 0x1.cc85ea5db4ed7p-45)
|
||||
A(0x1.d600000000000p-1, 0x1.5e95a4d97a000p-4, -0x1.c69063c5d1d1ep-45)
|
||||
A(0x1.d400000000000p-1, 0x1.700d30aeac000p-4, 0x1.c1e8da99ded32p-49)
|
||||
A(0x1.d000000000000p-1, 0x1.9335e5d594000p-4, 0x1.3115c3abd47dap-45)
|
||||
A(0x1.cc00000000000p-1, 0x1.b6ac88dad6000p-4, -0x1.390802bf768e5p-46)
|
||||
A(0x1.ca00000000000p-1, 0x1.c885801bc4000p-4, 0x1.646d1c65aacd3p-45)
|
||||
A(0x1.c600000000000p-1, 0x1.ec739830a2000p-4, -0x1.dc068afe645e0p-45)
|
||||
A(0x1.c400000000000p-1, 0x1.fe89139dbe000p-4, -0x1.534d64fa10afdp-45)
|
||||
A(0x1.c000000000000p-1, 0x1.1178e8227e000p-3, 0x1.1ef78ce2d07f2p-45)
|
||||
A(0x1.be00000000000p-1, 0x1.1aa2b7e23f000p-3, 0x1.ca78e44389934p-45)
|
||||
A(0x1.ba00000000000p-1, 0x1.2d1610c868000p-3, 0x1.39d6ccb81b4a1p-47)
|
||||
A(0x1.b800000000000p-1, 0x1.365fcb0159000p-3, 0x1.62fa8234b7289p-51)
|
||||
A(0x1.b400000000000p-1, 0x1.4913d8333b000p-3, 0x1.5837954fdb678p-45)
|
||||
A(0x1.b200000000000p-1, 0x1.527e5e4a1b000p-3, 0x1.633e8e5697dc7p-45)
|
||||
A(0x1.ae00000000000p-1, 0x1.6574ebe8c1000p-3, 0x1.9cf8b2c3c2e78p-46)
|
||||
A(0x1.ac00000000000p-1, 0x1.6f0128b757000p-3, -0x1.5118de59c21e1p-45)
|
||||
A(0x1.aa00000000000p-1, 0x1.7898d85445000p-3, -0x1.c661070914305p-46)
|
||||
A(0x1.a600000000000p-1, 0x1.8beafeb390000p-3, -0x1.73d54aae92cd1p-47)
|
||||
A(0x1.a400000000000p-1, 0x1.95a5adcf70000p-3, 0x1.7f22858a0ff6fp-47)
|
||||
A(0x1.a000000000000p-1, 0x1.a93ed3c8ae000p-3, -0x1.8724350562169p-45)
|
||||
A(0x1.9e00000000000p-1, 0x1.b31d8575bd000p-3, -0x1.c358d4eace1aap-47)
|
||||
A(0x1.9c00000000000p-1, 0x1.bd087383be000p-3, -0x1.d4bc4595412b6p-45)
|
||||
A(0x1.9a00000000000p-1, 0x1.c6ffbc6f01000p-3, -0x1.1ec72c5962bd2p-48)
|
||||
A(0x1.9600000000000p-1, 0x1.db13db0d49000p-3, -0x1.aff2af715b035p-45)
|
||||
A(0x1.9400000000000p-1, 0x1.e530effe71000p-3, 0x1.212276041f430p-51)
|
||||
A(0x1.9200000000000p-1, 0x1.ef5ade4dd0000p-3, -0x1.a211565bb8e11p-51)
|
||||
A(0x1.9000000000000p-1, 0x1.f991c6cb3b000p-3, 0x1.bcbecca0cdf30p-46)
|
||||
A(0x1.8c00000000000p-1, 0x1.07138604d5800p-2, 0x1.89cdb16ed4e91p-48)
|
||||
A(0x1.8a00000000000p-1, 0x1.0c42d67616000p-2, 0x1.7188b163ceae9p-45)
|
||||
A(0x1.8800000000000p-1, 0x1.1178e8227e800p-2, -0x1.c210e63a5f01cp-45)
|
||||
A(0x1.8600000000000p-1, 0x1.16b5ccbacf800p-2, 0x1.b9acdf7a51681p-45)
|
||||
A(0x1.8400000000000p-1, 0x1.1bf99635a6800p-2, 0x1.ca6ed5147bdb7p-45)
|
||||
A(0x1.8200000000000p-1, 0x1.214456d0eb800p-2, 0x1.a87deba46baeap-47)
|
||||
A(0x1.7e00000000000p-1, 0x1.2bef07cdc9000p-2, 0x1.a9cfa4a5004f4p-45)
|
||||
A(0x1.7c00000000000p-1, 0x1.314f1e1d36000p-2, -0x1.8e27ad3213cb8p-45)
|
||||
A(0x1.7a00000000000p-1, 0x1.36b6776be1000p-2, 0x1.16ecdb0f177c8p-46)
|
||||
A(0x1.7800000000000p-1, 0x1.3c25277333000p-2, 0x1.83b54b606bd5cp-46)
|
||||
A(0x1.7600000000000p-1, 0x1.419b423d5e800p-2, 0x1.8e436ec90e09dp-47)
|
||||
A(0x1.7400000000000p-1, 0x1.4718dc271c800p-2, -0x1.f27ce0967d675p-45)
|
||||
A(0x1.7200000000000p-1, 0x1.4c9e09e173000p-2, -0x1.e20891b0ad8a4p-45)
|
||||
A(0x1.7000000000000p-1, 0x1.522ae0738a000p-2, 0x1.ebe708164c759p-45)
|
||||
A(0x1.6e00000000000p-1, 0x1.57bf753c8d000p-2, 0x1.fadedee5d40efp-46)
|
||||
A(0x1.6c00000000000p-1, 0x1.5d5bddf596000p-2, -0x1.a0b2a08a465dcp-47)
|
||||
#endif
|
||||
},
|
||||
};
|
||||
@ -21,6 +21,7 @@
|
||||
|
||||
#include <math.h>
|
||||
#include <math_private.h>
|
||||
#include <nan-high-order-bit.h>
|
||||
#include <stdint.h>
|
||||
|
||||
#ifndef WANT_ROUNDING
|
||||
@ -78,6 +79,15 @@ asdouble (uint64_t i)
|
||||
return u.f;
|
||||
}
|
||||
|
||||
static inline int
|
||||
issignaling_inline (double x)
|
||||
{
|
||||
uint64_t ix = asuint64 (x);
|
||||
if (HIGH_ORDER_BIT_IS_SET_FOR_SNAN)
|
||||
return (ix & 0x7ff8000000000000) == 0x7ff8000000000000;
|
||||
return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL;
|
||||
}
|
||||
|
||||
#define NOINLINE __attribute__ ((noinline))
|
||||
|
||||
/* Error handling tail calls for special cases, with a sign argument.
|
||||
@ -165,4 +175,16 @@ extern const struct log2_data
|
||||
#endif
|
||||
} __log2_data attribute_hidden;
|
||||
|
||||
#define POW_LOG_TABLE_BITS 7
|
||||
#define POW_LOG_POLY_ORDER 8
|
||||
extern const struct pow_log_data
|
||||
{
|
||||
double ln2hi;
|
||||
double ln2lo;
|
||||
double poly[POW_LOG_POLY_ORDER - 1]; /* First coefficient is 1. */
|
||||
/* Note: the pad field is unused, but allows slightly faster indexing. */
|
||||
/* See e_pow_log_data.c for details. */
|
||||
struct {double invc, pad, logc, logctail;} tab[1 << POW_LOG_TABLE_BITS];
|
||||
} __pow_log_data attribute_hidden;
|
||||
|
||||
#endif
|
||||
|
||||
@ -1,76 +0,0 @@
|
||||
/*
|
||||
* IBM Accurate Mathematical Library
|
||||
* Written by International Business Machines Corp.
|
||||
* 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
|
||||
* 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:upow.h */
|
||||
/* */
|
||||
/* common data and variables prototype and definition */
|
||||
/******************************************************************/
|
||||
|
||||
#ifndef UPOW_H
|
||||
#define UPOW_H
|
||||
|
||||
#include "mydefs.h"
|
||||
|
||||
#ifdef BIG_ENDI
|
||||
const static mynumber
|
||||
/**/ nZERO = {{0x80000000, 0}}, /* -0.0 */
|
||||
/**/ INF = {{0x7ff00000, 0x00000000}}, /* INF */
|
||||
/**/ nINF = {{0xfff00000, 0x00000000}}, /* -INF */
|
||||
/**/ ln2a = {{0x3fe62e42, 0xfefa3800}}, /* ln(2) 43 bits */
|
||||
/**/ ln2b = {{0x3d2ef357, 0x93c76730}}, /* ln(2)-ln2a */
|
||||
/**/ bigu = {{0x4297ffff, 0xfffffd2c}}, /* 1.5*2**42 -724*2**-10 */
|
||||
/**/ bigv = {{0x4207ffff, 0xfff8016a}}, /* 1.5*2**33-1+362*2**-19 */
|
||||
/**/ t52 = {{0x43300000, 0x00000000}}, /* 2**52 */
|
||||
/**/ two52e = {{0x43300000, 0x000003ff}}; /* 2**52' */
|
||||
|
||||
#else
|
||||
#ifdef LITTLE_ENDI
|
||||
const static mynumber
|
||||
/**/ nZERO = {{0, 0x80000000}}, /* -0.0 */
|
||||
/**/ INF = {{0x00000000, 0x7ff00000}}, /* INF */
|
||||
/**/ nINF = {{0x00000000, 0xfff00000}}, /* -INF */
|
||||
/**/ ln2a = {{0xfefa3800, 0x3fe62e42}}, /* ln(2) 43 bits */
|
||||
/**/ ln2b = {{0x93c76730, 0x3d2ef357}}, /* ln(2)-ln2a */
|
||||
/**/ bigu = {{0xfffffd2c, 0x4297ffff}}, /* 1.5*2**42 -724*2**-10 */
|
||||
/**/ bigv = {{0xfff8016a, 0x4207ffff}}, /* 1.5*2**33-1+362*2**-19 */
|
||||
/**/ t52 = {{0x00000000, 0x43300000}}, /* 2**52 */
|
||||
/**/ two52e = {{0x000003ff, 0x43300000}}; /* 2**52' */
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
const static double p2=-0.5, p3 = 3.3333333333333333333e-1, p4 = -0.25,
|
||||
q2 = -0.5, q3 = 3.3333333333331404e-01, q4 = -2.4999999999996436e-01,
|
||||
q5 = 2.0000010500004459e-01, q6 = -1.6666678916688004e-01,
|
||||
r3 = 3.33333333333333333372884096563030E-01,
|
||||
r4 = -2.50000000000000000213574153875908E-01,
|
||||
r5 = 1.99999999999683593814072199830603E-01,
|
||||
r6 = -1.66666666666065494878165510225378E-01,
|
||||
r7 = 1.42857517857114380606360005067609E-01,
|
||||
r8 = -1.25000449999974370683775964001702E-01,
|
||||
s3 = 0.333251953125000000e0,
|
||||
ss3 = 8.138020833333333333e-05,
|
||||
s4 = -2.500000000000000000e-01,
|
||||
s5 = 1.999999999999960937e-01,
|
||||
s6 = -1.666666666666592447e-01,
|
||||
s7 = 1.428571845238194705e-01,
|
||||
s8 = -1.250000500000149097e-01;
|
||||
#endif
|
||||
File diff suppressed because it is too large
Load Diff
1
sysdeps/m68k/m680x0/fpu/e_pow_log_data.c
Normal file
1
sysdeps/m68k/m680x0/fpu/e_pow_log_data.c
Normal file
@ -0,0 +1 @@
|
||||
/* Not needed. */
|
||||
@ -20,7 +20,7 @@ CFLAGS-e_asin-fma.c = -mfma -mavx2
|
||||
CFLAGS-e_atan2-fma.c = -mfma -mavx2
|
||||
CFLAGS-e_exp-fma.c = -mfma -mavx2
|
||||
CFLAGS-e_log-fma.c = -mfma -mavx2
|
||||
CFLAGS-e_pow-fma.c = -mfma -mavx2 $(config-cflags-nofma)
|
||||
CFLAGS-e_pow-fma.c = -mfma -mavx2
|
||||
CFLAGS-mpa-fma.c = -mfma -mavx2
|
||||
CFLAGS-mpatan-fma.c = -mfma -mavx2
|
||||
CFLAGS-mpatan2-fma.c = -mfma -mavx2
|
||||
@ -57,7 +57,7 @@ CFLAGS-e_asin-fma4.c = -mfma4
|
||||
CFLAGS-e_atan2-fma4.c = -mfma4
|
||||
CFLAGS-e_exp-fma4.c = -mfma4
|
||||
CFLAGS-e_log-fma4.c = -mfma4
|
||||
CFLAGS-e_pow-fma4.c = -mfma4 $(config-cflags-nofma)
|
||||
CFLAGS-e_pow-fma4.c = -mfma4
|
||||
CFLAGS-mpa-fma4.c = -mfma4
|
||||
CFLAGS-mpatan-fma4.c = -mfma4
|
||||
CFLAGS-mpatan2-fma4.c = -mfma4
|
||||
|
||||
Loading…
Reference in New Issue
Block a user