simba/scalar/
complex.rs

1use num::{One, Signed, Zero};
2use num_traits::FromPrimitive;
3use std::any::Any;
4use std::fmt::{Debug, Display};
5use std::ops::Neg;
6use std::{f32, f64};
7
8use crate::scalar::{Field, RealField, SubsetOf, SupersetOf};
9#[cfg(all(not(feature = "std"), not(feature = "libm_force"), feature = "libm"))]
10use num::Float;
11//#[cfg(feature = "decimal")]
12//use decimal::d128;
13
14macro_rules! complex_trait_methods (
15    ($RealField: ident $(, $prefix: ident)*) => {
16        paste::item! {
17            /// Builds a pure-real complex number from the given value.
18            fn [<from_ $($prefix)* real>](re: Self::$RealField) -> Self;
19
20            /// The real part of this complex number.
21            fn [<$($prefix)* real>](self) -> Self::$RealField;
22
23            /// The imaginary part of this complex number.
24            fn [<$($prefix)* imaginary>](self) -> Self::$RealField;
25
26            /// The modulus of this complex number.
27            fn [<$($prefix)* modulus>](self) -> Self::$RealField;
28
29            /// The squared modulus of this complex number.
30            fn [<$($prefix)* modulus_squared>](self) -> Self::$RealField;
31
32            /// The argument of this complex number.
33            fn [<$($prefix)* argument>](self) -> Self::$RealField;
34
35            /// The sum of the absolute value of this complex number's real and imaginary part.
36            fn [<$($prefix)* norm1>](self) -> Self::$RealField;
37
38            /// Multiplies this complex number by `factor`.
39            fn [<$($prefix)* scale>](self, factor: Self::$RealField) -> Self;
40
41            /// Divides this complex number by `factor`.
42            fn [<$($prefix)* unscale>](self, factor: Self::$RealField) -> Self;
43
44            /// The polar form of this complex number: (modulus, arg)
45            fn [<$($prefix)* to_polar>](self) -> (Self::$RealField, Self::$RealField) {
46                (self.clone().[<$($prefix)* modulus>](), self.[<$($prefix)* argument>]())
47            }
48
49            /// The exponential form of this complex number: (modulus, e^{i arg})
50            fn [<$($prefix)* to_exp>](self) -> (Self::$RealField, Self) {
51                let m = self.clone().[<$($prefix)* modulus>]();
52
53                if !m.is_zero() {
54                    (m.clone(), self.[<$($prefix)* unscale>](m))
55                } else {
56                    (Self::$RealField::zero(), Self::one())
57                }
58            }
59
60            /// The exponential part of this complex number: `self / self.modulus()`
61            fn [<$($prefix)* signum>](self) -> Self {
62                self.[<$($prefix)* to_exp>]().1
63            }
64
65            fn [<$($prefix)* floor>](self) -> Self;
66            fn [<$($prefix)* ceil>](self) -> Self;
67            fn [<$($prefix)* round>](self) -> Self;
68            fn [<$($prefix)* trunc>](self) -> Self;
69            fn [<$($prefix)* fract>](self) -> Self;
70            fn [<$($prefix)* mul_add>](self, a: Self, b: Self) -> Self;
71
72            /// The absolute value of this complex number: `self / self.signum()`.
73            ///
74            /// This is equivalent to `self.modulus()`.
75            fn [<$($prefix)* abs>](self) -> Self::$RealField;
76
77            /// Computes (self.conjugate() * self + other.conjugate() * other).sqrt()
78            fn [<$($prefix)* hypot>](self, other: Self) -> Self::$RealField;
79
80            fn [<$($prefix)* recip>](self) -> Self;
81            fn [<$($prefix)* conjugate>](self) -> Self;
82            fn [<$($prefix)* sin>](self) -> Self;
83            fn [<$($prefix)* cos>](self) -> Self;
84            fn [<$($prefix)* sin_cos>](self) -> (Self, Self);
85            #[inline]
86            fn [<$($prefix)* sinh_cosh>](self) -> (Self, Self) {
87                (self.clone().[<$($prefix)* sinh>](), self.[<$($prefix)* cosh>]())
88            }
89            fn [<$($prefix)* tan>](self) -> Self;
90            fn [<$($prefix)* asin>](self) -> Self;
91            fn [<$($prefix)* acos>](self) -> Self;
92            fn [<$($prefix)* atan>](self) -> Self;
93            fn [<$($prefix)* sinh>](self) -> Self;
94            fn [<$($prefix)* cosh>](self) -> Self;
95            fn [<$($prefix)* tanh>](self) -> Self;
96            fn [<$($prefix)* asinh>](self) -> Self;
97            fn [<$($prefix)* acosh>](self) -> Self;
98            fn [<$($prefix)* atanh>](self) -> Self;
99
100            /// Cardinal sine
101            #[inline]
102            fn [<$($prefix)* sinc>](self) -> Self {
103                if self.is_zero() {
104                    Self::one()
105                } else {
106                    self.clone().[<$($prefix)* sin>]() / self
107                }
108            }
109
110            #[inline]
111            fn [<$($prefix)* sinhc>](self) -> Self {
112                if self.is_zero() {
113                    Self::one()
114                } else {
115                    self.clone().[<$($prefix)* sinh>]() / self
116                }
117            }
118
119            /// Cardinal cos
120            #[inline]
121            fn [<$($prefix)* cosc>](self) -> Self {
122                if self.is_zero() {
123                    Self::one()
124                } else {
125                    self.clone().[<$($prefix)* cos>]() / self
126                }
127            }
128
129            #[inline]
130            fn [<$($prefix)* coshc>](self) -> Self {
131                if self.is_zero() {
132                    Self::one()
133                } else {
134                    self.clone().[<$($prefix)* cosh>]() / self
135                }
136            }
137
138            fn [<$($prefix)* log>](self, base: Self::$RealField) -> Self;
139            fn [<$($prefix)* log2>](self) -> Self;
140            fn [<$($prefix)* log10>](self) -> Self;
141            fn [<$($prefix)* ln>](self) -> Self;
142            fn [<$($prefix)* ln_1p>](self) -> Self;
143            fn [<$($prefix)* sqrt>](self) -> Self;
144            fn [<$($prefix)* exp>](self) -> Self;
145            fn [<$($prefix)* exp2>](self) -> Self;
146            fn [<$($prefix)* exp_m1>](self) -> Self;
147            fn [<$($prefix)* powi>](self, n: i32) -> Self;
148            fn [<$($prefix)* powf>](self, n: Self::$RealField) -> Self;
149            fn [<$($prefix)* powc>](self, n: Self) -> Self;
150            fn [<$($prefix)* cbrt>](self) -> Self;
151        }
152    }
153);
154
155/// Trait shared by all complex fields and its subfields (like real numbers).
156///
157/// Complex numbers are equipped with functions that are commonly used on complex numbers and reals.
158/// The results of those functions only have to be approximately equal to the actual theoretical values.
159// TODO: SubsetOf should be removed when specialization will be supported by rustc. This will
160// allow a blanket impl: impl<T: Clone> SubsetOf<T> for T { ... }
161#[allow(missing_docs)]
162pub trait ComplexField:
163SubsetOf<Self>
164+ SupersetOf<f32>
165+ SupersetOf<f64>
166+ FromPrimitive
167+ Field<Element=Self, SimdBool=bool>
168+ Neg<Output=Self>
169+ Clone
170//    + MeetSemilattice
171//    + JoinSemilattice
172+ Send
173+ Sync
174+ Any
175+ 'static
176+ Debug
177+ Display
178{
179    type RealField: RealField;
180    complex_trait_methods!(RealField);
181
182    fn is_finite(&self) -> bool;
183    fn try_sqrt(self) -> Option<Self>;
184}
185
186#[cfg(not(feature = "libm_force"))]
187macro_rules! impl_complex (
188    ($($T:ty, $M:ident, $libm: ident);*) => ($(
189        impl ComplexField for $T {
190            type RealField = $T;
191
192            #[inline]
193            fn from_real(re: Self::RealField) -> Self {
194                re
195            }
196
197            #[inline]
198            fn real(self) -> Self::RealField {
199                self
200            }
201
202            #[inline]
203            fn imaginary(self) -> Self::RealField {
204                Self::zero()
205            }
206
207            #[inline]
208            fn norm1(self) -> Self::RealField {
209                $libm::abs(self)
210            }
211
212            #[inline]
213            fn modulus(self) -> Self::RealField {
214                $libm::abs(self)
215            }
216
217            #[inline]
218            fn modulus_squared(self) -> Self::RealField {
219                self * self
220            }
221
222            #[inline]
223            fn argument(self) -> Self::RealField {
224                if self >= Self::zero() {
225                    Self::zero()
226                } else {
227                    Self::pi()
228                }
229            }
230
231            #[inline]
232            fn to_exp(self) -> (Self, Self) {
233                if self >= Self::zero() {
234                    (self, Self::one())
235                } else {
236                    (-self, -Self::one())
237                }
238            }
239
240            #[inline]
241            fn recip(self) -> Self {
242                $M::recip(self)
243            }
244
245            #[inline]
246            fn conjugate(self) -> Self {
247                self
248            }
249
250            #[inline]
251            fn scale(self, factor: Self::RealField) -> Self {
252                self * factor
253            }
254
255            #[inline]
256            fn unscale(self, factor: Self::RealField) -> Self {
257                self / factor
258            }
259
260            #[inline]
261            fn floor(self) -> Self {
262                $libm::floor(self)
263            }
264
265            #[inline]
266            fn ceil(self) -> Self {
267                $libm::ceil(self)
268            }
269
270            #[inline]
271            fn round(self) -> Self {
272                $libm::round(self)
273            }
274
275            #[inline]
276            fn trunc(self) -> Self {
277                $libm::trunc(self)
278            }
279
280            #[inline]
281            fn fract(self) -> Self {
282                $libm::fract(self)
283            }
284
285            #[inline]
286            fn abs(self) -> Self {
287                $libm::abs(self)
288            }
289
290            #[inline]
291            fn signum(self) -> Self {
292                Signed::signum(&self)
293            }
294
295            #[inline]
296            fn mul_add(self, a: Self, b: Self) -> Self {
297                $libm::mul_add(self, a, b)
298            }
299
300            #[cfg(feature = "std")]
301            #[inline]
302            fn powi(self, n: i32) -> Self {
303                self.powi(n)
304            }
305
306            #[cfg(not(feature = "std"))]
307            #[inline]
308            fn powi(self, n: i32) -> Self {
309                // TODO: is there a more accurate solution?
310                $libm::powf(self, n as $T)
311            }
312
313            #[inline]
314            fn powf(self, n: Self) -> Self {
315                $libm::powf(self, n)
316            }
317
318            #[inline]
319            fn powc(self, n: Self) -> Self {
320                // Same as powf.
321                $libm::powf(self, n)
322            }
323
324            #[inline]
325            fn sqrt(self) -> Self {
326                $libm::sqrt(self)
327            }
328
329            #[inline]
330            fn try_sqrt(self) -> Option<Self> {
331                if self >= Self::zero() {
332                    Some($libm::sqrt(self))
333                } else {
334                    None
335                }
336            }
337
338            #[inline]
339            fn exp(self) -> Self {
340                $libm::exp(self)
341            }
342
343            #[inline]
344            fn exp2(self) -> Self {
345                $libm::exp2(self)
346            }
347
348
349            #[inline]
350            fn exp_m1(self) -> Self {
351                $libm::exp_m1(self)
352            }
353
354            #[inline]
355            fn ln_1p(self) -> Self {
356                $libm::ln_1p(self)
357            }
358
359            #[inline]
360            fn ln(self) -> Self {
361                $libm::ln(self)
362            }
363
364            #[inline]
365            fn log(self, base: Self) -> Self {
366                $libm::log(self, base)
367            }
368
369            #[inline]
370            fn log2(self) -> Self {
371                $libm::log2(self)
372            }
373
374            #[inline]
375            fn log10(self) -> Self {
376                $libm::log10(self)
377            }
378
379            #[inline]
380            fn cbrt(self) -> Self {
381                $libm::cbrt(self)
382            }
383
384            #[inline]
385            fn hypot(self, other: Self) -> Self::RealField {
386                $libm::hypot(self, other)
387            }
388
389            #[inline]
390            fn sin(self) -> Self {
391                $libm::sin(self)
392            }
393
394            #[inline]
395            fn cos(self) -> Self {
396                $libm::cos(self)
397            }
398
399            #[inline]
400            fn tan(self) -> Self {
401                $libm::tan(self)
402            }
403
404            #[inline]
405            fn asin(self) -> Self {
406                $libm::asin(self)
407            }
408
409            #[inline]
410            fn acos(self) -> Self {
411                $libm::acos(self)
412            }
413
414            #[inline]
415            fn atan(self) -> Self {
416                $libm::atan(self)
417            }
418
419            #[inline]
420            fn sin_cos(self) -> (Self, Self) {
421                $libm::sin_cos(self)
422            }
423
424//            #[inline]
425//            fn exp_m1(self) -> Self {
426//                $libm::exp_m1(self)
427//            }
428//
429//            #[inline]
430//            fn ln_1p(self) -> Self {
431//                $libm::ln_1p(self)
432//            }
433//
434            #[inline]
435            fn sinh(self) -> Self {
436                $libm::sinh(self)
437            }
438
439            #[inline]
440            fn cosh(self) -> Self {
441                $libm::cosh(self)
442            }
443
444            #[inline]
445            fn tanh(self) -> Self {
446                $libm::tanh(self)
447            }
448
449            #[inline]
450            fn asinh(self) -> Self {
451                $libm::asinh(self)
452            }
453
454            #[inline]
455            fn acosh(self) -> Self {
456                $libm::acosh(self)
457            }
458
459            #[inline]
460            fn atanh(self) -> Self {
461                $libm::atanh(self)
462            }
463
464            #[inline]
465            fn is_finite(&self) -> bool {
466                $M::is_finite(*self)
467            }
468        }
469    )*)
470);
471
472#[cfg(all(not(feature = "std"), not(feature = "libm_force"), feature = "libm"))]
473impl_complex!(
474    f32, f32, Float;
475    f64, f64, Float
476);
477
478#[cfg(all(feature = "std", not(feature = "libm_force")))]
479impl_complex!(
480    f32,f32,f32;
481    f64,f64,f64
482);
483
484#[cfg(feature = "libm_force")]
485impl ComplexField for f32 {
486    type RealField = f32;
487
488    #[inline]
489    fn from_real(re: Self::RealField) -> Self {
490        re
491    }
492
493    #[inline]
494    fn real(self) -> Self::RealField {
495        self
496    }
497
498    #[inline]
499    fn imaginary(self) -> Self::RealField {
500        Self::zero()
501    }
502
503    #[inline]
504    fn norm1(self) -> Self::RealField {
505        libm_force::fabsf(self)
506    }
507
508    #[inline]
509    fn modulus(self) -> Self::RealField {
510        libm_force::fabsf(self)
511    }
512
513    #[inline]
514    fn modulus_squared(self) -> Self::RealField {
515        self * self
516    }
517
518    #[inline]
519    fn argument(self) -> Self::RealField {
520        if self >= Self::zero() {
521            Self::zero()
522        } else {
523            Self::pi()
524        }
525    }
526
527    #[inline]
528    fn to_exp(self) -> (Self, Self) {
529        if self >= Self::zero() {
530            (self, Self::one())
531        } else {
532            (-self, -Self::one())
533        }
534    }
535
536    #[inline]
537    fn recip(self) -> Self {
538        f32::recip(self)
539    }
540
541    #[inline]
542    fn conjugate(self) -> Self {
543        self
544    }
545
546    #[inline]
547    fn scale(self, factor: Self::RealField) -> Self {
548        self * factor
549    }
550
551    #[inline]
552    fn unscale(self, factor: Self::RealField) -> Self {
553        self / factor
554    }
555
556    #[inline]
557    fn floor(self) -> Self {
558        libm_force::floorf(self)
559    }
560
561    #[inline]
562    fn ceil(self) -> Self {
563        libm_force::ceilf(self)
564    }
565
566    #[inline]
567    fn round(self) -> Self {
568        libm_force::roundf(self)
569    }
570
571    #[inline]
572    fn trunc(self) -> Self {
573        libm_force::truncf(self)
574    }
575
576    #[inline]
577    fn fract(self) -> Self {
578        self - libm_force::truncf(self)
579    }
580
581    #[inline]
582    fn abs(self) -> Self {
583        libm_force::fabsf(self)
584    }
585
586    #[inline]
587    fn signum(self) -> Self {
588        Signed::signum(&self)
589    }
590
591    #[inline]
592    fn mul_add(self, a: Self, b: Self) -> Self {
593        libm_force::fmaf(self, a, b)
594    }
595
596    #[inline]
597    fn powi(self, n: i32) -> Self {
598        // TODO: implement a more accurate/efficient solution?
599        libm_force::powf(self, n as f32)
600    }
601
602    #[inline]
603    fn powf(self, n: Self) -> Self {
604        libm_force::powf(self, n)
605    }
606
607    #[inline]
608    fn powc(self, n: Self) -> Self {
609        // Same as powf.
610        libm_force::powf(self, n)
611    }
612
613    #[inline]
614    fn sqrt(self) -> Self {
615        libm_force::sqrtf(self)
616    }
617
618    #[inline]
619    fn try_sqrt(self) -> Option<Self> {
620        if self >= Self::zero() {
621            Some(libm_force::sqrtf(self))
622        } else {
623            None
624        }
625    }
626
627    #[inline]
628    fn exp(self) -> Self {
629        libm_force::expf(self)
630    }
631
632    #[inline]
633    fn exp2(self) -> Self {
634        libm_force::exp2f(self)
635    }
636
637    #[inline]
638    fn exp_m1(self) -> Self {
639        libm_force::expm1f(self)
640    }
641
642    #[inline]
643    fn ln_1p(self) -> Self {
644        libm_force::log1pf(self)
645    }
646
647    #[inline]
648    fn ln(self) -> Self {
649        libm_force::logf(self)
650    }
651
652    #[inline]
653    fn log(self, base: Self) -> Self {
654        libm_force::logf(self) / libm_force::logf(base)
655    }
656
657    #[inline]
658    fn log2(self) -> Self {
659        libm_force::log2f(self)
660    }
661
662    #[inline]
663    fn log10(self) -> Self {
664        libm_force::log10f(self)
665    }
666
667    #[inline]
668    fn cbrt(self) -> Self {
669        libm_force::cbrtf(self)
670    }
671
672    #[inline]
673    fn hypot(self, other: Self) -> Self::RealField {
674        libm_force::hypotf(self, other)
675    }
676
677    #[inline]
678    fn sin(self) -> Self {
679        libm_force::sinf(self)
680    }
681
682    #[inline]
683    fn cos(self) -> Self {
684        libm_force::cosf(self)
685    }
686
687    #[inline]
688    fn tan(self) -> Self {
689        libm_force::tanf(self)
690    }
691
692    #[inline]
693    fn asin(self) -> Self {
694        libm_force::asinf(self)
695    }
696
697    #[inline]
698    fn acos(self) -> Self {
699        libm_force::acosf(self)
700    }
701
702    #[inline]
703    fn atan(self) -> Self {
704        libm_force::atanf(self)
705    }
706
707    #[inline]
708    fn sin_cos(self) -> (Self, Self) {
709        libm_force::sincosf(self)
710    }
711
712    //            #[inline]
713    //            fn exp_m1(self) -> Self {
714    //                libm_force::exp_m1(self)
715    //            }
716    //
717    //            #[inline]
718    //            fn ln_1p(self) -> Self {
719    //                libm_force::ln_1p(self)
720    //            }
721    //
722    #[inline]
723    fn sinh(self) -> Self {
724        libm_force::sinhf(self)
725    }
726
727    #[inline]
728    fn cosh(self) -> Self {
729        libm_force::coshf(self)
730    }
731
732    #[inline]
733    fn tanh(self) -> Self {
734        libm_force::tanhf(self)
735    }
736
737    #[inline]
738    fn asinh(self) -> Self {
739        libm_force::asinhf(self)
740    }
741
742    #[inline]
743    fn acosh(self) -> Self {
744        libm_force::acoshf(self)
745    }
746
747    #[inline]
748    fn atanh(self) -> Self {
749        libm_force::atanhf(self)
750    }
751
752    #[inline]
753    fn is_finite(&self) -> bool {
754        f32::is_finite(*self)
755    }
756}
757
758#[cfg(feature = "libm_force")]
759impl ComplexField for f64 {
760    type RealField = f64;
761
762    #[inline]
763    fn from_real(re: Self::RealField) -> Self {
764        re
765    }
766
767    #[inline]
768    fn real(self) -> Self::RealField {
769        self
770    }
771
772    #[inline]
773    fn imaginary(self) -> Self::RealField {
774        Self::zero()
775    }
776
777    #[inline]
778    fn norm1(self) -> Self::RealField {
779        libm_force::fabs(self)
780    }
781
782    #[inline]
783    fn modulus(self) -> Self::RealField {
784        libm_force::fabs(self)
785    }
786
787    #[inline]
788    fn modulus_squared(self) -> Self::RealField {
789        self * self
790    }
791
792    #[inline]
793    fn argument(self) -> Self::RealField {
794        if self >= Self::zero() {
795            Self::zero()
796        } else {
797            Self::pi()
798        }
799    }
800
801    #[inline]
802    fn to_exp(self) -> (Self, Self) {
803        if self >= Self::zero() {
804            (self, Self::one())
805        } else {
806            (-self, -Self::one())
807        }
808    }
809
810    #[inline]
811    fn recip(self) -> Self {
812        f64::recip(self)
813    }
814
815    #[inline]
816    fn conjugate(self) -> Self {
817        self
818    }
819
820    #[inline]
821    fn scale(self, factor: Self::RealField) -> Self {
822        self * factor
823    }
824
825    #[inline]
826    fn unscale(self, factor: Self::RealField) -> Self {
827        self / factor
828    }
829
830    #[inline]
831    fn floor(self) -> Self {
832        libm_force::floor(self)
833    }
834
835    #[inline]
836    fn ceil(self) -> Self {
837        libm_force::ceil(self)
838    }
839
840    #[inline]
841    fn round(self) -> Self {
842        libm_force::round(self)
843    }
844
845    #[inline]
846    fn trunc(self) -> Self {
847        libm_force::trunc(self)
848    }
849
850    #[inline]
851    fn fract(self) -> Self {
852        self - libm_force::trunc(self)
853    }
854
855    #[inline]
856    fn abs(self) -> Self {
857        libm_force::fabs(self)
858    }
859
860    #[inline]
861    fn signum(self) -> Self {
862        Signed::signum(&self)
863    }
864
865    #[inline]
866    fn mul_add(self, a: Self, b: Self) -> Self {
867        libm_force::fma(self, a, b)
868    }
869
870    #[inline]
871    fn powi(self, n: i32) -> Self {
872        // TODO: implement a more accurate solution?
873        libm_force::pow(self, n as f64)
874    }
875
876    #[inline]
877    fn powf(self, n: Self) -> Self {
878        libm_force::pow(self, n)
879    }
880
881    #[inline]
882    fn powc(self, n: Self) -> Self {
883        // Same as powf.
884        libm_force::pow(self, n)
885    }
886
887    #[inline]
888    fn sqrt(self) -> Self {
889        libm_force::sqrt(self)
890    }
891
892    #[inline]
893    fn try_sqrt(self) -> Option<Self> {
894        if self >= Self::zero() {
895            Some(libm_force::sqrt(self))
896        } else {
897            None
898        }
899    }
900
901    #[inline]
902    fn exp(self) -> Self {
903        libm_force::exp(self)
904    }
905
906    #[inline]
907    fn exp2(self) -> Self {
908        libm_force::exp2(self)
909    }
910
911    #[inline]
912    fn exp_m1(self) -> Self {
913        libm_force::expm1(self)
914    }
915
916    #[inline]
917    fn ln_1p(self) -> Self {
918        libm_force::log1p(self)
919    }
920
921    #[inline]
922    fn ln(self) -> Self {
923        libm_force::log(self)
924    }
925
926    #[inline]
927    fn log(self, base: Self) -> Self {
928        libm_force::log(self) / libm_force::log(base)
929    }
930
931    #[inline]
932    fn log2(self) -> Self {
933        libm_force::log2(self)
934    }
935
936    #[inline]
937    fn log10(self) -> Self {
938        libm_force::log10(self)
939    }
940
941    #[inline]
942    fn cbrt(self) -> Self {
943        libm_force::cbrt(self)
944    }
945
946    #[inline]
947    fn hypot(self, other: Self) -> Self::RealField {
948        libm_force::hypot(self, other)
949    }
950
951    #[inline]
952    fn sin(self) -> Self {
953        libm_force::sin(self)
954    }
955
956    #[inline]
957    fn cos(self) -> Self {
958        libm_force::cos(self)
959    }
960
961    #[inline]
962    fn tan(self) -> Self {
963        libm_force::tan(self)
964    }
965
966    #[inline]
967    fn asin(self) -> Self {
968        libm_force::asin(self)
969    }
970
971    #[inline]
972    fn acos(self) -> Self {
973        libm_force::acos(self)
974    }
975
976    #[inline]
977    fn atan(self) -> Self {
978        libm_force::atan(self)
979    }
980
981    #[inline]
982    fn sin_cos(self) -> (Self, Self) {
983        libm_force::sincos(self)
984    }
985
986    //            #[inline]
987    //            fn exp_m1(self) -> Self {
988    //                libm_force::exp_m1(self)
989    //            }
990    //
991    //            #[inline]
992    //            fn ln_1p(self) -> Self {
993    //                libm_force::ln_1p(self)
994    //            }
995    //
996    #[inline]
997    fn sinh(self) -> Self {
998        libm_force::sinh(self)
999    }
1000
1001    #[inline]
1002    fn cosh(self) -> Self {
1003        libm_force::cosh(self)
1004    }
1005
1006    #[inline]
1007    fn tanh(self) -> Self {
1008        libm_force::tanh(self)
1009    }
1010
1011    #[inline]
1012    fn asinh(self) -> Self {
1013        libm_force::asinh(self)
1014    }
1015
1016    #[inline]
1017    fn acosh(self) -> Self {
1018        libm_force::acosh(self)
1019    }
1020
1021    #[inline]
1022    fn atanh(self) -> Self {
1023        libm_force::atanh(self)
1024    }
1025
1026    #[inline]
1027    fn is_finite(&self) -> bool {
1028        f64::is_finite(*self)
1029    }
1030}
1031
1032//#[cfg(feature = "decimal")]
1033//impl_real!(d128, d128, d128);
1034
1035// NOTE: all those impls have been copied-pasted to `simd_impl.rs` to implement
1036// SimdComplexField for Complex.
1037impl<N: RealField + PartialOrd> ComplexField for num_complex::Complex<N> {
1038    type RealField = N;
1039
1040    #[inline]
1041    fn from_real(re: Self::RealField) -> Self {
1042        Self::new(re, Self::RealField::zero())
1043    }
1044
1045    #[inline]
1046    fn real(self) -> Self::RealField {
1047        self.re
1048    }
1049
1050    #[inline]
1051    fn imaginary(self) -> Self::RealField {
1052        self.im
1053    }
1054
1055    #[inline]
1056    fn argument(self) -> Self::RealField {
1057        self.im.atan2(self.re)
1058    }
1059
1060    #[inline]
1061    fn modulus(self) -> Self::RealField {
1062        self.re.hypot(self.im)
1063    }
1064
1065    #[inline]
1066    fn modulus_squared(self) -> Self::RealField {
1067        self.re.clone() * self.re + self.im.clone() * self.im
1068    }
1069
1070    #[inline]
1071    fn norm1(self) -> Self::RealField {
1072        self.re.abs() + self.im.abs()
1073    }
1074
1075    #[inline]
1076    fn recip(self) -> Self {
1077        Self::one() / self
1078    }
1079
1080    #[inline]
1081    fn conjugate(self) -> Self {
1082        self.conj()
1083    }
1084
1085    #[inline]
1086    fn scale(self, factor: Self::RealField) -> Self {
1087        self * factor
1088    }
1089
1090    #[inline]
1091    fn unscale(self, factor: Self::RealField) -> Self {
1092        self / factor
1093    }
1094
1095    #[inline]
1096    fn floor(self) -> Self {
1097        Self::new(self.re.floor(), self.im.floor())
1098    }
1099
1100    #[inline]
1101    fn ceil(self) -> Self {
1102        Self::new(self.re.ceil(), self.im.ceil())
1103    }
1104
1105    #[inline]
1106    fn round(self) -> Self {
1107        Self::new(self.re.round(), self.im.round())
1108    }
1109
1110    #[inline]
1111    fn trunc(self) -> Self {
1112        Self::new(self.re.trunc(), self.im.trunc())
1113    }
1114
1115    #[inline]
1116    fn fract(self) -> Self {
1117        Self::new(self.re.fract(), self.im.fract())
1118    }
1119
1120    #[inline]
1121    fn mul_add(self, a: Self, b: Self) -> Self {
1122        self * a + b
1123    }
1124
1125    #[inline]
1126    fn abs(self) -> Self::RealField {
1127        self.modulus()
1128    }
1129
1130    #[inline]
1131    fn exp2(self) -> Self {
1132        let _2 = N::one() + N::one();
1133        num_complex::Complex::new(_2, N::zero()).powc(self)
1134    }
1135
1136    #[inline]
1137    fn exp_m1(self) -> Self {
1138        self.exp() - Self::one()
1139    }
1140
1141    #[inline]
1142    fn ln_1p(self) -> Self {
1143        (Self::one() + self).ln()
1144    }
1145
1146    #[inline]
1147    fn log2(self) -> Self {
1148        let _2 = N::one() + N::one();
1149        self.log(_2)
1150    }
1151
1152    #[inline]
1153    fn log10(self) -> Self {
1154        let _10 = N::from_subset(&10.0f64);
1155        self.log(_10)
1156    }
1157
1158    #[inline]
1159    fn cbrt(self) -> Self {
1160        let one_third = N::from_subset(&(1.0 / 3.0));
1161        self.powf(one_third)
1162    }
1163
1164    #[inline]
1165    fn powi(self, n: i32) -> Self {
1166        // TODO: is there a more accurate solution?
1167        let n = N::from_subset(&(n as f64));
1168        self.powf(n)
1169    }
1170
1171    #[inline]
1172    fn is_finite(&self) -> bool {
1173        self.re.is_finite() && self.im.is_finite()
1174    }
1175
1176    /*
1177     *
1178     *
1179     * Unfortunately we are forced to copy-paste all
1180     * those impls from https://github.com/rust-num/num-complex/blob/master/src/lib.rs
1181     * to avoid requiring `std`.
1182     *
1183     *
1184     */
1185    /// Computes `e^(self)`, where `e` is the base of the natural logarithm.
1186    #[inline]
1187    fn exp(self) -> Self {
1188        // formula: e^(a + bi) = e^a (cos(b) + i*sin(b))
1189        // = from_polar(e^a, b)
1190        complex_from_polar(self.re.exp(), self.im)
1191    }
1192
1193    /// Computes the principal value of natural logarithm of `self`.
1194    ///
1195    /// This function has one branch cut:
1196    ///
1197    /// * `(-∞, 0]`, continuous from above.
1198    ///
1199    /// The branch satisfies `-π ≤ arg(ln(z)) ≤ π`.
1200    #[inline]
1201    fn ln(self) -> Self {
1202        // formula: ln(z) = ln|z| + i*arg(z)
1203        let (r, theta) = self.to_polar();
1204        Self::new(r.ln(), theta)
1205    }
1206
1207    /// Computes the principal value of the square root of `self`.
1208    ///
1209    /// This function has one branch cut:
1210    ///
1211    /// * `(-∞, 0)`, continuous from above.
1212    ///
1213    /// The branch satisfies `-π/2 ≤ arg(sqrt(z)) ≤ π/2`.
1214    #[inline]
1215    fn sqrt(self) -> Self {
1216        // formula: sqrt(r e^(it)) = sqrt(r) e^(it/2)
1217        let two = N::one() + N::one();
1218        let (r, theta) = self.to_polar();
1219        complex_from_polar(r.sqrt(), theta / two)
1220    }
1221
1222    #[inline]
1223    fn try_sqrt(self) -> Option<Self> {
1224        Some(self.sqrt())
1225    }
1226
1227    #[inline]
1228    fn hypot(self, b: Self) -> Self::RealField {
1229        (self.modulus_squared() + b.modulus_squared()).sqrt()
1230    }
1231
1232    /// Raises `self` to a floating point power.
1233    #[inline]
1234    fn powf(self, exp: Self::RealField) -> Self {
1235        // formula: x^y = (ρ e^(i θ))^y = ρ^y e^(i θ y)
1236        // = from_polar(ρ^y, θ y)
1237        let (r, theta) = self.to_polar();
1238        complex_from_polar(r.powf(exp.clone()), theta * exp)
1239    }
1240
1241    /// Returns the logarithm of `self` with respect to an arbitrary base.
1242    #[inline]
1243    fn log(self, base: N) -> Self {
1244        // formula: log_y(x) = log_y(ρ e^(i θ))
1245        // = log_y(ρ) + log_y(e^(i θ)) = log_y(ρ) + ln(e^(i θ)) / ln(y)
1246        // = log_y(ρ) + i θ / ln(y)
1247        let (r, theta) = self.to_polar();
1248        Self::new(r.log(base.clone()), theta / base.ln())
1249    }
1250
1251    /// Raises `self` to a complex power.
1252    #[inline]
1253    fn powc(self, exp: Self) -> Self {
1254        // formula: x^y = (a + i b)^(c + i d)
1255        // = (ρ e^(i θ))^c (ρ e^(i θ))^(i d)
1256        //    where ρ=|x| and θ=arg(x)
1257        // = ρ^c e^(−d θ) e^(i c θ) ρ^(i d)
1258        // = p^c e^(−d θ) (cos(c θ)
1259        //   + i sin(c θ)) (cos(d ln(ρ)) + i sin(d ln(ρ)))
1260        // = p^c e^(−d θ) (
1261        //   cos(c θ) cos(d ln(ρ)) − sin(c θ) sin(d ln(ρ))
1262        //   + i(cos(c θ) sin(d ln(ρ)) + sin(c θ) cos(d ln(ρ))))
1263        // = p^c e^(−d θ) (cos(c θ + d ln(ρ)) + i sin(c θ + d ln(ρ)))
1264        // = from_polar(p^c e^(−d θ), c θ + d ln(ρ))
1265        let (r, theta) = self.to_polar();
1266        complex_from_polar(
1267            r.clone().powf(exp.re.clone()) * (-exp.im.clone() * theta.clone()).exp(),
1268            exp.re * theta + exp.im * r.ln(),
1269        )
1270    }
1271
1272    /*
1273    /// Raises a floating point number to the complex power `self`.
1274    #[inline]
1275    fn expf(&self, base: T) -> Self {
1276        // formula: x^(a+bi) = x^a x^bi = x^a e^(b ln(x) i)
1277        // = from_polar(x^a, b ln(x))
1278        Self::from_polar(&base.powf(self.re), &(self.im * base.ln()))
1279    }
1280    */
1281
1282    /// Computes the sine of `self`.
1283    #[inline]
1284    fn sin(self) -> Self {
1285        // formula: sin(a + bi) = sin(a)cosh(b) + i*cos(a)sinh(b)
1286        Self::new(
1287            self.re.clone().sin() * self.im.clone().cosh(),
1288            self.re.cos() * self.im.sinh(),
1289        )
1290    }
1291
1292    /// Computes the cosine of `self`.
1293    #[inline]
1294    fn cos(self) -> Self {
1295        // formula: cos(a + bi) = cos(a)cosh(b) - i*sin(a)sinh(b)
1296        Self::new(
1297            self.re.clone().cos() * self.im.clone().cosh(),
1298            -self.re.sin() * self.im.sinh(),
1299        )
1300    }
1301
1302    #[inline]
1303    fn sin_cos(self) -> (Self, Self) {
1304        let (rsin, rcos) = self.re.sin_cos();
1305        let (isinh, icosh) = self.im.sinh_cosh();
1306        let sin = Self::new(rsin.clone() * icosh.clone(), rcos.clone() * isinh.clone());
1307        let cos = Self::new(rcos * icosh, -rsin * isinh);
1308
1309        (sin, cos)
1310    }
1311
1312    /// Computes the tangent of `self`.
1313    #[inline]
1314    fn tan(self) -> Self {
1315        // formula: tan(a + bi) = (sin(2a) + i*sinh(2b))/(cos(2a) + cosh(2b))
1316        let (two_re, two_im) = (self.re.clone() + self.re, self.im.clone() + self.im);
1317        Self::new(two_re.clone().sin(), two_im.clone().sinh()).unscale(two_re.cos() + two_im.cosh())
1318    }
1319
1320    /// Computes the principal value of the inverse sine of `self`.
1321    ///
1322    /// This function has two branch cuts:
1323    ///
1324    /// * `(-∞, -1)`, continuous from above.
1325    /// * `(1, ∞)`, continuous from below.
1326    ///
1327    /// The branch satisfies `-π/2 ≤ Re(asin(z)) ≤ π/2`.
1328    #[inline]
1329    fn asin(self) -> Self {
1330        // formula: arcsin(z) = -i ln(sqrt(1-z^2) + iz)
1331        let i = Self::i();
1332        -i.clone() * ((Self::one() - self.clone() * self.clone()).sqrt() + i * self).ln()
1333    }
1334
1335    /// Computes the principal value of the inverse cosine of `self`.
1336    ///
1337    /// This function has two branch cuts:
1338    ///
1339    /// * `(-∞, -1)`, continuous from above.
1340    /// * `(1, ∞)`, continuous from below.
1341    ///
1342    /// The branch satisfies `0 ≤ Re(acos(z)) ≤ π`.
1343    #[inline]
1344    fn acos(self) -> Self {
1345        // formula: arccos(z) = -i ln(i sqrt(1-z^2) + z)
1346        let i = Self::i();
1347        -i.clone() * (i * (Self::one() - self.clone() * self.clone()).sqrt() + self).ln()
1348    }
1349
1350    /// Computes the principal value of the inverse tangent of `self`.
1351    ///
1352    /// This function has two branch cuts:
1353    ///
1354    /// * `(-∞i, -i]`, continuous from the left.
1355    /// * `[i, ∞i)`, continuous from the right.
1356    ///
1357    /// The branch satisfies `-π/2 ≤ Re(atan(z)) ≤ π/2`.
1358    #[inline]
1359    fn atan(self) -> Self {
1360        // formula: arctan(z) = (ln(1+iz) - ln(1-iz))/(2i)
1361        let i = Self::i();
1362        let one = Self::one();
1363        let two = one.clone() + one.clone();
1364
1365        if self == i {
1366            return Self::new(N::zero(), N::one() / N::zero());
1367        } else if self == -i.clone() {
1368            return Self::new(N::zero(), -N::one() / N::zero());
1369        }
1370
1371        ((one.clone() + i.clone() * self.clone()).ln() - (one - i.clone() * self).ln()) / (two * i)
1372    }
1373
1374    /// Computes the hyperbolic sine of `self`.
1375    #[inline]
1376    fn sinh(self) -> Self {
1377        // formula: sinh(a + bi) = sinh(a)cos(b) + i*cosh(a)sin(b)
1378        Self::new(
1379            self.re.clone().sinh() * self.im.clone().cos(),
1380            self.re.cosh() * self.im.sin(),
1381        )
1382    }
1383
1384    /// Computes the hyperbolic cosine of `self`.
1385    #[inline]
1386    fn cosh(self) -> Self {
1387        // formula: cosh(a + bi) = cosh(a)cos(b) + i*sinh(a)sin(b)
1388        Self::new(
1389            self.re.clone().cosh() * self.im.clone().cos(),
1390            self.re.sinh() * self.im.sin(),
1391        )
1392    }
1393
1394    #[inline]
1395    fn sinh_cosh(self) -> (Self, Self) {
1396        let (rsinh, rcosh) = self.re.sinh_cosh();
1397        let (isin, icos) = self.im.sin_cos();
1398        let sin = Self::new(rsinh.clone() * icos.clone(), rcosh.clone() * isin.clone());
1399        let cos = Self::new(rcosh * icos, rsinh * isin);
1400
1401        (sin, cos)
1402    }
1403
1404    /// Computes the hyperbolic tangent of `self`.
1405    #[inline]
1406    fn tanh(self) -> Self {
1407        // formula: tanh(a + bi) = (sinh(2a) + i*sin(2b))/(cosh(2a) + cos(2b))
1408        let (two_re, two_im) = (self.re.clone() + self.re, self.im.clone() + self.im);
1409        Self::new(two_re.clone().sinh(), two_im.clone().sin()).unscale(two_re.cosh() + two_im.cos())
1410    }
1411
1412    /// Computes the principal value of inverse hyperbolic sine of `self`.
1413    ///
1414    /// This function has two branch cuts:
1415    ///
1416    /// * `(-∞i, -i)`, continuous from the left.
1417    /// * `(i, ∞i)`, continuous from the right.
1418    ///
1419    /// The branch satisfies `-π/2 ≤ Im(asinh(z)) ≤ π/2`.
1420    #[inline]
1421    fn asinh(self) -> Self {
1422        // formula: arcsinh(z) = ln(z + sqrt(1+z^2))
1423        let one = Self::one();
1424        (self.clone() + (one + self.clone() * self).sqrt()).ln()
1425    }
1426
1427    /// Computes the principal value of inverse hyperbolic cosine of `self`.
1428    ///
1429    /// This function has one branch cut:
1430    ///
1431    /// * `(-∞, 1)`, continuous from above.
1432    ///
1433    /// The branch satisfies `-π ≤ Im(acosh(z)) ≤ π` and `0 ≤ Re(acosh(z)) < ∞`.
1434    #[inline]
1435    fn acosh(self) -> Self {
1436        // formula: arccosh(z) = 2 ln(sqrt((z+1)/2) + sqrt((z-1)/2))
1437        let one = Self::one();
1438        let two = one.clone() + one.clone();
1439        two.clone()
1440            * (((self.clone() + one.clone()) / two.clone()).sqrt() + ((self - one) / two).sqrt())
1441                .ln()
1442    }
1443
1444    /// Computes the principal value of inverse hyperbolic tangent of `self`.
1445    ///
1446    /// This function has two branch cuts:
1447    ///
1448    /// * `(-∞, -1]`, continuous from above.
1449    /// * `[1, ∞)`, continuous from below.
1450    ///
1451    /// The branch satisfies `-π/2 ≤ Im(atanh(z)) ≤ π/2`.
1452    #[inline]
1453    fn atanh(self) -> Self {
1454        // formula: arctanh(z) = (ln(1+z) - ln(1-z))/2
1455        let one = Self::one();
1456        let two = one.clone() + one.clone();
1457        if self == one {
1458            return Self::new(N::one() / N::zero(), N::zero());
1459        } else if self == -one.clone() {
1460            return Self::new(-N::one() / N::zero(), N::zero());
1461        }
1462        ((one.clone() + self.clone()).ln() - (one - self).ln()) / two
1463    }
1464}
1465
1466#[inline]
1467fn complex_from_polar<N: RealField>(r: N, theta: N) -> num_complex::Complex<N> {
1468    num_complex::Complex::new(r.clone() * theta.clone().cos(), r * theta.sin())
1469}
simba/scalar/complex.rs

simba/scalar/
complex.rs
