Skip to content

Commit 0451165

Browse files
dpgeorgedpgeorge
committed
lib: Add libm_dbl, a double-precision math library, from musl-1.1.16.
1 parent 409fc8f commit 0451165

43 files changed

Lines changed: 3743 additions & 0 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎lib/libm_dbl/README‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
This directory contains source code for the standard double-precision math
2+
functions.
3+
4+
The files lgamma.c, log10.c and tanh.c are too small to have a meaningful
5+
copyright or license.
6+
7+
The rest of the files in this directory are copied from the musl library,
8+
v1.1.16, and, unless otherwise stated in the individual file, have the
9+
following copyright and MIT license:
10+
11+
----------------------------------------------------------------------
12+
Copyright © 2005-2014 Rich Felker, et al.
13+
14+
Permission is hereby granted, free of charge, to any person obtaining
15+
a copy of this software and associated documentation files (the
16+
"Software"), to deal in the Software without restriction, including
17+
without limitation the rights to use, copy, modify, merge, publish,
18+
distribute, sublicense, and/or sell copies of the Software, and to
19+
permit persons to whom the Software is furnished to do so, subject to
20+
the following conditions:
21+
22+
The above copyright notice and this permission notice shall be
23+
included in all copies or substantial portions of the Software.
24+
25+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
26+
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
27+
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
28+
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
29+
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
30+
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
31+
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
32+
----------------------------------------------------------------------

‎lib/libm_dbl/__cos.c‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
/* origin: FreeBSD /usr/src/lib/msun/src/k_cos.c */
2+
/*
3+
* ====================================================
4+
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5+
*
6+
* Developed at SunSoft, a Sun Microsystems, Inc. business.
7+
* Permission to use, copy, modify, and distribute this
8+
* software is freely granted, provided that this notice
9+
* is preserved.
10+
* ====================================================
11+
*/
12+
/*
13+
* __cos( x, y )
14+
* kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
15+
* Input x is assumed to be bounded by ~pi/4 in magnitude.
16+
* Input y is the tail of x.
17+
*
18+
* Algorithm
19+
* 1. Since cos(-x) = cos(x), we need only to consider positive x.
20+
* 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
21+
* 3. cos(x) is approximated by a polynomial of degree 14 on
22+
* [0,pi/4]
23+
* 4 14
24+
* cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
25+
* where the remez error is
26+
*
27+
* | 2 4 6 8 10 12 14 | -58
28+
* |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
29+
* | |
30+
*
31+
* 4 6 8 10 12 14
32+
* 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
33+
* cos(x) ~ 1 - x*x/2 + r
34+
* since cos(x+y) ~ cos(x) - sin(x)*y
35+
* ~ cos(x) - x*y,
36+
* a correction term is necessary in cos(x) and hence
37+
* cos(x+y) = 1 - (x*x/2 - (r - x*y))
38+
* For better accuracy, rearrange to
39+
* cos(x+y) ~ w + (tmp + (r-x*y))
40+
* where w = 1 - x*x/2 and tmp is a tiny correction term
41+
* (1 - x*x/2 == w + tmp exactly in infinite precision).
42+
* The exactness of w + tmp in infinite precision depends on w
43+
* and tmp having the same precision as x. If they have extra
44+
* precision due to compiler bugs, then the extra precision is
45+
* only good provided it is retained in all terms of the final
46+
* expression for cos(). Retention happens in all cases tested
47+
* under FreeBSD, so don't pessimize things by forcibly clipping
48+
* any extra precision in w.
49+
*/
50+
51+
#include "libm.h"
52+
53+
static const double
54+
C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
55+
C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
56+
C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
57+
C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
58+
C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
59+
C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
60+
61+
double __cos(double x, double y)
62+
{
63+
double_t hz,z,r,w;
64+
65+
z = x*x;
66+
w = z*z;
67+
r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
68+
hz = 0.5*z;
69+
w = 1.0-hz;
70+
return w + (((1.0-w)-hz) + (z*r-x*y));
71+
}

‎lib/libm_dbl/__expo2.c‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
#include "libm.h"
2+
3+
/* k is such that k*ln2 has minimal relative error and x - kln2 > log(DBL_MIN) */
4+
static const int k = 2043;
5+
static const double kln2 = 0x1.62066151add8bp+10;
6+
7+
/* exp(x)/2 for x >= log(DBL_MAX), slightly better than 0.5*exp(x/2)*exp(x/2) */
8+
double __expo2(double x)
9+
{
10+
double scale;
11+
12+
/* note that k is odd and scale*scale overflows */
13+
INSERT_WORDS(scale, (uint32_t)(0x3ff + k/2) << 20, 0);
14+
/* exp(x - k ln2) * 2**(k-1) */
15+
return exp(x - kln2) * scale * scale;
16+
}

‎lib/libm_dbl/__fpclassify.c‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
#include <math.h>
2+
#include <stdint.h>
3+
4+
int __fpclassifyd(double x)
5+
{
6+
union {double f; uint64_t i;} u = {x};
7+
int e = u.i>>52 & 0x7ff;
8+
if (!e) return u.i<<1 ? FP_SUBNORMAL : FP_ZERO;
9+
if (e==0x7ff) return u.i<<12 ? FP_NAN : FP_INFINITE;
10+
return FP_NORMAL;
11+
}

‎lib/libm_dbl/__rem_pio2.c‎

Lines changed: 177 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,177 @@
1+
/* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c */
2+
/*
3+
* ====================================================
4+
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5+
*
6+
* Developed at SunSoft, a Sun Microsystems, Inc. business.
7+
* Permission to use, copy, modify, and distribute this
8+
* software is freely granted, provided that this notice
9+
* is preserved.
10+
* ====================================================
11+
*
12+
* Optimized by Bruce D. Evans.
13+
*/
14+
/* __rem_pio2(x,y)
15+
*
16+
* return the remainder of x rem pi/2 in y[0]+y[1]
17+
* use __rem_pio2_large() for large x
18+
*/
19+
20+
#include "libm.h"
21+
22+
#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
23+
#define EPS DBL_EPSILON
24+
#elif FLT_EVAL_METHOD==2
25+
#define EPS LDBL_EPSILON
26+
#endif
27+
28+
/*
29+
* invpio2: 53 bits of 2/pi
30+
* pio2_1: first 33 bit of pi/2
31+
* pio2_1t: pi/2 - pio2_1
32+
* pio2_2: second 33 bit of pi/2
33+
* pio2_2t: pi/2 - (pio2_1+pio2_2)
34+
* pio2_3: third 33 bit of pi/2
35+
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
36+
*/
37+
static const double
38+
toint = 1.5/EPS,
39+
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
40+
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
41+
pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
42+
pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
43+
pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
44+
pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
45+
pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
46+
47+
/* caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
48+
int __rem_pio2(double x, double *y)
49+
{
50+
union {double f; uint64_t i;} u = {x};
51+
double_t z,w,t,r,fn;
52+
double tx[3],ty[2];
53+
uint32_t ix;
54+
int sign, n, ex, ey, i;
55+
56+
sign = u.i>>63;
57+
ix = u.i>>32 & 0x7fffffff;
58+
if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */
59+
if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */
60+
goto medium; /* cancellation -- use medium case */
61+
if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */
62+
if (!sign) {
63+
z = x - pio2_1; /* one round good to 85 bits */
64+
y[0] = z - pio2_1t;
65+
y[1] = (z-y[0]) - pio2_1t;
66+
return 1;
67+
} else {
68+
z = x + pio2_1;
69+
y[0] = z + pio2_1t;
70+
y[1] = (z-y[0]) + pio2_1t;
71+
return -1;
72+
}
73+
} else {
74+
if (!sign) {
75+
z = x - 2*pio2_1;
76+
y[0] = z - 2*pio2_1t;
77+
y[1] = (z-y[0]) - 2*pio2_1t;
78+
return 2;
79+
} else {
80+
z = x + 2*pio2_1;
81+
y[0] = z + 2*pio2_1t;
82+
y[1] = (z-y[0]) + 2*pio2_1t;
83+
return -2;
84+
}
85+
}
86+
}
87+
if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */
88+
if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */
89+
if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */
90+
goto medium;
91+
if (!sign) {
92+
z = x - 3*pio2_1;
93+
y[0] = z - 3*pio2_1t;
94+
y[1] = (z-y[0]) - 3*pio2_1t;
95+
return 3;
96+
} else {
97+
z = x + 3*pio2_1;
98+
y[0] = z + 3*pio2_1t;
99+
y[1] = (z-y[0]) + 3*pio2_1t;
100+
return -3;
101+
}
102+
} else {
103+
if (ix == 0x401921fb) /* |x| ~= 4pi/2 */
104+
goto medium;
105+
if (!sign) {
106+
z = x - 4*pio2_1;
107+
y[0] = z - 4*pio2_1t;
108+
y[1] = (z-y[0]) - 4*pio2_1t;
109+
return 4;
110+
} else {
111+
z = x + 4*pio2_1;
112+
y[0] = z + 4*pio2_1t;
113+
y[1] = (z-y[0]) + 4*pio2_1t;
114+
return -4;
115+
}
116+
}
117+
}
118+
if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */
119+
medium:
120+
/* rint(x/(pi/2)), Assume round-to-nearest. */
121+
fn = (double_t)x*invpio2 + toint - toint;
122+
n = (int32_t)fn;
123+
r = x - fn*pio2_1;
124+
w = fn*pio2_1t; /* 1st round, good to 85 bits */
125+
y[0] = r - w;
126+
u.f = y[0];
127+
ey = u.i>>52 & 0x7ff;
128+
ex = ix>>20;
129+
if (ex - ey > 16) { /* 2nd round, good to 118 bits */
130+
t = r;
131+
w = fn*pio2_2;
132+
r = t - w;
133+
w = fn*pio2_2t - ((t-r)-w);
134+
y[0] = r - w;
135+
u.f = y[0];
136+
ey = u.i>>52 & 0x7ff;
137+
if (ex - ey > 49) { /* 3rd round, good to 151 bits, covers all cases */
138+
t = r;
139+
w = fn*pio2_3;
140+
r = t - w;
141+
w = fn*pio2_3t - ((t-r)-w);
142+
y[0] = r - w;
143+
}
144+
}
145+
y[1] = (r - y[0]) - w;
146+
return n;
147+
}
148+
/*
149+
* all other (large) arguments
150+
*/
151+
if (ix >= 0x7ff00000) { /* x is inf or NaN */
152+
y[0] = y[1] = x - x;
153+
return 0;
154+
}
155+
/* set z = scalbn(|x|,-ilogb(x)+23) */
156+
u.f = x;
157+
u.i &= (uint64_t)-1>>12;
158+
u.i |= (uint64_t)(0x3ff + 23)<<52;
159+
z = u.f;
160+
for (i=0; i < 2; i++) {
161+
tx[i] = (double)(int32_t)z;
162+
z = (z-tx[i])*0x1p24;
163+
}
164+
tx[i] = z;
165+
/* skip zero terms, first term is non-zero */
166+
while (tx[i] == 0.0)
167+
i--;
168+
n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1);
169+
if (sign) {
170+
y[0] = -ty[0];
171+
y[1] = -ty[1];
172+
return -n;
173+
}
174+
y[0] = ty[0];
175+
y[1] = ty[1];
176+
return n;
177+
}

0 commit comments

Comments
 (0)