Skip to content

Commit 8a23643

Browse files
dosisoddosisod
authored
Add more math functions (#128)
1 parent ab64ef7 commit 8a23643

File tree

13 files changed

+1249
-0
lines changed

13 files changed

+1249
-0
lines changed

‎libc/tinymath/erff.c

Lines changed: 220 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,220 @@
1+
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
2+
│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│
3+
╚──────────────────────────────────────────────────────────────────────────────╝
4+
│ │
5+
│ Musl Libc │
6+
│ Copyright © 2005-2020 Rich Felker, et al. │
7+
│ │
8+
│ Permission is hereby granted, free of charge, to any person obtaining │
9+
│ a copy of this software and associated documentation files (the │
10+
│ "Software"), to deal in the Software without restriction, including │
11+
│ without limitation the rights to use, copy, modify, merge, publish, │
12+
│ distribute, sublicense, and/or sell copies of the Software, and to │
13+
│ permit persons to whom the Software is furnished to do so, subject to │
14+
│ the following conditions: │
15+
│ │
16+
│ The above copyright notice and this permission notice shall be │
17+
│ included in all copies or substantial portions of the Software. │
18+
│ │
19+
│ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, │
20+
│ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF │
21+
│ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. │
22+
│ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY │
23+
│ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, │
24+
│ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE │
25+
│ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │
26+
│ │
27+
╚─────────────────────────────────────────────────────────────────────────────*/
28+
29+
#include "libc/math.h"
30+
31+
asm(".ident\t\"\\n\\n\
32+
fdlibm (fdlibm license)\\n\
33+
Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
34+
asm(".ident\t\"\\n\\n\
35+
Musl libc (MIT License)\\n\
36+
Copyright 2005-2020 Rich Felker, et. al.\"");
37+
asm(".include \"libc/disclaimer.inc\"");
38+
39+
/* clang-format off */
40+
/* origin: FreeBSD /usr/src/lib/msun/src/s_erff.c */
41+
/*
42+
* ====================================================
43+
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
44+
*
45+
* Developed at SunPro, a Sun Microsystems, Inc. business.
46+
* Permission to use, copy, modify, and distribute this
47+
* software is freely granted, provided that this notice
48+
* is preserved.
49+
* ====================================================
50+
*/
51+
52+
#define asuint(f) ((union{float _f; uint32_t _i;}){f})._i
53+
#define asfloat(i) ((union{uint32_t _i; float _f;}){i})._f
54+
55+
static const float
56+
erx = 8.4506291151e-01, /* 0x3f58560b */
57+
/*
58+
* Coefficients for approximation to erf on [0,0.84375]
59+
*/
60+
efx8 = 1.0270333290e+00, /* 0x3f8375d4 */
61+
pp0 = 1.2837916613e-01, /* 0x3e0375d4 */
62+
pp1 = -3.2504209876e-01, /* 0xbea66beb */
63+
pp2 = -2.8481749818e-02, /* 0xbce9528f */
64+
pp3 = -5.7702702470e-03, /* 0xbbbd1489 */
65+
pp4 = -2.3763017452e-05, /* 0xb7c756b1 */
66+
qq1 = 3.9791721106e-01, /* 0x3ecbbbce */
67+
qq2 = 6.5022252500e-02, /* 0x3d852a63 */
68+
qq3 = 5.0813062117e-03, /* 0x3ba68116 */
69+
qq4 = 1.3249473704e-04, /* 0x390aee49 */
70+
qq5 = -3.9602282413e-06, /* 0xb684e21a */
71+
/*
72+
* Coefficients for approximation to erf in [0.84375,1.25]
73+
*/
74+
pa0 = -2.3621185683e-03, /* 0xbb1acdc6 */
75+
pa1 = 4.1485610604e-01, /* 0x3ed46805 */
76+
pa2 = -3.7220788002e-01, /* 0xbebe9208 */
77+
pa3 = 3.1834661961e-01, /* 0x3ea2fe54 */
78+
pa4 = -1.1089469492e-01, /* 0xbde31cc2 */
79+
pa5 = 3.5478305072e-02, /* 0x3d1151b3 */
80+
pa6 = -2.1663755178e-03, /* 0xbb0df9c0 */
81+
qa1 = 1.0642088205e-01, /* 0x3dd9f331 */
82+
qa2 = 5.4039794207e-01, /* 0x3f0a5785 */
83+
qa3 = 7.1828655899e-02, /* 0x3d931ae7 */
84+
qa4 = 1.2617121637e-01, /* 0x3e013307 */
85+
qa5 = 1.3637083583e-02, /* 0x3c5f6e13 */
86+
qa6 = 1.1984500103e-02, /* 0x3c445aa3 */
87+
/*
88+
* Coefficients for approximation to erfc in [1.25,1/0.35]
89+
*/
90+
ra0 = -9.8649440333e-03, /* 0xbc21a093 */
91+
ra1 = -6.9385856390e-01, /* 0xbf31a0b7 */
92+
ra2 = -1.0558626175e+01, /* 0xc128f022 */
93+
ra3 = -6.2375331879e+01, /* 0xc2798057 */
94+
ra4 = -1.6239666748e+02, /* 0xc322658c */
95+
ra5 = -1.8460508728e+02, /* 0xc3389ae7 */
96+
ra6 = -8.1287437439e+01, /* 0xc2a2932b */
97+
ra7 = -9.8143291473e+00, /* 0xc11d077e */
98+
sa1 = 1.9651271820e+01, /* 0x419d35ce */
99+
sa2 = 1.3765776062e+02, /* 0x4309a863 */
100+
sa3 = 4.3456588745e+02, /* 0x43d9486f */
101+
sa4 = 6.4538726807e+02, /* 0x442158c9 */
102+
sa5 = 4.2900814819e+02, /* 0x43d6810b */
103+
sa6 = 1.0863500214e+02, /* 0x42d9451f */
104+
sa7 = 6.5702495575e+00, /* 0x40d23f7c */
105+
sa8 = -6.0424413532e-02, /* 0xbd777f97 */
106+
/*
107+
* Coefficients for approximation to erfc in [1/.35,28]
108+
*/
109+
rb0 = -9.8649431020e-03, /* 0xbc21a092 */
110+
rb1 = -7.9928326607e-01, /* 0xbf4c9dd4 */
111+
rb2 = -1.7757955551e+01, /* 0xc18e104b */
112+
rb3 = -1.6063638306e+02, /* 0xc320a2ea */
113+
rb4 = -6.3756646729e+02, /* 0xc41f6441 */
114+
rb5 = -1.0250950928e+03, /* 0xc480230b */
115+
rb6 = -4.8351919556e+02, /* 0xc3f1c275 */
116+
sb1 = 3.0338060379e+01, /* 0x41f2b459 */
117+
sb2 = 3.2579251099e+02, /* 0x43a2e571 */
118+
sb3 = 1.5367296143e+03, /* 0x44c01759 */
119+
sb4 = 3.1998581543e+03, /* 0x4547fdbb */
120+
sb5 = 2.5530502930e+03, /* 0x451f90ce */
121+
sb6 = 4.7452853394e+02, /* 0x43ed43a7 */
122+
sb7 = -2.2440952301e+01; /* 0xc1b38712 */
123+
124+
static float erfc1(float x)
125+
{
126+
float_t s,P,Q;
127+
128+
s = fabsf(x) - 1;
129+
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
130+
Q = 1+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
131+
return 1 - erx - P/Q;
132+
}
133+
134+
static float erfc2(uint32_t ix, float x)
135+
{
136+
float_t s,R,S;
137+
float z;
138+
139+
if (ix < 0x3fa00000) /* |x| < 1.25 */
140+
return erfc1(x);
141+
142+
x = fabsf(x);
143+
s = 1/(x*x);
144+
if (ix < 0x4036db6d) { /* |x| < 1/0.35 */
145+
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
146+
ra5+s*(ra6+s*ra7))))));
147+
S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
148+
sa5+s*(sa6+s*(sa7+s*sa8)))))));
149+
} else { /* |x| >= 1/0.35 */
150+
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
151+
rb5+s*rb6)))));
152+
S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
153+
sb5+s*(sb6+s*sb7))))));
154+
}
155+
ix = asuint(x);
156+
z = asfloat(ix&0xffffe000);
157+
return expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S)/x;
158+
}
159+
160+
float erff(float x)
161+
{
162+
float r,s,z,y;
163+
uint32_t ix;
164+
int sign;
165+
166+
ix = asuint(x);
167+
sign = ix>>31;
168+
ix &= 0x7fffffff;
169+
if (ix >= 0x7f800000) {
170+
/* erf(nan)=nan, erf(+-inf)=+-1 */
171+
return 1-2*sign + 1/x;
172+
}
173+
if (ix < 0x3f580000) { /* |x| < 0.84375 */
174+
if (ix < 0x31800000) { /* |x| < 2**-28 */
175+
/*avoid underflow */
176+
return 0.125f*(8*x + efx8*x);
177+
}
178+
z = x*x;
179+
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
180+
s = 1+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
181+
y = r/s;
182+
return x + x*y;
183+
}
184+
if (ix < 0x40c00000) /* |x| < 6 */
185+
y = 1 - erfc2(ix,x);
186+
else
187+
y = 1 - 0x1p-120f;
188+
return sign ? -y : y;
189+
}
190+
191+
float erfcf(float x)
192+
{
193+
float r,s,z,y;
194+
uint32_t ix;
195+
int sign;
196+
197+
ix = asuint(x);
198+
sign = ix>>31;
199+
ix &= 0x7fffffff;
200+
if (ix >= 0x7f800000) {
201+
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
202+
return 2*sign + 1/x;
203+
}
204+
205+
if (ix < 0x3f580000) { /* |x| < 0.84375 */
206+
if (ix < 0x23800000) /* |x| < 2**-56 */
207+
return 1.0f - x;
208+
z = x*x;
209+
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
210+
s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
211+
y = r/s;
212+
if (sign || ix < 0x3e800000) /* x < 1/4 */
213+
return 1.0f - (x+x*y);
214+
return 0.5f - (x - 0.5f + x*y);
215+
}
216+
if (ix < 0x41e00000) { /* |x| < 28 */
217+
return sign ? 2 - erfc2(ix,x) : erfc2(ix,x);
218+
}
219+
return sign ? 2 - 0x1p-120f : 0x1p-120f*0x1p-120f;
220+
}

‎libc/tinymath/fdim.c

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
2+
│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│
3+
╚──────────────────────────────────────────────────────────────────────────────╝
4+
│ │
5+
│ Musl Libc │
6+
│ Copyright © 2005-2020 Rich Felker, et al. │
7+
│ │
8+
│ Permission is hereby granted, free of charge, to any person obtaining │
9+
│ a copy of this software and associated documentation files (the │
10+
│ "Software"), to deal in the Software without restriction, including │
11+
│ without limitation the rights to use, copy, modify, merge, publish, │
12+
│ distribute, sublicense, and/or sell copies of the Software, and to │
13+
│ permit persons to whom the Software is furnished to do so, subject to │
14+
│ the following conditions: │
15+
│ │
16+
│ The above copyright notice and this permission notice shall be │
17+
│ included in all copies or substantial portions of the Software. │
18+
│ │
19+
│ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, │
20+
│ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF │
21+
│ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. │
22+
│ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY │
23+
│ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, │
24+
│ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE │
25+
│ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │
26+
│ │
27+
╚─────────────────────────────────────────────────────────────────────────────*/
28+
29+
#include "libc/math.h"
30+
31+
asm(".ident\t\"\\n\\n\
32+
Musl libc (MIT License)\\n\
33+
Copyright 2005-2020 Rich Felker, et. al.\"");
34+
asm(".include \"libc/disclaimer.inc\"");
35+
36+
/* clang-format off */
37+
38+
double fdim(double x, double y)
39+
{
40+
if (isnan(x))
41+
return x;
42+
if (isnan(y))
43+
return y;
44+
return x > y ? x - y : 0;
45+
}

‎libc/tinymath/fdimf.c

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
2+
│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│
3+
╚──────────────────────────────────────────────────────────────────────────────╝
4+
│ │
5+
│ Musl Libc │
6+
│ Copyright © 2005-2020 Rich Felker, et al. │
7+
│ │
8+
│ Permission is hereby granted, free of charge, to any person obtaining │
9+
│ a copy of this software and associated documentation files (the │
10+
│ "Software"), to deal in the Software without restriction, including │
11+
│ without limitation the rights to use, copy, modify, merge, publish, │
12+
│ distribute, sublicense, and/or sell copies of the Software, and to │
13+
│ permit persons to whom the Software is furnished to do so, subject to │
14+
│ the following conditions: │
15+
│ │
16+
│ The above copyright notice and this permission notice shall be │
17+
│ included in all copies or substantial portions of the Software. │
18+
│ │
19+
│ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, │
20+
│ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF │
21+
│ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. │
22+
│ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY │
23+
│ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, │
24+
│ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE │
25+
│ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │
26+
│ │
27+
╚─────────────────────────────────────────────────────────────────────────────*/
28+
29+
#include "libc/math.h"
30+
31+
asm(".ident\t\"\\n\\n\
32+
Musl libc (MIT License)\\n\
33+
Copyright 2005-2020 Rich Felker, et. al.\"");
34+
asm(".include \"libc/disclaimer.inc\"");
35+
36+
/* clang-format off */
37+
38+
float fdimf(float x, float y)
39+
{
40+
if (isnan(x))
41+
return x;
42+
if (isnan(y))
43+
return y;
44+
return x > y ? x - y : 0;
45+
}

‎libc/tinymath/fdiml.c

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
2+
│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│
3+
╚──────────────────────────────────────────────────────────────────────────────╝
4+
│ │
5+
│ Musl Libc │
6+
│ Copyright © 2005-2020 Rich Felker, et al. │
7+
│ │
8+
│ Permission is hereby granted, free of charge, to any person obtaining │
9+
│ a copy of this software and associated documentation files (the │
10+
│ "Software"), to deal in the Software without restriction, including │
11+
│ without limitation the rights to use, copy, modify, merge, publish, │
12+
│ distribute, sublicense, and/or sell copies of the Software, and to │
13+
│ permit persons to whom the Software is furnished to do so, subject to │
14+
│ the following conditions: │
15+
│ │
16+
│ The above copyright notice and this permission notice shall be │
17+
│ included in all copies or substantial portions of the Software. │
18+
│ │
19+
│ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, │
20+
│ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF │
21+
│ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. │
22+
│ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY │
23+
│ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, │
24+
│ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE │
25+
│ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │
26+
│ │
27+
╚─────────────────────────────────────────────────────────────────────────────*/
28+
29+
#include "libc/math.h"
30+
31+
asm(".ident\t\"\\n\\n\
32+
Musl libc (MIT License)\\n\
33+
Copyright 2005-2020 Rich Felker, et. al.\"");
34+
asm(".include \"libc/disclaimer.inc\"");
35+
36+
/* clang-format off */
37+
38+
long double fdiml(long double x, long double y)
39+
{
40+
if (isnan(x))
41+
return x;
42+
if (isnan(y))
43+
return y;
44+
return x > y ? x - y : 0;
45+
}

0 commit comments

Comments
 (0)