jcs.org
Serenity Operating System
1/*
2 * Copyright (c) 2018-2020, Andreas Kling <kling@serenityos.org>
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice, this
9 * list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
16 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
19 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
21 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
22 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
23 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
24 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 */
26
27#include <LibC/assert.h>
28#include <LibM/math.h>
29#include <stdint.h>
30#include <stdlib.h>
31
32template<size_t>
33constexpr double e_to_power();
34template<>
35constexpr double e_to_power<0>() { return 1; }
36template<size_t exponent>
37constexpr double e_to_power() { return M_E * e_to_power<exponent - 1>(); }
38
39template<size_t>
40constexpr size_t factorial();
41template<>
42constexpr size_t factorial<0>() { return 1; }
43template<size_t value>
44constexpr size_t factorial() { return value * factorial<value - 1>(); }
45
46template<size_t>
47constexpr size_t product_even();
48template<>
49constexpr size_t product_even<2>() { return 2; }
50template<size_t value>
51constexpr size_t product_even() { return value * product_even<value - 2>(); }
52
53template<size_t>
54constexpr size_t product_odd();
55template<>
56constexpr size_t product_odd<1>() { return 1; }
57template<size_t value>
58constexpr size_t product_odd() { return value * product_odd<value - 2>(); }
59
60extern "C" {
61
62double trunc(double x)
63{
64 return (int64_t)x;
65}
66
67double cos(double angle)
68{
69 return sin(angle + M_PI_2);
70}
71
72// This can also be done with a taylor expansion, but for
73// now this works pretty well (and doesn't mess anything up
74// in quake in particular, which is very Floating-Point precision
75// heavy)
76double sin(double angle)
77{
78 double ret = 0.0;
79 __asm__(
80 "fsin"
81 : "=t"(ret)
82 : "0"(angle));
83
84 return ret;
85}
86
87double pow(double x, double y)
88{
89 //FIXME: Extremely unlikely to be standards compliant.
90 return exp(y * log(x));
91}
92
93double ldexp(double x, int exp)
94{
95 // FIXME: Please fix me. I am naive.
96 double val = pow(2, exp);
97 return x * val;
98}
99
100double tanh(double x)
101{
102 if (x > 0) {
103 double exponentiated = exp(2 * x);
104 return (exponentiated - 1) / (exponentiated + 1);
105 }
106 double plusX = exp(x);
107 double minusX = 1 / plusX;
108 return (plusX - minusX) / (plusX + minusX);
109}
110
111double ampsin(double angle)
112{
113 double looped_angle = fmod(M_PI + angle, M_TAU) - M_PI;
114 double looped_angle_squared = looped_angle * looped_angle;
115
116 double quadratic_term;
117 if (looped_angle > 0) {
118 quadratic_term = -looped_angle_squared;
119 } else {
120 quadratic_term = looped_angle_squared;
121 }
122
123 double linear_term = M_PI * looped_angle;
124
125 return quadratic_term + linear_term;
126}
127
128double tan(double angle)
129{
130 return ampsin(angle) / ampsin(M_PI_2 + angle);
131}
132
133double sqrt(double x)
134{
135 double res;
136 __asm__("fsqrt"
137 : "=t"(res)
138 : "0"(x));
139 return res;
140}
141
142double sinh(double x)
143{
144 double exponentiated = exp(x);
145 if (x > 0)
146 return (exponentiated * exponentiated - 1) / 2 / exponentiated;
147 return (exponentiated - 1 / exponentiated) / 2;
148}
149
150double log10(double x)
151{
152 return log(x) / M_LN10;
153}
154
155double log(double x)
156{
157 if (x < 0)
158 return __builtin_nan("");
159 if (x == 0)
160 return -__builtin_huge_val();
161 double y = 1 + 2 * (x - 1) / (x + 1);
162 double exponentiated = exp(y);
163 y = y + 2 * (x - exponentiated) / (x + exponentiated);
164 exponentiated = exp(y);
165 y = y + 2 * (x - exponentiated) / (x + exponentiated);
166 exponentiated = exp(y);
167 return y + 2 * (x - exponentiated) / (x + exponentiated);
168}
169
170double fmod(double index, double period)
171{
172 return index - trunc(index / period) * period;
173}
174
175double exp(double exponent)
176{
177 double result = 1;
178 if (exponent >= 1) {
179 size_t integer_part = (size_t)exponent;
180 if (integer_part & 1)
181 result *= e_to_power<1>();
182 if (integer_part & 2)
183 result *= e_to_power<2>();
184 if (integer_part > 3) {
185 if (integer_part & 4)
186 result *= e_to_power<4>();
187 if (integer_part & 8)
188 result *= e_to_power<8>();
189 if (integer_part & 16)
190 result *= e_to_power<16>();
191 if (integer_part & 32)
192 result *= e_to_power<32>();
193 if (integer_part >= 64)
194 return __builtin_huge_val();
195 }
196 exponent -= integer_part;
197 } else if (exponent < 0)
198 return 1 / exp(-exponent);
199 double taylor_series_result = 1 + exponent;
200 double taylor_series_numerator = exponent * exponent;
201 taylor_series_result += taylor_series_numerator / factorial<2>();
202 taylor_series_numerator *= exponent;
203 taylor_series_result += taylor_series_numerator / factorial<3>();
204 taylor_series_numerator *= exponent;
205 taylor_series_result += taylor_series_numerator / factorial<4>();
206 taylor_series_numerator *= exponent;
207 taylor_series_result += taylor_series_numerator / factorial<5>();
208 return result * taylor_series_result;
209}
210
211double cosh(double x)
212{
213 double exponentiated = exp(-x);
214 if (x < 0)
215 return (1 + exponentiated * exponentiated) / 2 / exponentiated;
216 return (1 / exponentiated + exponentiated) / 2;
217}
218
219double atan2(double y, double x)
220{
221 if (x > 0)
222 return atan(y / x);
223 if (x == 0) {
224 if (y > 0)
225 return M_PI_2;
226 if (y < 0)
227 return -M_PI_2;
228 return 0;
229 }
230 if (y >= 0)
231 return atan(y / x) + M_PI;
232 return atan(y / x) - M_PI;
233}
234
235double atan(double x)
236{
237 if (x < 0)
238 return -atan(-x);
239 if (x > 1)
240 return M_PI_2 - atan(1 / x);
241 double squared = x * x;
242 return x / (1 + 1 * 1 * squared / (3 + 2 * 2 * squared / (5 + 3 * 3 * squared / (7 + 4 * 4 * squared / (9 + 5 * 5 * squared / (11 + 6 * 6 * squared / (13 + 7 * 7 * squared)))))));
243}
244
245double asin(double x)
246{
247 if (x > 1 || x < -1)
248 return __builtin_nan("");
249 if (x > 0.5 || x < -0.5)
250 return 2 * atan(x / (1 + sqrt(1 - x * x)));
251 double squared = x * x;
252 double value = x;
253 double i = x * squared;
254 value += i * product_odd<1>() / product_even<2>() / 3;
255 i *= squared;
256 value += i * product_odd<3>() / product_even<4>() / 5;
257 i *= squared;
258 value += i * product_odd<5>() / product_even<6>() / 7;
259 i *= squared;
260 value += i * product_odd<7>() / product_even<8>() / 9;
261 i *= squared;
262 value += i * product_odd<9>() / product_even<10>() / 11;
263 i *= squared;
264 value += i * product_odd<11>() / product_even<12>() / 13;
265 return value;
266}
267
268double acos(double x)
269{
270 return M_PI_2 - asin(x);
271}
272
273double fabs(double value)
274{
275 return value < 0 ? -value : value;
276}
277
278double log2(double x)
279{
280 return log(x) / M_LN2;
281}
282
283float log2f(float x)
284{
285 return log2(x);
286}
287
288long double log2l(long double x)
289{
290 return log2(x);
291}
292
293double frexp(double, int*)
294{
295 ASSERT_NOT_REACHED();
296 return 0;
297}
298
299float frexpf(float, int*)
300{
301 ASSERT_NOT_REACHED();
302 return 0;
303}
304
305long double frexpl(long double, int*)
306{
307 ASSERT_NOT_REACHED();
308 return 0;
309}
310
311float roundf(float value)
312{
313 // FIXME: Please fix me. I am naive.
314 if (value >= 0.0f)
315 return (float)(int)(value + 0.5f);
316 return (float)(int)(value - 0.5f);
317}
318
319double floor(double value)
320{
321 return (int)value;
322}
323
324double rint(double value)
325{
326 return (int)roundf(value);
327}
328
329float ceilf(float value)
330{
331 // FIXME: Please fix me. I am naive.
332 int as_int = (int)value;
333 if (value == (float)as_int) {
334 return (float)as_int;
335 }
336 return as_int + 1;
337}
338
339double ceil(double value)
340{
341 // FIXME: Please fix me. I am naive.
342 int as_int = (int)value;
343 if (value == (double)as_int) {
344 return (double)as_int;
345 }
346 return as_int + 1;
347}
348
349double modf(double x, double* intpart)
350{
351 *intpart = (double)((int)(x));
352 return x - (int)x;
353}
354}