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
72float cosf(float angle)
73{
74 return sinf(angle + M_PI_2);
75}
76
77// This can also be done with a taylor expansion, but for
78// now this works pretty well (and doesn't mess anything up
79// in quake in particular, which is very Floating-Point precision
80// heavy)
81double sin(double angle)
82{
83 double ret = 0.0;
84 __asm__(
85 "fsin"
86 : "=t"(ret)
87 : "0"(angle));
88
89 return ret;
90}
91
92float sinf(float angle)
93{
94 float ret = 0.0f;
95 __asm__(
96 "fsin"
97 : "=t"(ret)
98 : "0"(angle));
99 return ret;
100}
101
102double pow(double x, double y)
103{
104 // FIXME: Please fix me. I am naive.
105 if (y == 0)
106 return 1;
107 if (y == 1)
108 return x;
109 int y_as_int = (int)y;
110 if (y == (double)y_as_int) {
111 double result = x;
112 for (int i = 0; i < abs(y) - 1; ++i)
113 result *= x;
114 if (y < 0)
115 result = 1.0 / result;
116 return result;
117 }
118 return exp(y * log(x));
119}
120
121float powf(float x, float y)
122{
123 // FIXME: Please fix me. I am naive.
124 if (y == 0)
125 return 1;
126 if (y == 1)
127 return x;
128 int y_as_int = (int)y;
129 if (y == (float)y_as_int) {
130 float result = x;
131 for (int i = 0; i < abs(y) - 1; ++i)
132 result *= x;
133 if (y < 0)
134 result = 1.0 / result;
135 return result;
136 }
137 return (float)exp((double)y * log((double)x));
138}
139
140double ldexp(double x, int exp)
141{
142 // FIXME: Please fix me. I am naive.
143 double val = pow(2, exp);
144 return x * val;
145}
146
147double tanh(double x)
148{
149 if (x > 0) {
150 double exponentiated = exp(2 * x);
151 return (exponentiated - 1) / (exponentiated + 1);
152 }
153 double plusX = exp(x);
154 double minusX = 1 / plusX;
155 return (plusX - minusX) / (plusX + minusX);
156}
157
158double ampsin(double angle)
159{
160 double looped_angle = fmod(M_PI + angle, M_TAU) - M_PI;
161 double looped_angle_squared = looped_angle * looped_angle;
162
163 double quadratic_term;
164 if (looped_angle > 0) {
165 quadratic_term = -looped_angle_squared;
166 } else {
167 quadratic_term = looped_angle_squared;
168 }
169
170 double linear_term = M_PI * looped_angle;
171
172 return quadratic_term + linear_term;
173}
174
175double tan(double angle)
176{
177 return ampsin(angle) / ampsin(M_PI_2 + angle);
178}
179
180double sqrt(double x)
181{
182 double res;
183 __asm__("fsqrt"
184 : "=t"(res)
185 : "0"(x));
186 return res;
187}
188
189float sqrtf(float x)
190{
191 float res;
192 __asm__("fsqrt"
193 : "=t"(res)
194 : "0"(x));
195 return res;
196}
197
198double sinh(double x)
199{
200 double exponentiated = exp(x);
201 if (x > 0)
202 return (exponentiated * exponentiated - 1) / 2 / exponentiated;
203 return (exponentiated - 1 / exponentiated) / 2;
204}
205
206double log10(double x)
207{
208 return log(x) / M_LN10;
209}
210
211double log(double x)
212{
213 if (x < 0)
214 return __builtin_nan("");
215 if (x == 0)
216 return -__builtin_huge_val();
217 double y = 1 + 2 * (x - 1) / (x + 1);
218 double exponentiated = exp(y);
219 y = y + 2 * (x - exponentiated) / (x + exponentiated);
220 exponentiated = exp(y);
221 y = y + 2 * (x - exponentiated) / (x + exponentiated);
222 exponentiated = exp(y);
223 return y + 2 * (x - exponentiated) / (x + exponentiated);
224}
225
226float logf(float x)
227{
228 return (float)log(x);
229}
230
231double fmod(double index, double period)
232{
233 return index - trunc(index / period) * period;
234}
235
236double exp(double exponent)
237{
238 double result = 1;
239 if (exponent >= 1) {
240 size_t integer_part = (size_t)exponent;
241 if (integer_part & 1)
242 result *= e_to_power<1>();
243 if (integer_part & 2)
244 result *= e_to_power<2>();
245 if (integer_part > 3) {
246 if (integer_part & 4)
247 result *= e_to_power<4>();
248 if (integer_part & 8)
249 result *= e_to_power<8>();
250 if (integer_part & 16)
251 result *= e_to_power<16>();
252 if (integer_part & 32)
253 result *= e_to_power<32>();
254 if (integer_part >= 64)
255 return __builtin_huge_val();
256 }
257 exponent -= integer_part;
258 } else if (exponent < 0)
259 return 1 / exp(-exponent);
260 double taylor_series_result = 1 + exponent;
261 double taylor_series_numerator = exponent * exponent;
262 taylor_series_result += taylor_series_numerator / factorial<2>();
263 taylor_series_numerator *= exponent;
264 taylor_series_result += taylor_series_numerator / factorial<3>();
265 taylor_series_numerator *= exponent;
266 taylor_series_result += taylor_series_numerator / factorial<4>();
267 taylor_series_numerator *= exponent;
268 taylor_series_result += taylor_series_numerator / factorial<5>();
269 return result * taylor_series_result;
270}
271
272float expf(float exponent)
273{
274 return (float)exp(exponent);
275}
276
277double cosh(double x)
278{
279 double exponentiated = exp(-x);
280 if (x < 0)
281 return (1 + exponentiated * exponentiated) / 2 / exponentiated;
282 return (1 / exponentiated + exponentiated) / 2;
283}
284
285double atan2(double y, double x)
286{
287 if (x > 0)
288 return atan(y / x);
289 if (x == 0) {
290 if (y > 0)
291 return M_PI_2;
292 if (y < 0)
293 return -M_PI_2;
294 return 0;
295 }
296 if (y >= 0)
297 return atan(y / x) + M_PI;
298 return atan(y / x) - M_PI;
299}
300
301float atan2f(float y, float x)
302{
303 return (float)atan2(y, x);
304}
305
306double atan(double x)
307{
308 if (x < 0)
309 return -atan(-x);
310 if (x > 1)
311 return M_PI_2 - atan(1 / x);
312 double squared = x * x;
313 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)))))));
314}
315
316double asin(double x)
317{
318 if (x > 1 || x < -1)
319 return __builtin_nan("");
320 if (x > 0.5 || x < -0.5)
321 return 2 * atan(x / (1 + sqrt(1 - x * x)));
322 double squared = x * x;
323 double value = x;
324 double i = x * squared;
325 value += i * product_odd<1>() / product_even<2>() / 3;
326 i *= squared;
327 value += i * product_odd<3>() / product_even<4>() / 5;
328 i *= squared;
329 value += i * product_odd<5>() / product_even<6>() / 7;
330 i *= squared;
331 value += i * product_odd<7>() / product_even<8>() / 9;
332 i *= squared;
333 value += i * product_odd<9>() / product_even<10>() / 11;
334 i *= squared;
335 value += i * product_odd<11>() / product_even<12>() / 13;
336 return value;
337}
338
339float asinf(float x)
340{
341 return (float)asin(x);
342}
343
344double acos(double x)
345{
346 return M_PI_2 - asin(x);
347}
348
349float acosf(float x)
350{
351 return M_PI_2 - asinf(x);
352}
353
354double fabs(double value)
355{
356 return value < 0 ? -value : value;
357}
358
359double log2(double x)
360{
361 return log(x) / M_LN2;
362}
363
364float log2f(float x)
365{
366 return log2(x);
367}
368
369long double log2l(long double x)
370{
371 return log2(x);
372}
373
374double frexp(double, int*)
375{
376 ASSERT_NOT_REACHED();
377 return 0;
378}
379
380float frexpf(float, int*)
381{
382 ASSERT_NOT_REACHED();
383 return 0;
384}
385
386long double frexpl(long double, int*)
387{
388 ASSERT_NOT_REACHED();
389 return 0;
390}
391
392float roundf(float value)
393{
394 // FIXME: Please fix me. I am naive.
395 if (value >= 0.0f)
396 return (float)(int)(value + 0.5f);
397 return (float)(int)(value - 0.5f);
398}
399
400double floor(double value)
401{
402 return (int)value;
403}
404
405double rint(double value)
406{
407 return (int)roundf(value);
408}
409
410float ceilf(float value)
411{
412 // FIXME: Please fix me. I am naive.
413 int as_int = (int)value;
414 if (value == (float)as_int)
415 return as_int;
416 if (value < 0) {
417 if (as_int == 0)
418 return -0;
419 return as_int;
420 }
421 return as_int + 1;
422}
423
424double ceil(double value)
425{
426 // FIXME: Please fix me. I am naive.
427 int as_int = (int)value;
428 if (value == (double)as_int)
429 return as_int;
430 if (value < 0) {
431 if (as_int == 0)
432 return -0;
433 return as_int;
434 }
435 return as_int + 1;
436}
437
438double modf(double x, double* intpart)
439{
440 *intpart = (double)((int)(x));
441 return x - (int)x;
442}
443}