jcs.org
Serenity Operating System
1/*
2 * Copyright (c) 2018-2020, Andreas Kling <kling@serenityos.org>
3 * Copyright (c) 2021, Mițca Dumitru <dumitru0mitca@gmail.com>
4 * Copyright (c) 2022, the SerenityOS developers.
5 * Copyright (c) 2022, Leon Albrecht <leon.a@serenityos.org>
6 *
7 * SPDX-License-Identifier: BSD-2-Clause
8 */
9
10#include <AK/BuiltinWrappers.h>
11#include <AK/ExtraMathConstants.h>
12#include <AK/FloatingPoint.h>
13#ifndef AK_ARCH_AARCH64
14# include <AK/FPControl.h>
15#endif
16#include <AK/Math.h>
17#include <AK/Platform.h>
18#include <AK/StdLibExtras.h>
19#include <assert.h>
20#include <fenv.h>
21#include <math.h>
22#include <stdint.h>
23
24#if defined(AK_COMPILER_CLANG)
25# pragma clang diagnostic push
26# pragma clang diagnostic ignored "-Wdouble-promotion"
27#endif
28
29template<size_t>
30constexpr double e_to_power();
31template<>
32constexpr double e_to_power<0>() { return 1; }
33template<size_t exponent>
34constexpr double e_to_power() { return M_E * e_to_power<exponent - 1>(); }
35
36template<size_t>
37constexpr size_t factorial();
38template<>
39constexpr size_t factorial<0>() { return 1; }
40template<size_t value>
41constexpr size_t factorial() { return value * factorial<value - 1>(); }
42
43template<size_t>
44constexpr size_t product_even();
45template<>
46constexpr size_t product_even<2>() { return 2; }
47template<size_t value>
48constexpr size_t product_even() { return value * product_even<value - 2>(); }
49
50template<size_t>
51constexpr size_t product_odd();
52template<>
53constexpr size_t product_odd<1>() { return 1; }
54template<size_t value>
55constexpr size_t product_odd() { return value * product_odd<value - 2>(); }
56
57enum class RoundingMode {
58 ToZero = FE_TOWARDZERO,
59 Up = FE_UPWARD,
60 Down = FE_DOWNWARD,
61 ToEven = FE_TONEAREST
62};
63
64// This is much branchier than it really needs to be
65template<typename FloatType>
66static FloatType internal_to_integer(FloatType x, RoundingMode rounding_mode)
67{
68 if (!isfinite(x))
69 return x;
70
71 using Extractor = FloatExtractor<decltype(x)>;
72 Extractor extractor;
73 extractor.d = x;
74
75 auto unbiased_exponent = extractor.exponent - Extractor::exponent_bias;
76
77 bool has_half_fraction = false;
78 bool has_nonhalf_fraction = false;
79 if (unbiased_exponent < 0) {
80 // it was easier to special case [0..1) as it saves us from
81 // handling subnormals, underflows, etc
82 if (unbiased_exponent == -1) {
83 has_half_fraction = true;
84 }
85
86 has_nonhalf_fraction = unbiased_exponent < -1 || extractor.mantissa != 0;
87 extractor.mantissa = 0;
88 extractor.exponent = 0;
89 } else {
90 if (unbiased_exponent >= Extractor::mantissa_bits)
91 return x;
92
93 auto dead_bitcount = Extractor::mantissa_bits - unbiased_exponent;
94 auto dead_mask = (1ull << dead_bitcount) - 1;
95 auto dead_bits = extractor.mantissa & dead_mask;
96 extractor.mantissa &= ~dead_mask;
97
98 auto nonhalf_fraction_mask = dead_mask >> 1;
99 has_nonhalf_fraction = (dead_bits & nonhalf_fraction_mask) != 0;
100 has_half_fraction = (dead_bits & ~nonhalf_fraction_mask) != 0;
101 }
102
103 bool should_round = false;
104 switch (rounding_mode) {
105 case RoundingMode::ToEven:
106 should_round = has_half_fraction;
107 break;
108 case RoundingMode::Up:
109 if (!extractor.sign)
110 should_round = has_nonhalf_fraction || has_half_fraction;
111 break;
112 case RoundingMode::Down:
113 if (extractor.sign)
114 should_round = has_nonhalf_fraction || has_half_fraction;
115 break;
116 case RoundingMode::ToZero:
117 break;
118 }
119
120 if (should_round) {
121 // We could do this ourselves, but this saves us from manually
122 // handling overflow.
123 if (extractor.sign)
124 extractor.d -= static_cast<FloatType>(1.0);
125 else
126 extractor.d += static_cast<FloatType>(1.0);
127 }
128
129 return extractor.d;
130}
131
132// This is much branchier than it really needs to be
133template<typename FloatType>
134static FloatType internal_nextafter(FloatType x, bool up)
135{
136 if (!isfinite(x))
137 return x;
138 using Extractor = FloatExtractor<decltype(x)>;
139 Extractor extractor;
140 extractor.d = x;
141 if (x == 0) {
142 if (!extractor.sign) {
143 extractor.mantissa = 1;
144 extractor.sign = !up;
145 return extractor.d;
146 }
147 if (up) {
148 extractor.sign = false;
149 extractor.mantissa = 1;
150 return extractor.d;
151 }
152 extractor.mantissa = 1;
153 extractor.sign = up != extractor.sign;
154 return extractor.d;
155 }
156 if (up != extractor.sign) {
157 extractor.mantissa++;
158 if (!extractor.mantissa) {
159 // no need to normalize the mantissa as we just hit a power
160 // of two.
161 extractor.exponent++;
162 if (extractor.exponent == Extractor::exponent_max) {
163 extractor.exponent = Extractor::exponent_max - 1;
164 extractor.mantissa = Extractor::mantissa_max;
165 }
166 }
167 return extractor.d;
168 }
169
170 if (!extractor.mantissa) {
171 if (extractor.exponent) {
172 extractor.exponent--;
173 extractor.mantissa = Extractor::mantissa_max;
174 } else {
175 extractor.d = 0;
176 }
177 return extractor.d;
178 }
179
180 extractor.mantissa--;
181 if (extractor.mantissa != Extractor::mantissa_max)
182 return extractor.d;
183 if (extractor.exponent) {
184 extractor.exponent--;
185 // normalize
186 extractor.mantissa <<= 1;
187 } else {
188 if (extractor.sign) {
189 // Negative infinity
190 extractor.mantissa = 0;
191 extractor.exponent = Extractor::exponent_max;
192 }
193 }
194 return extractor.d;
195}
196
197template<typename FloatT>
198static int internal_ilogb(FloatT x) NOEXCEPT
199{
200 if (x == 0)
201 return FP_ILOGB0;
202
203 if (isnan(x))
204 return FP_ILOGNAN;
205
206 if (!isfinite(x))
207 return INT_MAX;
208
209 using Extractor = FloatExtractor<FloatT>;
210
211 Extractor extractor;
212 extractor.d = x;
213
214 return (int)extractor.exponent - Extractor::exponent_bias;
215}
216
217template<typename FloatT>
218static FloatT internal_modf(FloatT x, FloatT* intpart) NOEXCEPT
219{
220 FloatT integer_part = internal_to_integer(x, RoundingMode::ToZero);
221 *intpart = integer_part;
222 auto fraction = x - integer_part;
223 if (signbit(fraction) != signbit(x))
224 fraction = -fraction;
225 return fraction;
226}
227
228template<typename FloatT>
229static FloatT internal_scalbn(FloatT x, int exponent) NOEXCEPT
230{
231 if (x == 0 || !isfinite(x) || isnan(x) || exponent == 0)
232 return x;
233
234 using Extractor = FloatExtractor<FloatT>;
235 Extractor extractor;
236 extractor.d = x;
237
238 if (extractor.exponent != 0) {
239 extractor.exponent = clamp((int)extractor.exponent + exponent, 0, (int)Extractor::exponent_max);
240 return extractor.d;
241 }
242
243 unsigned leading_mantissa_zeroes = extractor.mantissa == 0 ? 32 : count_leading_zeroes(extractor.mantissa);
244 int shift = min((int)leading_mantissa_zeroes, exponent);
245 exponent = max(exponent - shift, 0);
246
247 extractor.exponent <<= shift;
248 extractor.exponent = exponent + 1;
249
250 return extractor.d;
251}
252
253template<typename FloatT>
254static FloatT internal_copysign(FloatT x, FloatT y) NOEXCEPT
255{
256 using Extractor = FloatExtractor<FloatT>;
257 Extractor ex, ey;
258 ex.d = x;
259 ey.d = y;
260 ex.sign = ey.sign;
261 return ex.d;
262}
263
264template<typename FloatT>
265static FloatT internal_gamma(FloatT x) NOEXCEPT
266{
267 if (isnan(x))
268 return (FloatT)NAN;
269
270 if (x == (FloatT)0.0)
271 return signbit(x) ? (FloatT)-INFINITY : (FloatT)INFINITY;
272
273 if (x < (FloatT)0 && (rintl(x) == x || isinf(x)))
274 return (FloatT)NAN;
275
276 if (isinf(x))
277 return (FloatT)INFINITY;
278
279 using Extractor = FloatExtractor<FloatT>;
280 // These constants were obtained through use of WolframAlpha
281 constexpr long long max_integer_whose_factorial_fits = (Extractor::mantissa_bits == FloatExtractor<long double>::mantissa_bits ? 20 : (Extractor::mantissa_bits == FloatExtractor<double>::mantissa_bits ? 18 : (Extractor::mantissa_bits == FloatExtractor<float>::mantissa_bits ? 10 : 0)));
282 static_assert(max_integer_whose_factorial_fits != 0, "internal_gamma needs to be aware of the integer factorial that fits in this floating point type.");
283 if ((int)x == x && x <= max_integer_whose_factorial_fits + 1) {
284 long long result = 1;
285 for (long long cursor = 2; cursor < (long long)x; cursor++)
286 result *= cursor;
287 return (FloatT)result;
288 }
289
290 // Stirling approximation
291 return sqrtl(2.0 * M_PIl / static_cast<long double>(x)) * powl(static_cast<long double>(x) / M_El, static_cast<long double>(x));
292}
293
294extern "C" {
295
296float nanf(char const* s) NOEXCEPT
297{
298 return __builtin_nanf(s);
299}
300
301double nan(char const* s) NOEXCEPT
302{
303 return __builtin_nan(s);
304}
305
306long double nanl(char const* s) NOEXCEPT
307{
308 return __builtin_nanl(s);
309}
310
311#define MAKE_AK_BACKED1(name) \
312 long double name##l(long double arg) NOEXCEPT \
313 { \
314 return AK::name<long double>(arg); \
315 } \
316 double name(double arg) NOEXCEPT \
317 { \
318 return AK::name<double>(arg); \
319 } \
320 float name##f(float arg) NOEXCEPT \
321 { \
322 return AK::name<float>(arg); \
323 }
324#define MAKE_AK_BACKED2(name) \
325 long double name##l(long double arg1, long double arg2) NOEXCEPT \
326 { \
327 return AK::name<long double>(arg1, arg2); \
328 } \
329 double name(double arg1, double arg2) NOEXCEPT \
330 { \
331 return AK::name<double>(arg1, arg2); \
332 } \
333 float name##f(float arg1, float arg2) NOEXCEPT \
334 { \
335 return AK::name<float>(arg1, arg2); \
336 }
337
338MAKE_AK_BACKED1(sin);
339MAKE_AK_BACKED1(cos);
340MAKE_AK_BACKED1(tan);
341MAKE_AK_BACKED1(asin);
342MAKE_AK_BACKED1(acos);
343MAKE_AK_BACKED1(atan);
344MAKE_AK_BACKED1(sinh);
345MAKE_AK_BACKED1(cosh);
346MAKE_AK_BACKED1(tanh);
347MAKE_AK_BACKED1(asinh);
348MAKE_AK_BACKED1(acosh);
349MAKE_AK_BACKED1(atanh);
350MAKE_AK_BACKED1(sqrt);
351MAKE_AK_BACKED1(cbrt);
352MAKE_AK_BACKED1(log);
353MAKE_AK_BACKED1(log2);
354MAKE_AK_BACKED1(log10);
355MAKE_AK_BACKED1(exp);
356MAKE_AK_BACKED1(exp2);
357MAKE_AK_BACKED1(fabs);
358
359MAKE_AK_BACKED2(atan2);
360MAKE_AK_BACKED2(hypot);
361MAKE_AK_BACKED2(fmod);
362MAKE_AK_BACKED2(pow);
363MAKE_AK_BACKED2(remainder);
364
365long double truncl(long double x) NOEXCEPT
366{
367#ifndef AK_ARCH_AARCH64
368 if (fabsl(x) < LONG_LONG_MAX) {
369 // This is 1.6 times faster than the implementation using the "internal_to_integer"
370 // helper (on x86_64)
371 // https://quick-bench.com/q/xBmxuY8am9qibSYVna90Y6PIvqA
372 u64 temp;
373 asm(
374 "fisttpq %[temp]\n"
375 "fildq %[temp]"
376 : "+t"(x)
377 : [temp] "m"(temp));
378 return x;
379 }
380#endif
381
382 return internal_to_integer(x, RoundingMode::ToZero);
383}
384
385double trunc(double x) NOEXCEPT
386{
387#ifndef AK_ARCH_AARCH64
388 if (fabs(x) < LONG_LONG_MAX) {
389 u64 temp;
390 asm(
391 "fisttpq %[temp]\n"
392 "fildq %[temp]"
393 : "+t"(x)
394 : [temp] "m"(temp));
395 return x;
396 }
397#endif
398
399 return internal_to_integer(x, RoundingMode::ToZero);
400}
401
402float truncf(float x) NOEXCEPT
403{
404#ifndef AK_ARCH_AARCH64
405 if (fabsf(x) < LONG_LONG_MAX) {
406 u64 temp;
407 asm(
408 "fisttpq %[temp]\n"
409 "fildq %[temp]"
410 : "+t"(x)
411 : [temp] "m"(temp));
412 return x;
413 }
414#endif
415
416 return internal_to_integer(x, RoundingMode::ToZero);
417}
418
419long double rintl(long double value)
420{
421#ifdef AK_ARCH_AARCH64
422 (void)value;
423 TODO_AARCH64();
424#else
425 long double res;
426 asm(
427 "frndint\n"
428 : "=t"(res)
429 : "0"(value));
430 return res;
431#endif
432}
433double rint(double value)
434{
435#ifdef AK_ARCH_AARCH64
436 (void)value;
437 TODO_AARCH64();
438#else
439 double res;
440 asm(
441 "frndint\n"
442 : "=t"(res)
443 : "0"(value));
444 return res;
445#endif
446}
447float rintf(float value)
448{
449#ifdef AK_ARCH_AARCH64
450 (void)value;
451 TODO_AARCH64();
452#else
453 float res;
454 asm(
455 "frndint\n"
456 : "=t"(res)
457 : "0"(value));
458 return res;
459#endif
460}
461
462long lrintl(long double value)
463{
464#ifdef AK_ARCH_AARCH64
465 (void)value;
466 TODO_AARCH64();
467#else
468 long res;
469 asm(
470 "fistpl %0\n"
471 : "+m"(res)
472 : "t"(value)
473 : "st");
474 return res;
475#endif
476}
477long lrint(double value)
478{
479#ifdef AK_ARCH_AARCH64
480 (void)value;
481 TODO_AARCH64();
482#else
483 long res;
484 asm(
485 "fistpl %0\n"
486 : "+m"(res)
487 : "t"(value)
488 : "st");
489 return res;
490#endif
491}
492long lrintf(float value)
493{
494#ifdef AK_ARCH_AARCH64
495 (void)value;
496 TODO_AARCH64();
497#else
498 long res;
499 asm(
500 "fistpl %0\n"
501 : "+m"(res)
502 : "t"(value)
503 : "st");
504 return res;
505#endif
506}
507
508long long llrintl(long double value)
509{
510#ifdef AK_ARCH_AARCH64
511 (void)value;
512 TODO_AARCH64();
513#else
514 long long res;
515 asm(
516 "fistpq %0\n"
517 : "+m"(res)
518 : "t"(value)
519 : "st");
520 return res;
521#endif
522}
523long long llrint(double value)
524{
525#ifdef AK_ARCH_AARCH64
526 (void)value;
527 TODO_AARCH64();
528#else
529 long long res;
530 asm(
531 "fistpq %0\n"
532 : "+m"(res)
533 : "t"(value)
534 : "st");
535 return res;
536#endif
537}
538long long llrintf(float value)
539{
540#ifdef AK_ARCH_AARCH64
541 (void)value;
542 TODO_AARCH64();
543#else
544 long long res;
545 asm(
546 "fistpq %0\n"
547 : "+m"(res)
548 : "t"(value)
549 : "st");
550 return res;
551#endif
552}
553
554// On systems where FLT_RADIX == 2, ldexp is equivalent to scalbn
555long double ldexpl(long double x, int exp) NOEXCEPT
556{
557 return internal_scalbn(x, exp);
558}
559
560double ldexp(double x, int exp) NOEXCEPT
561{
562 return internal_scalbn(x, exp);
563}
564
565float ldexpf(float x, int exp) NOEXCEPT
566{
567 return internal_scalbn(x, exp);
568}
569
570[[maybe_unused]] static long double ampsin(long double angle) NOEXCEPT
571{
572 long double looped_angle = fmodl(M_PI + angle, M_TAU) - M_PI;
573 long double looped_angle_squared = looped_angle * looped_angle;
574
575 long double quadratic_term;
576 if (looped_angle > 0) {
577 quadratic_term = -looped_angle_squared;
578 } else {
579 quadratic_term = looped_angle_squared;
580 }
581
582 long double linear_term = M_PI * looped_angle;
583
584 return quadratic_term + linear_term;
585}
586
587int ilogbl(long double x) NOEXCEPT
588{
589 return internal_ilogb(x);
590}
591
592int ilogb(double x) NOEXCEPT
593{
594 return internal_ilogb(x);
595}
596
597int ilogbf(float x) NOEXCEPT
598{
599 return internal_ilogb(x);
600}
601
602long double logbl(long double x) NOEXCEPT
603{
604 return ilogbl(x);
605}
606
607double logb(double x) NOEXCEPT
608{
609 return ilogb(x);
610}
611
612float logbf(float x) NOEXCEPT
613{
614 return ilogbf(x);
615}
616
617double frexp(double x, int* exp) NOEXCEPT
618{
619 *exp = (x == 0) ? 0 : (1 + ilogb(x));
620 return scalbn(x, -(*exp));
621}
622
623float frexpf(float x, int* exp) NOEXCEPT
624{
625 *exp = (x == 0) ? 0 : (1 + ilogbf(x));
626 return scalbnf(x, -(*exp));
627}
628
629long double frexpl(long double x, int* exp) NOEXCEPT
630{
631 *exp = (x == 0) ? 0 : (1 + ilogbl(x));
632 return scalbnl(x, -(*exp));
633}
634
635#if !(ARCH(X86_64))
636
637double round(double value) NOEXCEPT
638{
639 return internal_to_integer(value, RoundingMode::ToEven);
640}
641
642float roundf(float value) NOEXCEPT
643{
644 return internal_to_integer(value, RoundingMode::ToEven);
645}
646
647long double roundl(long double value) NOEXCEPT
648{
649 return internal_to_integer(value, RoundingMode::ToEven);
650}
651
652long lroundf(float value) NOEXCEPT
653{
654 return internal_to_integer(value, RoundingMode::ToEven);
655}
656
657long lround(double value) NOEXCEPT
658{
659 return internal_to_integer(value, RoundingMode::ToEven);
660}
661
662long lroundl(long double value) NOEXCEPT
663{
664 return internal_to_integer(value, RoundingMode::ToEven);
665}
666
667long long llroundf(float value) NOEXCEPT
668{
669 return internal_to_integer(value, RoundingMode::ToEven);
670}
671
672long long llround(double value) NOEXCEPT
673{
674 return internal_to_integer(value, RoundingMode::ToEven);
675}
676
677long long llroundd(long double value) NOEXCEPT
678{
679 return internal_to_integer(value, RoundingMode::ToEven);
680}
681
682float floorf(float value) NOEXCEPT
683{
684 return internal_to_integer(value, RoundingMode::Down);
685}
686
687double floor(double value) NOEXCEPT
688{
689 return internal_to_integer(value, RoundingMode::Down);
690}
691
692long double floorl(long double value) NOEXCEPT
693{
694 return internal_to_integer(value, RoundingMode::Down);
695}
696
697float ceilf(float value) NOEXCEPT
698{
699 return internal_to_integer(value, RoundingMode::Up);
700}
701
702double ceil(double value) NOEXCEPT
703{
704 return internal_to_integer(value, RoundingMode::Up);
705}
706
707long double ceill(long double value) NOEXCEPT
708{
709 return internal_to_integer(value, RoundingMode::Up);
710}
711
712#else
713
714double round(double x) NOEXCEPT
715{
716 // Note: This is break-tie-away-from-zero, so not the hw's understanding of
717 // "nearest", which would be towards even.
718 if (x == 0.)
719 return x;
720 if (x > 0.)
721 return floor(x + .5);
722 return ceil(x - .5);
723}
724
725float roundf(float x) NOEXCEPT
726{
727 if (x == 0.f)
728 return x;
729 if (x > 0.f)
730 return floorf(x + .5f);
731 return ceilf(x - .5f);
732}
733
734long double roundl(long double x) NOEXCEPT
735{
736 if (x == 0.L)
737 return x;
738 if (x > 0.L)
739 return floorl(x + .5L);
740 return ceill(x - .5L);
741}
742
743long lroundf(float value) NOEXCEPT
744{
745 return static_cast<long>(roundf(value));
746}
747
748long lround(double value) NOEXCEPT
749{
750 return static_cast<long>(round(value));
751}
752
753long lroundl(long double value) NOEXCEPT
754{
755 return static_cast<long>(roundl(value));
756}
757
758long long llroundf(float value) NOEXCEPT
759{
760 return static_cast<long long>(roundf(value));
761}
762
763long long llround(double value) NOEXCEPT
764{
765 return static_cast<long long>(round(value));
766}
767
768long long llroundd(long double value) NOEXCEPT
769{
770 return static_cast<long long>(roundl(value));
771}
772
773float floorf(float value) NOEXCEPT
774{
775 AK::X87RoundingModeScope scope { AK::RoundingMode::DOWN };
776 asm("frndint"
777 : "+t"(value));
778 return value;
779}
780
781double floor(double value) NOEXCEPT
782{
783 AK::X87RoundingModeScope scope { AK::RoundingMode::DOWN };
784 asm("frndint"
785 : "+t"(value));
786 return value;
787}
788
789long double floorl(long double value) NOEXCEPT
790{
791 AK::X87RoundingModeScope scope { AK::RoundingMode::DOWN };
792 asm("frndint"
793 : "+t"(value));
794 return value;
795}
796
797float ceilf(float value) NOEXCEPT
798{
799 AK::X87RoundingModeScope scope { AK::RoundingMode::UP };
800 asm("frndint"
801 : "+t"(value));
802 return value;
803}
804
805double ceil(double value) NOEXCEPT
806{
807 AK::X87RoundingModeScope scope { AK::RoundingMode::UP };
808 asm("frndint"
809 : "+t"(value));
810 return value;
811}
812
813long double ceill(long double value) NOEXCEPT
814{
815 AK::X87RoundingModeScope scope { AK::RoundingMode::UP };
816 asm("frndint"
817 : "+t"(value));
818 return value;
819}
820
821#endif
822
823long double modfl(long double x, long double* intpart) NOEXCEPT
824{
825 return internal_modf(x, intpart);
826}
827
828double modf(double x, double* intpart) NOEXCEPT
829{
830 return internal_modf(x, intpart);
831}
832
833float modff(float x, float* intpart) NOEXCEPT
834{
835 return internal_modf(x, intpart);
836}
837
838double gamma(double x) NOEXCEPT
839{
840 // Stirling approximation
841 return sqrt(2.0 * M_PI / x) * pow(x / M_E, x);
842}
843
844long double tgammal(long double value) NOEXCEPT
845{
846 return internal_gamma(value);
847}
848
849double tgamma(double value) NOEXCEPT
850{
851 return internal_gamma(value);
852}
853
854float tgammaf(float value) NOEXCEPT
855{
856 return internal_gamma(value);
857}
858
859int signgam = 0;
860
861long double lgammal(long double value) NOEXCEPT
862{
863 return lgammal_r(value, &signgam);
864}
865
866double lgamma(double value) NOEXCEPT
867{
868 return lgamma_r(value, &signgam);
869}
870
871float lgammaf(float value) NOEXCEPT
872{
873 return lgammaf_r(value, &signgam);
874}
875
876long double lgammal_r(long double value, int* sign) NOEXCEPT
877{
878 if (value == 1.0 || value == 2.0)
879 return 0.0;
880 if (isinf(value) || value == 0.0)
881 return INFINITY;
882 long double result = logl(internal_gamma(value));
883 *sign = signbit(result) ? -1 : 1;
884 return result;
885}
886
887double lgamma_r(double value, int* sign) NOEXCEPT
888{
889 if (value == 1.0 || value == 2.0)
890 return 0.0;
891 if (isinf(value) || value == 0.0)
892 return INFINITY;
893 double result = log(internal_gamma(value));
894 *sign = signbit(result) ? -1 : 1;
895 return result;
896}
897
898float lgammaf_r(float value, int* sign) NOEXCEPT
899{
900 if (value == 1.0f || value == 2.0f)
901 return 0.0;
902 if (isinf(value) || value == 0.0f)
903 return INFINITY;
904 float result = logf(internal_gamma(value));
905 *sign = signbit(result) ? -1 : 1;
906 return result;
907}
908
909long double expm1l(long double x) NOEXCEPT
910{
911 return expl(x) - 1;
912}
913
914double expm1(double x) NOEXCEPT
915{
916 return exp(x) - 1;
917}
918
919float expm1f(float x) NOEXCEPT
920{
921 return expf(x) - 1;
922}
923
924long double log1pl(long double x) NOEXCEPT
925{
926 return logl(1 + x);
927}
928
929double log1p(double x) NOEXCEPT
930{
931 return log(1 + x);
932}
933
934float log1pf(float x) NOEXCEPT
935{
936 return logf(1 + x);
937}
938
939long double erfl(long double x) NOEXCEPT
940{
941 // algorithm taken from Abramowitz and Stegun (no. 26.2.17)
942 long double t = 1 / (1 + 0.47047l * fabsl(x));
943 long double poly = t * (0.3480242l + t * (-0.958798l + t * 0.7478556l));
944 long double answer = 1 - poly * expl(-x * x);
945 if (x < 0)
946 return -answer;
947
948 return answer;
949}
950
951double erf(double x) NOEXCEPT
952{
953 return (double)erfl(x);
954}
955
956float erff(float x) NOEXCEPT
957{
958 return (float)erf(x);
959}
960
961long double erfcl(long double x) NOEXCEPT
962{
963 return 1 - erfl(x);
964}
965
966double erfc(double x) NOEXCEPT
967{
968 return 1 - erf(x);
969}
970
971float erfcf(float x) NOEXCEPT
972{
973 return 1 - erff(x);
974}
975
976double nextafter(double x, double target) NOEXCEPT
977{
978 if (x == target)
979 return target;
980 return internal_nextafter(x, target >= x);
981}
982
983float nextafterf(float x, float target) NOEXCEPT
984{
985 if (x == target)
986 return target;
987 return internal_nextafter(x, target >= x);
988}
989
990long double nextafterl(long double x, long double target) NOEXCEPT
991{
992 return internal_nextafter(x, target >= x);
993}
994
995double nexttoward(double x, long double target) NOEXCEPT
996{
997 if (x == target)
998 return target;
999 return internal_nextafter(x, target >= x);
1000}
1001
1002float nexttowardf(float x, long double target) NOEXCEPT
1003{
1004 if (x == target)
1005 return target;
1006 return internal_nextafter(x, target >= x);
1007}
1008
1009long double nexttowardl(long double x, long double target) NOEXCEPT
1010{
1011 if (x == target)
1012 return target;
1013 return internal_nextafter(x, target >= x);
1014}
1015
1016float copysignf(float x, float y) NOEXCEPT
1017{
1018 return internal_copysign(x, y);
1019}
1020
1021double copysign(double x, double y) NOEXCEPT
1022{
1023 return internal_copysign(x, y);
1024}
1025
1026long double copysignl(long double x, long double y) NOEXCEPT
1027{
1028 return internal_copysign(x, y);
1029}
1030
1031float scalbnf(float x, int exponent) NOEXCEPT
1032{
1033 return internal_scalbn(x, exponent);
1034}
1035
1036double scalbn(double x, int exponent) NOEXCEPT
1037{
1038 return internal_scalbn(x, exponent);
1039}
1040
1041long double scalbnl(long double x, int exponent) NOEXCEPT
1042{
1043 return internal_scalbn(x, exponent);
1044}
1045
1046float scalbnlf(float x, long exponent) NOEXCEPT
1047{
1048 return internal_scalbn(x, exponent);
1049}
1050
1051double scalbln(double x, long exponent) NOEXCEPT
1052{
1053 return internal_scalbn(x, exponent);
1054}
1055
1056long double scalblnl(long double x, long exponent) NOEXCEPT
1057{
1058 return internal_scalbn(x, exponent);
1059}
1060
1061long double fmaxl(long double x, long double y) NOEXCEPT
1062{
1063 if (isnan(x))
1064 return y;
1065 if (isnan(y))
1066 return x;
1067
1068 return x > y ? x : y;
1069}
1070
1071double fmax(double x, double y) NOEXCEPT
1072{
1073 if (isnan(x))
1074 return y;
1075 if (isnan(y))
1076 return x;
1077
1078 return x > y ? x : y;
1079}
1080
1081float fmaxf(float x, float y) NOEXCEPT
1082{
1083 if (isnan(x))
1084 return y;
1085 if (isnan(y))
1086 return x;
1087
1088 return x > y ? x : y;
1089}
1090
1091long double fminl(long double x, long double y) NOEXCEPT
1092{
1093 if (isnan(x))
1094 return y;
1095 if (isnan(y))
1096 return x;
1097
1098 return x < y ? x : y;
1099}
1100
1101double fmin(double x, double y) NOEXCEPT
1102{
1103 if (isnan(x))
1104 return y;
1105 if (isnan(y))
1106 return x;
1107
1108 return x < y ? x : y;
1109}
1110
1111float fminf(float x, float y) NOEXCEPT
1112{
1113 if (isnan(x))
1114 return y;
1115 if (isnan(y))
1116 return x;
1117
1118 return x < y ? x : y;
1119}
1120
1121// https://pubs.opengroup.org/onlinepubs/9699919799/functions/fma.html
1122long double fmal(long double x, long double y, long double z) NOEXCEPT
1123{
1124 return (x * y) + z;
1125}
1126
1127double fma(double x, double y, double z) NOEXCEPT
1128{
1129 return (x * y) + z;
1130}
1131
1132float fmaf(float x, float y, float z) NOEXCEPT
1133{
1134 return (x * y) + z;
1135}
1136
1137long double nearbyintl(long double value) NOEXCEPT
1138{
1139 return internal_to_integer(value, RoundingMode { fegetround() });
1140}
1141
1142double nearbyint(double value) NOEXCEPT
1143{
1144 return internal_to_integer(value, RoundingMode { fegetround() });
1145}
1146
1147float nearbyintf(float value) NOEXCEPT
1148{
1149 return internal_to_integer(value, RoundingMode { fegetround() });
1150}
1151}
1152
1153#if defined(AK_COMPILER_CLANG)
1154# pragma clang diagnostic pop
1155#endif