simba/scalar/
complex.rs

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

simba/scalar/
complex.rs
