Bring over more math functions from musl
* asin * asinf * acosf * tan * tanf * atan * atan2 * atan2f * atanf * cbrt * cbrtf * log * logf * log2 * log2f
This commit is contained in:
parent
116e614ab6
commit
e85d774245
|
|
@ -2,13 +2,12 @@
|
|||
#define _MATH_H
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C"
|
||||
{
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
// NOTE(orca): not doing anything fancy for float_t and double_t
|
||||
typedef float float_t;
|
||||
typedef double double_t;
|
||||
// NOTE(orca): not doing anything fancy for float_t and double_t
|
||||
typedef float float_t;
|
||||
typedef double double_t;
|
||||
|
||||
#define NAN __builtin_nanf("")
|
||||
#define INFINITY __builtin_inff()
|
||||
|
|
@ -19,33 +18,33 @@ extern "C"
|
|||
#define FP_SUBNORMAL 3
|
||||
#define FP_NORMAL 4
|
||||
|
||||
int __fpclassify(double);
|
||||
int __fpclassifyf(float);
|
||||
int __fpclassifyl(long double);
|
||||
int __fpclassify(double);
|
||||
int __fpclassifyf(float);
|
||||
int __fpclassifyl(long double);
|
||||
|
||||
static __inline unsigned __FLOAT_BITS(float __f)
|
||||
static __inline unsigned __FLOAT_BITS(float __f)
|
||||
{
|
||||
union
|
||||
{
|
||||
union
|
||||
{
|
||||
float __f;
|
||||
unsigned __i;
|
||||
} __u;
|
||||
float __f;
|
||||
unsigned __i;
|
||||
} __u;
|
||||
|
||||
__u.__f = __f;
|
||||
return __u.__i;
|
||||
}
|
||||
__u.__f = __f;
|
||||
return __u.__i;
|
||||
}
|
||||
|
||||
static __inline unsigned long long __DOUBLE_BITS(double __f)
|
||||
static __inline unsigned long long __DOUBLE_BITS(double __f)
|
||||
{
|
||||
union
|
||||
{
|
||||
union
|
||||
{
|
||||
double __f;
|
||||
unsigned long long __i;
|
||||
} __u;
|
||||
double __f;
|
||||
unsigned long long __i;
|
||||
} __u;
|
||||
|
||||
__u.__f = __f;
|
||||
return __u.__i;
|
||||
}
|
||||
__u.__f = __f;
|
||||
return __u.__i;
|
||||
}
|
||||
|
||||
#define fpclassify(x) ( \
|
||||
sizeof(x) == sizeof(float) ? __fpclassifyf(x) : sizeof(x) == sizeof(double) ? __fpclassify(x) \
|
||||
|
|
@ -59,29 +58,49 @@ extern "C"
|
|||
sizeof(x) == sizeof(float) ? (__FLOAT_BITS(x) & 0x7fffffff) > 0x7f800000 : sizeof(x) == sizeof(double) ? (__DOUBLE_BITS(x) & -1ULL >> 1) > 0x7ffULL << 52 \
|
||||
: __fpclassifyl(x) == FP_NAN)
|
||||
|
||||
double acos(double);
|
||||
double fabs(double);
|
||||
float fabsf(float);
|
||||
|
||||
double ceil(double);
|
||||
double acos(double);
|
||||
float acosf(float);
|
||||
|
||||
double cos(double);
|
||||
float cosf(float);
|
||||
double cbrt(double);
|
||||
float cbrtf(float);
|
||||
|
||||
double fabs(double);
|
||||
float fabsf(float);
|
||||
double ceil(double);
|
||||
|
||||
double floor(double);
|
||||
double cos(double);
|
||||
float cosf(float);
|
||||
|
||||
double fmod(double, double);
|
||||
double floor(double);
|
||||
|
||||
double pow(double, double);
|
||||
double fmod(double, double);
|
||||
|
||||
double scalbn(double, int);
|
||||
double log(double);
|
||||
float logf(float);
|
||||
double log2(double);
|
||||
float log2f(float);
|
||||
|
||||
double sin(double);
|
||||
float sinf(float);
|
||||
double pow(double, double);
|
||||
|
||||
double sqrt(double);
|
||||
float sqrtf(float);
|
||||
double scalbn(double, int);
|
||||
|
||||
double sin(double);
|
||||
float sinf(float);
|
||||
|
||||
double asin(double);
|
||||
float asinf(float);
|
||||
|
||||
double tan(double);
|
||||
float tanf(float);
|
||||
|
||||
double atan(double);
|
||||
float atanf(float);
|
||||
double atan2(double, double);
|
||||
float atan2f(float, float);
|
||||
|
||||
double sqrt(double);
|
||||
float sqrtf(float);
|
||||
|
||||
#define M_E 2.7182818284590452354 /* e */
|
||||
#define M_LOG2E 1.4426950408889634074 /* log_2 e */
|
||||
|
|
@ -97,6 +116,11 @@ extern "C"
|
|||
#define M_SQRT2 1.41421356237309504880 /* sqrt(2) */
|
||||
#define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */
|
||||
|
||||
//NOTE(orca) - implementation details
|
||||
typedef unsigned uint32_t;
|
||||
double __math_divzero(uint32_t sign);
|
||||
float __math_divzerof(uint32_t sign);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -0,0 +1,3 @@
|
|||
---
|
||||
DisableFormat: true
|
||||
SortIncludes: Never
|
||||
|
|
@ -0,0 +1,7 @@
|
|||
#include "libm.h"
|
||||
|
||||
double __math_divzero(uint32_t sign)
|
||||
{
|
||||
// NOTE(orca): no fp barriers
|
||||
return (sign ? -1.0 : 1.0) / 0.0;
|
||||
}
|
||||
|
|
@ -0,0 +1,7 @@
|
|||
#include "libm.h"
|
||||
|
||||
float __math_divzerof(uint32_t sign)
|
||||
{
|
||||
// NOTE(orca): no fp barriers
|
||||
return (sign ? -1.0f : 1.0f) / 0.0f;
|
||||
}
|
||||
|
|
@ -0,0 +1,110 @@
|
|||
/* 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 odd indicates whether tan (if odd = 0) or -1/tan (if odd = 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 */
|
||||
},
|
||||
pio4 = 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
|
||||
pio4lo = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */
|
||||
|
||||
double __tan(double x, double y, int odd)
|
||||
{
|
||||
double_t z, r, v, w, s, a;
|
||||
double w0, a0;
|
||||
uint32_t hx;
|
||||
int big, sign;
|
||||
|
||||
GET_HIGH_WORD(hx,x);
|
||||
big = (hx&0x7fffffff) >= 0x3FE59428; /* |x| >= 0.6744 */
|
||||
if (big) {
|
||||
sign = hx>>31;
|
||||
if (sign) {
|
||||
x = -x;
|
||||
y = -y;
|
||||
}
|
||||
x = (pio4 - x) + (pio4lo - y);
|
||||
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) + s*T[0];
|
||||
w = x + r;
|
||||
if (big) {
|
||||
s = 1 - 2*odd;
|
||||
v = s - 2.0 * (x + (r - w*w/(w + s)));
|
||||
return sign ? -v : v;
|
||||
}
|
||||
if (!odd)
|
||||
return w;
|
||||
/* -1.0/(x+r) has up to 2ulp error, so compute it accurately */
|
||||
w0 = w;
|
||||
SET_LOW_WORD(w0, 0);
|
||||
v = r - (w0 - x); /* w0+v = r+x */
|
||||
a0 = a = -1.0 / w;
|
||||
SET_LOW_WORD(a0, 0);
|
||||
return a0 + a*(1.0 + a0*w0 + a0*v);
|
||||
}
|
||||
|
|
@ -0,0 +1,54 @@
|
|||
/* 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 odd)
|
||||
{
|
||||
double_t 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);
|
||||
return odd ? -1.0/r : r;
|
||||
}
|
||||
|
|
@ -0,0 +1,71 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_acosf.c */
|
||||
/*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
/*
|
||||
* ====================================================
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "libm.h"
|
||||
|
||||
static const float
|
||||
pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */
|
||||
pio2_lo = 7.5497894159e-08, /* 0x33a22168 */
|
||||
pS0 = 1.6666586697e-01,
|
||||
pS1 = -4.2743422091e-02,
|
||||
pS2 = -8.6563630030e-03,
|
||||
qS1 = -7.0662963390e-01;
|
||||
|
||||
static float R(float z)
|
||||
{
|
||||
float_t p, q;
|
||||
p = z*(pS0+z*(pS1+z*pS2));
|
||||
q = 1.0f+z*qS1;
|
||||
return p/q;
|
||||
}
|
||||
|
||||
float acosf(float x)
|
||||
{
|
||||
float z,w,s,c,df;
|
||||
uint32_t hx,ix;
|
||||
|
||||
GET_FLOAT_WORD(hx, x);
|
||||
ix = hx & 0x7fffffff;
|
||||
/* |x| >= 1 or nan */
|
||||
if (ix >= 0x3f800000) {
|
||||
if (ix == 0x3f800000) {
|
||||
if (hx >> 31)
|
||||
return 2*pio2_hi + 0x1p-120f;
|
||||
return 0;
|
||||
}
|
||||
return 0/(x-x);
|
||||
}
|
||||
/* |x| < 0.5 */
|
||||
if (ix < 0x3f000000) {
|
||||
if (ix <= 0x32800000) /* |x| < 2**-26 */
|
||||
return pio2_hi + 0x1p-120f;
|
||||
return pio2_hi - (x - (pio2_lo-x*R(x*x)));
|
||||
}
|
||||
/* x < -0.5 */
|
||||
if (hx >> 31) {
|
||||
z = (1+x)*0.5f;
|
||||
s = sqrtf(z);
|
||||
w = R(z)*s-pio2_lo;
|
||||
return 2*(pio2_hi - (s+w));
|
||||
}
|
||||
/* x > 0.5 */
|
||||
z = (1-x)*0.5f;
|
||||
s = sqrtf(z);
|
||||
GET_FLOAT_WORD(hx,s);
|
||||
SET_FLOAT_WORD(df,hx&0xfffff000);
|
||||
c = (z-df*df)/(s+df);
|
||||
w = R(z)*s+c;
|
||||
return 2*(df+w);
|
||||
}
|
||||
|
|
@ -0,0 +1,107 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_asin.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.
|
||||
* ====================================================
|
||||
*/
|
||||
/* asin(x)
|
||||
* Method :
|
||||
* Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
|
||||
* we approximate asin(x) on [0,0.5] by
|
||||
* asin(x) = x + x*x^2*R(x^2)
|
||||
* where
|
||||
* R(x^2) is a rational approximation of (asin(x)-x)/x^3
|
||||
* and its remez error is bounded by
|
||||
* |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
|
||||
*
|
||||
* For x in [0.5,1]
|
||||
* asin(x) = pi/2-2*asin(sqrt((1-x)/2))
|
||||
* Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
|
||||
* then for x>0.98
|
||||
* asin(x) = pi/2 - 2*(s+s*z*R(z))
|
||||
* = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
|
||||
* For x<=0.98, let pio4_hi = pio2_hi/2, then
|
||||
* f = hi part of s;
|
||||
* c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
|
||||
* and
|
||||
* asin(x) = pi/2 - 2*(s+s*z*R(z))
|
||||
* = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
|
||||
* = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
|
||||
*
|
||||
* Special cases:
|
||||
* if x is NaN, return x itself;
|
||||
* if |x|>1, return NaN with invalid signal.
|
||||
*
|
||||
*/
|
||||
|
||||
#include "libm.h"
|
||||
|
||||
static const double
|
||||
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
|
||||
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
|
||||
/* coefficients for R(x^2) */
|
||||
pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
|
||||
pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
|
||||
pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
|
||||
pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
|
||||
pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
|
||||
pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
|
||||
qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
|
||||
qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
|
||||
qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
|
||||
qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
|
||||
|
||||
static double R(double z)
|
||||
{
|
||||
double_t p, q;
|
||||
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
|
||||
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
|
||||
return p/q;
|
||||
}
|
||||
|
||||
double asin(double x)
|
||||
{
|
||||
double z,r,s;
|
||||
uint32_t hx,ix;
|
||||
|
||||
GET_HIGH_WORD(hx, x);
|
||||
ix = hx & 0x7fffffff;
|
||||
/* |x| >= 1 or nan */
|
||||
if (ix >= 0x3ff00000) {
|
||||
uint32_t lx;
|
||||
GET_LOW_WORD(lx, x);
|
||||
if ((ix-0x3ff00000 | lx) == 0)
|
||||
/* asin(1) = +-pi/2 with inexact */
|
||||
return x*pio2_hi + 0x1p-120f;
|
||||
return 0/(x-x);
|
||||
}
|
||||
/* |x| < 0.5 */
|
||||
if (ix < 0x3fe00000) {
|
||||
/* if 0x1p-1022 <= |x| < 0x1p-26, avoid raising underflow */
|
||||
if (ix < 0x3e500000 && ix >= 0x00100000)
|
||||
return x;
|
||||
return x + x*R(x*x);
|
||||
}
|
||||
/* 1 > |x| >= 0.5 */
|
||||
z = (1 - fabs(x))*0.5;
|
||||
s = sqrt(z);
|
||||
r = R(z);
|
||||
if (ix >= 0x3fef3333) { /* if |x| > 0.975 */
|
||||
x = pio2_hi-(2*(s+s*r)-pio2_lo);
|
||||
} else {
|
||||
double f,c;
|
||||
/* f+c = sqrt(z) */
|
||||
f = s;
|
||||
SET_LOW_WORD(f,0);
|
||||
c = (z-f*f)/(s+f);
|
||||
x = 0.5*pio2_hi - (2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));
|
||||
}
|
||||
if (hx >> 31)
|
||||
return -x;
|
||||
return x;
|
||||
}
|
||||
|
|
@ -0,0 +1,61 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.c */
|
||||
/*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
/*
|
||||
* ====================================================
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
#include "libm.h"
|
||||
|
||||
static const double
|
||||
pio2 = 1.570796326794896558e+00;
|
||||
|
||||
static const float
|
||||
/* coefficients for R(x^2) */
|
||||
pS0 = 1.6666586697e-01,
|
||||
pS1 = -4.2743422091e-02,
|
||||
pS2 = -8.6563630030e-03,
|
||||
qS1 = -7.0662963390e-01;
|
||||
|
||||
static float R(float z)
|
||||
{
|
||||
float_t p, q;
|
||||
p = z*(pS0+z*(pS1+z*pS2));
|
||||
q = 1.0f+z*qS1;
|
||||
return p/q;
|
||||
}
|
||||
|
||||
float asinf(float x)
|
||||
{
|
||||
double s;
|
||||
float z;
|
||||
uint32_t hx,ix;
|
||||
|
||||
GET_FLOAT_WORD(hx, x);
|
||||
ix = hx & 0x7fffffff;
|
||||
if (ix >= 0x3f800000) { /* |x| >= 1 */
|
||||
if (ix == 0x3f800000) /* |x| == 1 */
|
||||
return x*pio2 + 0x1p-120f; /* asin(+-1) = +-pi/2 with inexact */
|
||||
return 0/(x-x); /* asin(|x|>1) is NaN */
|
||||
}
|
||||
if (ix < 0x3f000000) { /* |x| < 0.5 */
|
||||
/* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */
|
||||
if (ix < 0x39800000 && ix >= 0x00800000)
|
||||
return x;
|
||||
return x + x*R(x*x);
|
||||
}
|
||||
/* 1 > |x| >= 0.5 */
|
||||
z = (1 - fabsf(x))*0.5f;
|
||||
s = sqrt(z);
|
||||
x = pio2 - 2*(s+s*R(z));
|
||||
if (hx >> 31)
|
||||
return -x;
|
||||
return x;
|
||||
}
|
||||
|
|
@ -0,0 +1,116 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_atan.c */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
/* atan(x)
|
||||
* Method
|
||||
* 1. Reduce x to positive by atan(x) = -atan(-x).
|
||||
* 2. According to the integer k=4t+0.25 chopped, t=x, the argument
|
||||
* is further reduced to one of the following intervals and the
|
||||
* arctangent of t is evaluated by the corresponding formula:
|
||||
*
|
||||
* [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
|
||||
* [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
|
||||
* [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
|
||||
* [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
|
||||
* [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
|
||||
*
|
||||
* 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 double atanhi[] = {
|
||||
4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
|
||||
7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
|
||||
9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
|
||||
1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
|
||||
};
|
||||
|
||||
static const double atanlo[] = {
|
||||
2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
|
||||
3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
|
||||
1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
|
||||
6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
|
||||
};
|
||||
|
||||
static const double aT[] = {
|
||||
3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
|
||||
-1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
|
||||
1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
|
||||
-1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
|
||||
9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
|
||||
-7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
|
||||
6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
|
||||
-5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
|
||||
4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
|
||||
-3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
|
||||
1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
|
||||
};
|
||||
|
||||
double atan(double x)
|
||||
{
|
||||
double_t w,s1,s2,z;
|
||||
uint32_t ix,sign;
|
||||
int id;
|
||||
|
||||
GET_HIGH_WORD(ix, x);
|
||||
sign = ix >> 31;
|
||||
ix &= 0x7fffffff;
|
||||
if (ix >= 0x44100000) { /* if |x| >= 2^66 */
|
||||
if (isnan(x))
|
||||
return x;
|
||||
z = atanhi[3] + 0x1p-120f;
|
||||
return sign ? -z : z;
|
||||
}
|
||||
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
|
||||
if (ix < 0x3e400000) { /* |x| < 2^-27 */
|
||||
if (ix < 0x00100000)
|
||||
/* raise underflow for subnormal x */
|
||||
FORCE_EVAL((float)x);
|
||||
return x;
|
||||
}
|
||||
id = -1;
|
||||
} else {
|
||||
x = fabs(x);
|
||||
if (ix < 0x3ff30000) { /* |x| < 1.1875 */
|
||||
if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */
|
||||
id = 0;
|
||||
x = (2.0*x-1.0)/(2.0+x);
|
||||
} else { /* 11/16 <= |x| < 19/16 */
|
||||
id = 1;
|
||||
x = (x-1.0)/(x+1.0);
|
||||
}
|
||||
} else {
|
||||
if (ix < 0x40038000) { /* |x| < 2.4375 */
|
||||
id = 2;
|
||||
x = (x-1.5)/(1.0+1.5*x);
|
||||
} else { /* 2.4375 <= |x| < 2^66 */
|
||||
id = 3;
|
||||
x = -1.0/x;
|
||||
}
|
||||
}
|
||||
}
|
||||
/* end of argument reduction */
|
||||
z = x*x;
|
||||
w = z*z;
|
||||
/* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
|
||||
s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
|
||||
s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
|
||||
if (id < 0)
|
||||
return x - x*(s1+s2);
|
||||
z = atanhi[id] - (x*(s1+s2) - atanlo[id] - x);
|
||||
return sign ? -z : z;
|
||||
}
|
||||
|
|
@ -0,0 +1,107 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_atan2.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.
|
||||
* ====================================================
|
||||
*
|
||||
*/
|
||||
/* atan2(y,x)
|
||||
* Method :
|
||||
* 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
|
||||
* 2. Reduce x to positive by (if x and y are unexceptional):
|
||||
* ARG (x+iy) = arctan(y/x) ... if x > 0,
|
||||
* ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
|
||||
*
|
||||
* Special cases:
|
||||
*
|
||||
* ATAN2((anything), NaN ) is NaN;
|
||||
* ATAN2(NAN , (anything) ) is NaN;
|
||||
* ATAN2(+-0, +(anything but NaN)) is +-0 ;
|
||||
* ATAN2(+-0, -(anything but NaN)) is +-pi ;
|
||||
* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
|
||||
* ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
|
||||
* ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
|
||||
* ATAN2(+-INF,+INF ) is +-pi/4 ;
|
||||
* ATAN2(+-INF,-INF ) is +-3pi/4;
|
||||
* ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
|
||||
*
|
||||
* 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 double
|
||||
pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
|
||||
pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
|
||||
|
||||
double atan2(double y, double x)
|
||||
{
|
||||
double z;
|
||||
uint32_t m,lx,ly,ix,iy;
|
||||
|
||||
if (isnan(x) || isnan(y))
|
||||
return x+y;
|
||||
EXTRACT_WORDS(ix, lx, x);
|
||||
EXTRACT_WORDS(iy, ly, y);
|
||||
if ((ix-0x3ff00000 | lx) == 0) /* x = 1.0 */
|
||||
return atan(y);
|
||||
m = ((iy>>31)&1) | ((ix>>30)&2); /* 2*sign(x)+sign(y) */
|
||||
ix = ix & 0x7fffffff;
|
||||
iy = iy & 0x7fffffff;
|
||||
|
||||
/* when y = 0 */
|
||||
if ((iy|ly) == 0) {
|
||||
switch(m) {
|
||||
case 0:
|
||||
case 1: return y; /* atan(+-0,+anything)=+-0 */
|
||||
case 2: return pi; /* atan(+0,-anything) = pi */
|
||||
case 3: return -pi; /* atan(-0,-anything) =-pi */
|
||||
}
|
||||
}
|
||||
/* when x = 0 */
|
||||
if ((ix|lx) == 0)
|
||||
return m&1 ? -pi/2 : pi/2;
|
||||
/* when x is INF */
|
||||
if (ix == 0x7ff00000) {
|
||||
if (iy == 0x7ff00000) {
|
||||
switch(m) {
|
||||
case 0: return pi/4; /* atan(+INF,+INF) */
|
||||
case 1: return -pi/4; /* atan(-INF,+INF) */
|
||||
case 2: return 3*pi/4; /* atan(+INF,-INF) */
|
||||
case 3: return -3*pi/4; /* atan(-INF,-INF) */
|
||||
}
|
||||
} else {
|
||||
switch(m) {
|
||||
case 0: return 0.0; /* atan(+...,+INF) */
|
||||
case 1: return -0.0; /* atan(-...,+INF) */
|
||||
case 2: return pi; /* atan(+...,-INF) */
|
||||
case 3: return -pi; /* atan(-...,-INF) */
|
||||
}
|
||||
}
|
||||
}
|
||||
/* |y/x| > 0x1p64 */
|
||||
if (ix+(64<<20) < iy || iy == 0x7ff00000)
|
||||
return m&1 ? -pi/2 : pi/2;
|
||||
|
||||
/* z = atan(|y/x|) without spurious underflow */
|
||||
if ((m&2) && iy+(64<<20) < ix) /* |y/x| < 0x1p-64, x<0 */
|
||||
z = 0;
|
||||
else
|
||||
z = atan(fabs(y/x));
|
||||
switch (m) {
|
||||
case 0: return z; /* atan(+,+) */
|
||||
case 1: return -z; /* atan(-,+) */
|
||||
case 2: return pi - (z-pi_lo); /* atan(+,-) */
|
||||
default: /* case 3 */
|
||||
return (z-pi_lo) - pi; /* atan(-,-) */
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,83 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_atan2f.c */
|
||||
/*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
/*
|
||||
* ====================================================
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "libm.h"
|
||||
|
||||
static const float
|
||||
pi = 3.1415927410e+00, /* 0x40490fdb */
|
||||
pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */
|
||||
|
||||
float atan2f(float y, float x)
|
||||
{
|
||||
float z;
|
||||
uint32_t m,ix,iy;
|
||||
|
||||
if (isnan(x) || isnan(y))
|
||||
return x+y;
|
||||
GET_FLOAT_WORD(ix, x);
|
||||
GET_FLOAT_WORD(iy, y);
|
||||
if (ix == 0x3f800000) /* x=1.0 */
|
||||
return atanf(y);
|
||||
m = ((iy>>31)&1) | ((ix>>30)&2); /* 2*sign(x)+sign(y) */
|
||||
ix &= 0x7fffffff;
|
||||
iy &= 0x7fffffff;
|
||||
|
||||
/* when y = 0 */
|
||||
if (iy == 0) {
|
||||
switch (m) {
|
||||
case 0:
|
||||
case 1: return y; /* atan(+-0,+anything)=+-0 */
|
||||
case 2: return pi; /* atan(+0,-anything) = pi */
|
||||
case 3: return -pi; /* atan(-0,-anything) =-pi */
|
||||
}
|
||||
}
|
||||
/* when x = 0 */
|
||||
if (ix == 0)
|
||||
return m&1 ? -pi/2 : pi/2;
|
||||
/* when x is INF */
|
||||
if (ix == 0x7f800000) {
|
||||
if (iy == 0x7f800000) {
|
||||
switch (m) {
|
||||
case 0: return pi/4; /* atan(+INF,+INF) */
|
||||
case 1: return -pi/4; /* atan(-INF,+INF) */
|
||||
case 2: return 3*pi/4; /*atan(+INF,-INF)*/
|
||||
case 3: return -3*pi/4; /*atan(-INF,-INF)*/
|
||||
}
|
||||
} else {
|
||||
switch (m) {
|
||||
case 0: return 0.0f; /* atan(+...,+INF) */
|
||||
case 1: return -0.0f; /* atan(-...,+INF) */
|
||||
case 2: return pi; /* atan(+...,-INF) */
|
||||
case 3: return -pi; /* atan(-...,-INF) */
|
||||
}
|
||||
}
|
||||
}
|
||||
/* |y/x| > 0x1p26 */
|
||||
if (ix+(26<<23) < iy || iy == 0x7f800000)
|
||||
return m&1 ? -pi/2 : pi/2;
|
||||
|
||||
/* z = atan(|y/x|) with correct underflow */
|
||||
if ((m&2) && iy+(26<<23) < ix) /*|y/x| < 0x1p-26, x < 0 */
|
||||
z = 0.0;
|
||||
else
|
||||
z = atanf(fabsf(y/x));
|
||||
switch (m) {
|
||||
case 0: return z; /* atan(+,+) */
|
||||
case 1: return -z; /* atan(-,+) */
|
||||
case 2: return pi - (z-pi_lo); /* atan(+,-) */
|
||||
default: /* case 3 */
|
||||
return (z-pi_lo) - pi; /* atan(-,-) */
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,94 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_atanf.c */
|
||||
/*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
/*
|
||||
* ====================================================
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
|
||||
#include "libm.h"
|
||||
|
||||
static const float atanhi[] = {
|
||||
4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */
|
||||
7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */
|
||||
9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */
|
||||
1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */
|
||||
};
|
||||
|
||||
static const float atanlo[] = {
|
||||
5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */
|
||||
3.7748947079e-08, /* atan(1.0)lo 0x33222168 */
|
||||
3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */
|
||||
7.5497894159e-08, /* atan(inf)lo 0x33a22168 */
|
||||
};
|
||||
|
||||
static const float aT[] = {
|
||||
3.3333328366e-01,
|
||||
-1.9999158382e-01,
|
||||
1.4253635705e-01,
|
||||
-1.0648017377e-01,
|
||||
6.1687607318e-02,
|
||||
};
|
||||
|
||||
float atanf(float x)
|
||||
{
|
||||
float_t w,s1,s2,z;
|
||||
uint32_t ix,sign;
|
||||
int id;
|
||||
|
||||
GET_FLOAT_WORD(ix, x);
|
||||
sign = ix>>31;
|
||||
ix &= 0x7fffffff;
|
||||
if (ix >= 0x4c800000) { /* if |x| >= 2**26 */
|
||||
if (isnan(x))
|
||||
return x;
|
||||
z = atanhi[3] + 0x1p-120f;
|
||||
return sign ? -z : z;
|
||||
}
|
||||
if (ix < 0x3ee00000) { /* |x| < 0.4375 */
|
||||
if (ix < 0x39800000) { /* |x| < 2**-12 */
|
||||
if (ix < 0x00800000)
|
||||
/* raise underflow for subnormal x */
|
||||
FORCE_EVAL(x*x);
|
||||
return x;
|
||||
}
|
||||
id = -1;
|
||||
} else {
|
||||
x = fabsf(x);
|
||||
if (ix < 0x3f980000) { /* |x| < 1.1875 */
|
||||
if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */
|
||||
id = 0;
|
||||
x = (2.0f*x - 1.0f)/(2.0f + x);
|
||||
} else { /* 11/16 <= |x| < 19/16 */
|
||||
id = 1;
|
||||
x = (x - 1.0f)/(x + 1.0f);
|
||||
}
|
||||
} else {
|
||||
if (ix < 0x401c0000) { /* |x| < 2.4375 */
|
||||
id = 2;
|
||||
x = (x - 1.5f)/(1.0f + 1.5f*x);
|
||||
} else { /* 2.4375 <= |x| < 2**26 */
|
||||
id = 3;
|
||||
x = -1.0f/x;
|
||||
}
|
||||
}
|
||||
}
|
||||
/* end of argument reduction */
|
||||
z = x*x;
|
||||
w = z*z;
|
||||
/* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
|
||||
s1 = z*(aT[0]+w*(aT[2]+w*aT[4]));
|
||||
s2 = w*(aT[1]+w*aT[3]);
|
||||
if (id < 0)
|
||||
return x - x*(s1+s2);
|
||||
z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
|
||||
return sign ? -z : z;
|
||||
}
|
||||
|
|
@ -0,0 +1,103 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrt.c */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*
|
||||
* Optimized by Bruce D. Evans.
|
||||
*/
|
||||
/* cbrt(x)
|
||||
* Return cube root of x
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include <stdint.h>
|
||||
|
||||
static const uint32_t
|
||||
B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */
|
||||
B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */
|
||||
|
||||
/* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */
|
||||
static const double
|
||||
P0 = 1.87595182427177009643, /* 0x3ffe03e6, 0x0f61e692 */
|
||||
P1 = -1.88497979543377169875, /* 0xbffe28e0, 0x92f02420 */
|
||||
P2 = 1.621429720105354466140, /* 0x3ff9f160, 0x4a49d6c2 */
|
||||
P3 = -0.758397934778766047437, /* 0xbfe844cb, 0xbee751d9 */
|
||||
P4 = 0.145996192886612446982; /* 0x3fc2b000, 0xd4e4edd7 */
|
||||
|
||||
double cbrt(double x)
|
||||
{
|
||||
union {double f; uint64_t i;} u = {x};
|
||||
double_t r,s,t,w;
|
||||
uint32_t hx = u.i>>32 & 0x7fffffff;
|
||||
|
||||
if (hx >= 0x7ff00000) /* cbrt(NaN,INF) is itself */
|
||||
return x+x;
|
||||
|
||||
/*
|
||||
* Rough cbrt to 5 bits:
|
||||
* cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3)
|
||||
* where e is integral and >= 0, m is real and in [0, 1), and "/" and
|
||||
* "%" are integer division and modulus with rounding towards minus
|
||||
* infinity. The RHS is always >= the LHS and has a maximum relative
|
||||
* error of about 1 in 16. Adding a bias of -0.03306235651 to the
|
||||
* (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE
|
||||
* floating point representation, for finite positive normal values,
|
||||
* ordinary integer divison of the value in bits magically gives
|
||||
* almost exactly the RHS of the above provided we first subtract the
|
||||
* exponent bias (1023 for doubles) and later add it back. We do the
|
||||
* subtraction virtually to keep e >= 0 so that ordinary integer
|
||||
* division rounds towards minus infinity; this is also efficient.
|
||||
*/
|
||||
if (hx < 0x00100000) { /* zero or subnormal? */
|
||||
u.f = x*0x1p54;
|
||||
hx = u.i>>32 & 0x7fffffff;
|
||||
if (hx == 0)
|
||||
return x; /* cbrt(0) is itself */
|
||||
hx = hx/3 + B2;
|
||||
} else
|
||||
hx = hx/3 + B1;
|
||||
u.i &= 1ULL<<63;
|
||||
u.i |= (uint64_t)hx << 32;
|
||||
t = u.f;
|
||||
|
||||
/*
|
||||
* New cbrt to 23 bits:
|
||||
* cbrt(x) = t*cbrt(x/t**3) ~= t*P(t**3/x)
|
||||
* where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r)
|
||||
* to within 2**-23.5 when |r - 1| < 1/10. The rough approximation
|
||||
* has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this
|
||||
* gives us bounds for r = t**3/x.
|
||||
*
|
||||
* Try to optimize for parallel evaluation as in __tanf.c.
|
||||
*/
|
||||
r = (t*t)*(t/x);
|
||||
t = t*((P0+r*(P1+r*P2))+((r*r)*r)*(P3+r*P4));
|
||||
|
||||
/*
|
||||
* Round t away from zero to 23 bits (sloppily except for ensuring that
|
||||
* the result is larger in magnitude than cbrt(x) but not much more than
|
||||
* 2 23-bit ulps larger). With rounding towards zero, the error bound
|
||||
* would be ~5/6 instead of ~4/6. With a maximum error of 2 23-bit ulps
|
||||
* in the rounded t, the infinite-precision error in the Newton
|
||||
* approximation barely affects third digit in the final error
|
||||
* 0.667; the error in the rounded t can be up to about 3 23-bit ulps
|
||||
* before the final error is larger than 0.667 ulps.
|
||||
*/
|
||||
u.f = t;
|
||||
u.i = (u.i + 0x80000000) & 0xffffffffc0000000ULL;
|
||||
t = u.f;
|
||||
|
||||
/* one step Newton iteration to 53 bits with error < 0.667 ulps */
|
||||
s = t*t; /* t*t is exact */
|
||||
r = x/s; /* error <= 0.5 ulps; |r| < |t| */
|
||||
w = t+t; /* t+t is exact */
|
||||
r = (r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
|
||||
t = t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */
|
||||
return t;
|
||||
}
|
||||
|
|
@ -0,0 +1,66 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtf.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.
|
||||
* ====================================================
|
||||
*/
|
||||
/* cbrtf(x)
|
||||
* Return cube root of x
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include <stdint.h>
|
||||
|
||||
static const unsigned
|
||||
B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */
|
||||
B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
|
||||
|
||||
float cbrtf(float x)
|
||||
{
|
||||
double_t r,T;
|
||||
union {float f; uint32_t i;} u = {x};
|
||||
uint32_t hx = u.i & 0x7fffffff;
|
||||
|
||||
if (hx >= 0x7f800000) /* cbrt(NaN,INF) is itself */
|
||||
return x + x;
|
||||
|
||||
/* rough cbrt to 5 bits */
|
||||
if (hx < 0x00800000) { /* zero or subnormal? */
|
||||
if (hx == 0)
|
||||
return x; /* cbrt(+-0) is itself */
|
||||
u.f = x*0x1p24f;
|
||||
hx = u.i & 0x7fffffff;
|
||||
hx = hx/3 + B2;
|
||||
} else
|
||||
hx = hx/3 + B1;
|
||||
u.i &= 0x80000000;
|
||||
u.i |= hx;
|
||||
|
||||
/*
|
||||
* First step Newton iteration (solving t*t-x/t == 0) to 16 bits. In
|
||||
* double precision so that its terms can be arranged for efficiency
|
||||
* without causing overflow or underflow.
|
||||
*/
|
||||
T = u.f;
|
||||
r = T*T*T;
|
||||
T = T*((double_t)x+x+r)/(x+r+r);
|
||||
|
||||
/*
|
||||
* Second step Newton iteration to 47 bits. In double precision for
|
||||
* efficiency and accuracy.
|
||||
*/
|
||||
r = T*T*T;
|
||||
T = T*((double_t)x+x+r)/(x+r+r);
|
||||
|
||||
/* rounding to 24 bits is perfect in round-to-nearest mode */
|
||||
return T;
|
||||
}
|
||||
|
|
@ -0,0 +1,112 @@
|
|||
/*
|
||||
* Double-precision log(x) function.
|
||||
*
|
||||
* Copyright (c) 2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include <stdint.h>
|
||||
#include "libm.h"
|
||||
#include "log_data.h"
|
||||
|
||||
#define T __log_data.tab
|
||||
#define T2 __log_data.tab2
|
||||
#define B __log_data.poly1
|
||||
#define A __log_data.poly
|
||||
#define Ln2hi __log_data.ln2hi
|
||||
#define Ln2lo __log_data.ln2lo
|
||||
#define N (1 << LOG_TABLE_BITS)
|
||||
#define OFF 0x3fe6000000000000
|
||||
|
||||
/* Top 16 bits of a double. */
|
||||
static inline uint32_t top16(double x)
|
||||
{
|
||||
return asuint64(x) >> 48;
|
||||
}
|
||||
|
||||
double log(double x)
|
||||
{
|
||||
double_t w, z, r, r2, r3, y, invc, logc, kd, hi, lo;
|
||||
uint64_t ix, iz, tmp;
|
||||
uint32_t top;
|
||||
int k, i;
|
||||
|
||||
ix = asuint64(x);
|
||||
top = top16(x);
|
||||
#define LO asuint64(1.0 - 0x1p-4)
|
||||
#define HI asuint64(1.0 + 0x1.09p-4)
|
||||
if (predict_false(ix - LO < HI - LO)) {
|
||||
/* Handle close to 1.0 inputs separately. */
|
||||
/* Fix sign of zero with downward rounding when x==1. */
|
||||
if (WANT_ROUNDING && predict_false(ix == asuint64(1.0)))
|
||||
return 0;
|
||||
r = x - 1.0;
|
||||
r2 = r * r;
|
||||
r3 = r * r2;
|
||||
y = r3 *
|
||||
(B[1] + r * B[2] + r2 * B[3] +
|
||||
r3 * (B[4] + r * B[5] + r2 * B[6] +
|
||||
r3 * (B[7] + r * B[8] + r2 * B[9] + r3 * B[10])));
|
||||
/* Worst-case error is around 0.507 ULP. */
|
||||
w = r * 0x1p27;
|
||||
double_t rhi = r + w - w;
|
||||
double_t rlo = r - rhi;
|
||||
w = rhi * rhi * B[0]; /* B[0] == -0.5. */
|
||||
hi = r + w;
|
||||
lo = r - hi + w;
|
||||
lo += B[0] * rlo * (rhi + r);
|
||||
y += lo;
|
||||
y += hi;
|
||||
return eval_as_double(y);
|
||||
}
|
||||
if (predict_false(top - 0x0010 >= 0x7ff0 - 0x0010)) {
|
||||
/* x < 0x1p-1022 or inf or nan. */
|
||||
if (ix * 2 == 0)
|
||||
return __math_divzero(1);
|
||||
if (ix == asuint64(INFINITY)) /* log(inf) == inf. */
|
||||
return x;
|
||||
if ((top & 0x8000) || (top & 0x7ff0) == 0x7ff0)
|
||||
return __math_invalid(x);
|
||||
/* x is subnormal, normalize it. */
|
||||
ix = asuint64(x * 0x1p52);
|
||||
ix -= 52ULL << 52;
|
||||
}
|
||||
|
||||
/* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
|
||||
The range is split into N subintervals.
|
||||
The ith subinterval contains z and c is near its center. */
|
||||
tmp = ix - OFF;
|
||||
i = (tmp >> (52 - LOG_TABLE_BITS)) % N;
|
||||
k = (int64_t)tmp >> 52; /* arithmetic shift */
|
||||
iz = ix - (tmp & 0xfffULL << 52);
|
||||
invc = T[i].invc;
|
||||
logc = T[i].logc;
|
||||
z = asdouble(iz);
|
||||
|
||||
/* log(x) = log1p(z/c-1) + log(c) + k*Ln2. */
|
||||
/* r ~= z/c - 1, |r| < 1/(2*N). */
|
||||
#if __FP_FAST_FMA
|
||||
/* rounding error: 0x1p-55/N. */
|
||||
r = __builtin_fma(z, invc, -1.0);
|
||||
#else
|
||||
/* rounding error: 0x1p-55/N + 0x1p-66. */
|
||||
r = (z - T2[i].chi - T2[i].clo) * invc;
|
||||
#endif
|
||||
kd = (double_t)k;
|
||||
|
||||
/* hi + lo = r + log(c) + k*Ln2. */
|
||||
w = kd * Ln2hi + logc;
|
||||
hi = w + r;
|
||||
lo = w - hi + r + kd * Ln2lo;
|
||||
|
||||
/* log(x) = lo + (log1p(r) - r) + hi. */
|
||||
r2 = r * r; /* rounding error: 0x1p-54/N^2. */
|
||||
/* Worst case error if |y| > 0x1p-5:
|
||||
0.5 + 4.13/N + abs-poly-error*2^57 ULP (+ 0.002 ULP without fma)
|
||||
Worst case error if |y| > 0x1p-4:
|
||||
0.5 + 2.06/N + abs-poly-error*2^56 ULP (+ 0.001 ULP without fma). */
|
||||
y = lo + r2 * A[0] +
|
||||
r * r2 * (A[1] + r * A[2] + r2 * (A[3] + r * A[4])) + hi;
|
||||
return eval_as_double(y);
|
||||
}
|
||||
|
|
@ -0,0 +1,122 @@
|
|||
/*
|
||||
* Double-precision log2(x) function.
|
||||
*
|
||||
* Copyright (c) 2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include <stdint.h>
|
||||
#include "libm.h"
|
||||
#include "log2_data.h"
|
||||
|
||||
#define T __log2_data.tab
|
||||
#define T2 __log2_data.tab2
|
||||
#define B __log2_data.poly1
|
||||
#define A __log2_data.poly
|
||||
#define InvLn2hi __log2_data.invln2hi
|
||||
#define InvLn2lo __log2_data.invln2lo
|
||||
#define N (1 << LOG2_TABLE_BITS)
|
||||
#define OFF 0x3fe6000000000000
|
||||
|
||||
/* Top 16 bits of a double. */
|
||||
static inline uint32_t top16(double x)
|
||||
{
|
||||
return asuint64(x) >> 48;
|
||||
}
|
||||
|
||||
double log2(double x)
|
||||
{
|
||||
double_t z, r, r2, r4, y, invc, logc, kd, hi, lo, t1, t2, t3, p;
|
||||
uint64_t ix, iz, tmp;
|
||||
uint32_t top;
|
||||
int k, i;
|
||||
|
||||
ix = asuint64(x);
|
||||
top = top16(x);
|
||||
#define LO asuint64(1.0 - 0x1.5b51p-5)
|
||||
#define HI asuint64(1.0 + 0x1.6ab2p-5)
|
||||
if (predict_false(ix - LO < HI - LO)) {
|
||||
/* Handle close to 1.0 inputs separately. */
|
||||
/* Fix sign of zero with downward rounding when x==1. */
|
||||
if (WANT_ROUNDING && predict_false(ix == asuint64(1.0)))
|
||||
return 0;
|
||||
r = x - 1.0;
|
||||
#if __FP_FAST_FMA
|
||||
hi = r * InvLn2hi;
|
||||
lo = r * InvLn2lo + __builtin_fma(r, InvLn2hi, -hi);
|
||||
#else
|
||||
double_t rhi, rlo;
|
||||
rhi = asdouble(asuint64(r) & -1ULL << 32);
|
||||
rlo = r - rhi;
|
||||
hi = rhi * InvLn2hi;
|
||||
lo = rlo * InvLn2hi + r * InvLn2lo;
|
||||
#endif
|
||||
r2 = r * r; /* rounding error: 0x1p-62. */
|
||||
r4 = r2 * r2;
|
||||
/* Worst-case error is less than 0.54 ULP (0.55 ULP without fma). */
|
||||
p = r2 * (B[0] + r * B[1]);
|
||||
y = hi + p;
|
||||
lo += hi - y + p;
|
||||
lo += r4 * (B[2] + r * B[3] + r2 * (B[4] + r * B[5]) +
|
||||
r4 * (B[6] + r * B[7] + r2 * (B[8] + r * B[9])));
|
||||
y += lo;
|
||||
return eval_as_double(y);
|
||||
}
|
||||
if (predict_false(top - 0x0010 >= 0x7ff0 - 0x0010)) {
|
||||
/* x < 0x1p-1022 or inf or nan. */
|
||||
if (ix * 2 == 0)
|
||||
return __math_divzero(1);
|
||||
if (ix == asuint64(INFINITY)) /* log(inf) == inf. */
|
||||
return x;
|
||||
if ((top & 0x8000) || (top & 0x7ff0) == 0x7ff0)
|
||||
return __math_invalid(x);
|
||||
/* x is subnormal, normalize it. */
|
||||
ix = asuint64(x * 0x1p52);
|
||||
ix -= 52ULL << 52;
|
||||
}
|
||||
|
||||
/* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
|
||||
The range is split into N subintervals.
|
||||
The ith subinterval contains z and c is near its center. */
|
||||
tmp = ix - OFF;
|
||||
i = (tmp >> (52 - LOG2_TABLE_BITS)) % N;
|
||||
k = (int64_t)tmp >> 52; /* arithmetic shift */
|
||||
iz = ix - (tmp & 0xfffULL << 52);
|
||||
invc = T[i].invc;
|
||||
logc = T[i].logc;
|
||||
z = asdouble(iz);
|
||||
kd = (double_t)k;
|
||||
|
||||
/* log2(x) = log2(z/c) + log2(c) + k. */
|
||||
/* r ~= z/c - 1, |r| < 1/(2*N). */
|
||||
#if __FP_FAST_FMA
|
||||
/* rounding error: 0x1p-55/N. */
|
||||
r = __builtin_fma(z, invc, -1.0);
|
||||
t1 = r * InvLn2hi;
|
||||
t2 = r * InvLn2lo + __builtin_fma(r, InvLn2hi, -t1);
|
||||
#else
|
||||
double_t rhi, rlo;
|
||||
/* rounding error: 0x1p-55/N + 0x1p-65. */
|
||||
r = (z - T2[i].chi - T2[i].clo) * invc;
|
||||
rhi = asdouble(asuint64(r) & -1ULL << 32);
|
||||
rlo = r - rhi;
|
||||
t1 = rhi * InvLn2hi;
|
||||
t2 = rlo * InvLn2hi + r * InvLn2lo;
|
||||
#endif
|
||||
|
||||
/* hi + lo = r/ln2 + log2(c) + k. */
|
||||
t3 = kd + logc;
|
||||
hi = t3 + t1;
|
||||
lo = t3 - hi + t1 + t2;
|
||||
|
||||
/* log2(r+1) = r/ln2 + r^2*poly(r). */
|
||||
/* Evaluation is optimized assuming superscalar pipelined execution. */
|
||||
r2 = r * r; /* rounding error: 0x1p-54/N^2. */
|
||||
r4 = r2 * r2;
|
||||
/* Worst-case error if |y| > 0x1p-4: 0.547 ULP (0.550 ULP without fma).
|
||||
~ 0.5 + 2/N/ln2 + abs-poly-error*0x1p56 ULP (+ 0.003 ULP without fma). */
|
||||
p = A[0] + r * A[1] + r2 * (A[2] + r * A[3]) + r4 * (A[4] + r * A[5]);
|
||||
y = lo + r2 * p + hi;
|
||||
return eval_as_double(y);
|
||||
}
|
||||
|
|
@ -0,0 +1,201 @@
|
|||
/*
|
||||
* Data for log2.
|
||||
*
|
||||
* Copyright (c) 2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
|
||||
#include "log2_data.h"
|
||||
|
||||
#define N (1 << LOG2_TABLE_BITS)
|
||||
|
||||
const struct log2_data __log2_data = {
|
||||
// First coefficient: 0x1.71547652b82fe1777d0ffda0d24p0
|
||||
.invln2hi = 0x1.7154765200000p+0,
|
||||
.invln2lo = 0x1.705fc2eefa200p-33,
|
||||
.poly1 = {
|
||||
// relative error: 0x1.2fad8188p-63
|
||||
// in -0x1.5b51p-5 0x1.6ab2p-5
|
||||
-0x1.71547652b82fep-1,
|
||||
0x1.ec709dc3a03f7p-2,
|
||||
-0x1.71547652b7c3fp-2,
|
||||
0x1.2776c50f05be4p-2,
|
||||
-0x1.ec709dd768fe5p-3,
|
||||
0x1.a61761ec4e736p-3,
|
||||
-0x1.7153fbc64a79bp-3,
|
||||
0x1.484d154f01b4ap-3,
|
||||
-0x1.289e4a72c383cp-3,
|
||||
0x1.0b32f285aee66p-3,
|
||||
},
|
||||
.poly = {
|
||||
// relative error: 0x1.a72c2bf8p-58
|
||||
// abs error: 0x1.67a552c8p-66
|
||||
// in -0x1.f45p-8 0x1.f45p-8
|
||||
-0x1.71547652b8339p-1,
|
||||
0x1.ec709dc3a04bep-2,
|
||||
-0x1.7154764702ffbp-2,
|
||||
0x1.2776c50034c48p-2,
|
||||
-0x1.ec7b328ea92bcp-3,
|
||||
0x1.a6225e117f92ep-3,
|
||||
},
|
||||
/* Algorithm:
|
||||
|
||||
x = 2^k z
|
||||
log2(x) = k + log2(c) + log2(z/c)
|
||||
log2(z/c) = poly(z/c - 1)
|
||||
|
||||
where z is in [1.6p-1; 1.6p0] which is split into N subintervals and z falls
|
||||
into the ith one, then table entries are computed as
|
||||
|
||||
tab[i].invc = 1/c
|
||||
tab[i].logc = (double)log2(c)
|
||||
tab2[i].chi = (double)c
|
||||
tab2[i].clo = (double)(c - (double)c)
|
||||
|
||||
where c is near the center of the subinterval and is chosen by trying +-2^29
|
||||
floating point invc candidates around 1/center and selecting one for which
|
||||
|
||||
1) the rounding error in 0x1.8p10 + logc is 0,
|
||||
2) the rounding error in z - chi - clo is < 0x1p-64 and
|
||||
3) the rounding error in (double)log2(c) is minimized (< 0x1p-68).
|
||||
|
||||
Note: 1) ensures that k + logc can be computed without rounding error, 2)
|
||||
ensures that z/c - 1 can be computed as (z - chi - clo)*invc with close to a
|
||||
single rounding error when there is no fast fma for z*invc - 1, 3) ensures
|
||||
that logc + poly(z/c - 1) has small error, however near x == 1 when
|
||||
|log2(x)| < 0x1p-4, this is not enough so that is special cased. */
|
||||
.tab = {
|
||||
{0x1.724286bb1acf8p+0, -0x1.1095feecdb000p-1},
|
||||
{0x1.6e1f766d2cca1p+0, -0x1.08494bd76d000p-1},
|
||||
{0x1.6a13d0e30d48ap+0, -0x1.00143aee8f800p-1},
|
||||
{0x1.661ec32d06c85p+0, -0x1.efec5360b4000p-2},
|
||||
{0x1.623fa951198f8p+0, -0x1.dfdd91ab7e000p-2},
|
||||
{0x1.5e75ba4cf026cp+0, -0x1.cffae0cc79000p-2},
|
||||
{0x1.5ac055a214fb8p+0, -0x1.c043811fda000p-2},
|
||||
{0x1.571ed0f166e1ep+0, -0x1.b0b67323ae000p-2},
|
||||
{0x1.53909590bf835p+0, -0x1.a152f5a2db000p-2},
|
||||
{0x1.5014fed61adddp+0, -0x1.9217f5af86000p-2},
|
||||
{0x1.4cab88e487bd0p+0, -0x1.8304db0719000p-2},
|
||||
{0x1.49539b4334feep+0, -0x1.74189f9a9e000p-2},
|
||||
{0x1.460cbdfafd569p+0, -0x1.6552bb5199000p-2},
|
||||
{0x1.42d664ee4b953p+0, -0x1.56b23a29b1000p-2},
|
||||
{0x1.3fb01111dd8a6p+0, -0x1.483650f5fa000p-2},
|
||||
{0x1.3c995b70c5836p+0, -0x1.39de937f6a000p-2},
|
||||
{0x1.3991c4ab6fd4ap+0, -0x1.2baa1538d6000p-2},
|
||||
{0x1.3698e0ce099b5p+0, -0x1.1d98340ca4000p-2},
|
||||
{0x1.33ae48213e7b2p+0, -0x1.0fa853a40e000p-2},
|
||||
{0x1.30d191985bdb1p+0, -0x1.01d9c32e73000p-2},
|
||||
{0x1.2e025cab271d7p+0, -0x1.e857da2fa6000p-3},
|
||||
{0x1.2b404cf13cd82p+0, -0x1.cd3c8633d8000p-3},
|
||||
{0x1.288b02c7ccb50p+0, -0x1.b26034c14a000p-3},
|
||||
{0x1.25e2263944de5p+0, -0x1.97c1c2f4fe000p-3},
|
||||
{0x1.234563d8615b1p+0, -0x1.7d6023f800000p-3},
|
||||
{0x1.20b46e33eaf38p+0, -0x1.633a71a05e000p-3},
|
||||
{0x1.1e2eefdcda3ddp+0, -0x1.494f5e9570000p-3},
|
||||
{0x1.1bb4a580b3930p+0, -0x1.2f9e424e0a000p-3},
|
||||
{0x1.19453847f2200p+0, -0x1.162595afdc000p-3},
|
||||
{0x1.16e06c0d5d73cp+0, -0x1.f9c9a75bd8000p-4},
|
||||
{0x1.1485f47b7e4c2p+0, -0x1.c7b575bf9c000p-4},
|
||||
{0x1.12358ad0085d1p+0, -0x1.960c60ff48000p-4},
|
||||
{0x1.0fef00f532227p+0, -0x1.64ce247b60000p-4},
|
||||
{0x1.0db2077d03a8fp+0, -0x1.33f78b2014000p-4},
|
||||
{0x1.0b7e6d65980d9p+0, -0x1.0387d1a42c000p-4},
|
||||
{0x1.0953efe7b408dp+0, -0x1.a6f9208b50000p-5},
|
||||
{0x1.07325cac53b83p+0, -0x1.47a954f770000p-5},
|
||||
{0x1.05197e40d1b5cp+0, -0x1.d23a8c50c0000p-6},
|
||||
{0x1.03091c1208ea2p+0, -0x1.16a2629780000p-6},
|
||||
{0x1.0101025b37e21p+0, -0x1.720f8d8e80000p-8},
|
||||
{0x1.fc07ef9caa76bp-1, 0x1.6fe53b1500000p-7},
|
||||
{0x1.f4465d3f6f184p-1, 0x1.11ccce10f8000p-5},
|
||||
{0x1.ecc079f84107fp-1, 0x1.c4dfc8c8b8000p-5},
|
||||
{0x1.e573a99975ae8p-1, 0x1.3aa321e574000p-4},
|
||||
{0x1.de5d6f0bd3de6p-1, 0x1.918a0d08b8000p-4},
|
||||
{0x1.d77b681ff38b3p-1, 0x1.e72e9da044000p-4},
|
||||
{0x1.d0cb5724de943p-1, 0x1.1dcd2507f6000p-3},
|
||||
{0x1.ca4b2dc0e7563p-1, 0x1.476ab03dea000p-3},
|
||||
{0x1.c3f8ee8d6cb51p-1, 0x1.7074377e22000p-3},
|
||||
{0x1.bdd2b4f020c4cp-1, 0x1.98ede8ba94000p-3},
|
||||
{0x1.b7d6c006015cap-1, 0x1.c0db86ad2e000p-3},
|
||||
{0x1.b20366e2e338fp-1, 0x1.e840aafcee000p-3},
|
||||
{0x1.ac57026295039p-1, 0x1.0790ab4678000p-2},
|
||||
{0x1.a6d01bc2731ddp-1, 0x1.1ac056801c000p-2},
|
||||
{0x1.a16d3bc3ff18bp-1, 0x1.2db11d4fee000p-2},
|
||||
{0x1.9c2d14967feadp-1, 0x1.406464ec58000p-2},
|
||||
{0x1.970e4f47c9902p-1, 0x1.52dbe093af000p-2},
|
||||
{0x1.920fb3982bcf2p-1, 0x1.651902050d000p-2},
|
||||
{0x1.8d30187f759f1p-1, 0x1.771d2cdeaf000p-2},
|
||||
{0x1.886e5ebb9f66dp-1, 0x1.88e9c857d9000p-2},
|
||||
{0x1.83c97b658b994p-1, 0x1.9a80155e16000p-2},
|
||||
{0x1.7f405ffc61022p-1, 0x1.abe186ed3d000p-2},
|
||||
{0x1.7ad22181415cap-1, 0x1.bd0f2aea0e000p-2},
|
||||
{0x1.767dcf99eff8cp-1, 0x1.ce0a43dbf4000p-2},
|
||||
},
|
||||
#if !__FP_FAST_FMA
|
||||
.tab2 = {
|
||||
{0x1.6200012b90a8ep-1, 0x1.904ab0644b605p-55},
|
||||
{0x1.66000045734a6p-1, 0x1.1ff9bea62f7a9p-57},
|
||||
{0x1.69fffc325f2c5p-1, 0x1.27ecfcb3c90bap-55},
|
||||
{0x1.6e00038b95a04p-1, 0x1.8ff8856739326p-55},
|
||||
{0x1.71fffe09994e3p-1, 0x1.afd40275f82b1p-55},
|
||||
{0x1.7600015590e1p-1, -0x1.2fd75b4238341p-56},
|
||||
{0x1.7a00012655bd5p-1, 0x1.808e67c242b76p-56},
|
||||
{0x1.7e0003259e9a6p-1, -0x1.208e426f622b7p-57},
|
||||
{0x1.81fffedb4b2d2p-1, -0x1.402461ea5c92fp-55},
|
||||
{0x1.860002dfafcc3p-1, 0x1.df7f4a2f29a1fp-57},
|
||||
{0x1.89ffff78c6b5p-1, -0x1.e0453094995fdp-55},
|
||||
{0x1.8e00039671566p-1, -0x1.a04f3bec77b45p-55},
|
||||
{0x1.91fffe2bf1745p-1, -0x1.7fa34400e203cp-56},
|
||||
{0x1.95fffcc5c9fd1p-1, -0x1.6ff8005a0695dp-56},
|
||||
{0x1.9a0003bba4767p-1, 0x1.0f8c4c4ec7e03p-56},
|
||||
{0x1.9dfffe7b92da5p-1, 0x1.e7fd9478c4602p-55},
|
||||
{0x1.a1fffd72efdafp-1, -0x1.a0c554dcdae7ep-57},
|
||||
{0x1.a5fffde04ff95p-1, 0x1.67da98ce9b26bp-55},
|
||||
{0x1.a9fffca5e8d2bp-1, -0x1.284c9b54c13dep-55},
|
||||
{0x1.adfffddad03eap-1, 0x1.812c8ea602e3cp-58},
|
||||
{0x1.b1ffff10d3d4dp-1, -0x1.efaddad27789cp-55},
|
||||
{0x1.b5fffce21165ap-1, 0x1.3cb1719c61237p-58},
|
||||
{0x1.b9fffd950e674p-1, 0x1.3f7d94194cep-56},
|
||||
{0x1.be000139ca8afp-1, 0x1.50ac4215d9bcp-56},
|
||||
{0x1.c20005b46df99p-1, 0x1.beea653e9c1c9p-57},
|
||||
{0x1.c600040b9f7aep-1, -0x1.c079f274a70d6p-56},
|
||||
{0x1.ca0006255fd8ap-1, -0x1.a0b4076e84c1fp-56},
|
||||
{0x1.cdfffd94c095dp-1, 0x1.8f933f99ab5d7p-55},
|
||||
{0x1.d1ffff975d6cfp-1, -0x1.82c08665fe1bep-58},
|
||||
{0x1.d5fffa2561c93p-1, -0x1.b04289bd295f3p-56},
|
||||
{0x1.d9fff9d228b0cp-1, 0x1.70251340fa236p-55},
|
||||
{0x1.de00065bc7e16p-1, -0x1.5011e16a4d80cp-56},
|
||||
{0x1.e200002f64791p-1, 0x1.9802f09ef62ep-55},
|
||||
{0x1.e600057d7a6d8p-1, -0x1.e0b75580cf7fap-56},
|
||||
{0x1.ea00027edc00cp-1, -0x1.c848309459811p-55},
|
||||
{0x1.ee0006cf5cb7cp-1, -0x1.f8027951576f4p-55},
|
||||
{0x1.f2000782b7dccp-1, -0x1.f81d97274538fp-55},
|
||||
{0x1.f6000260c450ap-1, -0x1.071002727ffdcp-59},
|
||||
{0x1.f9fffe88cd533p-1, -0x1.81bdce1fda8bp-58},
|
||||
{0x1.fdfffd50f8689p-1, 0x1.7f91acb918e6ep-55},
|
||||
{0x1.0200004292367p+0, 0x1.b7ff365324681p-54},
|
||||
{0x1.05fffe3e3d668p+0, 0x1.6fa08ddae957bp-55},
|
||||
{0x1.0a0000a85a757p+0, -0x1.7e2de80d3fb91p-58},
|
||||
{0x1.0e0001a5f3fccp+0, -0x1.1823305c5f014p-54},
|
||||
{0x1.11ffff8afbaf5p+0, -0x1.bfabb6680bac2p-55},
|
||||
{0x1.15fffe54d91adp+0, -0x1.d7f121737e7efp-54},
|
||||
{0x1.1a00011ac36e1p+0, 0x1.c000a0516f5ffp-54},
|
||||
{0x1.1e00019c84248p+0, -0x1.082fbe4da5dap-54},
|
||||
{0x1.220000ffe5e6ep+0, -0x1.8fdd04c9cfb43p-55},
|
||||
{0x1.26000269fd891p+0, 0x1.cfe2a7994d182p-55},
|
||||
{0x1.2a00029a6e6dap+0, -0x1.00273715e8bc5p-56},
|
||||
{0x1.2dfffe0293e39p+0, 0x1.b7c39dab2a6f9p-54},
|
||||
{0x1.31ffff7dcf082p+0, 0x1.df1336edc5254p-56},
|
||||
{0x1.35ffff05a8b6p+0, -0x1.e03564ccd31ebp-54},
|
||||
{0x1.3a0002e0eaeccp+0, 0x1.5f0e74bd3a477p-56},
|
||||
{0x1.3e000043bb236p+0, 0x1.c7dcb149d8833p-54},
|
||||
{0x1.4200002d187ffp+0, 0x1.e08afcf2d3d28p-56},
|
||||
{0x1.460000d387cb1p+0, 0x1.20837856599a6p-55},
|
||||
{0x1.4a00004569f89p+0, -0x1.9fa5c904fbcd2p-55},
|
||||
{0x1.4e000043543f3p+0, -0x1.81125ed175329p-56},
|
||||
{0x1.51fffcc027f0fp+0, 0x1.883d8847754dcp-54},
|
||||
{0x1.55ffffd87b36fp+0, -0x1.709e731d02807p-55},
|
||||
{0x1.59ffff21df7bap+0, 0x1.7f79f68727b02p-55},
|
||||
{0x1.5dfffebfc3481p+0, -0x1.180902e30e93ep-54},
|
||||
},
|
||||
#endif
|
||||
};
|
||||
|
|
@ -0,0 +1,28 @@
|
|||
/*
|
||||
* Copyright (c) 2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
#ifndef _LOG2_DATA_H
|
||||
#define _LOG2_DATA_H
|
||||
|
||||
#include <features.h>
|
||||
|
||||
#define LOG2_TABLE_BITS 6
|
||||
#define LOG2_POLY_ORDER 7
|
||||
#define LOG2_POLY1_ORDER 11
|
||||
extern const struct log2_data {
|
||||
double invln2hi;
|
||||
double invln2lo;
|
||||
double poly[LOG2_POLY_ORDER - 1];
|
||||
double poly1[LOG2_POLY1_ORDER - 1];
|
||||
struct {
|
||||
double invc, logc;
|
||||
} tab[1 << LOG2_TABLE_BITS];
|
||||
#if !__FP_FAST_FMA
|
||||
struct {
|
||||
double chi, clo;
|
||||
} tab2[1 << LOG2_TABLE_BITS];
|
||||
#endif
|
||||
} __log2_data;
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,72 @@
|
|||
/*
|
||||
* Single-precision log2 function.
|
||||
*
|
||||
* Copyright (c) 2017-2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include <stdint.h>
|
||||
#include "libm.h"
|
||||
#include "log2f_data.h"
|
||||
|
||||
/*
|
||||
LOG2F_TABLE_BITS = 4
|
||||
LOG2F_POLY_ORDER = 4
|
||||
|
||||
ULP error: 0.752 (nearest rounding.)
|
||||
Relative error: 1.9 * 2^-26 (before rounding.)
|
||||
*/
|
||||
|
||||
#define N (1 << LOG2F_TABLE_BITS)
|
||||
#define T __log2f_data.tab
|
||||
#define A __log2f_data.poly
|
||||
#define OFF 0x3f330000
|
||||
|
||||
float log2f(float x)
|
||||
{
|
||||
double_t z, r, r2, p, y, y0, invc, logc;
|
||||
uint32_t ix, iz, top, tmp;
|
||||
int k, i;
|
||||
|
||||
ix = asuint(x);
|
||||
/* Fix sign of zero with downward rounding when x==1. */
|
||||
if (WANT_ROUNDING && predict_false(ix == 0x3f800000))
|
||||
return 0;
|
||||
if (predict_false(ix - 0x00800000 >= 0x7f800000 - 0x00800000)) {
|
||||
/* x < 0x1p-126 or inf or nan. */
|
||||
if (ix * 2 == 0)
|
||||
return __math_divzerof(1);
|
||||
if (ix == 0x7f800000) /* log2(inf) == inf. */
|
||||
return x;
|
||||
if ((ix & 0x80000000) || ix * 2 >= 0xff000000)
|
||||
return __math_invalidf(x);
|
||||
/* x is subnormal, normalize it. */
|
||||
ix = asuint(x * 0x1p23f);
|
||||
ix -= 23 << 23;
|
||||
}
|
||||
|
||||
/* x = 2^k z; where z is in range [OFF,2*OFF] and exact.
|
||||
The range is split into N subintervals.
|
||||
The ith subinterval contains z and c is near its center. */
|
||||
tmp = ix - OFF;
|
||||
i = (tmp >> (23 - LOG2F_TABLE_BITS)) % N;
|
||||
top = tmp & 0xff800000;
|
||||
iz = ix - top;
|
||||
k = (int32_t)tmp >> 23; /* arithmetic shift */
|
||||
invc = T[i].invc;
|
||||
logc = T[i].logc;
|
||||
z = (double_t)asfloat(iz);
|
||||
|
||||
/* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k */
|
||||
r = z * invc - 1;
|
||||
y0 = logc + (double_t)k;
|
||||
|
||||
/* Pipelined polynomial evaluation to approximate log1p(r)/ln2. */
|
||||
r2 = r * r;
|
||||
y = A[1] * r + A[2];
|
||||
y = A[0] * r2 + y;
|
||||
p = A[3] * r + y0;
|
||||
y = y * r2 + p;
|
||||
return eval_as_float(y);
|
||||
}
|
||||
|
|
@ -0,0 +1,33 @@
|
|||
/*
|
||||
* Data definition for log2f.
|
||||
*
|
||||
* Copyright (c) 2017-2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
|
||||
#include "log2f_data.h"
|
||||
|
||||
const struct log2f_data __log2f_data = {
|
||||
.tab = {
|
||||
{ 0x1.661ec79f8f3bep+0, -0x1.efec65b963019p-2 },
|
||||
{ 0x1.571ed4aaf883dp+0, -0x1.b0b6832d4fca4p-2 },
|
||||
{ 0x1.49539f0f010bp+0, -0x1.7418b0a1fb77bp-2 },
|
||||
{ 0x1.3c995b0b80385p+0, -0x1.39de91a6dcf7bp-2 },
|
||||
{ 0x1.30d190c8864a5p+0, -0x1.01d9bf3f2b631p-2 },
|
||||
{ 0x1.25e227b0b8eap+0, -0x1.97c1d1b3b7afp-3 },
|
||||
{ 0x1.1bb4a4a1a343fp+0, -0x1.2f9e393af3c9fp-3 },
|
||||
{ 0x1.12358f08ae5bap+0, -0x1.960cbbf788d5cp-4 },
|
||||
{ 0x1.0953f419900a7p+0, -0x1.a6f9db6475fcep-5 },
|
||||
{ 0x1p+0, 0x0p+0 },
|
||||
{ 0x1.e608cfd9a47acp-1, 0x1.338ca9f24f53dp-4 },
|
||||
{ 0x1.ca4b31f026aap-1, 0x1.476a9543891bap-3 },
|
||||
{ 0x1.b2036576afce6p-1, 0x1.e840b4ac4e4d2p-3 },
|
||||
{ 0x1.9c2d163a1aa2dp-1, 0x1.40645f0c6651cp-2 },
|
||||
{ 0x1.886e6037841edp-1, 0x1.88e9c2c1b9ff8p-2 },
|
||||
{ 0x1.767dcf5534862p-1, 0x1.ce0a44eb17bccp-2 },
|
||||
},
|
||||
.poly = {
|
||||
-0x1.712b6f70a7e4dp-2, 0x1.ecabf496832ep-2, -0x1.715479ffae3dep-1,
|
||||
0x1.715475f35c8b8p0,
|
||||
}
|
||||
};
|
||||
|
|
@ -0,0 +1,19 @@
|
|||
/*
|
||||
* Copyright (c) 2017-2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
#ifndef _LOG2F_DATA_H
|
||||
#define _LOG2F_DATA_H
|
||||
|
||||
#include <features.h>
|
||||
|
||||
#define LOG2F_TABLE_BITS 4
|
||||
#define LOG2F_POLY_ORDER 4
|
||||
extern const struct log2f_data {
|
||||
struct {
|
||||
double invc, logc;
|
||||
} tab[1 << LOG2F_TABLE_BITS];
|
||||
double poly[LOG2F_POLY_ORDER];
|
||||
} __log2f_data;
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,328 @@
|
|||
/*
|
||||
* Data for log.
|
||||
*
|
||||
* Copyright (c) 2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
|
||||
#include "log_data.h"
|
||||
|
||||
#define N (1 << LOG_TABLE_BITS)
|
||||
|
||||
const struct log_data __log_data = {
|
||||
.ln2hi = 0x1.62e42fefa3800p-1,
|
||||
.ln2lo = 0x1.ef35793c76730p-45,
|
||||
.poly1 = {
|
||||
// relative error: 0x1.c04d76cp-63
|
||||
// in -0x1p-4 0x1.09p-4 (|log(1+x)| > 0x1p-4 outside the interval)
|
||||
-0x1p-1,
|
||||
0x1.5555555555577p-2,
|
||||
-0x1.ffffffffffdcbp-3,
|
||||
0x1.999999995dd0cp-3,
|
||||
-0x1.55555556745a7p-3,
|
||||
0x1.24924a344de3p-3,
|
||||
-0x1.fffffa4423d65p-4,
|
||||
0x1.c7184282ad6cap-4,
|
||||
-0x1.999eb43b068ffp-4,
|
||||
0x1.78182f7afd085p-4,
|
||||
-0x1.5521375d145cdp-4,
|
||||
},
|
||||
.poly = {
|
||||
// relative error: 0x1.926199e8p-56
|
||||
// abs error: 0x1.882ff33p-65
|
||||
// in -0x1.fp-9 0x1.fp-9
|
||||
-0x1.0000000000001p-1,
|
||||
0x1.555555551305bp-2,
|
||||
-0x1.fffffffeb459p-3,
|
||||
0x1.999b324f10111p-3,
|
||||
-0x1.55575e506c89fp-3,
|
||||
},
|
||||
/* Algorithm:
|
||||
|
||||
x = 2^k z
|
||||
log(x) = k ln2 + log(c) + log(z/c)
|
||||
log(z/c) = poly(z/c - 1)
|
||||
|
||||
where z is in [1.6p-1; 1.6p0] which is split into N subintervals and z falls
|
||||
into the ith one, then table entries are computed as
|
||||
|
||||
tab[i].invc = 1/c
|
||||
tab[i].logc = (double)log(c)
|
||||
tab2[i].chi = (double)c
|
||||
tab2[i].clo = (double)(c - (double)c)
|
||||
|
||||
where c is near the center of the subinterval and is chosen by trying +-2^29
|
||||
floating point invc candidates around 1/center and selecting one for which
|
||||
|
||||
1) the rounding error in 0x1.8p9 + logc is 0,
|
||||
2) the rounding error in z - chi - clo is < 0x1p-66 and
|
||||
3) the rounding error in (double)log(c) is minimized (< 0x1p-66).
|
||||
|
||||
Note: 1) ensures that k*ln2hi + logc can be computed without rounding error,
|
||||
2) ensures that z/c - 1 can be computed as (z - chi - clo)*invc with close to
|
||||
a single rounding error when there is no fast fma for z*invc - 1, 3) ensures
|
||||
that logc + poly(z/c - 1) has small error, however near x == 1 when
|
||||
|log(x)| < 0x1p-4, this is not enough so that is special cased. */
|
||||
.tab = {
|
||||
{0x1.734f0c3e0de9fp+0, -0x1.7cc7f79e69000p-2},
|
||||
{0x1.713786a2ce91fp+0, -0x1.76feec20d0000p-2},
|
||||
{0x1.6f26008fab5a0p+0, -0x1.713e31351e000p-2},
|
||||
{0x1.6d1a61f138c7dp+0, -0x1.6b85b38287800p-2},
|
||||
{0x1.6b1490bc5b4d1p+0, -0x1.65d5590807800p-2},
|
||||
{0x1.69147332f0cbap+0, -0x1.602d076180000p-2},
|
||||
{0x1.6719f18224223p+0, -0x1.5a8ca86909000p-2},
|
||||
{0x1.6524f99a51ed9p+0, -0x1.54f4356035000p-2},
|
||||
{0x1.63356aa8f24c4p+0, -0x1.4f637c36b4000p-2},
|
||||
{0x1.614b36b9ddc14p+0, -0x1.49da7fda85000p-2},
|
||||
{0x1.5f66452c65c4cp+0, -0x1.445923989a800p-2},
|
||||
{0x1.5d867b5912c4fp+0, -0x1.3edf439b0b800p-2},
|
||||
{0x1.5babccb5b90dep+0, -0x1.396ce448f7000p-2},
|
||||
{0x1.59d61f2d91a78p+0, -0x1.3401e17bda000p-2},
|
||||
{0x1.5805612465687p+0, -0x1.2e9e2ef468000p-2},
|
||||
{0x1.56397cee76bd3p+0, -0x1.2941b3830e000p-2},
|
||||
{0x1.54725e2a77f93p+0, -0x1.23ec58cda8800p-2},
|
||||
{0x1.52aff42064583p+0, -0x1.1e9e129279000p-2},
|
||||
{0x1.50f22dbb2bddfp+0, -0x1.1956d2b48f800p-2},
|
||||
{0x1.4f38f4734ded7p+0, -0x1.141679ab9f800p-2},
|
||||
{0x1.4d843cfde2840p+0, -0x1.0edd094ef9800p-2},
|
||||
{0x1.4bd3ec078a3c8p+0, -0x1.09aa518db1000p-2},
|
||||
{0x1.4a27fc3e0258ap+0, -0x1.047e65263b800p-2},
|
||||
{0x1.4880524d48434p+0, -0x1.feb224586f000p-3},
|
||||
{0x1.46dce1b192d0bp+0, -0x1.f474a7517b000p-3},
|
||||
{0x1.453d9d3391854p+0, -0x1.ea4443d103000p-3},
|
||||
{0x1.43a2744b4845ap+0, -0x1.e020d44e9b000p-3},
|
||||
{0x1.420b54115f8fbp+0, -0x1.d60a22977f000p-3},
|
||||
{0x1.40782da3ef4b1p+0, -0x1.cc00104959000p-3},
|
||||
{0x1.3ee8f5d57fe8fp+0, -0x1.c202956891000p-3},
|
||||
{0x1.3d5d9a00b4ce9p+0, -0x1.b81178d811000p-3},
|
||||
{0x1.3bd60c010c12bp+0, -0x1.ae2c9ccd3d000p-3},
|
||||
{0x1.3a5242b75dab8p+0, -0x1.a45402e129000p-3},
|
||||
{0x1.38d22cd9fd002p+0, -0x1.9a877681df000p-3},
|
||||
{0x1.3755bc5847a1cp+0, -0x1.90c6d69483000p-3},
|
||||
{0x1.35dce49ad36e2p+0, -0x1.87120a645c000p-3},
|
||||
{0x1.34679984dd440p+0, -0x1.7d68fb4143000p-3},
|
||||
{0x1.32f5cceffcb24p+0, -0x1.73cb83c627000p-3},
|
||||
{0x1.3187775a10d49p+0, -0x1.6a39a9b376000p-3},
|
||||
{0x1.301c8373e3990p+0, -0x1.60b3154b7a000p-3},
|
||||
{0x1.2eb4ebb95f841p+0, -0x1.5737d76243000p-3},
|
||||
{0x1.2d50a0219a9d1p+0, -0x1.4dc7b8fc23000p-3},
|
||||
{0x1.2bef9a8b7fd2ap+0, -0x1.4462c51d20000p-3},
|
||||
{0x1.2a91c7a0c1babp+0, -0x1.3b08abc830000p-3},
|
||||
{0x1.293726014b530p+0, -0x1.31b996b490000p-3},
|
||||
{0x1.27dfa5757a1f5p+0, -0x1.2875490a44000p-3},
|
||||
{0x1.268b39b1d3bbfp+0, -0x1.1f3b9f879a000p-3},
|
||||
{0x1.2539d838ff5bdp+0, -0x1.160c8252ca000p-3},
|
||||
{0x1.23eb7aac9083bp+0, -0x1.0ce7f57f72000p-3},
|
||||
{0x1.22a012ba940b6p+0, -0x1.03cdc49fea000p-3},
|
||||
{0x1.2157996cc4132p+0, -0x1.f57bdbc4b8000p-4},
|
||||
{0x1.201201dd2fc9bp+0, -0x1.e370896404000p-4},
|
||||
{0x1.1ecf4494d480bp+0, -0x1.d17983ef94000p-4},
|
||||
{0x1.1d8f5528f6569p+0, -0x1.bf9674ed8a000p-4},
|
||||
{0x1.1c52311577e7cp+0, -0x1.adc79202f6000p-4},
|
||||
{0x1.1b17c74cb26e9p+0, -0x1.9c0c3e7288000p-4},
|
||||
{0x1.19e010c2c1ab6p+0, -0x1.8a646b372c000p-4},
|
||||
{0x1.18ab07bb670bdp+0, -0x1.78d01b3ac0000p-4},
|
||||
{0x1.1778a25efbcb6p+0, -0x1.674f145380000p-4},
|
||||
{0x1.1648d354c31dap+0, -0x1.55e0e6d878000p-4},
|
||||
{0x1.151b990275fddp+0, -0x1.4485cdea1e000p-4},
|
||||
{0x1.13f0ea432d24cp+0, -0x1.333d94d6aa000p-4},
|
||||
{0x1.12c8b7210f9dap+0, -0x1.22079f8c56000p-4},
|
||||
{0x1.11a3028ecb531p+0, -0x1.10e4698622000p-4},
|
||||
{0x1.107fbda8434afp+0, -0x1.ffa6c6ad20000p-5},
|
||||
{0x1.0f5ee0f4e6bb3p+0, -0x1.dda8d4a774000p-5},
|
||||
{0x1.0e4065d2a9fcep+0, -0x1.bbcece4850000p-5},
|
||||
{0x1.0d244632ca521p+0, -0x1.9a1894012c000p-5},
|
||||
{0x1.0c0a77ce2981ap+0, -0x1.788583302c000p-5},
|
||||
{0x1.0af2f83c636d1p+0, -0x1.5715e67d68000p-5},
|
||||
{0x1.09ddb98a01339p+0, -0x1.35c8a49658000p-5},
|
||||
{0x1.08cabaf52e7dfp+0, -0x1.149e364154000p-5},
|
||||
{0x1.07b9f2f4e28fbp+0, -0x1.e72c082eb8000p-6},
|
||||
{0x1.06ab58c358f19p+0, -0x1.a55f152528000p-6},
|
||||
{0x1.059eea5ecf92cp+0, -0x1.63d62cf818000p-6},
|
||||
{0x1.04949cdd12c90p+0, -0x1.228fb8caa0000p-6},
|
||||
{0x1.038c6c6f0ada9p+0, -0x1.c317b20f90000p-7},
|
||||
{0x1.02865137932a9p+0, -0x1.419355daa0000p-7},
|
||||
{0x1.0182427ea7348p+0, -0x1.81203c2ec0000p-8},
|
||||
{0x1.008040614b195p+0, -0x1.0040979240000p-9},
|
||||
{0x1.fe01ff726fa1ap-1, 0x1.feff384900000p-9},
|
||||
{0x1.fa11cc261ea74p-1, 0x1.7dc41353d0000p-7},
|
||||
{0x1.f6310b081992ep-1, 0x1.3cea3c4c28000p-6},
|
||||
{0x1.f25f63ceeadcdp-1, 0x1.b9fc114890000p-6},
|
||||
{0x1.ee9c8039113e7p-1, 0x1.1b0d8ce110000p-5},
|
||||
{0x1.eae8078cbb1abp-1, 0x1.58a5bd001c000p-5},
|
||||
{0x1.e741aa29d0c9bp-1, 0x1.95c8340d88000p-5},
|
||||
{0x1.e3a91830a99b5p-1, 0x1.d276aef578000p-5},
|
||||
{0x1.e01e009609a56p-1, 0x1.07598e598c000p-4},
|
||||
{0x1.dca01e577bb98p-1, 0x1.253f5e30d2000p-4},
|
||||
{0x1.d92f20b7c9103p-1, 0x1.42edd8b380000p-4},
|
||||
{0x1.d5cac66fb5ccep-1, 0x1.606598757c000p-4},
|
||||
{0x1.d272caa5ede9dp-1, 0x1.7da76356a0000p-4},
|
||||
{0x1.cf26e3e6b2ccdp-1, 0x1.9ab434e1c6000p-4},
|
||||
{0x1.cbe6da2a77902p-1, 0x1.b78c7bb0d6000p-4},
|
||||
{0x1.c8b266d37086dp-1, 0x1.d431332e72000p-4},
|
||||
{0x1.c5894bd5d5804p-1, 0x1.f0a3171de6000p-4},
|
||||
{0x1.c26b533bb9f8cp-1, 0x1.067152b914000p-3},
|
||||
{0x1.bf583eeece73fp-1, 0x1.147858292b000p-3},
|
||||
{0x1.bc4fd75db96c1p-1, 0x1.2266ecdca3000p-3},
|
||||
{0x1.b951e0c864a28p-1, 0x1.303d7a6c55000p-3},
|
||||
{0x1.b65e2c5ef3e2cp-1, 0x1.3dfc33c331000p-3},
|
||||
{0x1.b374867c9888bp-1, 0x1.4ba366b7a8000p-3},
|
||||
{0x1.b094b211d304ap-1, 0x1.5933928d1f000p-3},
|
||||
{0x1.adbe885f2ef7ep-1, 0x1.66acd2418f000p-3},
|
||||
{0x1.aaf1d31603da2p-1, 0x1.740f8ec669000p-3},
|
||||
{0x1.a82e63fd358a7p-1, 0x1.815c0f51af000p-3},
|
||||
{0x1.a5740ef09738bp-1, 0x1.8e92954f68000p-3},
|
||||
{0x1.a2c2a90ab4b27p-1, 0x1.9bb3602f84000p-3},
|
||||
{0x1.a01a01393f2d1p-1, 0x1.a8bed1c2c0000p-3},
|
||||
{0x1.9d79f24db3c1bp-1, 0x1.b5b515c01d000p-3},
|
||||
{0x1.9ae2505c7b190p-1, 0x1.c2967ccbcc000p-3},
|
||||
{0x1.9852ef297ce2fp-1, 0x1.cf635d5486000p-3},
|
||||
{0x1.95cbaeea44b75p-1, 0x1.dc1bd3446c000p-3},
|
||||
{0x1.934c69de74838p-1, 0x1.e8c01b8cfe000p-3},
|
||||
{0x1.90d4f2f6752e6p-1, 0x1.f5509c0179000p-3},
|
||||
{0x1.8e6528effd79dp-1, 0x1.00e6c121fb800p-2},
|
||||
{0x1.8bfce9fcc007cp-1, 0x1.071b80e93d000p-2},
|
||||
{0x1.899c0dabec30ep-1, 0x1.0d46b9e867000p-2},
|
||||
{0x1.87427aa2317fbp-1, 0x1.13687334bd000p-2},
|
||||
{0x1.84f00acb39a08p-1, 0x1.1980d67234800p-2},
|
||||
{0x1.82a49e8653e55p-1, 0x1.1f8ffe0cc8000p-2},
|
||||
{0x1.8060195f40260p-1, 0x1.2595fd7636800p-2},
|
||||
{0x1.7e22563e0a329p-1, 0x1.2b9300914a800p-2},
|
||||
{0x1.7beb377dcb5adp-1, 0x1.3187210436000p-2},
|
||||
{0x1.79baa679725c2p-1, 0x1.377266dec1800p-2},
|
||||
{0x1.77907f2170657p-1, 0x1.3d54ffbaf3000p-2},
|
||||
{0x1.756cadbd6130cp-1, 0x1.432eee32fe000p-2},
|
||||
},
|
||||
#if !__FP_FAST_FMA
|
||||
.tab2 = {
|
||||
{0x1.61000014fb66bp-1, 0x1.e026c91425b3cp-56},
|
||||
{0x1.63000034db495p-1, 0x1.dbfea48005d41p-55},
|
||||
{0x1.650000d94d478p-1, 0x1.e7fa786d6a5b7p-55},
|
||||
{0x1.67000074e6fadp-1, 0x1.1fcea6b54254cp-57},
|
||||
{0x1.68ffffedf0faep-1, -0x1.c7e274c590efdp-56},
|
||||
{0x1.6b0000763c5bcp-1, -0x1.ac16848dcda01p-55},
|
||||
{0x1.6d0001e5cc1f6p-1, 0x1.33f1c9d499311p-55},
|
||||
{0x1.6efffeb05f63ep-1, -0x1.e80041ae22d53p-56},
|
||||
{0x1.710000e86978p-1, 0x1.bff6671097952p-56},
|
||||
{0x1.72ffffc67e912p-1, 0x1.c00e226bd8724p-55},
|
||||
{0x1.74fffdf81116ap-1, -0x1.e02916ef101d2p-57},
|
||||
{0x1.770000f679c9p-1, -0x1.7fc71cd549c74p-57},
|
||||
{0x1.78ffffa7ec835p-1, 0x1.1bec19ef50483p-55},
|
||||
{0x1.7affffe20c2e6p-1, -0x1.07e1729cc6465p-56},
|
||||
{0x1.7cfffed3fc9p-1, -0x1.08072087b8b1cp-55},
|
||||
{0x1.7efffe9261a76p-1, 0x1.dc0286d9df9aep-55},
|
||||
{0x1.81000049ca3e8p-1, 0x1.97fd251e54c33p-55},
|
||||
{0x1.8300017932c8fp-1, -0x1.afee9b630f381p-55},
|
||||
{0x1.850000633739cp-1, 0x1.9bfbf6b6535bcp-55},
|
||||
{0x1.87000204289c6p-1, -0x1.bbf65f3117b75p-55},
|
||||
{0x1.88fffebf57904p-1, -0x1.9006ea23dcb57p-55},
|
||||
{0x1.8b00022bc04dfp-1, -0x1.d00df38e04b0ap-56},
|
||||
{0x1.8cfffe50c1b8ap-1, -0x1.8007146ff9f05p-55},
|
||||
{0x1.8effffc918e43p-1, 0x1.3817bd07a7038p-55},
|
||||
{0x1.910001efa5fc7p-1, 0x1.93e9176dfb403p-55},
|
||||
{0x1.9300013467bb9p-1, 0x1.f804e4b980276p-56},
|
||||
{0x1.94fffe6ee076fp-1, -0x1.f7ef0d9ff622ep-55},
|
||||
{0x1.96fffde3c12d1p-1, -0x1.082aa962638bap-56},
|
||||
{0x1.98ffff4458a0dp-1, -0x1.7801b9164a8efp-55},
|
||||
{0x1.9afffdd982e3ep-1, -0x1.740e08a5a9337p-55},
|
||||
{0x1.9cfffed49fb66p-1, 0x1.fce08c19bep-60},
|
||||
{0x1.9f00020f19c51p-1, -0x1.a3faa27885b0ap-55},
|
||||
{0x1.a10001145b006p-1, 0x1.4ff489958da56p-56},
|
||||
{0x1.a300007bbf6fap-1, 0x1.cbeab8a2b6d18p-55},
|
||||
{0x1.a500010971d79p-1, 0x1.8fecadd78793p-55},
|
||||
{0x1.a70001df52e48p-1, -0x1.f41763dd8abdbp-55},
|
||||
{0x1.a90001c593352p-1, -0x1.ebf0284c27612p-55},
|
||||
{0x1.ab0002a4f3e4bp-1, -0x1.9fd043cff3f5fp-57},
|
||||
{0x1.acfffd7ae1ed1p-1, -0x1.23ee7129070b4p-55},
|
||||
{0x1.aefffee510478p-1, 0x1.a063ee00edea3p-57},
|
||||
{0x1.b0fffdb650d5bp-1, 0x1.a06c8381f0ab9p-58},
|
||||
{0x1.b2ffffeaaca57p-1, -0x1.9011e74233c1dp-56},
|
||||
{0x1.b4fffd995badcp-1, -0x1.9ff1068862a9fp-56},
|
||||
{0x1.b7000249e659cp-1, 0x1.aff45d0864f3ep-55},
|
||||
{0x1.b8ffff987164p-1, 0x1.cfe7796c2c3f9p-56},
|
||||
{0x1.bafffd204cb4fp-1, -0x1.3ff27eef22bc4p-57},
|
||||
{0x1.bcfffd2415c45p-1, -0x1.cffb7ee3bea21p-57},
|
||||
{0x1.beffff86309dfp-1, -0x1.14103972e0b5cp-55},
|
||||
{0x1.c0fffe1b57653p-1, 0x1.bc16494b76a19p-55},
|
||||
{0x1.c2ffff1fa57e3p-1, -0x1.4feef8d30c6edp-57},
|
||||
{0x1.c4fffdcbfe424p-1, -0x1.43f68bcec4775p-55},
|
||||
{0x1.c6fffed54b9f7p-1, 0x1.47ea3f053e0ecp-55},
|
||||
{0x1.c8fffeb998fd5p-1, 0x1.383068df992f1p-56},
|
||||
{0x1.cb0002125219ap-1, -0x1.8fd8e64180e04p-57},
|
||||
{0x1.ccfffdd94469cp-1, 0x1.e7ebe1cc7ea72p-55},
|
||||
{0x1.cefffeafdc476p-1, 0x1.ebe39ad9f88fep-55},
|
||||
{0x1.d1000169af82bp-1, 0x1.57d91a8b95a71p-56},
|
||||
{0x1.d30000d0ff71dp-1, 0x1.9c1906970c7dap-55},
|
||||
{0x1.d4fffea790fc4p-1, -0x1.80e37c558fe0cp-58},
|
||||
{0x1.d70002edc87e5p-1, -0x1.f80d64dc10f44p-56},
|
||||
{0x1.d900021dc82aap-1, -0x1.47c8f94fd5c5cp-56},
|
||||
{0x1.dafffd86b0283p-1, 0x1.c7f1dc521617ep-55},
|
||||
{0x1.dd000296c4739p-1, 0x1.8019eb2ffb153p-55},
|
||||
{0x1.defffe54490f5p-1, 0x1.e00d2c652cc89p-57},
|
||||
{0x1.e0fffcdabf694p-1, -0x1.f8340202d69d2p-56},
|
||||
{0x1.e2fffdb52c8ddp-1, 0x1.b00c1ca1b0864p-56},
|
||||
{0x1.e4ffff24216efp-1, 0x1.2ffa8b094ab51p-56},
|
||||
{0x1.e6fffe88a5e11p-1, -0x1.7f673b1efbe59p-58},
|
||||
{0x1.e9000119eff0dp-1, -0x1.4808d5e0bc801p-55},
|
||||
{0x1.eafffdfa51744p-1, 0x1.80006d54320b5p-56},
|
||||
{0x1.ed0001a127fa1p-1, -0x1.002f860565c92p-58},
|
||||
{0x1.ef00007babcc4p-1, -0x1.540445d35e611p-55},
|
||||
{0x1.f0ffff57a8d02p-1, -0x1.ffb3139ef9105p-59},
|
||||
{0x1.f30001ee58ac7p-1, 0x1.a81acf2731155p-55},
|
||||
{0x1.f4ffff5823494p-1, 0x1.a3f41d4d7c743p-55},
|
||||
{0x1.f6ffffca94c6bp-1, -0x1.202f41c987875p-57},
|
||||
{0x1.f8fffe1f9c441p-1, 0x1.77dd1f477e74bp-56},
|
||||
{0x1.fafffd2e0e37ep-1, -0x1.f01199a7ca331p-57},
|
||||
{0x1.fd0001c77e49ep-1, 0x1.181ee4bceacb1p-56},
|
||||
{0x1.feffff7e0c331p-1, -0x1.e05370170875ap-57},
|
||||
{0x1.00ffff465606ep+0, -0x1.a7ead491c0adap-55},
|
||||
{0x1.02ffff3867a58p+0, -0x1.77f69c3fcb2ep-54},
|
||||
{0x1.04ffffdfc0d17p+0, 0x1.7bffe34cb945bp-54},
|
||||
{0x1.0700003cd4d82p+0, 0x1.20083c0e456cbp-55},
|
||||
{0x1.08ffff9f2cbe8p+0, -0x1.dffdfbe37751ap-57},
|
||||
{0x1.0b000010cda65p+0, -0x1.13f7faee626ebp-54},
|
||||
{0x1.0d00001a4d338p+0, 0x1.07dfa79489ff7p-55},
|
||||
{0x1.0effffadafdfdp+0, -0x1.7040570d66bcp-56},
|
||||
{0x1.110000bbafd96p+0, 0x1.e80d4846d0b62p-55},
|
||||
{0x1.12ffffae5f45dp+0, 0x1.dbffa64fd36efp-54},
|
||||
{0x1.150000dd59ad9p+0, 0x1.a0077701250aep-54},
|
||||
{0x1.170000f21559ap+0, 0x1.dfdf9e2e3deeep-55},
|
||||
{0x1.18ffffc275426p+0, 0x1.10030dc3b7273p-54},
|
||||
{0x1.1b000123d3c59p+0, 0x1.97f7980030188p-54},
|
||||
{0x1.1cffff8299eb7p+0, -0x1.5f932ab9f8c67p-57},
|
||||
{0x1.1effff48ad4p+0, 0x1.37fbf9da75bebp-54},
|
||||
{0x1.210000c8b86a4p+0, 0x1.f806b91fd5b22p-54},
|
||||
{0x1.2300003854303p+0, 0x1.3ffc2eb9fbf33p-54},
|
||||
{0x1.24fffffbcf684p+0, 0x1.601e77e2e2e72p-56},
|
||||
{0x1.26ffff52921d9p+0, 0x1.ffcbb767f0c61p-56},
|
||||
{0x1.2900014933a3cp+0, -0x1.202ca3c02412bp-56},
|
||||
{0x1.2b00014556313p+0, -0x1.2808233f21f02p-54},
|
||||
{0x1.2cfffebfe523bp+0, -0x1.8ff7e384fdcf2p-55},
|
||||
{0x1.2f0000bb8ad96p+0, -0x1.5ff51503041c5p-55},
|
||||
{0x1.30ffffb7ae2afp+0, -0x1.10071885e289dp-55},
|
||||
{0x1.32ffffeac5f7fp+0, -0x1.1ff5d3fb7b715p-54},
|
||||
{0x1.350000ca66756p+0, 0x1.57f82228b82bdp-54},
|
||||
{0x1.3700011fbf721p+0, 0x1.000bac40dd5ccp-55},
|
||||
{0x1.38ffff9592fb9p+0, -0x1.43f9d2db2a751p-54},
|
||||
{0x1.3b00004ddd242p+0, 0x1.57f6b707638e1p-55},
|
||||
{0x1.3cffff5b2c957p+0, 0x1.a023a10bf1231p-56},
|
||||
{0x1.3efffeab0b418p+0, 0x1.87f6d66b152bp-54},
|
||||
{0x1.410001532aff4p+0, 0x1.7f8375f198524p-57},
|
||||
{0x1.4300017478b29p+0, 0x1.301e672dc5143p-55},
|
||||
{0x1.44fffe795b463p+0, 0x1.9ff69b8b2895ap-55},
|
||||
{0x1.46fffe80475ep+0, -0x1.5c0b19bc2f254p-54},
|
||||
{0x1.48fffef6fc1e7p+0, 0x1.b4009f23a2a72p-54},
|
||||
{0x1.4afffe5bea704p+0, -0x1.4ffb7bf0d7d45p-54},
|
||||
{0x1.4d000171027dep+0, -0x1.9c06471dc6a3dp-54},
|
||||
{0x1.4f0000ff03ee2p+0, 0x1.77f890b85531cp-54},
|
||||
{0x1.5100012dc4bd1p+0, 0x1.004657166a436p-57},
|
||||
{0x1.530001605277ap+0, -0x1.6bfcece233209p-54},
|
||||
{0x1.54fffecdb704cp+0, -0x1.902720505a1d7p-55},
|
||||
{0x1.56fffef5f54a9p+0, 0x1.bbfe60ec96412p-54},
|
||||
{0x1.5900017e61012p+0, 0x1.87ec581afef9p-55},
|
||||
{0x1.5b00003c93e92p+0, -0x1.f41080abf0ccp-54},
|
||||
{0x1.5d0001d4919bcp+0, -0x1.8812afb254729p-54},
|
||||
{0x1.5efffe7b87a89p+0, -0x1.47eb780ed6904p-54},
|
||||
},
|
||||
#endif
|
||||
};
|
||||
|
|
@ -0,0 +1,28 @@
|
|||
/*
|
||||
* Copyright (c) 2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
#ifndef _LOG_DATA_H
|
||||
#define _LOG_DATA_H
|
||||
|
||||
#include <features.h>
|
||||
|
||||
#define LOG_TABLE_BITS 7
|
||||
#define LOG_POLY_ORDER 6
|
||||
#define LOG_POLY1_ORDER 12
|
||||
extern const struct log_data {
|
||||
double ln2hi;
|
||||
double ln2lo;
|
||||
double poly[LOG_POLY_ORDER - 1]; /* First coefficient is 1. */
|
||||
double poly1[LOG_POLY1_ORDER - 1];
|
||||
struct {
|
||||
double invc, logc;
|
||||
} tab[1 << LOG_TABLE_BITS];
|
||||
#if !__FP_FAST_FMA
|
||||
struct {
|
||||
double chi, clo;
|
||||
} tab2[1 << LOG_TABLE_BITS];
|
||||
#endif
|
||||
} __log_data;
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,71 @@
|
|||
/*
|
||||
* Single-precision log function.
|
||||
*
|
||||
* Copyright (c) 2017-2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include <stdint.h>
|
||||
#include "libm.h"
|
||||
#include "logf_data.h"
|
||||
|
||||
/*
|
||||
LOGF_TABLE_BITS = 4
|
||||
LOGF_POLY_ORDER = 4
|
||||
|
||||
ULP error: 0.818 (nearest rounding.)
|
||||
Relative error: 1.957 * 2^-26 (before rounding.)
|
||||
*/
|
||||
|
||||
#define T __logf_data.tab
|
||||
#define A __logf_data.poly
|
||||
#define Ln2 __logf_data.ln2
|
||||
#define N (1 << LOGF_TABLE_BITS)
|
||||
#define OFF 0x3f330000
|
||||
|
||||
float logf(float x)
|
||||
{
|
||||
double_t z, r, r2, y, y0, invc, logc;
|
||||
uint32_t ix, iz, tmp;
|
||||
int k, i;
|
||||
|
||||
ix = asuint(x);
|
||||
/* Fix sign of zero with downward rounding when x==1. */
|
||||
if (WANT_ROUNDING && predict_false(ix == 0x3f800000))
|
||||
return 0;
|
||||
if (predict_false(ix - 0x00800000 >= 0x7f800000 - 0x00800000)) {
|
||||
/* x < 0x1p-126 or inf or nan. */
|
||||
if (ix * 2 == 0)
|
||||
return __math_divzerof(1);
|
||||
if (ix == 0x7f800000) /* log(inf) == inf. */
|
||||
return x;
|
||||
if ((ix & 0x80000000) || ix * 2 >= 0xff000000)
|
||||
return __math_invalidf(x);
|
||||
/* x is subnormal, normalize it. */
|
||||
ix = asuint(x * 0x1p23f);
|
||||
ix -= 23 << 23;
|
||||
}
|
||||
|
||||
/* x = 2^k z; where z is in range [OFF,2*OFF] and exact.
|
||||
The range is split into N subintervals.
|
||||
The ith subinterval contains z and c is near its center. */
|
||||
tmp = ix - OFF;
|
||||
i = (tmp >> (23 - LOGF_TABLE_BITS)) % N;
|
||||
k = (int32_t)tmp >> 23; /* arithmetic shift */
|
||||
iz = ix - (tmp & 0xff800000);
|
||||
invc = T[i].invc;
|
||||
logc = T[i].logc;
|
||||
z = (double_t)asfloat(iz);
|
||||
|
||||
/* log(x) = log1p(z/c-1) + log(c) + k*Ln2 */
|
||||
r = z * invc - 1;
|
||||
y0 = logc + (double_t)k * Ln2;
|
||||
|
||||
/* Pipelined polynomial evaluation to approximate log1p(r). */
|
||||
r2 = r * r;
|
||||
y = A[1] * r + A[2];
|
||||
y = A[0] * r2 + y;
|
||||
y = y * r2 + (y0 + r);
|
||||
return eval_as_float(y);
|
||||
}
|
||||
|
|
@ -0,0 +1,33 @@
|
|||
/*
|
||||
* Data definition for logf.
|
||||
*
|
||||
* Copyright (c) 2017-2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
|
||||
#include "logf_data.h"
|
||||
|
||||
const struct logf_data __logf_data = {
|
||||
.tab = {
|
||||
{ 0x1.661ec79f8f3bep+0, -0x1.57bf7808caadep-2 },
|
||||
{ 0x1.571ed4aaf883dp+0, -0x1.2bef0a7c06ddbp-2 },
|
||||
{ 0x1.49539f0f010bp+0, -0x1.01eae7f513a67p-2 },
|
||||
{ 0x1.3c995b0b80385p+0, -0x1.b31d8a68224e9p-3 },
|
||||
{ 0x1.30d190c8864a5p+0, -0x1.6574f0ac07758p-3 },
|
||||
{ 0x1.25e227b0b8eap+0, -0x1.1aa2bc79c81p-3 },
|
||||
{ 0x1.1bb4a4a1a343fp+0, -0x1.a4e76ce8c0e5ep-4 },
|
||||
{ 0x1.12358f08ae5bap+0, -0x1.1973c5a611cccp-4 },
|
||||
{ 0x1.0953f419900a7p+0, -0x1.252f438e10c1ep-5 },
|
||||
{ 0x1p+0, 0x0p+0 },
|
||||
{ 0x1.e608cfd9a47acp-1, 0x1.aa5aa5df25984p-5 },
|
||||
{ 0x1.ca4b31f026aap-1, 0x1.c5e53aa362eb4p-4 },
|
||||
{ 0x1.b2036576afce6p-1, 0x1.526e57720db08p-3 },
|
||||
{ 0x1.9c2d163a1aa2dp-1, 0x1.bc2860d22477p-3 },
|
||||
{ 0x1.886e6037841edp-1, 0x1.1058bc8a07ee1p-2 },
|
||||
{ 0x1.767dcf5534862p-1, 0x1.4043057b6ee09p-2 },
|
||||
},
|
||||
.ln2 = 0x1.62e42fefa39efp-1,
|
||||
.poly = {
|
||||
-0x1.00ea348b88334p-2, 0x1.5575b0be00b6ap-2, -0x1.ffffef20a4123p-2,
|
||||
}
|
||||
};
|
||||
|
|
@ -0,0 +1,20 @@
|
|||
/*
|
||||
* Copyright (c) 2017-2018, Arm Limited.
|
||||
* SPDX-License-Identifier: MIT
|
||||
*/
|
||||
#ifndef _LOGF_DATA_H
|
||||
#define _LOGF_DATA_H
|
||||
|
||||
#include <features.h>
|
||||
|
||||
#define LOGF_TABLE_BITS 4
|
||||
#define LOGF_POLY_ORDER 4
|
||||
extern const struct logf_data {
|
||||
struct {
|
||||
double invc, logc;
|
||||
} tab[1 << LOGF_TABLE_BITS];
|
||||
double ln2;
|
||||
double poly[LOGF_POLY_ORDER - 1]; /* First order coefficient is 1. */
|
||||
} __logf_data;
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,72 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_tan.c */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
/* tan(x)
|
||||
* Return tangent function of x.
|
||||
*
|
||||
* kernel function:
|
||||
* __tan ... tangent function on [-pi/4,pi/4]
|
||||
* __rem_pio2 ... argument reduction routine
|
||||
*
|
||||
* Method.
|
||||
* Let S,C and T denote the sin, cos and tan respectively on
|
||||
* [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
|
||||
* in [-pi/4 , +pi/4], and let n = k mod 4.
|
||||
* We have
|
||||
*
|
||||
* n sin(x) cos(x) tan(x)
|
||||
* ----------------------------------------------------------
|
||||
* 0 S C T
|
||||
* 1 C -S -1/T
|
||||
* 2 -S -C T
|
||||
* 3 -C S -1/T
|
||||
* ----------------------------------------------------------
|
||||
*
|
||||
* Special cases:
|
||||
* Let trig be any of sin, cos, or tan.
|
||||
* trig(+-INF) is NaN, with signals;
|
||||
* trig(NaN) is that NaN;
|
||||
*
|
||||
* Accuracy:
|
||||
* TRIG(x) returns trig(x) nearly rounded
|
||||
*/
|
||||
|
||||
#include "libm.h"
|
||||
|
||||
double __tan(double x, double y, int odd);
|
||||
|
||||
double tan(double x)
|
||||
{
|
||||
double y[2];
|
||||
uint32_t ix;
|
||||
unsigned n;
|
||||
|
||||
GET_HIGH_WORD(ix, x);
|
||||
ix &= 0x7fffffff;
|
||||
|
||||
/* |x| ~< pi/4 */
|
||||
if (ix <= 0x3fe921fb) {
|
||||
if (ix < 0x3e400000) { /* |x| < 2**-27 */
|
||||
/* raise inexact if x!=0 and underflow if subnormal */
|
||||
FORCE_EVAL(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f);
|
||||
return x;
|
||||
}
|
||||
return __tan(x, 0.0, 0);
|
||||
}
|
||||
|
||||
/* tan(Inf or NaN) is NaN */
|
||||
if (ix >= 0x7ff00000)
|
||||
return x - x;
|
||||
|
||||
/* argument reduction */
|
||||
n = __rem_pio2(x, y);
|
||||
return __tan(y[0], y[1], n&1);
|
||||
}
|
||||
|
|
@ -0,0 +1,66 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_tanf.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.
|
||||
*
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "libm.h"
|
||||
|
||||
/* Small multiples of pi/2 rounded to double precision. */
|
||||
static const double
|
||||
t1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
|
||||
t2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
|
||||
t3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
|
||||
t4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */
|
||||
|
||||
float __tandf(double x, int odd);
|
||||
|
||||
float tanf(float x)
|
||||
{
|
||||
double y;
|
||||
uint32_t ix;
|
||||
unsigned n, sign;
|
||||
|
||||
GET_FLOAT_WORD(ix, x);
|
||||
sign = ix >> 31;
|
||||
ix &= 0x7fffffff;
|
||||
|
||||
if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */
|
||||
if (ix < 0x39800000) { /* |x| < 2**-12 */
|
||||
/* raise inexact if x!=0 and underflow if subnormal */
|
||||
FORCE_EVAL(ix < 0x00800000 ? x/0x1p120f : x+0x1p120f);
|
||||
return x;
|
||||
}
|
||||
return __tandf(x, 0);
|
||||
}
|
||||
if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */
|
||||
if (ix <= 0x4016cbe3) /* |x| ~<= 3pi/4 */
|
||||
return __tandf((sign ? x+t1pio2 : x-t1pio2), 1);
|
||||
else
|
||||
return __tandf((sign ? x+t2pio2 : x-t2pio2), 0);
|
||||
}
|
||||
if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */
|
||||
if (ix <= 0x40afeddf) /* |x| ~<= 7*pi/4 */
|
||||
return __tandf((sign ? x+t3pio2 : x-t3pio2), 1);
|
||||
else
|
||||
return __tandf((sign ? x+t4pio2 : x-t4pio2), 0);
|
||||
}
|
||||
|
||||
/* tan(Inf or NaN) is NaN */
|
||||
if (ix >= 0x7f800000)
|
||||
return x - x;
|
||||
|
||||
/* argument reduction */
|
||||
n = __rem_pio2f(x, &y);
|
||||
return __tandf(y, n&1);
|
||||
}
|
||||
Loading…
Reference in New Issue