nulib.math.intrinsics source code

1 /**
2     Bindings to compiler-specific intrinsics wrapped in a nicer interface.
3 
4     Copyright:
5         Copyright © 2023-2025, Kitsunebi Games
6         Copyright © 2023-2025, Inochi2D Project
7     
8     License:   $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
9     Authors:
10         Luna Nielsen
11 */
12 module nulib.math.intrinsics;
13 import nulib.math.fixed;
14 
15 version (GNU) import gcc.builtins;
16 else version (LDC) import ldc.intrinsics;
17 else {
18     import nmath = nulib.math.floating;
19     import cmath = nulib.c.math;
20 }
21 
22 @safe @nogc nothrow pure:
23 
24 /**
25     Computes the square root of the given value.
26 
27     Params:
28         x = The value
29     
30     Returns:
31         The square root of $(D x).
32 */
33 pragma(inline, true)
34 T sqrt(T)(T x) if (__traits(isFloating, T)) {
35     version(LDC) {
36         return x < 0 ? T.nan : llvm_sqrt(x);
37     }  else version(GNU) {
38         static if (is(T == float))
39             return __builtin_sqrtf(x);
40         else static if (is(T == double))
41             return __builtin_sqrt(x);
42         else static if (is(T == real))
43             return __builtin_sqrtl(x);
44     } else {
45         return cast(T)cmath.sqrt(cast(double)x);
46     }
47 }
48 
49 @("sqrt")
50 unittest {
51     assert(sqrt(9.0) == 3.0);
52     assert(sqrt(10.0) == 3.1622776601683795);
53 }
54 
55 /**
56     Computes sine of the given value.
57 
58     Params:
59         x = The value
60     
61     Returns:
62         The sine of $(D x).
63 */
64 pragma(inline, true)
65 T sin(T)(T x) if (__traits(isFloating, T)) {
66     version(LDC) {
67         return llvm_sin!T(x);
68     } else version(GNU) {
69         static if (is(T == float))
70             return __builtin_sinf(x);
71         else static if (is(T == double))
72             return __builtin_sin(x);
73         else static if (is(T == real))
74             return __builtin_sinl(x);
75     } else {
76         return cast(T)cmath.sin(cast(double)x);
77     }
78 }
79 
80 /**
81     Computes cosine of the given value.
82 
83     Params:
84         x = The value
85     
86     Returns:
87         The cosine of $(D x).
88 */
89 pragma(inline, true)
90 T cos(T)(T x) if (__traits(isFloating, T)) {
91     version(LDC) {
92         return llvm_cos!T(x);
93     } else version(GNU) {
94         static if (is(T == float))
95             return __builtin_cosf(x);
96         else static if (is(T == double))
97             return __builtin_cos(x);
98         else static if (is(T == real))
99             return __builtin_cosl(x);
100     } else {
101         return cast(T)cmath.cos(cast(double)x);
102     }
103 }
104 
105 /**
106     Computes tangent of the given value.
107 
108     Params:
109         x = The value
110     
111     Returns:
112         The tangent of $(D x).
113 */
114 pragma(inline, true)
115 T tan(T)(T x) if (__traits(isFloating, T)) {
116     version(LDC) {
117         return llvm_tan!T(x);
118     } else version(GNU) {
119         static if (is(T == float))
120             return __builtin_tanf(x);
121         else static if (is(T == double))
122             return __builtin_tan(x);
123         else static if (is(T == real))
124             return __builtin_tanl(x);
125     } else {
126         return cast(T)cmath.tan(cast(double)x);
127     }
128 }
129 
130 /**
131     Computes arc-sine of the given value.
132 
133     Params:
134         x = The value
135     
136     Returns:
137         The arc-sine of $(D x).
138 */
139 pragma(inline, true)
140 T asin(T)(T x) if (__traits(isFloating, T)) {
141     version(LDC) {
142         return llvm_asin!T(x);
143     } else version(GNU) {
144         static if (is(T == float))
145             return __builtin_asinf(x);
146         else static if (is(T == double))
147             return __builtin_asin(x);
148         else static if (is(T == real))
149             return __builtin_asinl(x);
150     } else {
151         return cast(T)cmath.asin(cast(double)x);
152     }
153 }
154 
155 /**
156     Computes arc-cosine of the given value.
157 
158     Params:
159         x = The value
160     
161     Returns:
162         The arc-cosine of $(D x).
163 */
164 pragma(inline, true)
165 T acos(T)(T x) if (__traits(isFloating, T)) {
166     version(LDC) {
167         return llvm_acos!T(x);
168     } else version(GNU) {
169         static if (is(T == float))
170             return __builtin_acosf(x);
171         else static if (is(T == double))
172             return __builtin_acos(x);
173         else static if (is(T == real))
174             return __builtin_acosl(x);
175     } else {
176         return cast(T)cmath.acos(cast(double)x);
177     }
178 }
179 
180 /**
181     Computes arc-tangent of the given value.
182 
183     Params:
184         x = The value
185     
186     Returns:
187         The arc-tangent of $(D x).
188 */
189 pragma(inline, true)
190 T atan(T)(T x) if (__traits(isFloating, T)) {
191     version(LDC) {
192         return llvm_atan!T(x);
193     } else version(GNU) {
194         static if (is(T == float))
195             return __builtin_atanf(x);
196         else static if (is(T == double))
197             return __builtin_atan(x);
198         else static if (is(T == real))
199             return __builtin_atanl(x);
200     } else version(DigitalMars) {
201         return nmath.atan(x); 
202     } else {
203         return cast(T)cmath.atan(cast(double)x);
204     }
205 }
206 
207 /**
208     Computes arc-tangent of the given value, using signs to determine quadrant.
209 
210     Params:
211         y = value
212         x = value
213     
214     Returns:
215         The arc-tangent of $(D y / x).
216 */
217 pragma(inline, true)
218 T atan2(T)(T y, T x) if (__traits(isFloating, T)) {
219     version(LDC) {
220         return llvm_atan2!T(y, x);
221     } else version(GNU) {
222         static if (is(T == float))
223             return __builtin_atan2f(y, x);
224         else static if (is(T == double))
225             return __builtin_atan2(y, x);
226         else static if (is(T == real))
227             return __builtin_atan2l(y, x);
228     } else version(DigitalMars) {
229         return nmath.atan2(y, x); 
230     } else {
231         return cast(T)cmath.atan2(cast(double)x, cast(double)y);
232     }
233 }
234 
235 /**
236     Computes hyperbolic sine of the given value.
237 
238     Params:
239         x = The value
240     
241     Returns:
242         The hyperbolic sine of $(D x).
243 */
244 pragma(inline, true)
245 T sinh(T)(T x) if (__traits(isFloating, T)) {
246     version(LDC) {
247         return llvm_sinh!T(x);
248     } else version(GNU) {
249         static if (is(T == float))
250             return __builtin_sinhf(x);
251         else static if (is(T == double))
252             return __builtin_sinh(x);
253         else static if (is(T == real))
254             return __builtin_sinhl(x);
255     } else {
256         return cast(T)cmath.sinh(cast(double)x);
257     }
258 }
259 
260 /**
261     Computes hyperbolic cosine of the given value.
262 
263     Params:
264         x = The value
265     
266     Returns:
267         The hyperbolic cosine of $(D x).
268 */
269 pragma(inline, true)
270 T cosh(T)(T x) if (__traits(isFloating, T)) {
271     version(LDC) {
272         return llvm_cosh!T(x);
273     } else version(GNU) {
274         static if (is(T == float))
275             return __builtin_coshf(x);
276         else static if (is(T == double))
277             return __builtin_cohs(x);
278         else static if (is(T == real))
279             return __builtin_coshl(x);
280     } else {
281         return cast(T)cmath.cosh(cast(double)x);
282     }
283 }
284 
285 /**
286     Computes hyperbolic tangent of the given value.
287 
288     Params:
289         x = The value
290     
291     Returns:
292         The hyperbolic tangent of $(D x).
293 */
294 pragma(inline, true)
295 T tanh(T)(T x) if (__traits(isFloating, T)) {
296     version(LDC) {
297         return llvm_tanh!T(x);
298     } else version(GNU) {
299         static if (is(T == float))
300             return __builtin_tanhf(x);
301         else static if (is(T == double))
302             return __builtin_tanh(x);
303         else static if (is(T == real))
304             return __builtin_tanhl(x);
305     } else {
306         return cast(T)cmath.tanh(cast(double)x);
307     }
308 }
309 
310 /**
311     Computes hyperbolic arc-sine of the given value.
312 
313     Params:
314         x = The value
315     
316     Returns:
317         The hyperbolic arc-sine of $(D x).
318 */
319 pragma(inline, true)
320 T asinh(T)(T x) if (__traits(isFloating, T)) {
321     version(LDC) {
322         return llvm_asin!T(x);
323     } else version(GNU) {
324         static if (is(T == float))
325             return __builtin_asinhf(x);
326         else static if (is(T == double))
327             return __builtin_asinh(x);
328         else static if (is(T == real))
329             return __builtin_asinhl(x);
330     } else {
331         return cast(T)cmath.asinh(cast(double)x);
332     }
333 }
334 
335 /**
336     Computes hyperbolic arc-cosine of the given value.
337 
338     Params:
339         x = The value
340     
341     Returns:
342         The hyperbolic arc-cosine of $(D x).
343 */
344 pragma(inline, true)
345 T acosh(T)(T x) if (__traits(isFloating, T)) {
346     version(LDC) {
347         return llvm_acosh!T(x);
348     } else version(GNU) {
349         static if (is(T == float))
350             return __builtin_acoshf(x);
351         else static if (is(T == double))
352             return __builtin_acosh(x);
353         else static if (is(T == real))
354             return __builtin_acoshl(x);
355     } else {
356         return cast(T)cmath.acosh(cast(double)x);
357     }
358 }
359 
360 /**
361     Computes hyperbolic arc-tangent of the given value.
362 
363     Params:
364         x = The value
365     
366     Returns:
367         The hyperbolic arc-tangent of $(D x).
368 */
369 pragma(inline, true)
370 T atanh(T)(T x) if (__traits(isFloating, T)) {
371     version(LDC) {
372         return llvm_atanh!T(x);
373     } else version(GNU) {
374         static if (is(T == float))
375             return __builtin_atanhf(x);
376         else static if (is(T == double))
377             return __builtin_atanh(x);
378         else static if (is(T == real))
379             return __builtin_atanhl(x);
380     } else {
381         return cast(T)cmath.atanh(cast(double)x);
382     }
383 }
384 
385 /**
386     Computes the nearest integer value lower in magnitude than
387     the given value.
388 
389     Params:
390         x = The value
391     
392     Returns:
393         The nearest integer value lower in magnitude than $(D x).
394 */
395 pragma(inline, true)
396 T trunc(T)(T x) if (__traits(isFloating, T)) {
397     version(LDC) {
398         return llvm_trunc!T(x);
399     } else version(GNU) {
400         static if (is(T == float))
401             return __builtin_truncf(x);
402         else static if (is(T == double))
403             return __builtin_trunc(x);
404         else static if (is(T == real))
405             return __builtin_truncl(x);
406     } else {
407         return cast(T)cmath.trunc(cast(double)x);
408     }
409 }
410 
411 @("trunc")
412 unittest {
413     assert(trunc(1.5) == 1);
414     assert(trunc(1.9999991) == 1);
415     assert(trunc(1.0000001) == 1);
416     assert(trunc(0.9999991) == 0);
417 
418     assert(trunc(cast(float)1.0) == 1);
419     assert(trunc(cast(double)1.0) == 1);
420     assert(trunc(cast(real)1.0) == 1);
421 }
422 
423 /**
424     Computes the nearest integer value, rounded away from 0.
425 
426     Params:
427         value = Input value
428     
429     Returns:
430         The nearest integer value to $(D x).
431 */
432 T round(T)(T value) if (__traits(isFloating, T)) {
433     version(LDC) {
434         return llvm_round!T(value);
435     } else version(GNU) {
436         static if (is(T == float))
437             return __builtin_roundf(value);
438         else static if (is(T == double))
439             return __builtin_round(value);
440         else static if (is(T == real))
441             return __builtin_roundl(value);
442     } else {
443         return cast(T)cmath.round(cast(double)value);
444     }
445 }
446 
447 @("round")
448 unittest {
449     assert(round(0.0) == 0);
450     assert(round(0.25) == 0);
451     assert(round(0.5) == 1);
452     assert(round(0.95) == 1);
453 
454     assert(round(cast(float)1.0) == 1);
455     assert(round(cast(double)1.0) == 1);
456     assert(round(cast(real)1.0) == 1);
457 }
458 
459 /**
460     Computes the nearest integer value lower than the given value.
461 
462     Params:
463         value = Input value
464     
465     Returns:
466         The nearest integer value lower than $(D x).
467 */
468 pragma(inline, true)
469 T floor(T)(T value) if (__traits(isFloating, T)) {
470     version(LDC) {
471         return llvm_floor!T(value);
472     } else version(GNU) {
473         static if (is(T == float))
474             return __builtin_floorf(value);
475         else static if (is(T == double))
476             return __builtin_floor(value);
477         else static if (is(T == real))
478             return __builtin_floorl(value);
479     } else {
480         return cast(T)cmath.floor(cast(double)value);
481     }
482 }
483 
484 @("floor")
485 unittest {
486     assert(floor(0.0) == 0);
487     assert(floor(0.25) == 0);
488     assert(floor(0.5) == 0);
489     assert(floor(0.95) == 0);
490 
491     assert(floor(cast(float)1.0) == 1);
492     assert(floor(cast(double)1.0) == 1);
493     assert(floor(cast(real)1.0) == 1);
494 }
495 
496 /**
497     Computes the nearest integer value lower than the given value.
498 
499     Params:
500         value = Input value
501     
502     Returns:
503         The nearest integer value lower than $(D x).
504 */
505 pragma(inline, true)
506 T ceil(T)(T value) if (__traits(isFloating, T)) {
507     version(LDC) {
508         return llvm_ceil!T(value);
509     } else version(GNU) {
510         static if (is(T == float))
511             return __builtin_ceilf(value);
512         else static if (is(T == double))
513             return __builtin_ceil(value);
514         else static if (is(T == real))
515             return __builtin_ceill(value);
516     } else {
517         return cast(T)cmath.ceil(cast(double)value);
518     }
519 }
520 
521 @("ceil")
522 unittest {
523     assert(ceil(0.0) == 0);
524     assert(ceil(0.25) == 1);
525     assert(ceil(0.5) == 1);
526     assert(ceil(0.95) == 1);
527     
528     assert(ceil(cast(float)1.0) == 1);
529     assert(ceil(cast(double)1.0) == 1);
530     assert(ceil(cast(real)1.0) == 1);
531 }
532 
533 /**
534     Computes the absolute value for the given value.
535 
536     Params:
537         value = the value
538     
539     Returns:
540         The absolute value of $(D value)
541 */
542 pragma(inline, true)
543 T abs(T)(T value) if (__traits(isScalar, T)) {
544     static if (__traits(isFloating, T)) {
545         version(LDC) {
546             return llvm_fabs!T(value);
547         } else version(GNU) {
548             static if (is(T == float))
549                 return __builtin_fabsf(value);
550             else static if (is(T == double))
551                 return __builtin_fabs(value);
552             else static if (is(T == real))
553                 return __builtin_fabsl(value);
554         } else {
555             return cast(T)cmath.fabs(cast(double)value);
556         }
557     } else {
558         return value < 0 ? -value : value;
559     }
560 }
561 
562 @("abs")
563 unittest {
564     foreach(i; 0..100) {
565         assert(abs(cast(float)-i) == cast(float)i);
566     }
567     
568     assert(abs(cast(float)1.0) == 1);
569     assert(abs(cast(double)1.0) == 1);
570     assert(abs(cast(real)1.0) == 1);
571 }
572 
573 /**
574     Rounds the given value to the nearest integer,
575     using the current rounding mode.
576 
577     Params:
578         x = The value
579     
580     Returns:
581         $(D x) rounded to the nearest integer 
582         in respect to rounding mode.
583 */
584 pragma(inline, true)
585 T rint(T)(T x) if (__traits(isFloating, T)) {
586     version(LDC) {
587         return llvm_rint!T(value);
588     } else version(GNU) {
589         static if (is(T == float))
590             return __builtin_rintf(value);
591         else static if (is(T == double))
592             return __builtin_rint(value);
593         else static if (is(T == real))
594             return __builtin_rintl(value);
595     } else {
596         return cast(T)cmath.rint(cast(double)value);
597     }
598 }
599 
600 /**
601     Computes n * 2$(SUPERSCRIPT exp).
602 
603     Params:
604         n =     The value
605         exp =   The exponent
606     
607     Returns:
608         $(D value * pow(exp, 2)).
609 */
610 pragma(inline, true)
611 T ldexp(T)(T n, int exp) if (__traits(isFloating, T)) {
612     version(LDC) {
613         enum RealFormat { ieeeSingle, ieeeDouble, ieeeExtended, ieeeQuadruple }
614 
615              static if (T.mant_dig ==  24) enum realFormat = RealFormat.ieeeSingle;
616         else static if (T.mant_dig ==  53) enum realFormat = RealFormat.ieeeDouble;
617         else static if (T.mant_dig ==  64) enum realFormat = RealFormat.ieeeExtended;
618         else static if (T.mant_dig == 113) enum realFormat = RealFormat.ieeeQuadruple;
619         else static assert(false, "Unsupported format for " ~ T.stringof);
620 
621         version (LittleEndian)
622         {
623             enum MANTISSA_LSB = 0;
624             enum MANTISSA_MSB = 1;
625         }
626         else
627         {
628             enum MANTISSA_LSB = 1;
629             enum MANTISSA_MSB = 0;
630         }
631 
632         static if (realFormat == RealFormat.ieeeExtended)
633         {
634             alias S = int;
635             alias U = ushort;
636             enum sig_mask = U(1) << (U.sizeof * 8 - 1);
637             enum exp_shft = 0;
638             enum man_mask = 0;
639             version (LittleEndian)
640                 enum idx = 4;
641             else
642                 enum idx = 0;
643         }
644         else
645         {
646             static if (realFormat == RealFormat.ieeeQuadruple || realFormat == RealFormat.ieeeDouble && double.sizeof == size_t.sizeof)
647             {
648                 alias S = long;
649                 alias U = ulong;
650             }
651             else
652             {
653                 alias S = int;
654                 alias U = uint;
655             }
656             static if (realFormat == RealFormat.ieeeQuadruple)
657                 alias M = ulong;
658             else
659                 alias M = U;
660             enum sig_mask = U(1) << (U.sizeof * 8 - 1);
661             enum uint exp_shft = T.mant_dig - 1 - (T.sizeof > U.sizeof ? U.sizeof * 8 : 0);
662             enum man_mask = (U(1) << exp_shft) - 1;
663             enum idx = T.sizeof > U.sizeof ? MANTISSA_MSB : 0;
664         }
665         enum exp_mask = (U.max >> (exp_shft + 1)) << exp_shft;
666         enum int exp_msh = exp_mask >> exp_shft;
667         enum intPartMask = man_mask + 1;
668 
669         import core.checkedint : adds;
670         alias _expect = llvm_expect;
671 
672         enum norm_factor = 1 / T.epsilon;
673         T vf = n;
674 
675         auto u = (cast(U*)&vf)[idx];
676         int e = (u & exp_mask) >> exp_shft;
677         if (_expect(e != exp_msh, true))
678         {
679             if (_expect(e == 0, false)) // subnormals input
680             {
681                 bool overflow;
682                 vf *= norm_factor;
683                 u = (cast(U*)&vf)[idx];
684                 e = int((u & exp_mask) >> exp_shft) - (T.mant_dig - 1);
685             }
686             bool overflow;
687             exp = adds(exp, e, overflow);
688             if (_expect(overflow || exp >= exp_msh, false)) // infs
689             {
690                 static if (realFormat == RealFormat.ieeeExtended)
691                 {
692                     return vf * T.infinity;
693                 }
694                 else
695                 {
696                     u &= sig_mask;
697                     u ^= exp_mask;
698                     static if (realFormat == RealFormat.ieeeExtended)
699                     {
700                         version (LittleEndian)
701                             auto mp = cast(ulong*)&vf;
702                         else
703                             auto mp = cast(ulong*)((cast(ushort*)&vf) + 1);
704                         *mp = 0;
705                     }
706                     else
707                     static if (T.sizeof > U.sizeof)
708                     {
709                         (cast(U*)&vf)[MANTISSA_LSB] = 0;
710                     }
711                 }
712             }
713             else
714             if (_expect(exp > 0, true)) // normal
715             {
716                 u = cast(U)((u & ~exp_mask) ^ (cast(typeof(U.init + 0))exp << exp_shft));
717             }
718             else // subnormal output
719             {
720                 exp = 1 - exp;
721                 static if (realFormat != RealFormat.ieeeExtended)
722                 {
723                     auto m = u & man_mask;
724                     if (exp > T.mant_dig)
725                     {
726                         exp = T.mant_dig;
727                         static if (T.sizeof > U.sizeof)
728                             (cast(U*)&vf)[MANTISSA_LSB] = 0;
729                     }
730                 }
731                 u &= sig_mask;
732                 static if (realFormat == RealFormat.ieeeExtended)
733                 {
734                     version (LittleEndian)
735                         auto mp = cast(ulong*)&vf;
736                     else
737                         auto mp = cast(ulong*)((cast(ushort*)&vf) + 1);
738                     if (exp >= ulong.sizeof * 8)
739                         *mp = 0;
740                     else
741                         *mp >>>= exp;
742                 }
743                 else
744                 {
745                     m ^= intPartMask;
746                     static if (T.sizeof > U.sizeof)
747                     {
748                         int exp2 = exp - int(U.sizeof) * 8;
749                         if (exp2 < 0)
750                         {
751                             (cast(U*)&vf)[MANTISSA_LSB] = ((cast(U*)&vf)[MANTISSA_LSB] >> exp) ^ (m << (U.sizeof * 8 - exp));
752                             m >>>= exp;
753                             u ^= cast(U) m;
754                         }
755                         else
756                         {
757                             exp = exp2;
758                             (cast(U*)&vf)[MANTISSA_LSB] = (exp < U.sizeof * 8) ? m >> exp : 0;
759                         }
760                     }
761                     else
762                     {
763                         m >>>= exp;
764                         u ^= cast(U) m;
765                     }
766                 }
767             }
768             (cast(U*)&vf)[idx] = u;
769         }
770         return vf;
771     } else version(GNU) {
772         static if (is(T == float))
773             return __builtin_ldexpf(n, exp);
774         else static if (is(T == double))
775             return __builtin_ldexp(n, exp);
776         else static if (is(T == real))
777             return __builtin_ldexpl(n, exp);
778     } else {
779         return cast(T)cmath.ldexp(cast(double)n, exp);
780     }
781 }
782 
783 
784 
785 //
786 //          EXTRAS
787 //
788 private:
789 
790 version(LDC) {
791 
792     // Trigonometry
793     pragma(LDC_intrinsic, "llvm.tan.f#")
794     T llvm_tan(T)(T val) if (__traits(isFloating, T));
795     
796     pragma(LDC_intrinsic, "llvm.asin.f#")
797     T llvm_asin(T)(T val) if (__traits(isFloating, T));
798     
799     pragma(LDC_intrinsic, "llvm.acos.f#")
800     T llvm_acos(T)(T val) if (__traits(isFloating, T));
801     
802     pragma(LDC_intrinsic, "llvm.atan.f#")
803     T llvm_atan(T)(T val) if (__traits(isFloating, T));
804 
805     pragma(LDC_intrinsic, "llvm.atan2.f#")
806     T llvm_atan2(T)(T y, T x) if (__traits(isFloating, T));
807     
808     
809     pragma(LDC_intrinsic, "llvm.sinh.f#")
810     T llvm_sinh(T)(T val) if (__traits(isFloating, T));
811     
812     pragma(LDC_intrinsic, "llvm.cosh.f#")
813     T llvm_cosh(T)(T val) if (__traits(isFloating, T));
814     
815     pragma(LDC_intrinsic, "llvm.tanh.f#")
816     T llvm_tanh(T)(T val) if (__traits(isFloating, T));
817     
818     
819 
820     pragma(LDC_intrinsic, "llvm.asinh.f#")
821     T llvm_asinh(T)(T val) if (__traits(isFloating, T));
822     
823     pragma(LDC_intrinsic, "llvm.acosh.f#")
824     T llvm_acosh(T)(T val) if (__traits(isFloating, T));
825     
826     pragma(LDC_intrinsic, "llvm.atanh.f#")
827     T llvm_atanh(T)(T val) if (__traits(isFloating, T));
828 }