mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-26 15:00:06 +00:00
424c4f60ed
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.
241 lines
5.8 KiB
C
241 lines
5.8 KiB
C
/*
|
|
* ====================================================
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
*
|
|
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
* Permission to use, copy, modify, and distribute this
|
|
* software is freely granted, provided that this notice
|
|
* is preserved.
|
|
* ====================================================
|
|
*/
|
|
|
|
/*
|
|
* from: @(#)fdlibm.h 5.1 93/09/24
|
|
*/
|
|
|
|
#ifndef _MATH_PRIVATE_H_
|
|
#define _MATH_PRIVATE_H_
|
|
|
|
#include <endian.h>
|
|
#include <stdbool.h>
|
|
#include <stdint.h>
|
|
#include <sys/types.h>
|
|
|
|
/* Gather machine dependent _Floatn support. */
|
|
#include <bits/floatn.h>
|
|
|
|
/* The original fdlibm code used statements like:
|
|
n0 = ((*(int*)&one)>>29)^1; * index of high word *
|
|
ix0 = *(n0+(int*)&x); * high word of x *
|
|
ix1 = *((1-n0)+(int*)&x); * low word of x *
|
|
to dig two 32 bit words out of the 64 bit IEEE floating point
|
|
value. That is non-ANSI, and, moreover, the gcc instruction
|
|
scheduler gets it wrong. We instead use the following macros.
|
|
Unlike the original code, we determine the endianness at compile
|
|
time, not at run time; I don't see much benefit to selecting
|
|
endianness at run time. */
|
|
|
|
/* A union which permits us to convert between a double and two 32 bit
|
|
ints. */
|
|
|
|
#if __FLOAT_WORD_ORDER == __BIG_ENDIAN
|
|
|
|
typedef union
|
|
{
|
|
double value;
|
|
struct
|
|
{
|
|
uint32_t msw;
|
|
uint32_t lsw;
|
|
} parts;
|
|
uint64_t word;
|
|
} ieee_double_shape_type;
|
|
|
|
#endif
|
|
|
|
#if __FLOAT_WORD_ORDER == __LITTLE_ENDIAN
|
|
|
|
typedef union
|
|
{
|
|
double value;
|
|
struct
|
|
{
|
|
uint32_t lsw;
|
|
uint32_t msw;
|
|
} parts;
|
|
uint64_t word;
|
|
} ieee_double_shape_type;
|
|
|
|
#endif
|
|
|
|
/* Get two 32 bit ints from a double. */
|
|
|
|
#define EXTRACT_WORDS(ix0,ix1,d) \
|
|
do { \
|
|
ieee_double_shape_type ew_u; \
|
|
ew_u.value = (d); \
|
|
(ix0) = ew_u.parts.msw; \
|
|
(ix1) = ew_u.parts.lsw; \
|
|
} while (0)
|
|
|
|
/* Get the more significant 32 bit int from a double. */
|
|
|
|
#ifndef GET_HIGH_WORD
|
|
# define GET_HIGH_WORD(i,d) \
|
|
do { \
|
|
ieee_double_shape_type gh_u; \
|
|
gh_u.value = (d); \
|
|
(i) = gh_u.parts.msw; \
|
|
} while (0)
|
|
#endif
|
|
|
|
/* Get the less significant 32 bit int from a double. */
|
|
|
|
#ifndef GET_LOW_WORD
|
|
# define GET_LOW_WORD(i,d) \
|
|
do { \
|
|
ieee_double_shape_type gl_u; \
|
|
gl_u.value = (d); \
|
|
(i) = gl_u.parts.lsw; \
|
|
} while (0)
|
|
#endif
|
|
|
|
/* Get all in one, efficient on 64-bit machines. */
|
|
#ifndef EXTRACT_WORDS64
|
|
# define EXTRACT_WORDS64(i,d) \
|
|
do { \
|
|
ieee_double_shape_type gh_u; \
|
|
gh_u.value = (d); \
|
|
(i) = gh_u.word; \
|
|
} while (0)
|
|
#endif
|
|
|
|
/* Set a double from two 32 bit ints. */
|
|
#ifndef INSERT_WORDS
|
|
# define INSERT_WORDS(d,ix0,ix1) \
|
|
do { \
|
|
ieee_double_shape_type iw_u; \
|
|
iw_u.parts.msw = (ix0); \
|
|
iw_u.parts.lsw = (ix1); \
|
|
(d) = iw_u.value; \
|
|
} while (0)
|
|
#endif
|
|
|
|
/* Get all in one, efficient on 64-bit machines. */
|
|
#ifndef INSERT_WORDS64
|
|
# define INSERT_WORDS64(d,i) \
|
|
do { \
|
|
ieee_double_shape_type iw_u; \
|
|
iw_u.word = (i); \
|
|
(d) = iw_u.value; \
|
|
} while (0)
|
|
#endif
|
|
|
|
/* Set the more significant 32 bits of a double from an int. */
|
|
#ifndef SET_HIGH_WORD
|
|
#define SET_HIGH_WORD(d,v) \
|
|
do { \
|
|
ieee_double_shape_type sh_u; \
|
|
sh_u.value = (d); \
|
|
sh_u.parts.msw = (v); \
|
|
(d) = sh_u.value; \
|
|
} while (0)
|
|
#endif
|
|
|
|
/* Set the less significant 32 bits of a double from an int. */
|
|
#ifndef SET_LOW_WORD
|
|
# define SET_LOW_WORD(d,v) \
|
|
do { \
|
|
ieee_double_shape_type sl_u; \
|
|
sl_u.value = (d); \
|
|
sl_u.parts.lsw = (v); \
|
|
(d) = sl_u.value; \
|
|
} while (0)
|
|
#endif
|
|
|
|
/* A union which permits us to convert between a float and a 32 bit
|
|
int. */
|
|
|
|
typedef union
|
|
{
|
|
float value;
|
|
uint32_t word;
|
|
} ieee_float_shape_type;
|
|
|
|
/* Get a 32 bit int from a float. */
|
|
#ifndef GET_FLOAT_WORD
|
|
# define GET_FLOAT_WORD(i,d) \
|
|
do { \
|
|
ieee_float_shape_type gf_u; \
|
|
gf_u.value = (d); \
|
|
(i) = gf_u.word; \
|
|
} while (0)
|
|
#endif
|
|
|
|
/* Set a float from a 32 bit int. */
|
|
#ifndef SET_FLOAT_WORD
|
|
# define SET_FLOAT_WORD(d,i) \
|
|
do { \
|
|
ieee_float_shape_type sf_u; \
|
|
sf_u.word = (i); \
|
|
(d) = sf_u.value; \
|
|
} while (0)
|
|
#endif
|
|
|
|
/* We need to guarantee an expansion of name when building
|
|
ldbl-128 files as another type (e.g _Float128). */
|
|
#define mathx_hidden_def(name) hidden_def(name)
|
|
|
|
/* Get long double macros from a separate header. */
|
|
#include <math_ldbl.h>
|
|
|
|
/* Include function declarations for each floating-point. */
|
|
#define _Mdouble_ double
|
|
#define _MSUF_
|
|
#include <math_private_calls.h>
|
|
#undef _MSUF_
|
|
#undef _Mdouble_
|
|
|
|
#define _Mdouble_ float
|
|
#define _MSUF_ f
|
|
#define __MATH_DECLARING_FLOAT
|
|
#include <math_private_calls.h>
|
|
#undef __MATH_DECLARING_FLOAT
|
|
#undef _MSUF_
|
|
#undef _Mdouble_
|
|
|
|
#define _Mdouble_ long double
|
|
#define _MSUF_ l
|
|
#define __MATH_DECLARING_LONG_DOUBLE
|
|
#include <math_private_calls.h>
|
|
#undef __MATH_DECLARING_LONG_DOUBLE
|
|
#undef _MSUF_
|
|
#undef _Mdouble_
|
|
|
|
#if __HAVE_DISTINCT_FLOAT128
|
|
# define _Mdouble_ _Float128
|
|
# define _MSUF_ f128
|
|
# define __MATH_DECLARING_FLOATN
|
|
# include <math_private_calls.h>
|
|
# undef __MATH_DECLARING_FLOATN
|
|
# undef _MSUF_
|
|
# undef _Mdouble_
|
|
#endif
|
|
|
|
|
|
|
|
/* Prototypes for functions of the IBM Accurate Mathematical Library. */
|
|
extern double __sin (double __x);
|
|
extern double __cos (double __x);
|
|
extern int __branred (double __x, double *__a, double *__aa);
|
|
extern void __doasin (double __x, double __dx, double __v[]);
|
|
extern void __dubsin (double __x, double __dx, double __v[]);
|
|
extern void __dubcos (double __x, double __dx, double __v[]);
|
|
extern double __sin32 (double __x, double __res, double __res1);
|
|
extern double __cos32 (double __x, double __res, double __res1);
|
|
extern double __mpsin (double __x, double __dx, bool __range_reduce);
|
|
extern double __mpcos (double __x, double __dx, bool __range_reduce);
|
|
extern void __docos (double __x, double __dx, double __v[]);
|
|
|
|
#endif /* _MATH_PRIVATE_H_ */
|