Skip to content

T
Third Party Musl

OpenHarmony / Third Party Musl
1 年多 前同步成功

通知 37
Star 125
Fork 0
T
Third Party Musl
提交 b69f695a 编写于 作者: R Rich Felker
浏览文件
操作

first commit of the new libm! 

thanks to the hard work of Szabolcs Nagy (nsz), identifying the best
(from correctness and license standpoint) implementations from freebsd
and openbsd and cleaning them up! musl should now fully support c99
float and long double math functions, and has near-complete complex
math support. tgmath should also work (fully on gcc-compatible
compilers, and mostly on any c99 compiler).

based largely on commit 0376d44a890fea261506f1fc63833e7a686dca19 from
nsz's libm git repo, with some additions (dummy versions of a few
missing long double complex functions, etc.) by me.

various cleanups still need to be made, including re-adding (if
they're correct) some asm functions that were dropped.
上级 d46cf2e1
展开全部 隐藏空白更改
内联 并排
Showing with 20547 addition and 5917 deletion
include/complex.h 0 → 100644
浏览文件 @ b69f695a
#ifndef _COMPLEX_H
#define _COMPLEX_H
#ifdef __cplusplus
extern "C" {
#endif
#define complex _Complex
#define _Complex_I 1.0fi
#define I _Complex_I
double complex cacos(double complex);
float complex cacosf(float complex);
long double complex cacosl(long double complex);
double complex casin(double complex);
float complex casinf(float complex);
long double complex casinl(long double complex);
double complex catan(double complex);
float complex catanf(float complex);
long double complex catanl(long double complex);
double complex ccos(double complex);
float complex ccosf(float complex);
long double complex ccosl(long double complex);
double complex csin(double complex);
float complex csinf(float complex);
long double complex csinl(long double complex);
double complex ctan(double complex);
float complex ctanf(float complex);
long double complex ctanl(long double complex);
double complex cacosh(double complex);
float complex cacoshf(float complex);
long double complex cacoshl(long double complex);
double complex casinh(double complex);
float complex casinhf(float complex);
long double complex casinhl(long double complex);
double complex catanh(double complex);
float complex catanhf(float complex);
long double complex catanhl(long double complex);
double complex ccosh(double complex);
float complex ccoshf(float complex);
long double complex ccoshl(long double complex);
double complex csinh(double complex);
float complex csinhf(float complex);
long double complex csinhl(long double complex);
double complex ctanh(double complex);
float complex ctanhf(float complex);
long double complex ctanhl(long double complex);
double complex cexp(double complex);
float complex cexpf(float complex);
long double complex cexpl(long double complex);
double complex clog(double complex);
float complex clogf(float complex);
long double complex clogl(long double complex);
double cabs(double complex);
float cabsf(float complex);
long double cabsl(long double complex);
double complex cpow(double complex, double complex);
float complex cpowf(float complex, float complex);
long double complex cpowl(long double complex, long double complex);
double complex csqrt(double complex);
float complex csqrtf(float complex);
long double complex csqrtl(long double complex);
double carg(double complex);
float cargf(float complex);
long double cargl(long double complex);
double cimag(double complex);
float cimagf(float complex);
long double cimagl(long double complex);
double complex conj(double complex);
float complex conjf(float complex);
long double complex conjl(long double complex);
double complex cproj(double complex);
float complex cprojf(float complex);
long double complex cprojl(long double complex);
double creal(double complex);
float crealf(float complex);
long double creall(long double complex);
#ifdef __cplusplus
}
#endif
#endif
include/math.h
浏览文件 @ b69f695a
......@@ -37,27 +37,53 @@ extern "C" {
#define FP_SUBNORMAL 3
#define FP_NORMAL 4
int __fpclassifyf(float);
int __fpclassify(double);
int __fpclassifyf(float);
int __fpclassifyl(long double);
#define __FLOAT_BITS(f) (((union { float __f; __uint32_t __i; }){ (f) }).__i)
#define __DOUBLE_BITS(f) (((union { double __f; __uint64_t __i; }){ (f) }).__i)
#define fpclassify(x) ( \
sizeof(x) == sizeof(float) ? __fpclassifyf(x) : \
sizeof(x) == sizeof(double) ? __fpclassify(x) : \
__fpclassifyl(x) )
#define isinf(x) (fpclassify(x) == FP_INFINITE)
#define isnan(x) (fpclassify(x) == FP_NAN)
#define isnormal(x) (fpclassify(x) == FP_NORMAL)
#define isfinite(x) (fpclassify(x) > FP_INFINITE)
#define isinf(x) ( \
sizeof(x) == sizeof(float) ? (__FLOAT_BITS(x) & 0x7fffffff) == 0x7f800000 : \
sizeof(x) == sizeof(double) ? (__DOUBLE_BITS(x) & (__uint64_t)-1>>1) == (__uint64_t)0x7ff<<52 : \
__fpclassifyl(x) == FP_INFINITE)
#define isnan(x) ( \
sizeof(x) == sizeof(float) ? (__FLOAT_BITS(x) & 0x7fffffff) > 0x7f800000 : \
sizeof(x) == sizeof(double) ? (__DOUBLE_BITS(x) & (__uint64_t)-1>>1) > (__uint64_t)0x7ff<<52 : \
__fpclassifyl(x) == FP_NAN)
#define isnormal(x) ( \
sizeof(x) == sizeof(float) ? ((__FLOAT_BITS(x)+0x00800000) & 0x7fffffff) >= 0x01000000 : \
sizeof(x) == sizeof(double) ? ((__DOUBLE_BITS(x)+((__uint64_t)1<<52)) & (__uint64_t)-1>>1) >= (__uint64_t)1<<53 : \
__fpclassifyl(x) == FP_NORMAL)
#define isunordered(x,y) (isnan((x)) ? ((y),1) : isnan((y)))
#define isfinite(x) ( \
sizeof(x) == sizeof(float) ? (__FLOAT_BITS(x) & 0x7fffffff) < 0x7f800000 : \
sizeof(x) == sizeof(double) ? (__DOUBLE_BITS(x) & (__uint64_t)-1>>1) < (__uint64_t)0x7ff<<52 : \
__fpclassifyl(x) > FP_INFINITE)
int __signbit(double);
int __signbitf(float);
int __signbitl(long double);
#define signbit(x) ( \
sizeof(x) == sizeof(float) ? !!(__FLOAT_BITS(x) & 0x80000000) : \
sizeof(x) == sizeof(double) ? !!(__DOUBLE_BITS(x) & (__uint64_t)1<<63) : \
__signbitl(x) )
#define isunordered(x,y) (isnan((x)) ? ((void)(y),1) : isnan((y)))
static
#if __STDC_VERSION__ >= 199901L
inline
#endif
int __isrel(long double __x, long double __y, int __rel)
static int __isrel(long double __x, long double __y, int __rel)
{
if (isunordered(__x, __y)) return 0;
if (__rel==-2) return __x < __y;
......@@ -316,17 +342,46 @@ long double truncl(long double);
#define M_2_SQRTPI 1.12837916709551257390 /* 2/sqrt(pi) */
#define M_SQRT2 1.41421356237309504880 /* sqrt(2) */
#define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */
extern int signgam;
double gamma(double);
float gammaf(float);
long double gammal(long double);
double lgamma_r(double, int*);
float lgammaf_r(float, int*);
long double lgammal_r(long double, int*);
double j0(double);
float j0f(float);
long double j0l(long double);
double j1(double);
float j1f(float);
long double j1l(long double);
double jn(int, double);
float jnf(int, float);
long double jnl(int, long double);
double y0(double);
float y0f(float);
long double y0l(long double);
double y1(double);
float y1f(float);
long double y1l(long double);
double yn(int, double);
extern int signgam;
float ynf(int, float);
long double ynl(int, long double);
#endif
#ifdef _GNU_SOURCE
double scalb(double, double);
float scalbf(float, float);
long double scalbl(long double, long double);
#endif
#ifdef __cplusplus
......
include/tgmath.h 0 → 100644
浏览文件 @ b69f695a
#ifndef _TGMATH_H
#define _TGMATH_H
/*
the return types are only correct with gcc (__GNUC__)
otherwise they are long double or long double complex
the long double version of a function is never chosen when
sizeof(double) == sizeof(long double)
(but the return type is set correctly with gcc)
*/
#include <math.h>
#include <complex.h>
#define __IS_FP(x) !!((1?1:(x))/2)
#define __IS_CX(x) (__IS_FP(x) && sizeof(x) == sizeof((x)+I))
#define __IS_REAL(x) (__IS_FP(x) && 2*sizeof(x) == sizeof((x)+I))
#define __FLT(x) (__IS_REAL(x) && sizeof(x) == sizeof(float))
#define __LDBL(x) (__IS_REAL(x) && sizeof(x) == sizeof(long double) && sizeof(long double) != sizeof(double))
#define __FLTCX(x) (__IS_CX(x) && sizeof(x) == sizeof(float complex))
#define __DBLCX(x) (__IS_CX(x) && sizeof(x) == sizeof(double complex))
#define __LDBLCX(x) (__IS_CX(x) && sizeof(x) == sizeof(long double complex) && sizeof(long double) != sizeof(double))
/* return type */
#ifdef __GNUC__
/* cast to double when x is integral, otherwise use typeof(x) */
#define __RETCAST(x) (__typeof__(*( \
0 ? (__typeof__(0 ? (double *)0 : (void *)__IS_FP(x)))0 : \
(__typeof__(0 ? (__typeof__(x) *)0 : (void *)!__IS_FP(x)))0 )))
/* 2 args case, consider complex types (for cpow) */
#define __RETCAST_2(x, y) (__typeof__(*( \
0 ? (__typeof__(0 ? (double *)0 : \
(void *)!((!__IS_FP(x) || !__IS_FP(y)) && __FLT((x)+(y)+1.0f))))0 : \
0 ? (__typeof__(0 ? (double complex *)0 : \
(void *)!((!__IS_FP(x) || !__IS_FP(y)) && __FLTCX((x)+(y)))))0 : \
(__typeof__(0 ? (__typeof__((x)+(y)) *)0 : \
(void *)((!__IS_FP(x) || !__IS_FP(y)) && (__FLT((x)+(y)+1.0f) || __FLTCX((x)+(y))))))0 )))
/* 3 args case, don't consider complex types (fma only) */
#define __RETCAST_3(x, y, z) (__typeof__(*( \
0 ? (__typeof__(0 ? (double *)0 : \
(void *)!((!__IS_FP(x) || !__IS_FP(y) || !__IS_FP(z)) && __FLT((x)+(y)+(z)+1.0f))))0 : \
(__typeof__(0 ? (__typeof__((x)+(y)) *)0 : \
(void *)((!__IS_FP(x) || !__IS_FP(y) || !__IS_FP(z)) && __FLT((x)+(y)+(z)+1.0f))))0 )))
/* drop complex from the type of x */
#define __TO_REAL(x) *( \
0 ? (__typeof__(0 ? (double *)0 : (void *)!__DBLCX(x)))0 : \
0 ? (__typeof__(0 ? (float *)0 : (void *)!__FLTCX(x)))0 : \
0 ? (__typeof__(0 ? (long double *)0 : (void *)!__LDBLCX(x)))0 : \
(__typeof__(0 ? (__typeof__(x) *)0 : (void *)__IS_CX(x)))0 )
#else
#define __RETCAST(x)
#define __RETCAST_2(x, y)
#define __RETCAST_3(x, y, z)
#endif
/* function selection */
#define __tg_real(fun, x) (__RETCAST(x)( \
__FLT(x) ? fun ## f (x) : \
__LDBL(x) ? fun ## l (x) : \
fun(x) ))
#define __tg_real_2_1(fun, x, y) (__RETCAST(x)( \
__FLT(x) ? fun ## f (x, y) : \
__LDBL(x) ? fun ## l (x, y) : \
fun(x, y) ))
#define __tg_real_2(fun, x, y) (__RETCAST_2(x, y)( \
__FLT(x) && __FLT(y) ? fun ## f (x, y) : \
__LDBL((x)+(y)) ? fun ## l (x, y) : \
fun(x, y) ))
#define __tg_complex(fun, x) (__RETCAST((x)+I)( \
__FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
__LDBLCX((x)+I) ? fun ## l (x) : \
fun(x) ))
#define __tg_complex_retreal(fun, x) (__RETCAST(__TO_REAL(x))( \
__FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
__LDBLCX((x)+I) ? fun ## l (x) : \
fun(x) ))
#define __tg_real_complex(fun, x) (__RETCAST(x)( \
__FLTCX(x) ? c ## fun ## f (x) : \
__DBLCX(x) ? c ## fun (x) : \
__LDBLCX(x) ? c ## fun ## l (x) : \
__FLT(x) ? fun ## f (x) : \
__LDBL(x) ? fun ## l (x) : \
fun(x) ))
/* special cases */
#define __tg_real_remquo(x, y, z) (__RETCAST_2(x, y)( \
__FLT(x) && __FLT(y) ? remquof(x, y, z) : \
__LDBL((x)+(y)) ? remquol(x, y, z) : \
remquo(x, y, z) ))
#define __tg_real_fma(x, y, z) (__RETCAST_3(x, y, z)( \
__FLT(x) && __FLT(y) && __FLT(z) ? fmaf(x, y, z) : \
__LDBL((x)+(y)+(z)) ? fmal(x, y, z) : \
fma(x, y, z) ))
#define __tg_real_complex_pow(x, y) (__RETCAST_2(x, y)( \
__FLTCX((x)+(y)) && __IS_FP(x) && __IS_FP(y) ? cpowf(x, y) : \
__FLTCX((x)+(y)) ? cpow(x, y) : \
__DBLCX((x)+(y)) ? cpow(x, y) : \
__LDBLCX((x)+(y)) ? cpowl(x, y) : \
__FLT(x) && __FLT(y) ? powf(x, y) : \
__LDBL((x)+(y)) ? powl(x, y) : \
pow(x, y) ))
#define __tg_real_complex_fabs(x) (__RETCAST(__TO_REAL(x))( \
__FLTCX(x) ? cabsf(x) : \
__DBLCX(x) ? cabs(x) : \
__LDBLCX(x) ? cabsl(x) : \
__FLT(x) ? fabsf(x) : \
__LDBL(x) ? fabsl(x) : \
fabs(x) ))
/* tg functions */
#define acos(x) __tg_real_complex(acos, (x))
#define acosh(x) __tg_real_complex(acosh, (x))
#define asin(x) __tg_real_complex(asin, (x))
#define asinh(x) __tg_real_complex(asinh, (x))
#define atan(x) __tg_real_complex(atan, (x))
#define atan2(x,y) __tg_real_2(atan2, (x), (y))
#define atanh(x) __tg_real_complex(atanh, (x))
#define carg(x) __tg_complex_retreal(carg, (x))
#define cbrt(x) __tg_real(cbrt, (x))
#define ceil(x) __tg_real(ceil, (x))
#define cimag(x) __tg_complex_retreal(cimag, (x))
#define conj(x) __tg_complex(conj, (x))
#define copysign(x,y) __tg_real_2(copysign, (x), (y))
#define cos(x) __tg_real_complex(cos, (x))
#define cosh(x) __tg_real_complex(cosh, (x))
#define cproj(x) __tg_complex(cproj, (x))
#define creal(x) __tg_complex_retreal(creal, (x))
#define erf(x) __tg_real(erf, (x))
#define erfc(x) __tg_real(erfc, (x))
#define exp(x) __tg_real_complex(exp, (x))
#define exp2(x) __tg_real(exp2, (x))
#define expm1(x) __tg_real(expm1, (x))
#define fabs(x) __tg_real_complex_fabs(x)
#define fdim(x,y) __tg_real_2(fdim, (x), (y))
#define floor(x) __tg_real(floor, (x))
#define fma(x,y,z) __tg_real_fma((x), (y), (z))
#define fmax(x,y) __tg_real_2(fmax, (x), (y))
#define fmin(x,y) __tg_real_2(fmin, (x), (y))
#define fmod(x,y) __tg_real_2(fmod, (x), (y))
#define frexp(x,y) __tg_real_2_1(frexp, (x), (y))
#define hypot(x,y) __tg_real_2(hypot, (x), (y))
#define ilogb(x) __tg_real(ilogb, (x))
#define ldexp(x,y) __tg_real_2_1(ldexp, (x), (y))
#define lgamma(x) __tg_real(lgamma, (x))
#define llrint(x) __tg_real(llrint, (x))
#define llround(x) __tg_real(llround, (x))
#define log(x) __tg_real_complex(log, (x))
#define log10(x) __tg_real(log10, (x))
#define log1p(x) __tg_real(log1p, (x))
#define log2(x) __tg_real(log2, (x))
#define logb(x) __tg_real(logb, (x))
#define lrint(x) __tg_real(lrint, (x))
#define lround(x) __tg_real(lround, (x))
#define nearbyint(x) __tg_real(nearbyint, (x))
#define nextafter(x,y) __tg_real_2(nextafter, (x), (y)
#define nexttoward(x,y) __tg_real_2(nexttoward, (x), (y))
#define pow(x,y) __tg_real_complex_pow((x), (y))
#define remainder(x,y) __tg_real_2(remainder, (x), (y))
#define remquo(x,y,z) __tg_real_remquo((x), (y), (z))
#define rint(x) __tg_real(rint, (x))
#define round(x) __tg_real(round, (x))
#define scalbln(x,y) __tg_real_2_1(scalbln, (x), (y))
#define scalbn(x,y) __tg_real_2_1(scalbn, (x), (y))
#define sin(x) __tg_real_complex(sin, (x))
#define sinh(x) __tg_real_complex(sinh, (x))
#define sqrt(x) __tg_real_complex(sqrt, (x))
#define tan(x) __tg_real_complex(tan, (x))
#define tanh(x) __tg_real_complex(tanh, (x))
#define tgamma(x) __tg_real(tgamma, (x))
#define trunc(x) __tg_real(trunc, (x))
#endif
src/complex/__cexp.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */
/*-
* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "libm.h"
static const uint32_t k = 1799; /* constant for reduction */
static const double kln2 = 1246.97177782734161156; /* k * ln2 */
/*
* Compute exp(x), scaled to avoid spurious overflow. An exponent is
* returned separately in 'expt'.
*
* Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
* Output: 2**1023 <= y < 2**1024
*/
static double __frexp_exp(double x, int *expt)
{
double exp_x;
uint32_t hx;
/*
* We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
* minimize |exp(kln2) - 2**k|. We also scale the exponent of
* exp_x to MAX_EXP so that the result can be multiplied by
* a tiny number without losing accuracy due to denormalization.
*/
exp_x = exp(x - kln2);
GET_HIGH_WORD(hx, exp_x);
*expt = (hx >> 20) - (0x3ff + 1023) + k;
SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20));
return exp_x;
}
/*
* __ldexp_cexp(x, expt) compute exp(x) * 2**expt.
* It is intended for large arguments (real part >= ln(DBL_MAX))
* where care is needed to avoid overflow.
*
* The present implementation is narrowly tailored for our hyperbolic and
* exponential functions. We assume expt is small (0 or -1), and the caller
* has filtered out very large x, for which overflow would be inevitable.
*/
double complex __ldexp_cexp(double complex z, int expt)
{
double x, y, exp_x, scale1, scale2;
int ex_expt, half_expt;
x = creal(z);
y = cimag(z);
exp_x = __frexp_exp(x, &ex_expt);
expt += ex_expt;
/*
* Arrange so that scale1 * scale2 == 2**expt. We use this to
* compensate for scalbn being horrendously slow.
*/
half_expt = expt / 2;
INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0);
half_expt = expt - half_expt;
INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
return cpack(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2);
}
src/complex/__cexpf.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/k_expf.c */
/*-
* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "libm.h"
static const uint32_t k = 235; /* constant for reduction */
static const float kln2 = 162.88958740F; /* k * ln2 */
/*
* See __cexp.c for details.
*
* Input: ln(FLT_MAX) <= x < ln(2 * FLT_MAX / FLT_MIN_DENORM) ~= 192.7
* Output: 2**127 <= y < 2**128
*/
static float __frexp_expf(float x, int *expt)
{
float exp_x;
uint32_t hx;
exp_x = expf(x - kln2);
GET_FLOAT_WORD(hx, exp_x);
*expt = (hx >> 23) - (0x7f + 127) + k;
SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 127) << 23));
return exp_x;
}
float complex __ldexp_cexpf(float complex z, int expt)
{
float x, y, exp_x, scale1, scale2;
int ex_expt, half_expt;
x = crealf(z);
y = cimagf(z);
exp_x = __frexp_expf(x, &ex_expt);
expt += ex_expt;
half_expt = expt / 2;
SET_FLOAT_WORD(scale1, (0x7f + half_expt) << 23);
half_expt = expt - half_expt;
SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);
return cpackf(cosf(y) * exp_x * scale1 * scale2,
sinf(y) * exp_x * scale1 * scale2);
}
src/complex/cabs.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
double cabs(double complex z)
{
return hypot(creal(z), cimag(z));
}
src/complex/cabsf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float cabsf(float complex z)
{
return hypotf(crealf(z), cimagf(z));
}
src/complex/cabsl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double cabsl(long double complex z)
{
return cabs(z);
}
#else
long double cabsl(long double complex z)
{
return hypotl(creall(z), cimagl(z));
}
#endif
src/complex/cacos.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
// FIXME: Hull et al. "Implementing the complex arcsine and arccosine functions using exception handling" 1997
/* acos(z) = pi/2 - asin(z) */
double complex cacos(double complex z)
{
z = casin(z);
return cpack(M_PI_2 - creal(z), -cimag(z));
}
src/complex/cacosf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
// FIXME
float complex cacosf(float complex z)
{
z = casinf(z);
return cpackf((float)M_PI_2 - crealf(z), -cimagf(z));
}
src/complex/cacosh.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
/* acosh(z) = i acos(z) */
double complex cacosh(double complex z)
{
z = cacos(z);
return cpack(-cimag(z), creal(z));
}
src/complex/cacoshf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float complex cacoshf(float complex z)
{
z = cacosf(z);
return cpackf(-cimagf(z), crealf(z));
}
src/complex/cacoshl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex cacoshl(long double complex z)
{
return cacosh(z);
}
#else
long double complex cacoshl(long double complex z)
{
z = cacosl(z);
return cpackl(-cimagl(z), creall(z));
}
#endif
src/complex/cacosl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex cacosl(long double complex z)
{
return cacos(z);
}
#else
// FIXME
#define PI_2 1.57079632679489661923132169163975144L
long double complex cacosl(long double complex z)
{
z = casinl(z);
return cpackl(PI_2 - creall(z), -cimagl(z));
}
#endif
src/complex/carg.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
double carg(double complex z)
{
return atan2(cimag(z), creal(z));
}
src/complex/cargf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float cargf(float complex z)
{
return atan2f(cimagf(z), crealf(z));
}
src/complex/cargl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double cargl(long double complex z)
{
return carg(z);
}
#else
long double cargl(long double complex z)
{
return atan2l(cimagl(z), creall(z));
}
#endif
src/complex/casin.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
// FIXME
/* asin(z) = -i log(i z + sqrt(1 - z*z)) */
double complex casin(double complex z)
{
double complex w;
double x, y;
x = creal(z);
y = cimag(z);
w = cpack(1.0 - (x - y)*(x + y), -2.0*x*y);
return clog(cpack(-y, x) + csqrt(w));
}
src/complex/casinf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
// FIXME
float complex casinf(float complex z)
{
float complex w;
float x, y;
x = crealf(z);
y = cimagf(z);
w = cpackf(1.0 - (x - y)*(x + y), -2.0*x*y);
return clogf(cpackf(-y, x) + csqrtf(w));
}
src/complex/casinh.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
/* asinh(z) = -i asin(i z) */
double complex casinh(double complex z)
{
z = casin(cpack(-cimag(z), creal(z)));
return cpack(cimag(z), -creal(z));
}
src/complex/casinhf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float complex casinhf(float complex z)
{
z = casinf(cpackf(-cimagf(z), crealf(z)));
return cpackf(cimagf(z), -crealf(z));
}
src/complex/casinhl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex casinhl(long double complex z)
{
return casinh(z);
}
#else
long double complex casinhl(long double complex z)
{
z = casinl(cpackl(-cimagl(z), creall(z)));
return cpackl(cimagl(z), -creall(z));
}
#endif
src/complex/casinl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex casinl(long double complex z)
{
return casin(z);
}
#else
// FIXME
long double complex casinl(long double complex z)
{
long double complex w;
long double x, y;
x = creall(z);
y = cimagl(z);
w = cpackl(1.0 - (x - y)*(x + y), -2.0*x*y);
return clogl(cpackl(-y, x) + csqrtl(w));
}
#endif
src/complex/catan.c 0 → 100644
浏览文件 @ b69f695a
/* origin: OpenBSD /usr/src/lib/libm/src/s_catan.c */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/*
* Complex circular arc tangent
*
*
* SYNOPSIS:
*
* double complex catan();
* double complex z, w;
*
* w = catan (z);
*
*
* DESCRIPTION:
*
* If
* z = x + iy,
*
* then
* 1 ( 2x )
* Re w = - arctan(-----------) + k PI
* 2 ( 2 2)
* (1 - x - y )
*
* ( 2 2)
* 1 (x + (y+1) )
* Im w = - log(------------)
* 4 ( 2 2)
* (x + (y-1) )
*
* Where k is an arbitrary integer.
*
* catan(z) = -i catanh(iz).
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -10,+10 5900 1.3e-16 7.8e-18
* IEEE -10,+10 30000 2.3e-15 8.5e-17
* The check catan( ctan(z) ) = z, with |x| and |y| < PI/2,
* had peak relative error 1.5e-16, rms relative error
* 2.9e-17. See also clog().
*/
#include "libm.h"
#define MAXNUM 1.0e308
static const double DP1 = 3.14159265160560607910E0;
static const double DP2 = 1.98418714791870343106E-9;
static const double DP3 = 1.14423774522196636802E-17;
static double _redupi(double x)
{
double t;
long i;
t = x/M_PI;
if (t >= 0.0)
t += 0.5;
else
t -= 0.5;
i = t; /* the multiple */
t = i;
t = ((x - t * DP1) - t * DP2) - t * DP3;
return t;
}
double complex catan(double complex z)
{
double complex w;
double a, t, x, x2, y;
x = creal(z);
y = cimag(z);
if (x == 0.0 && y > 1.0)
goto ovrf;
x2 = x * x;
a = 1.0 - x2 - (y * y);
if (a == 0.0)
goto ovrf;
t = 0.5 * atan2(2.0 * x, a);
w = _redupi(t);
t = y - 1.0;
a = x2 + (t * t);
if (a == 0.0)
goto ovrf;
t = y + 1.0;
a = (x2 + t * t)/a;
w = w + (0.25 * log(a)) * I;
return w;
ovrf:
// FIXME
w = MAXNUM + MAXNUM * I;
return w;
}
src/complex/catanf.c 0 → 100644
浏览文件 @ b69f695a
/* origin: OpenBSD /usr/src/lib/libm/src/s_catanf.c */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/*
* Complex circular arc tangent
*
*
* SYNOPSIS:
*
* float complex catanf();
* float complex z, w;
*
* w = catanf( z );
*
*
* DESCRIPTION:
*
* If
* z = x + iy,
*
* then
* 1 ( 2x )
* Re w = - arctan(-----------) + k PI
* 2 ( 2 2)
* (1 - x - y )
*
* ( 2 2)
* 1 (x + (y+1) )
* Im w = - log(------------)
* 4 ( 2 2)
* (x + (y-1) )
*
* Where k is an arbitrary integer.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 2.3e-6 5.2e-8
*/
#include "libm.h"
#define MAXNUMF 1.0e38F
static const double DP1 = 3.140625;
static const double DP2 = 9.67502593994140625E-4;
static const double DP3 = 1.509957990978376432E-7;
static float _redupif(float xx)
{
float x, t;
long i;
x = xx;
t = x/(float)M_PI;
if (t >= 0.0f)
t += 0.5f;
else
t -= 0.5f;
i = t; /* the multiple */
t = i;
t = ((x - t * DP1) - t * DP2) - t * DP3;
return t;
}
float complex catanf(float complex z)
{
float complex w;
float a, t, x, x2, y;
x = crealf(z);
y = cimagf(z);
if ((x == 0.0f) && (y > 1.0f))
goto ovrf;
x2 = x * x;
a = 1.0f - x2 - (y * y);
if (a == 0.0f)
goto ovrf;
t = 0.5f * atan2f(2.0f * x, a);
w = _redupif(t);
t = y - 1.0f;
a = x2 + (t * t);
if (a == 0.0f)
goto ovrf;
t = y + 1.0f;
a = (x2 + (t * t))/a;
w = w + (0.25f * logf (a)) * I;
return w;
ovrf:
// FIXME
w = MAXNUMF + MAXNUMF * I;
return w;
}
src/complex/catanh.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
/* atanh = -i atan(i z) */
double complex catanh(double complex z)
{
z = catan(cpack(-cimag(z), creal(z)));
return cpack(cimag(z), -creal(z));
}
src/complex/catanhf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float complex catanhf(float complex z)
{
z = catanf(cpackf(-cimagf(z), crealf(z)));
return cpackf(cimagf(z), -crealf(z));
}
src/complex/catanhl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex catanhl(long double complex z)
{
return catanh(z);
}
#else
long double complex catanhl(long double complex z)
{
z = catanl(cpackl(-cimagl(z), creall(z)));
return cpackl(cimagl(z), -creall(z));
}
#endif
src/complex/catanl.c 0 → 100644
浏览文件 @ b69f695a
/* origin: OpenBSD /usr/src/lib/libm/src/s_catanl.c */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/*
* Complex circular arc tangent
*
*
* SYNOPSIS:
*
* long double complex catanl();
* long double complex z, w;
*
* w = catanl( z );
*
*
* DESCRIPTION:
*
* If
* z = x + iy,
*
* then
* 1 ( 2x )
* Re w = - arctan(-----------) + k PI
* 2 ( 2 2)
* (1 - x - y )
*
* ( 2 2)
* 1 (x + (y+1) )
* Im w = - log(------------)
* 4 ( 2 2)
* (x + (y-1) )
*
* Where k is an arbitrary integer.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -10,+10 5900 1.3e-16 7.8e-18
* IEEE -10,+10 30000 2.3e-15 8.5e-17
* The check catan( ctan(z) ) = z, with |x| and |y| < PI/2,
* had peak relative error 1.5e-16, rms relative error
* 2.9e-17. See also clog().
*/
#include <complex.h>
#include <float.h>
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex catanl(long double complex z)
{
return catan(z);
}
#else
static const long double PIL = 3.141592653589793238462643383279502884197169L;
static const long double DP1 = 3.14159265358979323829596852490908531763125L;
static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L;
static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L;
static long double redupil(long double x)
{
long double t;
long i;
t = x / PIL;
if (t >= 0.0L)
t += 0.5L;
else
t -= 0.5L;
i = t; /* the multiple */
t = i;
t = ((x - t * DP1) - t * DP2) - t * DP3;
return t;
}
long double complex catanl(long double complex z)
{
long double complex w;
long double a, t, x, x2, y;
x = creall(z);
y = cimagl(z);
if ((x == 0.0L) && (y > 1.0L))
goto ovrf;
x2 = x * x;
a = 1.0L - x2 - (y * y);
if (a == 0.0L)
goto ovrf;
t = atan2l(2.0L * x, a) * 0.5L;
w = redupil(t);
t = y - 1.0L;
a = x2 + (t * t);
if (a == 0.0L)
goto ovrf;
t = y + 1.0L;
a = (x2 + (t * t)) / a;
w = w + (0.25L * logl(a)) * I;
return w;
ovrf:
// FIXME
w = LDBL_MAX + LDBL_MAX * I;
return w;
}
#endif
src/complex/ccos.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
/* cos(z) = cosh(i z) */
double complex ccos(double complex z)
{
return ccosh(cpack(-cimag(z), creal(z)));
}
src/complex/ccosf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float complex ccosf(float complex z)
{
return ccoshf(cpackf(-cimagf(z), crealf(z)));
}
src/complex/ccosh.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_ccosh.c */
/*-
* Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/*
* Hyperbolic cosine of a complex argument z = x + i y.
*
* cosh(z) = cosh(x+iy)
* = cosh(x) cos(y) + i sinh(x) sin(y).
*
* Exceptional values are noted in the comments within the source code.
* These values and the return value were taken from n1124.pdf.
*/
#include "libm.h"
static const double huge = 0x1p1023;
double complex ccosh(double complex z)
{
double x, y, h;
int32_t hx, hy, ix, iy, lx, ly;
x = creal(z);
y = cimag(z);
EXTRACT_WORDS(hx, lx, x);
EXTRACT_WORDS(hy, ly, y);
ix = 0x7fffffff & hx;
iy = 0x7fffffff & hy;
/* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0)
return cpack(cosh(x), x * y);
if (ix < 0x40360000) /* small x: normal case */
return cpack(cosh(x) * cos(y), sinh(x) * sin(y));
/* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5;
return cpack(h * cos(y), copysign(h, x) * sin(y));
} else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */
z = __ldexp_cexp(cpack(fabs(x), y), -1);
return cpack(creal(z), cimag(z) * copysign(1, x));
} else {
/* x >= 1455: the result always overflows */
h = huge * x;
return cpack(h * h * cos(y), h * sin(y));
}
}
/*
* cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as dNaN. Raise the invalid floating-point exception.
*
* cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as d(NaN).
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
return cpack(y - y, copysign(0, x * (y - y)));
/*
* cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
*
* cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0.
* The sign of 0 in the result is unspecified.
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
return cpack(x * x, copysign(0, x) * y);
return cpack(x * x, copysign(0, (x + x) * y));
}
/*
* cosh(x +- I Inf) = dNaN + I dNaN.
* Raise the invalid floating-point exception for finite nonzero x.
*
* cosh(x + I NaN) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception for finite
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
return cpack(y - y, x * (y - y));
/*
* cosh(+-Inf + I NaN) = +Inf + I d(NaN).
*
* cosh(+-Inf +- I Inf) = +Inf + I dNaN.
* The sign of Inf in the result is unspecified. Choice = always +.
* Raise the invalid floating-point exception.
*
* cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (iy >= 0x7ff00000)
return cpack(x * x, x * (y - y));
return cpack((x * x) * cos(y), x * sin(y));
}
/*
* cosh(NaN + I NaN) = d(NaN) + I d(NaN).
*
* cosh(NaN +- I Inf) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception.
* Choice = raise.
*
* cosh(NaN + I y) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
return cpack((x * x) * (y - y), (x + x) * (y - y));
}
src/complex/ccoshf.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_ccoshf.c */
/*-
* Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/*
* Hyperbolic cosine of a complex argument. See s_ccosh.c for details.
*/
#include "libm.h"
static const float huge = 0x1p127;
float complex ccoshf(float complex z)
{
float x, y, h;
int32_t hx, hy, ix, iy;
x = crealf(z);
y = cimagf(z);
GET_FLOAT_WORD(hx, x);
GET_FLOAT_WORD(hy, y);
ix = 0x7fffffff & hx;
iy = 0x7fffffff & hy;
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
return cpackf(coshf(x), x * y);
if (ix < 0x41100000) /* small x: normal case */
return cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y));
/* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f;
return cpackf(h * cosf(y), copysignf(h, x) * sinf(y));
} else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */
z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
return cpackf(crealf(z), cimagf(z) * copysignf(1, x));
} else {
/* x >= 192.7: the result always overflows */
h = huge * x;
return cpackf(h * h * cosf(y), h * sinf(y));
}
}
if (ix == 0 && iy >= 0x7f800000)
return cpackf(y - y, copysignf(0, x * (y - y)));
if (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0)
return cpackf(x * x, copysignf(0, x) * y);
return cpackf(x * x, copysignf(0, (x + x) * y));
}
if (ix < 0x7f800000 && iy >= 0x7f800000)
return cpackf(y - y, x * (y - y));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (iy >= 0x7f800000)
return cpackf(x * x, x * (y - y));
return cpackf((x * x) * cosf(y), x * sinf(y));
}
return cpackf((x * x) * (y - y), (x + x) * (y - y));
}
src/complex/ccoshl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
//FIXME
long double complex ccoshl(long double complex z)
{
return ccosh(z);
}
src/complex/ccosl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex ccosl(long double complex z)
{
return ccos(z);
}
#else
long double complex ccosl(long double complex z)
{
return ccoshl(cpackl(-cimagl(z), creall(z)));
}
#endif
src/complex/cexp.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_cexp.c */
/*-
* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "libm.h"
static const uint32_t
exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */
cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
double complex cexp(double complex z)
{
double x, y, exp_x;
uint32_t hx, hy, lx, ly;
x = creal(z);
y = cimag(z);
EXTRACT_WORDS(hy, ly, y);
hy &= 0x7fffffff;
/* cexp(x + I 0) = exp(x) + I 0 */
if ((hy | ly) == 0)
return cpack(exp(x), y);
EXTRACT_WORDS(hx, lx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if (((hx & 0x7fffffff) | lx) == 0)
return cpack(cos(y), sin(y));
if (hy >= 0x7ff00000) {
if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
return cpack(y - y, y - y);
} else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
return cpack(0.0, 0.0);
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
return cpack(x, y - y);
}
}
if (hx >= exp_ovfl && hx <= cexp_ovfl) {
/*
* x is between 709.7 and 1454.3, so we must scale to avoid
* overflow in exp(x).
*/
return __ldexp_cexp(z, 0);
} else {
/*
* Cases covered here:
* - x < exp_ovfl and exp(x) won't overflow (common case)
* - x > cexp_ovfl, so exp(x) * s overflows for all s > 0
* - x = +-Inf (generated by exp())
* - x = NaN (spurious inexact exception from y)
*/
exp_x = exp(x);
return cpack(exp_x * cos(y), exp_x * sin(y));
}
}
src/complex/cexpf.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_cexpf.c */
/*-
* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "libm.h"
static const uint32_t
exp_ovfl = 0x42b17218, /* MAX_EXP * ln2 ~= 88.722839355 */
cexp_ovfl = 0x43400074; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
float complex cexpf(float complex z)
{
float x, y, exp_x;
uint32_t hx, hy;
x = crealf(z);
y = cimagf(z);
GET_FLOAT_WORD(hy, y);
hy &= 0x7fffffff;
/* cexp(x + I 0) = exp(x) + I 0 */
if (hy == 0)
return cpackf(expf(x), y);
GET_FLOAT_WORD(hx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if ((hx & 0x7fffffff) == 0)
return cpackf(cosf(y), sinf(y));
if (hy >= 0x7f800000) {
if ((hx & 0x7fffffff) != 0x7f800000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
return cpackf(y - y, y - y);
} else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
return cpackf(0.0, 0.0);
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
return cpackf(x, y - y);
}
}
if (hx >= exp_ovfl && hx <= cexp_ovfl) {
/*
* x is between 88.7 and 192, so we must scale to avoid
* overflow in expf(x).
*/
return __ldexp_cexpf(z, 0);
} else {
/*
* Cases covered here:
* - x < exp_ovfl and exp(x) won't overflow (common case)
* - x > cexp_ovfl, so exp(x) * s overflows for all s > 0
* - x = +-Inf (generated by exp())
* - x = NaN (spurious inexact exception from y)
*/
exp_x = expf(x);
return cpackf(exp_x * cosf(y), exp_x * sinf(y));
}
}
src/complex/cexpl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
//FIXME
long double complex cexpl(long double complex z)
{
return cexp(z);
}
src/complex/cimag.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
double (cimag)(double complex z)
{
union dcomplex u = {z};
return u.a[1];
}
src/complex/cimagf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float (cimagf)(float complex z)
{
union fcomplex u = {z};
return u.a[1];
}
src/complex/cimagl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
long double (cimagl)(long double complex z)
{
union lcomplex u = {z};
return u.a[1];
}
src/complex/clog.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
// FIXME
/* log(z) = log(|z|) + i arg(z) */
double complex clog(double complex z)
{
double r, phi;
r = cabs(z);
phi = carg(z);
return cpack(log(r), phi);
}
src/complex/clogf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
// FIXME
float complex clogf(float complex z)
{
float r, phi;
r = cabsf(z);
phi = cargf(z);
return cpackf(logf(r), phi);
}
src/complex/clogl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex clogl(long double complex z)
{
return clog(z);
}
#else
// FIXME
long double complex clogl(long double complex z)
{
long double r, phi;
r = cabsl(z);
phi = cargl(z);
return cpackl(logl(r), phi);
}
#endif
src/complex/conj.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
double complex conj(double complex z)
{
return cpack(creal(z), -cimag(z));
}
src/complex/conjf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float complex conjf(float complex z)
{
return cpackf(crealf(z), -cimagf(z));
}
src/complex/conjl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
long double complex conjl(long double complex z)
{
return cpackl(creall(z), -cimagl(z));
}
src/complex/cpow.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
/* pow(z, c) = exp(c log(z)), See C99 G.6.4.1 */
double complex cpow(double complex z, double complex c)
{
return cexp(c * clog(z));
}
src/complex/cpowf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float complex cpowf(float complex z, float complex c)
{
return cexpf(c * clogf(z));
}
src/complex/cpowl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex cpowl(long double complex z, long double complex c)
{
return cpow(z, c);
}
#else
long double complex cpowl(long double complex z, long double complex c)
{
return cexpl(c * clogl(z));
}
#endif
src/complex/cproj.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
double complex cproj(double complex z)
{
if (isinf(creal(z)) || isinf(cimag(z)))
return cpack(INFINITY, copysign(0.0, creal(z)));
return z;
}
src/complex/cprojf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float complex cprojf(float complex z)
{
if (isinf(crealf(z)) || isinf(cimagf(z)))
return cpackf(INFINITY, copysignf(0.0, crealf(z)));
return z;
}
src/complex/cprojl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex cprojl(long double complex z)
{
return cproj(z);
}
#else
long double complex cprojl(long double complex z)
{
if (isinf(creall(z)) || isinf(cimagl(z)))
return cpackl(INFINITY, copysignl(0.0, creall(z)));
return z;
}
#endif
src/complex/creal.c 0 → 100644
浏览文件 @ b69f695a
#include <complex.h>
double creal(double complex z)
{
return z;
}
src/complex/crealf.c 0 → 100644
浏览文件 @ b69f695a
#include <complex.h>
float crealf(float complex z)
{
return z;
}
src/complex/creall.c 0 → 100644
浏览文件 @ b69f695a
#include <complex.h>
long double creall(long double complex z)
{
return z;
}
src/complex/csin.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
/* sin(z) = -i sinh(i z) */
double complex csin(double complex z)
{
z = csinh(cpack(-cimag(z), creal(z)));
return cpack(cimag(z), -creal(z));
}
src/complex/csinf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float complex csinf(float complex z)
{
z = csinhf(cpackf(-cimagf(z), crealf(z)));
return cpackf(cimagf(z), -crealf(z));
}
src/complex/csinh.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_csinh.c */
/*-
* Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/*
* Hyperbolic sine of a complex argument z = x + i y.
*
* sinh(z) = sinh(x+iy)
* = sinh(x) cos(y) + i cosh(x) sin(y).
*
* Exceptional values are noted in the comments within the source code.
* These values and the return value were taken from n1124.pdf.
*/
#include "libm.h"
static const double huge = 0x1p1023;
double complex csinh(double complex z)
{
double x, y, h;
int32_t hx, hy, ix, iy, lx, ly;
x = creal(z);
y = cimag(z);
EXTRACT_WORDS(hx, lx, x);
EXTRACT_WORDS(hy, ly, y);
ix = 0x7fffffff & hx;
iy = 0x7fffffff & hy;
/* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0)
return cpack(sinh(x), y);
if (ix < 0x40360000) /* small x: normal case */
return cpack(sinh(x) * cos(y), cosh(x) * sin(y));
/* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5;
return cpack(copysign(h, x) * cos(y), h * sin(y));
} else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */
z = __ldexp_cexp(cpack(fabs(x), y), -1);
return cpack(creal(z) * copysign(1, x), cimag(z));
} else {
/* x >= 1455: the result always overflows */
h = huge * x;
return cpack(h * cos(y), h * h * sin(y));
}
}
/*
* sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as dNaN. Raise the invalid floating-point exception.
*
* sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN).
* The sign of 0 in the result is unspecified. Choice = normally
* the same as d(NaN).
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
return cpack(copysign(0, x * (y - y)), y - y);
/*
* sinh(+-Inf +- I 0) = +-Inf + I +-0.
*
* sinh(NaN +- I 0) = d(NaN) + I +-0.
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
return cpack(x, y);
return cpack(x, copysign(0, y));
}
/*
* sinh(x +- I Inf) = dNaN + I dNaN.
* Raise the invalid floating-point exception for finite nonzero x.
*
* sinh(x + I NaN) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception for finite
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
return cpack(y - y, x * (y - y));
/*
* sinh(+-Inf + I NaN) = +-Inf + I d(NaN).
* The sign of Inf in the result is unspecified. Choice = normally
* the same as d(NaN).
*
* sinh(+-Inf +- I Inf) = +Inf + I dNaN.
* The sign of Inf in the result is unspecified. Choice = always +.
* Raise the invalid floating-point exception.
*
* sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (iy >= 0x7ff00000)
return cpack(x * x, x * (y - y));
return cpack(x * cos(y), INFINITY * sin(y));
}
/*
* sinh(NaN + I NaN) = d(NaN) + I d(NaN).
*
* sinh(NaN +- I Inf) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception.
* Choice = raise.
*
* sinh(NaN + I y) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
return cpack((x * x) * (y - y), (x + x) * (y - y));
}
src/complex/csinhf.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_csinhf.c */
/*-
* Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/*
* Hyperbolic sine of a complex argument z. See s_csinh.c for details.
*/
#include "libm.h"
static const float huge = 0x1p127;
float complex csinhf(float complex z)
{
float x, y, h;
int32_t hx, hy, ix, iy;
x = crealf(z);
y = cimagf(z);
GET_FLOAT_WORD(hx, x);
GET_FLOAT_WORD(hy, y);
ix = 0x7fffffff & hx;
iy = 0x7fffffff & hy;
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
return cpackf(sinhf(x), y);
if (ix < 0x41100000) /* small x: normal case */
return cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y));
/* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f;
return cpackf(copysignf(h, x) * cosf(y), h * sinf(y));
} else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */
z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
return cpackf(crealf(z) * copysignf(1, x), cimagf(z));
} else {
/* x >= 192.7: the result always overflows */
h = huge * x;
return cpackf(h * cosf(y), h * h * sinf(y));
}
}
if (ix == 0 && iy >= 0x7f800000)
return cpackf(copysignf(0, x * (y - y)), y - y);
if (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0)
return cpackf(x, y);
return cpackf(x, copysignf(0, y));
}
if (ix < 0x7f800000 && iy >= 0x7f800000)
return cpackf(y - y, x * (y - y));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (iy >= 0x7f800000)
return cpackf(x * x, x * (y - y));
return cpackf(x * cosf(y), INFINITY * sinf(y));
}
return cpackf((x * x) * (y - y), (x + x) * (y - y));
}
src/complex/csinhl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
//FIXME
long double complex csinhl(long double complex z)
{
return csinh(z);
}
src/complex/csinl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex csinl(long double complex z)
{
return csin(z);
}
#else
long double complex csinl(long double complex z)
{
z = csinhl(cpackl(-cimagl(z), creall(z)));
return cpackl(cimagl(z), -creall(z));
}
#endif
src/complex/csqrt.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_csqrt.c */
/*-
* Copyright (c) 2007 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "libm.h"
/*
* gcc doesn't implement complex multiplication or division correctly,
* so we need to handle infinities specially. We turn on this pragma to
* notify conforming c99 compilers that the fast-but-incorrect code that
* gcc generates is acceptable, since the special cases have already been
* handled.
*/
#pragma STDC CX_LIMITED_RANGE ON
/* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */
#define THRESH 0x1.a827999fcef32p+1022
double complex csqrt(double complex z)
{
double complex result;
double a, b;
double t;
int scale;
a = creal(z);
b = cimag(z);
/* Handle special cases. */
if (z == 0)
return cpack(0, b);
if (isinf(b))
return cpack(INFINITY, b);
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return cpack(a, t); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
* csqrt(inf + NaN i) = inf + NaN i
* csqrt(inf + y i) = inf + 0 i
* csqrt(-inf + NaN i) = NaN +- inf i
* csqrt(-inf + y i) = 0 + inf i
*/
if (signbit(a))
return cpack(fabs(b - b), copysign(a, b));
else
return cpack(a, copysign(b - b, b));
}
/*
* The remaining special case (b is NaN) is handled just fine by
* the normal code path below.
*/
/* Scale to avoid overflow. */
if (fabs(a) >= THRESH || fabs(b) >= THRESH) {
a *= 0.25;
b *= 0.25;
scale = 1;
} else {
scale = 0;
}
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
result = cpack(t, b / (2 * t));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
result = cpack(fabs(b) / (2 * t), copysign(t, b));
}
/* Rescale. */
if (scale)
result *= 2;
return result;
}
src/complex/csqrtf.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_csqrtf.c */
/*-
* Copyright (c) 2007 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "libm.h"
/*
* gcc doesn't implement complex multiplication or division correctly,
* so we need to handle infinities specially. We turn on this pragma to
* notify conforming c99 compilers that the fast-but-incorrect code that
* gcc generates is acceptable, since the special cases have already been
* handled.
*/
#pragma STDC CX_LIMITED_RANGE ON
float complex csqrtf(float complex z)
{
float a = crealf(z), b = cimagf(z);
double t;
/* Handle special cases. */
if (z == 0)
return cpackf(0, b);
if (isinf(b))
return cpackf(INFINITY, b);
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return cpackf(a, t); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
* csqrtf(inf + NaN i) = inf + NaN i
* csqrtf(inf + y i) = inf + 0 i
* csqrtf(-inf + NaN i) = NaN +- inf i
* csqrtf(-inf + y i) = 0 + inf i
*/
if (signbit(a))
return cpackf(fabsf(b - b), copysignf(a, b));
else
return cpackf(a, copysignf(b - b, b));
}
/*
* The remaining special case (b is NaN) is handled just fine by
* the normal code path below.
*/
/*
* We compute t in double precision to avoid overflow and to
* provide correct rounding in nearly all cases.
* This is Algorithm 312, CACM vol 10, Oct 1967.
*/
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
return cpackf(t, b / (2.0 * t));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
return cpackf(fabsf(b) / (2.0 * t), copysignf(t, b));
}
}
src/complex/csqrtl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
//FIXME
long double complex csqrtl(long double complex z)
{
return csqrt(z);
}
src/complex/ctan.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
/* tan(z) = -i tanh(i z) */
double complex ctan(double complex z)
{
z = ctanh(cpack(-cimag(z), creal(z)));
return cpack(cimag(z), -creal(z));
}
src/complex/ctanf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
float complex ctanf(float complex z)
{
z = ctanhf(cpackf(-cimagf(z), crealf(z)));
return cpackf(cimagf(z), -crealf(z));
}
src/complex/ctanh.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_ctanh.c */
/*-
* Copyright (c) 2011 David Schultz
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/*
* Hyperbolic tangent of a complex argument z = x + i y.
*
* The algorithm is from:
*
* W. Kahan. Branch Cuts for Complex Elementary Functions or Much
* Ado About Nothing's Sign Bit. In The State of the Art in
* Numerical Analysis, pp. 165 ff. Iserles and Powell, eds., 1987.
*
* Method:
*
* Let t = tan(x)
* beta = 1/cos^2(y)
* s = sinh(x)
* rho = cosh(x)
*
* We have:
*
* tanh(z) = sinh(z) / cosh(z)
*
* sinh(x) cos(y) + i cosh(x) sin(y)
* = ---------------------------------
* cosh(x) cos(y) + i sinh(x) sin(y)
*
* cosh(x) sinh(x) / cos^2(y) + i tan(y)
* = -------------------------------------
* 1 + sinh^2(x) / cos^2(y)
*
* beta rho s + i t
* = ----------------
* 1 + beta s^2
*
* Modifications:
*
* I omitted the original algorithm's handling of overflow in tan(x) after
* verifying with nearpi.c that this can't happen in IEEE single or double
* precision. I also handle large x differently.
*/
#include "libm.h"
double complex ctanh(double complex z)
{
double x, y;
double t, beta, s, rho, denom;
uint32_t hx, ix, lx;
x = creal(z);
y = cimag(z);
EXTRACT_WORDS(hx, lx, x);
ix = hx & 0x7fffffff;
/*
* ctanh(NaN + i 0) = NaN + i 0
*
* ctanh(NaN + i y) = NaN + i NaN for y != 0
*
* The imaginary part has the sign of x*sin(2*y), but there's no
* special effort to get this right.
*
* ctanh(+-Inf +- i Inf) = +-1 +- 0
*
* ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite
*
* The imaginary part of the sign is unspecified. This special
* case is only needed to avoid a spurious invalid exception when
* y is infinite.
*/
if (ix >= 0x7ff00000) {
if ((ix & 0xfffff) | lx) /* x is NaN */
return cpack(x, (y == 0 ? y : x * y));
SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */
return cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)));
}
/*
* ctanh(x + i NAN) = NaN + i NaN
* ctanh(x +- i Inf) = NaN + i NaN
*/
if (!isfinite(y))
return cpack(y - y, y - y);
/*
* ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
* approximation sinh^2(huge) ~= exp(2*huge) / 4.
* We use a modified formula to avoid spurious overflow.
*/
if (ix >= 0x40360000) { /* x >= 22 */
double exp_mx = exp(-fabs(x));
return cpack(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx);
}
/* Kahan's algorithm */
t = tan(y);
beta = 1.0 + t * t; /* = 1 / cos^2(y) */
s = sinh(x);
rho = sqrt(1 + s * s); /* = cosh(x) */
denom = 1 + beta * s * s;
return cpack((beta * rho * s) / denom, t / denom);
}
src/complex/ctanhf.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/s_ctanhf.c */
/*-
* Copyright (c) 2011 David Schultz
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/*
* Hyperbolic tangent of a complex argument z. See s_ctanh.c for details.
*/
#include "libm.h"
float complex ctanhf(float complex z)
{
float x, y;
float t, beta, s, rho, denom;
uint32_t hx, ix;
x = crealf(z);
y = cimagf(z);
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000) {
if (ix & 0x7fffff)
return cpackf(x, (y == 0 ? y : x * y));
SET_FLOAT_WORD(x, hx - 0x40000000);
return cpackf(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)));
}
if (!isfinite(y))
return cpackf(y - y, y - y);
if (ix >= 0x41300000) { /* x >= 11 */
float exp_mx = expf(-fabsf(x));
return cpackf(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx);
}
t = tanf(y);
beta = 1.0 + t * t;
s = sinhf(x);
rho = sqrtf(1 + s * s);
denom = 1 + beta * s * s;
return cpackf((beta * rho * s) / denom, t / denom);
}
src/complex/ctanhl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
//FIXME
long double complex ctanhl(long double complex z)
{
return ctanh(z);
}
src/complex/ctanl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double complex ctanl(long double complex z)
{
return ctan(z);
}
#else
long double complex ctanl(long double complex z)
{
z = ctanhl(cpackl(-cimagl(z), creall(z)));
return cpackl(cimagl(z), -creall(z));
}
#endif
src/internal/libm.h 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/math_private.h */
/*
* ====================================================
* 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.
* ====================================================
*/
#ifndef _LIBM_H
#define _LIBM_H
#include <stdint.h>
#include <float.h>
#include <math.h>
#include <complex.h>
#include "longdbl.h"
union fshape {
float value;
uint32_t bits;
};
union dshape {
double value;
uint64_t bits;
};
/* Get two 32 bit ints from a double. */
#define EXTRACT_WORDS(hi,lo,d) \
do { \
union dshape __u; \
__u.value = (d); \
(hi) = __u.bits >> 32; \
(lo) = (uint32_t)__u.bits; \
} while (0)
/* Get a 64 bit int from a double. */
#define EXTRACT_WORD64(i,d) \
do { \
union dshape __u; \
__u.value = (d); \
(i) = __u.bits; \
} while (0)
/* Get the more significant 32 bit int from a double. */
#define GET_HIGH_WORD(i,d) \
do { \
union dshape __u; \
__u.value = (d); \
(i) = __u.bits >> 32; \
} while (0)
/* Get the less significant 32 bit int from a double. */
#define GET_LOW_WORD(i,d) \
do { \
union dshape __u; \
__u.value = (d); \
(i) = (uint32_t)__u.bits; \
} while (0)
/* Set a double from two 32 bit ints. */
#define INSERT_WORDS(d,hi,lo) \
do { \
union dshape __u; \
__u.bits = ((uint64_t)(hi) << 32) | (uint32_t)(lo); \
(d) = __u.value; \
} while (0)
/* Set a double from a 64 bit int. */
#define INSERT_WORD64(d,i) \
do { \
union dshape __u; \
__u.bits = (i); \
(d) = __u.value; \
} while (0)
/* Set the more significant 32 bits of a double from an int. */
#define SET_HIGH_WORD(d,hi) \
do { \
union dshape __u; \
__u.value = (d); \
__u.bits &= 0xffffffff; \
__u.bits |= (uint64_t)(hi) << 32; \
(d) = __u.value; \
} while (0)
/* Set the less significant 32 bits of a double from an int. */
#define SET_LOW_WORD(d,lo) \
do { \
union dshape __u; \
__u.value = (d); \
__u.bits &= 0xffffffff00000000ull; \
__u.bits |= (uint32_t)(lo); \
(d) = __u.value; \
} while (0)
/* Get a 32 bit int from a float. */
#define GET_FLOAT_WORD(i,d) \
do { \
union fshape __u; \
__u.value = (d); \
(i) = __u.bits; \
} while (0)
/* Set a float from a 32 bit int. */
#define SET_FLOAT_WORD(d,i) \
do { \
union fshape __u; \
__u.bits = (i); \
(d) = __u.value; \
} while (0)
/* fdlibm kernel functions */
int __rem_pio2_large(double*,double*,int,int,int);
int __rem_pio2(double,double*);
double __sin(double,double,int);
double __cos(double,double);
double __tan(double,double,int);
double __expo2(double);
double complex __ldexp_cexp(double complex,int);
int __rem_pio2f(float,double*);
float __sindf(double);
float __cosdf(double);
float __tandf(double,int);
float __expo2f(float);
float complex __ldexp_cexpf(float complex,int);
long double __sinl(long double, long double, int);
long double __cosl(long double, long double);
long double __tanl(long double, long double, int);
/* polynomial evaluation */
long double __polevll(long double, long double *, int);
long double __p1evll(long double, long double *, int);
// FIXME: not needed when -fexcess-precision=standard is supported (>=gcc4.5)
/*
* Attempt to get strict C99 semantics for assignment with non-C99 compilers.
*/
#if 1
#define STRICT_ASSIGN(type, lval, rval) do { \
volatile type __v = (rval); \
(lval) = __v; \
} while (0)
#else
#define STRICT_ASSIGN(type, lval, rval) ((lval) = (type)(rval))
#endif
/* complex */
union dcomplex {
double complex z;
double a[2];
};
union fcomplex {
float complex z;
float a[2];
};
union lcomplex {
long double complex z;
long double a[2];
};
// FIXME: move to complex.h ?
#define creal(z) ((double)(z))
#define crealf(z) ((float)(z))
#define creall(z) ((long double)(z))
#define cimag(z) ((union dcomplex){(z)}.a[1])
#define cimagf(z) ((union fcomplex){(z)}.a[1])
#define cimagl(z) ((union lcomplex){(z)}.a[1])
/* x + y*I is not supported properly by gcc */
#define cpack(x,y) ((union dcomplex){.a={(x),(y)}}.z)
#define cpackf(x,y) ((union fcomplex){.a={(x),(y)}}.z)
#define cpackl(x,y) ((union lcomplex){.a={(x),(y)}}.z)
#endif
src/internal/longdbl.h 0 → 100644
浏览文件 @ b69f695a
#ifndef _LDHACK_H
#define _LDHACK_H
#include <float.h>
#include <stdint.h>
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
union ldshape {
long double value;
struct {
uint64_t m;
uint16_t exp:15;
uint16_t sign:1;
uint16_t pad;
} bits;
};
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
union ldshape {
long double value;
struct {
uint64_t mlo;
uint64_t mhi:48;
uint16_t exp:15;
uint16_t sign:1;
} bits;
};
#else
#error Unsupported long double representation
#endif
// FIXME: hacks to make freebsd+openbsd long double code happy
// union and macros for freebsd
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
union IEEEl2bits {
long double e;
struct {
uint32_t manl:32;
uint32_t manh:32;
uint32_t exp:15;
uint32_t sign:1;
uint32_t pad:16;
} bits;
struct {
uint64_t man:64;
uint32_t expsign:16;
uint32_t pad:16;
} xbits;
};
#define LDBL_MANL_SIZE 32
#define LDBL_MANH_SIZE 32
#define LDBL_NBIT (1ull << LDBL_MANH_SIZE-1)
#undef LDBL_IMPLICIT_NBIT
#define mask_nbit_l(u) ((u).bits.manh &= ~LDBL_NBIT)
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
/*
// ld128 float.h
//#define LDBL_MAX 1.189731495357231765085759326628007016E+4932L
#define LDBL_MAX 0x1.ffffffffffffffffffffffffffffp+16383
#define LDBL_MAX_EXP 16384
#define LDBL_HAS_INFINITY 1
//#define LDBL_MIN 3.362103143112093506262677817321752603E-4932L
#define LDBL_MIN 0x1p-16382
#define LDBL_HAS_QUIET_NAN 1
#define LDBL_HAS_DENORM 1
//#define LDBL_EPSILON 1.925929944387235853055977942584927319E-34L
#define LDBL_EPSILON 0x1p-112
#define LDBL_MANT_DIG 113
#define LDBL_MIN_EXP (-16381)
#define LDBL_MAX_10_EXP 4932
#define LDBL_DENORM_MIN 0x0.0000000000000000000000000001p-16381
#define LDBL_MIN_10_EXP (-4931)
#define LDBL_DIG 33
*/
union IEEEl2bits {
long double e;
struct {
uint64_t manl:64;
uint64_t manh:48;
uint32_t exp:15;
uint32_t sign:1;
} bits;
struct {
uint64_t unused0:64;
uint64_t unused1:48;
uint32_t expsign:16;
} xbits;
};
#define LDBL_MANL_SIZE 64
#define LDBL_MANH_SIZE 48
#define LDBL_NBIT (1ull << LDBL_MANH_SIZE)
#define LDBL_IMPLICIT_NBIT 1
#define mask_nbit_l(u)
#endif
// macros for openbsd
#define GET_LDOUBLE_WORDS(se,mh,ml, f) do{ \
union IEEEl2bits u; \
u.e = (f); \
(se) = u.xbits.expsign; \
(mh) = u.bits.manh; \
(ml) = u.bits.manl; \
}while(0)
#define SET_LDOUBLE_WORDS(f, se,mh,ml) do{ \
union IEEEl2bits u; \
u.xbits.expsign = (se); \
u.bits.manh = (mh); \
u.bits.manl = (ml); \
(f) = u.e; \
}while(0)
#define GET_LDOUBLE_EXP(se, f) do{ \
union IEEEl2bits u; \
u.e = (f); \
(se) = u.xbits.expsign; \
}while(0)
#define SET_LDOUBLE_EXP(f, se) do{ \
union IEEEl2bits u; \
u.e = (f); \
u.xbits.expsign = (se); \
(f) = u.e; \
}while(0)
#endif
src/math/k_cos.c → src/math/__cos.c
浏览文件 @ b69f695a
/* @(#)k_cos.c 1.3 95/01/18 */
/* origin: FreeBSD /usr/src/lib/msun/src/k_cos.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* __kernel_cos( x, y )
* __cos( x, y )
* kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
* Input y is the tail of x.
*
* Algorithm
* 1. Since cos(-x) = cos(x), we need only to consider positive x.
......@@ -25,29 +23,32 @@
* 4 14
* cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
* where the remez error is
*
*
* | 2 4 6 8 10 12 14 | -58
* |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
* | |
*
* 4 6 8 10 12 14
* | |
*
* 4 6 8 10 12 14
* 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
* cos(x) = 1 - x*x/2 + r
* since cos(x+y) ~ cos(x) - sin(x)*y
* cos(x) ~ 1 - x*x/2 + r
* since cos(x+y) ~ cos(x) - sin(x)*y
* ~ cos(x) - x*y,
* a correction term is necessary in cos(x) and hence
* cos(x+y) = 1 - (x*x/2 - (r - x*y))
* For better accuracy when x > 0.3, let qx = |x|/4 with
* the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
* Then
* cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
* Note that 1-qx and (x*x/2-qx) is EXACT here, and the
* magnitude of the latter is at least a quarter of x*x/2,
* thus, reducing the rounding error in the subtraction.
* For better accuracy, rearrange to
* cos(x+y) ~ w + (tmp + (r-x*y))
* where w = 1 - x*x/2 and tmp is a tiny correction term
* (1 - x*x/2 == w + tmp exactly in infinite precision).
* The exactness of w + tmp in infinite precision depends on w
* and tmp having the same precision as x. If they have extra
* precision due to compiler bugs, then the extra precision is
* only good provided it is retained in all terms of the final
* expression for cos(). Retention happens in all cases tested
* under FreeBSD, so don't pessimize things by forcibly clipping
* any extra precision in w.
*/
#include <math.h>
#include "math_private.h"
#include "libm.h"
static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
......@@ -58,28 +59,14 @@ C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
double
__kernel_cos(double x, double y)
double __cos(double x, double y)
{
double a,hz,z,r,qx;
int32_t ix;
GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff; /* ix = |x|'s high word*/
if(ix<0x3e400000) { /* if x < 2**27 */
if(((int)x)==0) return one; /* generate inexact */
}
z = x*x;
r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
if(ix < 0x3FD33333) /* if |x| < 0.3 */
return one - (0.5*z - (z*r - x*y));
else {
if(ix > 0x3fe90000) { /* x > 0.78125 */
qx = 0.28125;
} else {
INSERT_WORDS(qx,ix-0x00200000,0); /* x/4 */
}
hz = 0.5*z-qx;
a = one-qx;
return a - (hz - (z*r-x*y));
}
double hz,z,r,w;
z = x*x;
w = z*z;
r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
hz = 0.5*z;
w = one-hz;
return w + (((one-w)-hz) + (z*r-x*y));
}
src/math/s_tanf.c → src/math/__cosdf.c
浏览文件 @ b69f695a
/* s_tanf.c -- float version of s_tan.c.
/* origin: FreeBSD /usr/src/lib/msun/src/k_cosf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
* Debugged and optimized by Bruce D. Evans.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
......@@ -13,28 +14,23 @@
* ====================================================
*/
#include <math.h>
#include "math_private.h"
float
tanf(float x)
{
float y[2],z=0.0;
int32_t n, ix;
GET_FLOAT_WORD(ix,x);
#include "libm.h"
/* |x| ~< pi/4 */
ix &= 0x7fffffff;
if(ix <= 0x3f490fda) return __kernel_tanf(x,z,1);
/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
static const double
one = 1.0,
C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */
C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
C3 = 0x199342e0ee5069.0p-68; /* 0.0000243904487962774090654 */
/* tan(Inf or NaN) is NaN */
else if (ix>=0x7f800000) return x-x; /* NaN */
float __cosdf(double x)
{
double r, w, z;
/* argument reduction needed */
else {
n = __ieee754_rem_pio2f(x,y);
return __kernel_tanf(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even
-1 -- n odd */
}
/* Try to optimize for parallel evaluation as in __tandf.c. */
z = x*x;
w = z*z;
r = C2+z*C3;
return ((one+z*C0) + w*C1) + (w*z)*r;
}
src/math/__cosl.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/ld80/k_cosl.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
* Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "libm.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/*
* ld80 version of __cos.c. See __cos.c for most comments.
*/
/*
* Domain [-0.7854, 0.7854], range ~[-2.43e-23, 2.425e-23]:
* |cos(x) - c(x)| < 2**-75.1
*
* The coefficients of c(x) were generated by a pari-gp script using
* a Remez algorithm that searches for the best higher coefficients
* after rounding leading coefficients to a specified precision.
*
* Simpler methods like Chebyshev or basic Remez barely suffice for
* cos() in 64-bit precision, because we want the coefficient of x^2
* to be precisely -0.5 so that multiplying by it is exact, and plain
* rounding of the coefficients of a good polynomial approximation only
* gives this up to about 64-bit precision. Plain rounding also gives
* a mediocre approximation for the coefficient of x^4, but a rounding
* error of 0.5 ulps for this coefficient would only contribute ~0.01
* ulps to the final error, so this is unimportant. Rounding errors in
* higher coefficients are even less important.
*
* In fact, coefficients above the x^4 one only need to have 53-bit
* precision, and this is more efficient. We get this optimization
* almost for free from the complications needed to search for the best
* higher coefficients.
*/
static const double one = 1.0;
// FIXME
/* Long double constants are slow on these arches, and broken on i386. */
static const volatile double
C1hi = 0.041666666666666664, /* 0x15555555555555.0p-57 */
C1lo = 2.2598839032744733e-18; /* 0x14d80000000000.0p-111 */
#define C1 ((long double)C1hi + C1lo)
#if 0
static const long double
C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */
#endif
static const double
C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */
C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */
C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */
C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */
C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */
C7 = 4.7383039476436467e-14; /* 0x1aac9d9af5c43e.0p-97 */
long double __cosl(long double x, long double y)
{
long double hz,z,r,w;
z = x*x;
r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))));
hz = 0.5*z;
w = one-hz;
return w + (((one-w)-hz) + (z*r-x*y));
}
#endif
src/math/__expo2.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */
/*-
* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "libm.h"
/*
* We use exp(x) = exp(x - kln2) * 2**k,
* k is carefully chosen to minimize |exp(kln2) - 2**k|
*/
static const uint32_t k = 1799;
static const double kln2 = 1246.97177782734161156;
/* exp(x)/2 when x is huge */
double __expo2(double x)
{
double scale;
int n;
/*
* efficient scalbn(y, k-1):
* 2**(k-1) cannot be represented
* so we use that k-1 is even and scale in two steps
*/
n = (k - 1)/2;
INSERT_WORDS(scale, (0x3ff + n) << 20, 0);
return exp(x - kln2) * scale * scale;
}
src/math/__expo2f.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/k_expf.c */
/*-
* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "libm.h"
/*
* We use exp(x) = exp(x - kln2) * 2**k,
* k is carefully chosen to minimize |exp(kln2) - 2**k|
*/
static const uint32_t k = 235;
static const float kln2 = 162.88958740f;
/* expf(x)/2 when x is huge */
float __expo2f(float x)
{
float scale;
int n;
/*
* efficient scalbnf(y, k-1):
* 2**(k-1) cannot be represented
* so we use that k-1 is even and scale in two steps
*/
n = (k - 1)/2;
SET_FLOAT_WORD(scale, (0x7f + n) << 23);
return expf(x - kln2) * scale * scale;
}
src/math/__fpclassify.c
浏览文件 @ b69f695a
#include <stdint.h>
#include <math.h>
#include "libm.h"
int __fpclassify(double __x)
int __fpclassify(double x)
{
union {
double __d;
__uint64_t __i;
} __y = { __x };
int __ee = __y.__i>>52 & 0x7ff;
if (!__ee) return __y.__i<<1 ? FP_SUBNORMAL : FP_ZERO;
if (__ee==0x7ff) return __y.__i<<12 ? FP_NAN : FP_INFINITE;
union dshape u = { x };
int e = u.bits>>52 & 0x7ff;
if (!e) return u.bits<<1 ? FP_SUBNORMAL : FP_ZERO;
if (e==0x7ff) return u.bits<<12 ? FP_NAN : FP_INFINITE;
return FP_NORMAL;
}
src/math/__fpclassifyf.c
浏览文件 @ b69f695a
#include <stdint.h>
#include <math.h>
#include "libm.h"
int __fpclassifyf(float __x)
int __fpclassifyf(float x)
{
union {
float __f;
__uint32_t __i;
} __y = { __x };
int __ee = __y.__i>>23 & 0xff;
if (!__ee) return __y.__i<<1 ? FP_SUBNORMAL : FP_ZERO;
if (__ee==0xff) return __y.__i<<9 ? FP_NAN : FP_INFINITE;
union fshape u = { x };
int e = u.bits>>23 & 0xff;
if (!e) return u.bits<<1 ? FP_SUBNORMAL : FP_ZERO;
if (e==0xff) return u.bits<<9 ? FP_NAN : FP_INFINITE;
return FP_NORMAL;
}
src/math/__fpclassifyl.c
浏览文件 @ b69f695a
#include <stdint.h>
#include <math.h>
#include "libm.h"
/* FIXME: move this to arch-specific file */
int __fpclassifyl(long double __x)
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
int __fpclassifyl(long double x)
{
union ldshape u = { x };
int e = u.bits.exp;
if (!e)
return u.bits.m ? FP_SUBNORMAL : FP_ZERO;
if (e == 0x7fff)
return u.bits.m & (uint64_t)-1>>1 ? FP_NAN : FP_INFINITE;
return u.bits.m & (uint64_t)1<<63 ? FP_NORMAL : FP_NAN;
}
#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
int __fpclassifyl(long double x)
{
union {
long double __ld;
__uint16_t __hw[5];
int64_t __m;
} __y = { __x };
int __ee = __y.__hw[4]&0x7fff;
if (!__ee) return __y.__m ? FP_SUBNORMAL : FP_ZERO;
if (__ee==0x7fff) return __y.__m ? FP_NAN : FP_INFINITE;
return __y.__m < 0 ? FP_NORMAL : FP_NAN;
union ldshape u = { x };
int e = u.bits.exp;
if (!e)
return u.bits.mlo | u.bits.mhi ? FP_SUBNORMAL : FP_ZERO;
if (e == 0x7fff)
return u.bits.mlo | u.bits.mhi ? FP_NAN : FP_INFINITE;
return FP_NORMAL;
}
#endif
src/math/__invtrigl.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/ld80/invtrig.c */
/*-
* Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "__invtrigl.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/*
* asinl() and acosl()
*/
const long double
pS0 = 1.66666666666666666631e-01L,
pS1 = -4.16313987993683104320e-01L,
pS2 = 3.69068046323246813704e-01L,
pS3 = -1.36213932016738603108e-01L,
pS4 = 1.78324189708471965733e-02L,
pS5 = -2.19216428382605211588e-04L,
pS6 = -7.10526623669075243183e-06L,
qS1 = -2.94788392796209867269e+00L,
qS2 = 3.27309890266528636716e+00L,
qS3 = -1.68285799854822427013e+00L,
qS4 = 3.90699412641738801874e-01L,
qS5 = -3.14365703596053263322e-02L;
/*
* atanl()
*/
const long double atanhi[] = {
4.63647609000806116202e-01L,
7.85398163397448309628e-01L,
9.82793723247329067960e-01L,
1.57079632679489661926e+00L,
};
const long double atanlo[] = {
1.18469937025062860669e-20L,
-1.25413940316708300586e-20L,
2.55232234165405176172e-20L,
-2.50827880633416601173e-20L,
};
const long double aT[] = {
3.33333333333333333017e-01L,
-1.99999999999999632011e-01L,
1.42857142857046531280e-01L,
-1.11111111100562372733e-01L,
9.09090902935647302252e-02L,
-7.69230552476207730353e-02L,
6.66661718042406260546e-02L,
-5.88158892835030888692e-02L,
5.25499891539726639379e-02L,
-4.70119845393155721494e-02L,
4.03539201366454414072e-02L,
-2.91303858419364158725e-02L,
1.24822046299269234080e-02L,
};
const long double pi_lo = -5.01655761266833202345e-20L;
#endif
src/math/__invtrigl.h 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/ld80/invtrig.h */
/*-
* Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "libm.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
#define BIAS (LDBL_MAX_EXP - 1)
#define MANH_SIZE LDBL_MANH_SIZE
/* Approximation thresholds. */
#define ASIN_LINEAR (BIAS - 32) /* 2**-32 */
#define ACOS_CONST (BIAS - 65) /* 2**-65 */
#define ATAN_CONST (BIAS + 65) /* 2**65 */
#define ATAN_LINEAR (BIAS - 32) /* 2**-32 */
/* 0.95 */
#define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)
/* Constants shared by the long double inverse trig functions. */
#define pS0 __pS0
#define pS1 __pS1
#define pS2 __pS2
#define pS3 __pS3
#define pS4 __pS4
#define pS5 __pS5
#define pS6 __pS6
#define qS1 __qS1
#define qS2 __qS2
#define qS3 __qS3
#define qS4 __qS4
#define qS5 __qS5
#define atanhi __atanhi
#define atanlo __atanlo
#define aT __aT
#define pi_lo __pi_lo
#define pio2_hi atanhi[3]
#define pio2_lo atanlo[3]
#define pio4_hi atanhi[1]
#ifdef STRUCT_DECLS
typedef struct longdouble {
uint64_t mant;
uint16_t expsign;
} LONGDOUBLE;
#else
typedef long double LONGDOUBLE;
#endif
extern const LONGDOUBLE pS0, pS1, pS2, pS3, pS4, pS5, pS6;
extern const LONGDOUBLE qS1, qS2, qS3, qS4, qS5;
extern const LONGDOUBLE atanhi[], atanlo[], aT[];
extern const LONGDOUBLE pi_lo;
#ifndef STRUCT_DECLS
static inline long double
P(long double x)
{
return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \
(pS4 + x * (pS5 + x * pS6)))))));
}
static inline long double
Q(long double x)
{
return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5)))));
}
static inline long double
T_even(long double x)
{
return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \
(aT[8] + x * (aT[10] + x * aT[12]))))));
}
static inline long double
T_odd(long double x)
{
return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \
(aT[9] + x * aT[11])))));
}
#endif
#endif
src/math/__log1p.h 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/k_log.h */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* __log1p(f):
* Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)].
*
* The following describes the overall strategy for computing
* logarithms in base e. The argument reduction and adding the final
* term of the polynomial are done by the caller for increased accuracy
* when different bases are used.
*
* Method :
* 1. Argument Reduction: find k and f such that
* x = 2^k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* 2. Approximation of log(1+f).
* Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* = 2s + s*R
* We use a special Reme algorithm on [0,0.1716] to generate
* a polynomial of degree 14 to approximate R The maximum error
* of this polynomial approximation is bounded by 2**-58.45. In
* other words,
* 2 4 6 8 10 12 14
* R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
* (the values of Lg1 to Lg7 are listed in the program)
* and
* | 2 14 | -58.45
* | Lg1*s +...+Lg7*s - R(z) | <= 2
* | |
* Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
* In order to guarantee error in log below 1ulp, we compute log
* by
* log(1+f) = f - s*(f - R) (if f is not too large)
* log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
*
* 3. Finally, log(x) = k*ln2 + log(1+f).
* = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
* Here ln2 is split into two floating point number:
* ln2_hi + ln2_lo,
* where n*ln2_hi is always exact for |n| < 2000.
*
* Special cases:
* log(x) is NaN with signal if x < 0 (including -INF) ;
* log(+INF) is +INF; log(0) is -INF with signal;
* log(NaN) is that NaN with no signal.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
static const double
Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
/*
* We always inline __log1p(), since doing so produces a
* substantial performance improvement (~40% on amd64).
*/
static inline double __log1p(double f)
{
double hfsq,s,z,R,w,t1,t2;
s = f/(2.0+f);
z = s*s;
w = z*z;
t1= w*(Lg2+w*(Lg4+w*Lg6));
t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
R = t2+t1;
hfsq = 0.5*f*f;
return s*(hfsq+R);
}
src/math/s_fabs.c → src/math/__log1pf.h
浏览文件 @ b69f695a
/* @(#)s_fabs.c 5.1 93/09/24 */
/* origin: FreeBSD /usr/src/lib/msun/src/k_logf.h */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
......@@ -9,19 +9,27 @@
* is preserved.
* ====================================================
*/
/*
* fabs(x) returns the absolute value of x.
* See comments in __log1p.h.
*/
#include <math.h>
#include "math_private.h"
/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
static const float
Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
double
fabs(double x)
static inline float __log1pf(float f)
{
uint32_t high;
GET_HIGH_WORD(high,x);
SET_HIGH_WORD(x,high&0x7fffffff);
return x;
float hfsq,s,z,R,w,t1,t2;
s = f/((float)2.0+f);
z = s*s;
w = z*z;
t1 = w*(Lg2+w*Lg4);
t2 = z*(Lg1+w*Lg3);
R = t2+t1;
hfsq = (float)0.5*f*f;
return s*(hfsq+R);
}
src/math/__polevll.c 0 → 100644
浏览文件 @ b69f695a
/* origin: OpenBSD /usr/src/lib/libm/src/polevll.c */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/*
* Evaluate polynomial
*
*
* SYNOPSIS:
*
* int N;
* long double x, y, coef[N+1], polevl[];
*
* y = polevll( x, coef, N );
*
*
* DESCRIPTION:
*
* Evaluates polynomial of degree N:
*
* 2 N
* y = C + C x + C x +...+ C x
* 0 1 2 N
*
* Coefficients are stored in reverse order:
*
* coef[0] = C , ..., coef[N] = C .
* N 0
*
* The function p1evll() assumes that coef[N] = 1.0 and is
* omitted from the array. Its calling arguments are
* otherwise the same as polevll().
*
*
* SPEED:
*
* In the interest of speed, there are no checks for out
* of bounds arithmetic. This routine is used by most of
* the functions in the library. Depending on available
* equipment features, the user may wish to rewrite the
* program in microcode or assembly language.
*
*/
#include "libm.h"
/*
* Polynomial evaluator:
* P[0] x^n + P[1] x^(n-1) + ... + P[n]
*/
long double __polevll(long double x, long double *P, int n)
{
long double y;
y = *P++;
do {
y = y * x + *P++;
} while (--n);
return y;
}
/*
* Polynomial evaluator:
* x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n]
*/
long double __p1evll(long double x, long double *P, int n)
{
long double y;
n -= 1;
y = x + *P++;
do {
y = y * x + *P++;
} while (--n);
return y;
}
src/math/__rem_pio2.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
* Optimized by Bruce D. Evans.
*/
/* __rem_pio2(x,y)
*
* return the remainder of x rem pi/2 in y[0]+y[1]
* use __rem_pio2_large() for large x
*/
#include "libm.h"
/*
* invpio2: 53 bits of 2/pi
* pio2_1: first 33 bit of pi/2
* pio2_1t: pi/2 - pio2_1
* pio2_2: second 33 bit of pi/2
* pio2_2t: pi/2 - (pio2_1+pio2_2)
* pio2_3: third 33 bit of pi/2
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const double
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
/* caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
int __rem_pio2(double x, double *y)
{
double z,w,t,r,fn;
double tx[3],ty[2];
int32_t e0,i,j,nx,n,ix,hx;
uint32_t low;
GET_HIGH_WORD(hx,x);
ix = hx & 0x7fffffff;
if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */
if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */
goto medium; /* cancellation -- use medium case */
if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */
if (hx > 0) {
z = x - pio2_1; /* one round good to 85 bits */
y[0] = z - pio2_1t;
y[1] = (z-y[0]) - pio2_1t;
return 1;
} else {
z = x + pio2_1;
y[0] = z + pio2_1t;
y[1] = (z-y[0]) + pio2_1t;
return -1;
}
} else {
if (hx > 0) {
z = x - 2*pio2_1;
y[0] = z - 2*pio2_1t;
y[1] = (z-y[0]) - 2*pio2_1t;
return 2;
} else {
z = x + 2*pio2_1;
y[0] = z + 2*pio2_1t;
y[1] = (z-y[0]) + 2*pio2_1t;
return -2;
}
}
}
if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */
if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */
if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */
goto medium;
if (hx > 0) {
z = x - 3*pio2_1;
y[0] = z - 3*pio2_1t;
y[1] = (z-y[0]) - 3*pio2_1t;
return 3;
} else {
z = x + 3*pio2_1;
y[0] = z + 3*pio2_1t;
y[1] = (z-y[0]) + 3*pio2_1t;
return -3;
}
} else {
if (ix == 0x401921fb) /* |x| ~= 4pi/2 */
goto medium;
if (hx > 0) {
z = x - 4*pio2_1;
y[0] = z - 4*pio2_1t;
y[1] = (z-y[0]) - 4*pio2_1t;
return 4;
} else {
z = x + 4*pio2_1;
y[0] = z + 4*pio2_1t;
y[1] = (z-y[0]) + 4*pio2_1t;
return -4;
}
}
}
if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */
uint32_t high;
medium:
/* Use a specialized rint() to get fn. Assume round-to-nearest. */
STRICT_ASSIGN(double, fn, x*invpio2 + 0x1.8p52);
fn = fn - 0x1.8p52;
// FIXME
#ifdef HAVE_EFFICIENT_IRINT
n = irint(fn);
#else
n = (int32_t)fn;
#endif
r = x - fn*pio2_1;
w = fn*pio2_1t; /* 1st round, good to 85 bits */
j = ix>>20;
y[0] = r - w;
GET_HIGH_WORD(high,y[0]);
i = j - ((high>>20)&0x7ff);
if (i > 16) { /* 2nd round, good to 118 bits */
t = r;
w = fn*pio2_2;
r = t - w;
w = fn*pio2_2t - ((t-r)-w);
y[0] = r - w;
GET_HIGH_WORD(high,y[0]);
i = j - ((high>>20)&0x7ff);
if (i > 49) { /* 3rd round, good to 151 bits, covers all cases */
t = r;
w = fn*pio2_3;
r = t - w;
w = fn*pio2_3t - ((t-r)-w);
y[0] = r - w;
}
}
y[1] = (r-y[0]) - w;
return n;
}
/*
* all other (large) arguments
*/
if (ix >= 0x7ff00000) { /* x is inf or NaN */
y[0] = y[1] = x - x;
return 0;
}
/* set z = scalbn(|x|,ilogb(x)-23) */
GET_LOW_WORD(low,x);
e0 = (ix>>20) - 1046; /* e0 = ilogb(z)-23; */
INSERT_WORDS(z, ix - ((int32_t)(e0<<20)), low);
for (i=0; i<2; i++) {
tx[i] = (double)((int32_t)(z));
z = (z-tx[i])*two24;
}
tx[2] = z;
nx = 3;
while (tx[nx-1] == zero) nx--; /* skip zero term */
n = __rem_pio2_large(tx,ty,e0,nx,1);
if (hx < 0) {
y[0] = -ty[0];
y[1] = -ty[1];
return -n;
}
y[0] = ty[0];
y[1] = ty[1];
return n;
}
src/math/__rem_pio2_large.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/k_rem_pio2.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* __rem_pio2_large(x,y,e0,nx,prec)
* double x[],y[]; int e0,nx,prec;
*
* __rem_pio2_large return the last three digits of N with
* y = x - N*pi/2
* so that |y| < pi/2.
*
* The method is to compute the integer (mod 8) and fraction parts of
* (2/pi)*x without doing the full multiplication. In general we
* skip the part of the product that are known to be a huge integer (
* more accurately, = 0 mod 8 ). Thus the number of operations are
* independent of the exponent of the input.
*
* (2/pi) is represented by an array of 24-bit integers in ipio2[].
*
* Input parameters:
* x[] The input value (must be positive) is broken into nx
* pieces of 24-bit integers in double precision format.
* x[i] will be the i-th 24 bit of x. The scaled exponent
* of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
* match x's up to 24 bits.
*
* Example of breaking a double positive z into x[0]+x[1]+x[2]:
* e0 = ilogb(z)-23
* z = scalbn(z,-e0)
* for i = 0,1,2
* x[i] = floor(z)
* z = (z-x[i])*2**24
*
*
* y[] ouput result in an array of double precision numbers.
* The dimension of y[] is:
* 24-bit precision 1
* 53-bit precision 2
* 64-bit precision 2
* 113-bit precision 3
* The actual value is the sum of them. Thus for 113-bit
* precison, one may have to do something like:
*
* long double t,w,r_head, r_tail;
* t = (long double)y[2] + (long double)y[1];
* w = (long double)y[0];
* r_head = t+w;
* r_tail = w - (r_head - t);
*
* e0 The exponent of x[0]. Must be <= 16360 or you need to
* expand the ipio2 table.
*
* nx dimension of x[]
*
* prec an integer indicating the precision:
* 0 24 bits (single)
* 1 53 bits (double)
* 2 64 bits (extended)
* 3 113 bits (quad)
*
* External function:
* double scalbn(), floor();
*
*
* Here is the description of some local variables:
*
* jk jk+1 is the initial number of terms of ipio2[] needed
* in the computation. The minimum and recommended value
* for jk is 3,4,4,6 for single, double, extended, and quad.
* jk+1 must be 2 larger than you might expect so that our
* recomputation test works. (Up to 24 bits in the integer
* part (the 24 bits of it that we compute) and 23 bits in
* the fraction part may be lost to cancelation before we
* recompute.)
*
* jz local integer variable indicating the number of
* terms of ipio2[] used.
*
* jx nx - 1
*
* jv index for pointing to the suitable ipio2[] for the
* computation. In general, we want
* ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
* is an integer. Thus
* e0-3-24*jv >= 0 or (e0-3)/24 >= jv
* Hence jv = max(0,(e0-3)/24).
*
* jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
*
* q[] double array with integral value, representing the
* 24-bits chunk of the product of x and 2/pi.
*
* q0 the corresponding exponent of q[0]. Note that the
* exponent for q[i] would be q0-24*i.
*
* PIo2[] double precision array, obtained by cutting pi/2
* into 24 bits chunks.
*
* f[] ipio2[] in floating point
*
* iq[] integer array by breaking up q[] in 24-bits chunk.
*
* fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
*
* ih integer. If >0 it indicates q[] is >= 0.5, hence
* it also indicates the *sign* of the result.
*
*/
/*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
#include "libm.h"
static const int init_jk[] = {3,4,4,6}; /* initial value for jk */
/*
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
*
* integer array, contains the (24*i)-th to (24*i+23)-th
* bit of 2/pi after binary point. The corresponding
* floating value is
*
* ipio2[i] * 2^(-24(i+1)).
*
* NB: This table must have at least (e0-3)/24 + jk terms.
* For quad precision (e0 <= 16360, jk = 6), this is 686.
*/
static const int32_t ipio2[] = {
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
#if LDBL_MAX_EXP > 1024
0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
#endif
};
static const double PIo2[] = {
1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
};
static const double
zero = 0.0,
one = 1.0,
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec)
{
int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
double z,fw,f[20],fq[20],q[20];
/* initialize jk*/
jk = init_jk[prec];
jp = jk;
/* determine jx,jv,q0, note that 3>q0 */
jx = nx-1;
jv = (e0-3)/24; if(jv<0) jv=0;
q0 = e0-24*(jv+1);
/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
j = jv-jx; m = jx+jk;
for (i=0; i<=m; i++,j++)
f[i] = j<0 ? zero : (double)ipio2[j];
/* compute q[0],q[1],...q[jk] */
for (i=0; i<=jk; i++) {
for (j=0,fw=0.0; j<=jx; j++)
fw += x[j]*f[jx+i-j];
q[i] = fw;
}
jz = jk;
recompute:
/* distill q[] into iq[] reversingly */
for (i=0,j=jz,z=q[jz]; j>0; i++,j--) {
fw = (double)((int32_t)(twon24* z));
iq[i] = (int32_t)(z-two24*fw);
z = q[j-1]+fw;
}
/* compute n */
z = scalbn(z,q0); /* actual value of z */
z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
n = (int32_t)z;
z -= (double)n;
ih = 0;
if (q0 > 0) { /* need iq[jz-1] to determine n */
i = iq[jz-1]>>(24-q0); n += i;
iq[jz-1] -= i<<(24-q0);
ih = iq[jz-1]>>(23-q0);
}
else if (q0 == 0) ih = iq[jz-1]>>23;
else if (z >= 0.5) ih = 2;
if (ih > 0) { /* q > 0.5 */
n += 1; carry = 0;
for (i=0; i<jz; i++) { /* compute 1-q */
j = iq[i];
if (carry == 0) {
if (j != 0) {
carry = 1;
iq[i] = 0x1000000- j;
}
} else
iq[i] = 0xffffff - j;
}
if (q0 > 0) { /* rare case: chance is 1 in 12 */
switch(q0) {
case 1:
iq[jz-1] &= 0x7fffff; break;
case 2:
iq[jz-1] &= 0x3fffff; break;
}
}
if (ih == 2) {
z = one - z;
if (carry != 0)
z -= scalbn(one,q0);
}
}
/* check if recomputation is needed */
if (z == zero) {
j = 0;
for (i=jz-1; i>=jk; i--) j |= iq[i];
if (j == 0) { /* need recomputation */
for (k=1; iq[jk-k]==0; k++); /* k = no. of terms needed */
for (i=jz+1; i<=jz+k; i++) { /* add q[jz+1] to q[jz+k] */
f[jx+i] = (double)ipio2[jv+i];
for (j=0,fw=0.0; j<=jx; j++)
fw += x[j]*f[jx+i-j];
q[i] = fw;
}
jz += k;
goto recompute;
}
}
/* chop off zero terms */
if (z == 0.0) {
jz -= 1;
q0 -= 24;
while (iq[jz] == 0) {
jz--;
q0 -= 24;
}
} else { /* break z into 24-bit if necessary */
z = scalbn(z,-q0);
if (z >= two24) {
fw = (double)((int32_t)(twon24*z));
iq[jz] = (int32_t)(z-two24*fw);
jz += 1;
q0 += 24;
iq[jz] = (int32_t)fw;
} else
iq[jz] = (int32_t)z;
}
/* convert integer "bit" chunk to floating-point value */
fw = scalbn(one,q0);
for (i=jz; i>=0; i--) {
q[i] = fw*(double)iq[i];
fw *= twon24;
}
/* compute PIo2[0,...,jp]*q[jz,...,0] */
for(i=jz; i>=0; i--) {
for (fw=0.0,k=0; k<=jp && k<=jz-i; k++)
fw += PIo2[k]*q[i+k];
fq[jz-i] = fw;
}
/* compress fq[] into y[] */
switch(prec) {
case 0:
fw = 0.0;
for (i=jz; i>=0; i--)
fw += fq[i];
y[0] = ih==0 ? fw : -fw;
break;
case 1:
case 2:
fw = 0.0;
for (i=jz; i>=0; i--)
fw += fq[i];
STRICT_ASSIGN(double,fw,fw);
y[0] = ih==0 ? fw : -fw;
fw = fq[0]-fw;
for (i=1; i<=jz; i++)
fw += fq[i];
y[1] = ih==0 ? fw : -fw;
break;
case 3: /* painful */
for (i=jz; i>0; i--) {
fw = fq[i-1]+fq[i];
fq[i] += fq[i-1]-fw;
fq[i-1] = fw;
}
for (i=jz; i>1; i--) {
fw = fq[i-1]+fq[i];
fq[i] += fq[i-1]-fw;
fq[i-1] = fw;
}
for (fw=0.0,i=jz; i>=2; i--)
fw += fq[i];
if (ih==0) {
y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
} else {
y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
}
}
return n&7;
}
src/math/__rem_pio2f.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2f.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
* Debugged and optimized by Bruce D. Evans.
*/
/*
* ====================================================
* 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.
* ====================================================
*/
/* __rem_pio2f(x,y)
*
* return the remainder of x rem pi/2 in *y
* use double precision for everything except passing x
* use __rem_pio2_large() for large x
*/
#include "libm.h"
/*
* invpio2: 53 bits of 2/pi
* pio2_1: first 33 bit of pi/2
* pio2_1t: pi/2 - pio2_1
*/
static const double
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */
pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
int __rem_pio2f(float x, double *y)
{
double w,r,fn;
double tx[1],ty[1];
float z;
int32_t e0,n,ix,hx;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
/* 33+53 bit pi is good enough for medium size */
if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */
/* Use a specialized rint() to get fn. Assume round-to-nearest. */
STRICT_ASSIGN(double, fn, x*invpio2 + 0x1.8p52);
fn = fn - 0x1.8p52;
// FIXME
#ifdef HAVE_EFFICIENT_IRINT
n = irint(fn);
#else
n = (int32_t)fn;
#endif
r = x - fn*pio2_1;
w = fn*pio2_1t;
*y = r - w;
return n;
}
/*
* all other (large) arguments
*/
if(ix>=0x7f800000) { /* x is inf or NaN */
*y = x-x;
return 0;
}
/* set z = scalbn(|x|,ilogb(|x|)-23) */
e0 = (ix>>23) - 150; /* e0 = ilogb(|x|)-23; */
SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
tx[0] = z;
n = __rem_pio2_large(tx,ty,e0,1,0);
if (hx < 0) {
*y = -ty[0];
return -n;
}
*y = ty[0];
return n;
}
src/math/__rem_pio2l.h 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/ld80/e_rem_pio2.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
* Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
* Optimized by Bruce D. Evans.
*/
#include "libm.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* ld80 version of __rem_pio2(x,y)
*
* return the remainder of x rem pi/2 in y[0]+y[1]
* use __rem_pio2_large() for large x
*/
#define BIAS (LDBL_MAX_EXP - 1)
/*
* invpio2: 64 bits of 2/pi
* pio2_1: first 39 bits of pi/2
* pio2_1t: pi/2 - pio2_1
* pio2_2: second 39 bits of pi/2
* pio2_2t: pi/2 - (pio2_1+pio2_2)
* pio2_3: third 39 bits of pi/2
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const double
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */
pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */
pio2_3 = 6.36831716351370313614e-25; /* 0x18a2e037074000.0p-133 */
// FIXME: this should be verified (maybe old gcc specific hack)
//#if defined(__amd64__) || defined(__i386__)
/* Long double constants are slow on these arches, and broken on i386. */
static const volatile double
invpio2hi = 6.3661977236758138e-01, /* 0x145f306dc9c883.0p-53 */
invpio2lo = -3.9356538861223811e-17, /* -0x16b00000000000.0p-107 */
pio2_1thi = -1.0746346554971943e-12, /* -0x12e7b9676733af.0p-92 */
pio2_1tlo = 8.8451028997905949e-29, /* 0x1c080000000000.0p-146 */
pio2_2thi = 6.3683171635109499e-25, /* 0x18a2e03707344a.0p-133 */
pio2_2tlo = 2.3183081793789774e-41, /* 0x10280000000000.0p-187 */
pio2_3thi = -2.7529965190440717e-37, /* -0x176b7ed8fbbacc.0p-174 */
pio2_3tlo = -4.2006647512740502e-54; /* -0x19c00000000000.0p-230 */
#define invpio2 ((long double)invpio2hi + invpio2lo)
#define pio2_1t ((long double)pio2_1thi + pio2_1tlo)
#define pio2_2t ((long double)pio2_2thi + pio2_2tlo)
#define pio2_3t ((long double)pio2_3thi + pio2_3tlo)
//#else
#if 0
static const long double
invpio2 = 6.36619772367581343076e-01L, /* 0xa2f9836e4e44152a.0p-64 */
pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */
pio2_2t = 6.36831716351095013979e-25L, /* 0xc51701b839a25205.0p-144 */
pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */
#endif
static inline int __rem_pio2l(long double x, long double *y)
{
union IEEEl2bits u,u1;
long double z,w,t,r,fn;
double tx[3],ty[2];
int e0,ex,i,j,nx,n;
int16_t expsign;
u.e = x;
expsign = u.xbits.expsign;
ex = expsign & 0x7fff;
if (ex < BIAS + 25 || (ex == BIAS + 25 && u.bits.manh < 0xc90fdaa2)) {
union IEEEl2bits u2;
int ex1;
/* |x| ~< 2^25*(pi/2), medium size */
/* Use a specialized rint() to get fn. Assume round-to-nearest. */
fn = x*invpio2 + 0x1.8p63;
fn = fn - 0x1.8p63;
// FIXME
//#ifdef HAVE_EFFICIENT_IRINT
// n = irint(fn);
//#else
n = fn;
//#endif
r = x-fn*pio2_1;
w = fn*pio2_1t; /* 1st round good to 102 bit */
j = ex;
y[0] = r-w;
u2.e = y[0];
ex1 = u2.xbits.expsign & 0x7fff;
i = j-ex1;
if (i > 22) { /* 2nd iteration needed, good to 141 */
t = r;
w = fn*pio2_2;
r = t-w;
w = fn*pio2_2t-((t-r)-w);
y[0] = r-w;
u2.e = y[0];
ex1 = u2.xbits.expsign & 0x7fff;
i = j-ex1;
if (i > 61) { /* 3rd iteration need, 180 bits acc */
t = r; /* will cover all possible cases */
w = fn*pio2_3;
r = t-w;
w = fn*pio2_3t-((t-r)-w);
y[0] = r-w;
}
}
y[1] = (r - y[0]) - w;
return n;
}
/*
* all other (large) arguments
*/
if (ex == 0x7fff) { /* x is inf or NaN */
y[0] = y[1] = x - x;
return 0;
}
/* set z = scalbn(|x|,ilogb(x)-23) */
u1.e = x;
e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */
u1.xbits.expsign = ex - e0;
z = u1.e;
for (i=0; i<2; i++) {
tx[i] = (double)(int32_t)z;
z = (z-tx[i])*two24;
}
tx[2] = z;
nx = 3;
while (tx[nx-1] == zero)
nx--; /* skip zero term */
n = __rem_pio2_large(tx,ty,e0,nx,2);
r = (long double)ty[0] + ty[1];
w = ty[1] - (r - ty[0]);
if (expsign < 0) {
y[0] = -r;
y[1] = -w;
return -n;
}
y[0] = r;
y[1] = w;
return n;
}
#endif
src/math/__signbit.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
// FIXME: macro
int __signbit(double x)
{
union {
double d;
uint64_t i;
} y = { x };
return y.i>>63;
}
src/math/__signbitf.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
// FIXME
int __signbitf(float x)
{
union {
float f;
uint32_t i;
} y = { x };
return y.i>>31;
}
src/math/__signbitl.c 0 → 100644
浏览文件 @ b69f695a
#include "libm.h"
// FIXME: should be a macro
#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
int __signbitl(long double x)
{
union ldshape u = {x};
return u.bits.sign;
}
#endif
src/math/k_sin.c → src/math/__sin.c
浏览文件 @ b69f695a
/* @(#)k_sin.c 1.3 95/01/18 */
/* origin: FreeBSD /usr/src/lib/msun/src/k_sin.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __kernel_sin( x, y, iy)
* kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
/* __sin( x, y, iy)
* kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
* Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
* Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
*
* Algorithm
* 1. Since sin(-x) = -sin(x), we need only to consider positive x.
* 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
* 1. Since sin(-x) = -sin(x), we need only to consider positive x.
* 2. Callers must return sin(-0) = -0 without calling here since our
* odd polynomial is not evaluated in a way that preserves -0.
* Callers may do the optimization sin(x) ~ x for tiny x.
* 3. sin(x) is approximated by a polynomial of degree 13 on
* [0,pi/4]
* 3 13
* sin(x) ~ x + S1*x + ... + S6*x
* where
*
*
* |sin(x) 2 4 6 8 10 12 | -58
* |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2
* | x |
*
* | x |
*
* 4. sin(x+y) = sin(x) + sin'(x')*y
* ~ sin(x) + (1-x*x/2)*y
* For better accuracy, let
* For better accuracy, let
* 3 2 2 2 2
* r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
* then 3 2
* sin(x) = x + (S1*x + (x *(r-y/2)+y))
*/
#include <math.h>
#include "math_private.h"
#include "libm.h"
static const double
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
......@@ -51,18 +50,16 @@ S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
double
__kernel_sin(double x, double y, int iy)
double __sin(double x, double y, int iy)
{
double z,r,v;
int32_t ix;
GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff; /* high word of x */
if(ix<0x3e400000) /* |x| < 2**-27 */
{if((int)x==0) return x;} /* generate inexact */
z = x*x;
v = z*x;
r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
if(iy==0) return x+v*(S1+z*r);
else return x-((z*(half*y-v*r)-y)-v*S1);
double z,r,v,w;
z = x*x;
w = z*z;
r = S2 + z*(S3 + z*S4) + z*w*(S5 + z*S6);
v = z*x;
if (iy == 0)
return x + v*(S1 + z*r);
else
return x - ((z*(half*y - v*r) - y) - v*S1);
}
src/math/s_fabsf.c → src/math/__sindf.c
浏览文件 @ b69f695a
/* s_fabsf.c -- float version of s_fabs.c.
/* origin: FreeBSD /usr/src/lib/msun/src/k_sinf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
* Optimized by Bruce D. Evans.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
......@@ -13,18 +14,23 @@
* ====================================================
*/
/*
* fabsf(x) returns the absolute value of x.
*/
#include "libm.h"
#include <math.h>
#include "math_private.h"
/* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */
static const double
S1 = -0x15555554cbac77.0p-55, /* -0.166666666416265235595 */
S2 = 0x111110896efbb2.0p-59, /* 0.0083333293858894631756 */
S3 = -0x1a00f9e2cae774.0p-65, /* -0.000198393348360966317347 */
S4 = 0x16cd878c3b46a7.0p-71; /* 0.0000027183114939898219064 */
float
fabsf(float x)
float __sindf(double x)
{
uint32_t ix;
GET_FLOAT_WORD(ix,x);
SET_FLOAT_WORD(x,ix&0x7fffffff);
return x;
double r, s, w, z;
/* Try to optimize for parallel evaluation as in __tandf.c. */
z = x*x;
w = z*z;
r = S3 + z*S4;
s = z*x;
return (x + s*(S1 + z*S2)) + s*w*r;
}
src/math/__sinl.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/ld80/k_sinl.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
* Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "libm.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/*
* ld80 version of __sin.c. See __sin.c for most comments.
*/
/*
* Domain [-0.7854, 0.7854], range ~[-1.89e-22, 1.915e-22]
* |sin(x)/x - s(x)| < 2**-72.1
*
* See __cosl.c for more details about the polynomial.
*/
static const double half = 0.5;
// FIXME
/* Long double constants are slow on these arches, and broken on i386. */
static const volatile double
S1hi = -0.16666666666666666, /* -0x15555555555555.0p-55 */
S1lo = -9.2563760475949941e-18; /* -0x15580000000000.0p-109 */
#define S1 ((long double)S1hi + S1lo)
#if 0
static const long double
S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */
#endif
static const double
S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */
S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */
S4 = 0.0000027557319223597490, /* 0x171de3a55560f7.0p-71 */
S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */
S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */
S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */
S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */
long double __sinl(long double x, long double y, int iy)
{
long double z,r,v;
z = x*x;
v = z*x;
r = S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8)))));
if (iy == 0)
return x+v*(S1+z*r);
return x-((z*(half*y-v*r)-y)-v*S1);
}
#endif
src/math/__tan.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/k_tan.c */
/*
* ====================================================
* Copyright 2004 Sun Microsystems, Inc. All Rights Reserved.
*
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __tan( x, y, k )
* kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
* Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.
*
* Algorithm
* 1. Since tan(-x) = -tan(x), we need only to consider positive x.
* 2. Callers must return tan(-0) = -0 without calling here since our
* odd polynomial is not evaluated in a way that preserves -0.
* Callers may do the optimization tan(x) ~ x for tiny x.
* 3. tan(x) is approximated by a odd polynomial of degree 27 on
* [0,0.67434]
* 3 27
* tan(x) ~ x + T1*x + ... + T13*x
* where
*
* |tan(x) 2 4 26 | -59.2
* |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2
* | x |
*
* Note: tan(x+y) = tan(x) + tan'(x)*y
* ~ tan(x) + (1+x*x)*y
* Therefore, for better accuracy in computing tan(x+y), let
* 3 2 2 2 2
* r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
* then
* 3 2
* tan(x+y) = x + (T1*x + (x *(r+y)+y))
*
* 4. For x in [0.67434,pi/4], let y = pi/4 - x, then
* tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
* = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
*/
#include "libm.h"
static const double T[] = {
3.33333333333334091986e-01, /* 3FD55555, 55555563 */
1.33333333333201242699e-01, /* 3FC11111, 1110FE7A */
5.39682539762260521377e-02, /* 3FABA1BA, 1BB341FE */
2.18694882948595424599e-02, /* 3F9664F4, 8406D637 */
8.86323982359930005737e-03, /* 3F8226E3, E96E8493 */
3.59207910759131235356e-03, /* 3F6D6D22, C9560328 */
1.45620945432529025516e-03, /* 3F57DBC8, FEE08315 */
5.88041240820264096874e-04, /* 3F4344D8, F2F26501 */
2.46463134818469906812e-04, /* 3F3026F7, 1A8D1068 */
7.81794442939557092300e-05, /* 3F147E88, A03792A6 */
7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */
-1.85586374855275456654e-05, /* BEF375CB, DB605373 */
2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */
/* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */
/* pio4 */ 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
/* pio4lo */ 3.06161699786838301793e-17 /* 3C81A626, 33145C07 */
};
#define one T[13]
#define pio4 T[14]
#define pio4lo T[15]
double __tan(double x, double y, int iy)
{
double z, r, v, w, s, sign;
int32_t ix, hx;
GET_HIGH_WORD(hx,x);
ix = hx & 0x7fffffff; /* high word of |x| */
if (ix >= 0x3FE59428) { /* |x| >= 0.6744 */
if (hx < 0) {
x = -x;
y = -y;
}
z = pio4 - x;
w = pio4lo - y;
x = z + w;
y = 0.0;
}
z = x * x;
w = z * z;
/*
* Break x^5*(T[1]+x^2*T[2]+...) into
* x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
* x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
*/
r = T[1] + w*(T[3] + w*(T[5] + w*(T[7] + w*(T[9] + w*T[11]))));
v = z*(T[2] + w*(T[4] + w*(T[6] + w*(T[8] + w*(T[10] + w*T[12])))));
s = z * x;
r = y + z * (s * (r + v) + y);
r += T[0] * s;
w = x + r;
if (ix >= 0x3FE59428) {
v = iy;
sign = 1 - ((hx >> 30) & 2);
return sign * (v - 2.0 * (x - (w * w / (w + v) - r)));
}
if (iy == 1)
return w;
else {
/*
* if allow error up to 2 ulp, simply return
* -1.0 / (x+r) here
*/
/* compute -1.0 / (x+r) accurately */
double a, t;
z = w;
SET_LOW_WORD(z,0);
v = r - (z - x); /* z+v = r+x */
t = a = -1.0 / w; /* a = -1.0/w */
SET_LOW_WORD(t,0);
s = 1.0 + t * z;
return t + a * (s + t * v);
}
}
src/math/__tandf.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/src/k_tanf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
* Optimized by Bruce D. Evans.
*/
/*
* ====================================================
* Copyright 2004 Sun Microsystems, Inc. All Rights Reserved.
*
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "libm.h"
/* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */
static const double T[] = {
0x15554d3418c99f.0p-54, /* 0.333331395030791399758 */
0x1112fd38999f72.0p-55, /* 0.133392002712976742718 */
0x1b54c91d865afe.0p-57, /* 0.0533812378445670393523 */
0x191df3908c33ce.0p-58, /* 0.0245283181166547278873 */
0x185dadfcecf44e.0p-61, /* 0.00297435743359967304927 */
0x1362b9bf971bcd.0p-59, /* 0.00946564784943673166728 */
};
float __tandf(double x, int iy)
{
double z,r,w,s,t,u;
z = x*x;
/*
* Split up the polynomial into small independent terms to give
* opportunities for parallel evaluation. The chosen splitting is
* micro-optimized for Athlons (XP, X64). It costs 2 multiplications
* relative to Horner's method on sequential machines.
*
* We add the small terms from lowest degree up for efficiency on
* non-sequential machines (the lowest degree terms tend to be ready
* earlier). Apart from this, we don't care about order of
* operations, and don't need to to care since we have precision to
* spare. However, the chosen splitting is good for accuracy too,
* and would give results as accurate as Horner's method if the
* small terms were added from highest degree down.
*/
r = T[4] + z*T[5];
t = T[2] + z*T[3];
w = z*z;
s = z*x;
u = T[0] + z*T[1];
r = (x + s*u) + (s*w)*(t + w*r);
if(iy==1) return r;
else return -1.0/r;
}
src/math/__tanl.c 0 → 100644
浏览文件 @ b69f695a
/* origin: FreeBSD /usr/src/lib/msun/ld80/k_tanl.c */
/*
* ====================================================
* Copyright 2004 Sun Microsystems, Inc. All Rights Reserved.
* Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
*
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "libm.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/*
* ld80 version of __tan.c. See __tan.c for most comments.
*/
/*
* Domain [-0.67434, 0.67434], range ~[-2.25e-22, 1.921e-22]
* |tan(x)/x - t(x)| < 2**-71.9
*
* See __cosl.c for more details about the polynomial.
*/
/* Long double constants are slow on these arches, and broken on i386. */
static const volatile double
T3hi = 0.33333333333333331, /* 0x15555555555555.0p-54 */
T3lo = 1.8350121769317163e-17, /* 0x15280000000000.0p-108 */
T5hi = 0.13333333333333336, /* 0x11111111111112.0p-55 */
T5lo = 1.3051083651294260e-17, /* 0x1e180000000000.0p-109 */
T7hi = 0.053968253968250494, /* 0x1ba1ba1ba1b827.0p-57 */
T7lo = 3.1509625637859973e-18, /* 0x1d100000000000.0p-111 */
pio4_hi = 0.78539816339744828, /* 0x1921fb54442d18.0p-53 */
pio4_lo = 3.0628711372715500e-17, /* 0x11a80000000000.0p-107 */
pio4lo_hi = -1.2541394031670831e-20, /* -0x1d9cceba3f91f2.0p-119 */
pio4lo_lo = 6.1493048227390915e-37; /* 0x1a280000000000.0p-173 */
#define T3 ((long double)T3hi + T3lo)
#define T5 ((long double)T5hi + T5lo)
#define T7 ((long double)T7hi + T7lo)
#define pio4 ((long double)pio4_hi + pio4_lo)
#define pio4lo ((long double)pio4lo_hi + pio4lo_lo)
#if 0
static const long double
T3 = 0.333333333333333333180L, /* 0xaaaaaaaaaaaaaaa5.0p-65 */
T5 = 0.133333333333333372290L, /* 0x88888888888893c3.0p-66 */
T7 = 0.0539682539682504975744L, /* 0xdd0dd0dd0dc13ba2.0p-68 */
pio4 = 0.785398163397448309628L, /* 0xc90fdaa22168c235.0p-64 */
pio4lo = -1.25413940316708300586e-20L; /* -0xece675d1fc8f8cbb.0p-130 */
#endif
static const double
T9 = 0.021869488536312216, /* 0x1664f4882cc1c2.0p-58 */
T11 = 0.0088632355256619590, /* 0x1226e355c17612.0p-59 */
T13 = 0.0035921281113786528, /* 0x1d6d3d185d7ff8.0p-61 */
T15 = 0.0014558334756312418, /* 0x17da354aa3f96b.0p-62 */
T17 = 0.00059003538700862256, /* 0x13559358685b83.0p-63 */
T19 = 0.00023907843576635544, /* 0x1f56242026b5be.0p-65 */
T21 = 0.000097154625656538905, /* 0x1977efc26806f4.0p-66 */
T23 = 0.000038440165747303162, /* 0x14275a09b3ceac.0p-67 */
T25 = 0.000018082171885432524, /* 0x12f5e563e5487e.0p-68 */
T27 = 0.0000024196006108814377, /* 0x144c0d80cc6896.0p-71 */
T29 = 0.0000078293456938132840, /* 0x106b59141a6cb3.0p-69 */
T31 = -0.0000032609076735050182, /* -0x1b5abef3ba4b59.0p-71 */
T33 = 0.0000023261313142559411; /* 0x13835436c0c87f.0p-71 */
long double __tanl(long double x, long double y, int iy) {
long double z, r, v, w, s, a, t;
long double osign;
int i;
iy = iy == 1 ? -1 : 1; /* XXX recover original interface */
// FIXME: this is wrong, use copysign, signbit or union bithack
osign = x >= 0 ? 1.0 : -1.0; /* XXX slow, probably wrong for -0 */
if (fabsl(x) >= 0.67434) {
if (x < 0) {
x = -x;
y = -y;
}
z = pio4 - x;
w = pio4lo - y;
x = z + w;
y = 0.0;
i = 1;
} else
i = 0;
z = x * x;
w = z * z;
r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
w * (T25 + w * (T29 + w * T33))))));
v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
w * (T27 + w * T31))))));
s = z * x;
r = y + z * (s * (r + v) + y);
r += T3 * s;
w = x + r;
if (i == 1) {
v = (long double)iy;
return osign * (v - 2.0 * (x - (w * w / (w + v) - r)));
}
if (iy == 1)
return w;
/*
* if allow error up to 2 ulp, simply return
* -1.0 / (x+r) here
*/
/* compute -1.0 / (x+r) accurately */
z = w;
z = z + 0x1p32 - 0x1p32;
v = r - (z - x); /* z+v = r+x */
t = a = -1.0 / w; /* a = -1.0/w */
t = t + 0x1p32 - 0x1p32;
s = 1.0 + t * z;
return t + a * (s + t * v);
}
#endif
src/math/acosf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/acoshl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/acosl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/asinhf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/asinhl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/asinl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/atan2f.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/atan2l.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/atanhl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/atanl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/cbrt.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/cbrtf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/cbrtl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/ceil.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/ceilf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/ceill.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/copysign.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/copysignf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/copysignl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/s_cos.c → src/math/cos.c
浏览文件 @ b69f695a
此差异已折叠。
src/math/cosf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/coshl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/cosl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/e_acosf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_acoshf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_asinf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_atan2f.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_coshf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_expf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_fmod.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_fmodf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_hypot.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_hypotf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_log10.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_logf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_pow.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_powf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_rem_pio2.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_rem_pio2f.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_remainder.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_remainderf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_scalbf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_sinhf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_sqrt.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/e_sqrtf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_erf.c → src/math/erf.c
浏览文件 @ b69f695a
此差异已折叠。
src/math/erfl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/e_exp.c → src/math/exp.c
浏览文件 @ b69f695a
此差异已折叠。
src/math/exp2.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/exp2f.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/exp2l.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/expf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/expl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/expm1f.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/expm1l.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fabs.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fabsf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fabsl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fdim.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fdimf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fdiml.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/floor.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/floorl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fma.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fmaf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fmal.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fmax.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fmaxf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fmaxl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fmin.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fminf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fminl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fmod.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fmodf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/fmodl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/hypot.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/hypotf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/hypotl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/i386/e_exp.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/e_expf.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/e_log.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/e_log10.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/e_log10f.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/e_logf.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/e_remainder.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/s_fabs.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/s_fabsf.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/s_rint.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/s_rintf.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/s_scalbln.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/s_scalblnf.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/s_trunc.s 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/i386/sqrtl.s 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/ilogb.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/ilogbf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/ilogbl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/j0.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/j0f.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/j1.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/j1f.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/jn.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/jnf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/k_cosf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/k_rem_pio2.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/k_rem_pio2f.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/k_sinf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/k_tan.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/k_tanf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/ldexpl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lgamma.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lgamma_r.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lgammaf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lgammaf_r.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lgammal.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/llrint.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/llrintf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/llrintl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/llround.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/llroundf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/llroundl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/e_log.c → src/math/log.c
浏览文件 @ b69f695a
此差异已折叠。
src/math/log10.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/log10f.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/log10l.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/log1pf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/log1pl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/log2.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/log2f.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/log2l.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/logb.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/logbf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/logbl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/logf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/logl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lrint.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lrintf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lrintl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lround.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lroundf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/lroundl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/math_private.h 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/modf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/modff.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/modfl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/nearbyint.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/nearbyintf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/nearbyintl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/nextafter.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/nextafterf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/nextafterl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/nexttoward.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/nexttowardf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/nexttowardl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/pow.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/powf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/powl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/remainder.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/remainderf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/remainderl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/remquo.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/remquof.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/remquol.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/rint.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/rintl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/roundl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/s_cbrt.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_cbrtf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_ceil.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_ceilf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_copysign.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_expm1f.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_floor.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_ilogb.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_llrint.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_log1pf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_logb.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_logbf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_lrint.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_lrintf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_modf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_modff.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_nextafter.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_nextafterf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_remquo.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_remquof.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_rint.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_scalblnf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/s_tanhf.c 已删除 100644 → 0
浏览文件 @ d46cf2e1
此差异已折叠。
src/math/scalbln.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/scalblnf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/scalblnl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/scalbnl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/signgam.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/s_sin.c → src/math/sin.c
浏览文件 @ b69f695a
此差异已折叠。
src/math/sinf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/sinhf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/sinhl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/sinl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/sqrt.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/sqrtf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/s_tan.c → src/math/tan.c
浏览文件 @ b69f695a
此差异已折叠。
src/math/tanf.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/tanhl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/tanl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/tgammal.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/truncl.c 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
src/math/x86_64/sqrtl.s 0 → 100644
浏览文件 @ b69f695a
此差异已折叠。
Markdown is supported
0% 或 .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑！
想要评论请 注册 或