numdiff/automatic_differentiation/
dual.rs

1use num_traits::{Float, Num, NumCast, One, ToPrimitive, Zero};
2use std::cmp::Ordering;
3use std::num::FpCategory;
4use std::ops::{
5    Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign,
6};
7
8#[cfg(feature = "trig")]
9use trig::Trig;
10
11/// First-order dual number.
12#[derive(Debug, Clone, Copy)]
13pub struct Dual {
14    /// Real part of the dual number.
15    real: f64,
16
17    /// Dual part of the dual number.
18    dual: f64,
19}
20
21impl Dual {
22    /// Constructor.
23    ///
24    /// # Arguments
25    ///
26    /// * `real` - Real part.
27    /// * `dual` - Dual part.
28    ///
29    /// # Returns
30    ///
31    /// Dual number.
32    ///
33    /// # Example
34    ///
35    /// ```
36    /// use numdiff::Dual;
37    ///
38    /// let num = Dual::new(1.0, 2.0);
39    /// ```
40    pub fn new(real: f64, dual: f64) -> Dual {
41        Dual { real, dual }
42    }
43
44    /// Get the real part of the dual number.
45    ///
46    /// # Returns
47    ///
48    /// Real part of the dual number.
49    ///
50    /// # Example
51    ///
52    /// ```
53    /// use numdiff::Dual;
54    ///
55    /// let num = Dual::new(1.0, 2.0);
56    /// assert_eq!(num.get_real(), 1.0);
57    /// ```
58    pub fn get_real(self) -> f64 {
59        self.real
60    }
61
62    /// Get the dual part of the dual number.
63    ///
64    /// # Returns
65    ///
66    /// Dual part of the dual number.
67    ///
68    /// # Example
69    ///
70    /// ```
71    /// use numdiff::Dual;
72    ///
73    /// let num = Dual::new(1.0, 2.0);
74    /// assert_eq!(num.get_dual(), 2.0);
75    /// ```
76    pub fn get_dual(self) -> f64 {
77        self.dual
78    }
79}
80
81// --------------------------------
82// Implementing num_traits::NumOps.
83// --------------------------------
84// https://docs.rs/num-traits/latest/num_traits/trait.NumOps.html
85//
86// pub trait NumOps<Rhs = Self, Output = Self>:
87//     Add<Rhs, Output = Output>
88//     + Sub<Rhs, Output = Output>
89//     + Mul<Rhs, Output = Output>
90//     + Div<Rhs, Output = Output>
91//     + Rem<Rhs, Output = Output>
92// {
93// }
94
95// Dual + Dual.
96impl Add for Dual {
97    type Output = Dual;
98    fn add(self, rhs: Dual) -> Dual {
99        Dual::new(self.real + rhs.real, self.dual + rhs.dual)
100    }
101}
102
103// Dual - Dual.
104impl Sub for Dual {
105    type Output = Dual;
106    fn sub(self, rhs: Dual) -> Dual {
107        Dual::new(self.real - rhs.real, self.dual - rhs.dual)
108    }
109}
110
111// Dual * Dual.
112impl Mul for Dual {
113    type Output = Dual;
114    fn mul(self, rhs: Dual) -> Dual {
115        Dual::new(
116            self.real * rhs.real,
117            self.dual * rhs.real + self.real * rhs.dual,
118        )
119    }
120}
121
122// Dual / Dual.
123impl Div for Dual {
124    type Output = Dual;
125    fn div(self, rhs: Dual) -> Dual {
126        Dual::new(
127            self.real / rhs.real,
128            (self.dual * rhs.real - self.real * rhs.dual) / rhs.real.powi(2),
129        )
130    }
131}
132
133// Remainder of Dual / Dual.
134impl Rem for Dual {
135    type Output = Self;
136    fn rem(self, rhs: Self) -> Self::Output {
137        Dual::new(
138            self.real % rhs.real,
139            self.dual - (self.real / rhs.real).trunc() * rhs.dual,
140        )
141    }
142}
143
144// -----------------------------
145// Implementing num_traits::Num.
146// -----------------------------
147// https://docs.rs/num-traits/latest/num_traits/trait.Num.html
148//
149// pub trait Num: PartialEq + Zero + One + NumOps {
150//     type FromStrRadixErr;
151//
152//     // Required method
153//     fn from_str_radix(
154//         str: &str,
155//         radix: u32
156//     ) -> Result<Self, Self::FromStrRadixErr>;
157// }
158
159// https://iel.ucdavis.edu/publication/journal/j_EC1.pdf (p. 11)
160impl PartialEq for Dual {
161    fn eq(&self, other: &Self) -> bool {
162        self.real == other.real && self.dual == other.dual
163    }
164}
165
166impl Zero for Dual {
167    fn zero() -> Self {
168        Dual::new(0.0, 0.0)
169    }
170    fn is_zero(&self) -> bool {
171        self.real.is_zero() && self.dual.is_zero()
172    }
173}
174
175impl One for Dual {
176    fn one() -> Self {
177        Dual::new(1.0, 0.0)
178    }
179}
180
181impl Num for Dual {
182    type FromStrRadixErr = <f64 as Num>::FromStrRadixErr;
183    fn from_str_radix(str: &str, radix: u32) -> Result<Self, Self::FromStrRadixErr> {
184        f64::from_str_radix(str, radix).map(|re| Dual::new(re, 0.0))
185    }
186}
187
188// -------------------------------------
189// Implementing num_traits::ToPrimitive.
190// -------------------------------------
191// https://docs.rs/num-traits/latest/num_traits/cast/trait.ToPrimitive.html
192//
193// pub trait ToPrimitive {
194//     // Required methods
195//     fn to_i64(&self) -> Option<i64>;
196//     fn to_u64(&self) -> Option<u64>;
197//
198//     // Provided methods
199//     ...
200//     fn to_f64(&self) -> Option<f64> { ... }
201// }
202
203impl ToPrimitive for Dual {
204    fn to_i64(&self) -> Option<i64> {
205        self.real.to_i64()
206    }
207    fn to_u64(&self) -> Option<u64> {
208        self.real.to_u64()
209    }
210    fn to_f64(&self) -> Option<f64> {
211        Some(self.real)
212    }
213}
214
215// ---------------------------------
216// Implementing num_traits::NumCast.
217// ---------------------------------
218// https://docs.rs/num-traits/latest/num_traits/cast/trait.NumCast.html
219//
220// pub trait NumCast: Sized + ToPrimitive {
221//     // Required method
222//     fn from<T: ToPrimitive>(n: T) -> Option<Self>;
223// }
224
225impl NumCast for Dual {
226    fn from<T: ToPrimitive>(n: T) -> Option<Self> {
227        n.to_f64().map(|re| Dual::new(re, 0.0))
228    }
229}
230
231// -------------------------------
232// Implementing num_traits::Float.
233// -------------------------------
234// https://docs.rs/num-traits/latest/num_traits/float/trait.Float.html
235//
236// pub trait Float: Num + Copy + NumCast + PartialOrd + Neg<Output = Self> {
237// [+] Show 60 methods
238// }
239
240// Only perform comparisons on the real part.
241//  --> This is primarily to support numerical methods where we want to check convergence on the
242//      actual function evaluation, and NOT its derivative.
243impl PartialOrd for Dual {
244    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
245        self.real.partial_cmp(&other.real)
246    }
247}
248
249impl Neg for Dual {
250    type Output = Self;
251    fn neg(self) -> Self::Output {
252        Dual::new(-self.real, -self.dual)
253    }
254}
255
256impl Float for Dual {
257    fn nan() -> Self {
258        Dual::new(f64::nan(), f64::nan())
259    }
260
261    fn infinity() -> Self {
262        Dual::new(f64::infinity(), f64::infinity())
263    }
264
265    fn neg_infinity() -> Self {
266        Dual::new(f64::neg_infinity(), f64::neg_infinity())
267    }
268
269    fn neg_zero() -> Self {
270        Dual::new(f64::neg_zero(), f64::neg_zero())
271    }
272
273    fn min_value() -> Self {
274        Dual::new(f64::min_value(), f64::min_value())
275    }
276
277    fn min_positive_value() -> Self {
278        Dual::new(f64::min_positive_value(), f64::min_positive_value())
279    }
280
281    fn max_value() -> Self {
282        Dual::new(f64::max_value(), f64::max_value())
283    }
284
285    fn is_nan(self) -> bool {
286        self.real.is_nan()
287    }
288
289    fn is_infinite(self) -> bool {
290        self.real.is_infinite()
291    }
292
293    fn is_finite(self) -> bool {
294        self.real.is_finite()
295    }
296
297    fn is_normal(self) -> bool {
298        self.real.is_normal()
299    }
300
301    fn classify(self) -> FpCategory {
302        self.real.classify()
303    }
304
305    fn floor(self) -> Self {
306        Dual::new(self.real.floor(), 0.0)
307    }
308
309    fn ceil(self) -> Self {
310        Dual::new(self.real.ceil(), 0.0)
311    }
312
313    fn round(self) -> Self {
314        Dual::new(self.real.round(), 0.0)
315    }
316
317    fn trunc(self) -> Self {
318        Dual::new(self.real.trunc(), 0.0)
319    }
320
321    fn fract(self) -> Self {
322        Dual::new(self.real.fract(), self.dual)
323    }
324
325    fn abs(self) -> Self {
326        Dual::new(self.real.abs(), self.dual * self.real.signum())
327    }
328
329    fn signum(self) -> Self {
330        Dual::new(self.real.signum(), f64::zero())
331    }
332
333    fn is_sign_positive(self) -> bool {
334        self.real.is_sign_positive()
335    }
336
337    fn is_sign_negative(self) -> bool {
338        self.real.is_sign_negative()
339    }
340
341    fn mul_add(self, a: Self, b: Self) -> Self {
342        Dual::new(
343            self.real.mul_add(a.real, b.real),
344            self.real * a.dual + self.dual * a.real + b.dual,
345        )
346    }
347
348    fn recip(self) -> Self {
349        Dual::new(self.real.recip(), -self.dual / self.real.powi(2))
350    }
351
352    fn powi(self, n: i32) -> Self {
353        if n == 0 {
354            Dual::one()
355        } else {
356            Dual::new(
357                self.real.powi(n),
358                (n as f64) * self.dual * self.real.powi(n - 1),
359            )
360        }
361    }
362
363    // Numerically-stable version.
364    fn powf(self, n: Self) -> Self {
365        (self.ln() * n).exp()
366    }
367
368    fn sqrt(self) -> Self {
369        let sqrt_re = self.real.sqrt();
370        Dual::new(sqrt_re, self.dual / (2.0 * sqrt_re))
371    }
372
373    fn exp(self) -> Self {
374        let exp_re = self.real.exp();
375        Dual::new(exp_re, exp_re * self.dual)
376    }
377
378    fn exp2(self) -> Self {
379        let exp2_re = self.real.exp2();
380        Dual::new(exp2_re, exp2_re * std::f64::consts::LN_2 * self.dual)
381    }
382
383    fn ln(self) -> Self {
384        Dual::new(self.real.ln(), self.dual / self.real)
385    }
386
387    fn log(self, base: Self) -> Self {
388        Dual::new(
389            self.real.log(base.real),
390            self.dual / (self.real * base.real.ln()),
391        )
392    }
393
394    fn log2(self) -> Self {
395        Dual::new(
396            self.real.ln() / std::f64::consts::LN_2,
397            self.dual / (self.real * std::f64::consts::LN_2),
398        )
399    }
400
401    fn log10(self) -> Self {
402        Dual::new(
403            self.real.ln() / std::f64::consts::LN_10,
404            self.dual / (self.real * std::f64::consts::LN_10),
405        )
406    }
407
408    fn max(self, other: Self) -> Self {
409        if self.real > other.real {
410            self
411        } else {
412            other
413        }
414    }
415
416    fn min(self, other: Self) -> Self {
417        if self.real < other.real {
418            self
419        } else {
420            other
421        }
422    }
423
424    #[allow(deprecated)]
425    fn abs_sub(self, other: Self) -> Self {
426        if self.real > other.real {
427            self - other
428        } else {
429            Self::zero()
430        }
431    }
432
433    fn cbrt(self) -> Self {
434        let cbrt_re = self.real.cbrt();
435        Dual::new(cbrt_re, self.dual / (3.0 * cbrt_re.powi(2)))
436    }
437
438    fn hypot(self, other: Self) -> Self {
439        let hypot_real = (self.real.powi(2) + other.real.powi(2)).sqrt();
440        Dual::new(
441            hypot_real,
442            (self.real * self.dual + other.real * other.dual) / hypot_real,
443        )
444    }
445
446    fn sin(self) -> Self {
447        Dual::new(self.real.sin(), self.real.cos() * self.dual)
448    }
449
450    fn cos(self) -> Self {
451        Dual::new(self.real.cos(), -self.real.sin() * self.dual)
452    }
453
454    fn tan(self) -> Self {
455        let re_tan = self.real.tan();
456        Dual::new(re_tan, self.dual / (self.real.cos().powi(2)))
457    }
458
459    fn asin(self) -> Self {
460        Dual::new(
461            self.real.asin(),
462            self.dual / (1.0 - self.real.powi(2)).sqrt(),
463        )
464    }
465
466    fn acos(self) -> Self {
467        Dual::new(
468            self.real.acos(),
469            -self.dual / (1.0 - self.real.powi(2)).sqrt(),
470        )
471    }
472
473    fn atan(self) -> Self {
474        Dual::new(self.real.atan(), self.dual / (1.0 + self.real.powi(2)))
475    }
476
477    fn atan2(self, other: Self) -> Self {
478        Dual::new(
479            self.real.atan2(other.real),
480            (self.dual * other.real - self.real * other.dual)
481                / (self.real.powi(2) + other.real.powi(2)),
482        )
483    }
484
485    fn sin_cos(self) -> (Self, Self) {
486        (
487            Dual::new(self.real.sin(), self.real.cos() * self.dual),
488            Dual::new(self.real.cos(), -self.real.sin() * self.dual),
489        )
490    }
491
492    fn exp_m1(self) -> Self {
493        let exp_re = self.real.exp();
494        Dual::new(exp_re - 1.0, self.dual * exp_re)
495    }
496
497    fn ln_1p(self) -> Self {
498        Dual::new(self.real.ln_1p(), self.dual / (1.0 + self.real))
499    }
500
501    fn sinh(self) -> Self {
502        Dual::new(self.real.sinh(), self.dual * self.real.cosh())
503    }
504
505    fn cosh(self) -> Self {
506        Dual::new(self.real.cosh(), self.dual * self.real.sinh())
507    }
508
509    fn tanh(self) -> Self {
510        let tanh_re = self.real.tanh();
511        Dual::new(tanh_re, self.dual * (1.0 - tanh_re.powi(2)))
512    }
513
514    fn asinh(self) -> Self {
515        Dual::new(
516            self.real.asinh(),
517            self.dual / (self.real.powi(2) + 1.0).sqrt(),
518        )
519    }
520
521    fn acosh(self) -> Self {
522        Dual::new(
523            self.real.acosh(),
524            self.dual / (self.real.powi(2) - 1.0).sqrt(),
525        )
526    }
527
528    fn atanh(self) -> Self {
529        Dual::new(self.real.atanh(), self.dual / (1.0 - self.real.powi(2)))
530    }
531
532    // This method is really irrelevant, but we need to implement it anyway to satisfy the Float
533    // trait.
534    fn integer_decode(self) -> (u64, i16, i8) {
535        self.real.integer_decode()
536    }
537}
538
539// -----------------------------------
540// Implementing linalg_traits::Scalar.
541// -----------------------------------
542// https://docs.rs/linalg-traits/latest/linalg_traits/trait.Scalar.html
543//
544// pub trait Scalar:
545//     Float
546//     + AddAssign<Self>
547//     + SubAssign<Self>
548//     + MulAssign<Self>
549//     + DivAssign<Self>
550//     + RemAssign<Self>
551//     + Add<f64, Output = Self>
552//     + Sub<f64, Output = Self>
553//     + Mul<f64, Output = Self>
554//     + Div<f64, Output = Self>
555//     + Rem<f64, Output = Self>
556//     + AddAssign<f64>
557//     + SubAssign<f64>
558//     + MulAssign<f64>
559//     + DivAssign<f64>
560//     + RemAssign<f64>
561//     + Debug
562//     + 'static { }
563
564// Dual += Dual.
565impl AddAssign for Dual {
566    fn add_assign(&mut self, other: Dual) {
567        self.real += other.real;
568        self.dual += other.dual;
569    }
570}
571
572// Dual -= Dual.
573impl SubAssign for Dual {
574    fn sub_assign(&mut self, other: Dual) {
575        self.real -= other.real;
576        self.dual -= other.dual;
577    }
578}
579
580// Dual *= Dual.
581impl MulAssign for Dual {
582    fn mul_assign(&mut self, other: Dual) {
583        self.dual = self.dual * other.real + self.real * other.dual;
584        self.real *= other.real;
585    }
586}
587
588// Dual /= Dual.
589impl DivAssign for Dual {
590    fn div_assign(&mut self, other: Dual) {
591        self.dual = (self.dual * other.real - self.real * other.dual) / other.real.powi(2);
592        self.real /= other.real;
593    }
594}
595
596// Dual %= Dual.
597impl RemAssign for Dual {
598    fn rem_assign(&mut self, rhs: Self) {
599        self.dual -= (self.real / rhs.real).floor() * rhs.dual;
600        self.real %= rhs.real;
601    }
602}
603
604// Dual + f64.
605impl Add<f64> for Dual {
606    type Output = Dual;
607    fn add(self, rhs: f64) -> Dual {
608        Dual::new(self.real + rhs, self.dual)
609    }
610}
611
612// Dual - f64.
613impl Sub<f64> for Dual {
614    type Output = Dual;
615    fn sub(self, rhs: f64) -> Dual {
616        Dual::new(self.real - rhs, self.dual)
617    }
618}
619
620// Dual * f64.
621impl Mul<f64> for Dual {
622    type Output = Dual;
623    fn mul(self, rhs: f64) -> Dual {
624        Dual::new(self.real * rhs, self.dual * rhs)
625    }
626}
627
628// Dual / f64.
629impl Div<f64> for Dual {
630    type Output = Dual;
631    fn div(self, rhs: f64) -> Dual {
632        Dual::new(self.real / rhs, self.dual / rhs)
633    }
634}
635
636// Dual % f64.
637impl Rem<f64> for Dual {
638    type Output = Dual;
639    fn rem(self, rhs: f64) -> Self::Output {
640        let rem_real = self.real % rhs;
641        let rem_dual = if self.real % rhs == 0.0 {
642            0.0
643        } else {
644            self.dual
645        };
646        Dual::new(rem_real, rem_dual)
647    }
648}
649
650// Dual += f64.
651impl AddAssign<f64> for Dual {
652    fn add_assign(&mut self, rhs: f64) {
653        self.real += rhs;
654    }
655}
656
657// Dual -= f64.
658impl SubAssign<f64> for Dual {
659    fn sub_assign(&mut self, rhs: f64) {
660        self.real -= rhs;
661    }
662}
663
664// Dual *= f64.
665impl MulAssign<f64> for Dual {
666    fn mul_assign(&mut self, rhs: f64) {
667        self.real *= rhs;
668        self.dual *= rhs;
669    }
670}
671
672// Dual /= f64.
673impl DivAssign<f64> for Dual {
674    fn div_assign(&mut self, rhs: f64) {
675        self.real /= rhs;
676        self.dual /= rhs;
677    }
678}
679
680// Dual %= f64.
681impl RemAssign<f64> for Dual {
682    fn rem_assign(&mut self, rhs: f64) {
683        self.real = self.real % rhs;
684        if self.real == 0.0 {
685            self.dual = 0.0;
686        }
687    }
688}
689
690// ---------------------------
691// Interoperability with f64s.
692// ---------------------------
693
694// f64 + Dual.
695impl Add<Dual> for f64 {
696    type Output = Dual;
697    fn add(self, rhs: Dual) -> Dual {
698        Dual::new(self + rhs.real, rhs.dual)
699    }
700}
701
702// f64 - Dual.
703impl Sub<Dual> for f64 {
704    type Output = Dual;
705    fn sub(self, rhs: Dual) -> Dual {
706        Dual::new(self - rhs.real, rhs.dual)
707    }
708}
709
710// f64 * Dual.
711impl Mul<Dual> for f64 {
712    type Output = Dual;
713    fn mul(self, rhs: Dual) -> Dual {
714        Dual::new(self * rhs.real, self * rhs.dual)
715    }
716}
717
718// f64 / Dual.
719impl Div<Dual> for f64 {
720    type Output = Dual;
721    fn div(self, rhs: Dual) -> Dual {
722        Dual::new(self / rhs.real, -self * rhs.dual / rhs.real.powi(2))
723    }
724}
725
726// f64 % Dual.
727impl Rem<Dual> for f64 {
728    type Output = Dual;
729    fn rem(self, rhs: Dual) -> Dual {
730        Dual::new(self % rhs.real, -(self / rhs.real).floor() * rhs.dual)
731    }
732}
733
734// ------------------------
735// Implementing trig::Trig.
736// ------------------------
737
738#[cfg(feature = "trig")]
739impl Trig for Dual {
740    fn sin(&self) -> Dual {
741        <Dual as Float>::sin(*self)
742    }
743    fn cos(&self) -> Dual {
744        <Dual as Float>::cos(*self)
745    }
746    fn tan(&self) -> Dual {
747        <Dual as Float>::tan(*self)
748    }
749    fn csc(&self) -> Dual {
750        1.0 / self.sin()
751    }
752    fn sec(&self) -> Dual {
753        1.0 / self.cos()
754    }
755    fn cot(&self) -> Dual {
756        1.0 / self.tan()
757    }
758    fn asin(&self) -> Dual {
759        <Dual as Float>::asin(*self)
760    }
761    fn acos(&self) -> Dual {
762        <Dual as Float>::acos(*self)
763    }
764    fn atan(&self) -> Dual {
765        <Dual as Float>::atan(*self)
766    }
767    fn atan2(&self, other: &Dual) -> Dual {
768        <Dual as Float>::atan2(*self, *other)
769    }
770    fn acsc(&self) -> Dual {
771        (Dual::new(1.0, 0.0) / *self).asin()
772    }
773    fn asec(&self) -> Dual {
774        (Dual::new(1.0, 0.0) / *self).acos()
775    }
776    fn acot(&self) -> Dual {
777        (Dual::new(1.0, 0.0) / *self).atan()
778    }
779    fn deg2rad(&self) -> Dual {
780        *self * Dual::new(std::f64::consts::PI / 180.0, 0.0)
781    }
782    fn rad2deg(&self) -> Dual {
783        *self * Dual::new(180.0 / std::f64::consts::PI, 0.0)
784    }
785    fn sind(&self) -> Dual {
786        self.deg2rad().sin()
787    }
788    fn cosd(&self) -> Dual {
789        self.deg2rad().cos()
790    }
791    fn tand(&self) -> Dual {
792        self.deg2rad().tan()
793    }
794    fn cscd(&self) -> Dual {
795        self.deg2rad().csc()
796    }
797    fn secd(&self) -> Dual {
798        self.deg2rad().sec()
799    }
800    fn cotd(&self) -> Dual {
801        self.deg2rad().cot()
802    }
803    fn asind(&self) -> Dual {
804        self.asin().rad2deg()
805    }
806    fn acosd(&self) -> Dual {
807        self.acos().rad2deg()
808    }
809    fn atand(&self) -> Dual {
810        self.atan().rad2deg()
811    }
812    fn atan2d(&self, other: &Dual) -> Dual {
813        self.atan2(other).rad2deg()
814    }
815    fn acscd(&self) -> Dual {
816        self.acsc().rad2deg()
817    }
818    fn asecd(&self) -> Dual {
819        self.asec().rad2deg()
820    }
821    fn acotd(&self) -> Dual {
822        self.acot().rad2deg()
823    }
824    fn sinh(&self) -> Dual {
825        <Dual as Float>::sinh(*self)
826    }
827    fn cosh(&self) -> Dual {
828        <Dual as Float>::cosh(*self)
829    }
830    fn tanh(&self) -> Dual {
831        <Dual as Float>::tanh(*self)
832    }
833    fn csch(&self) -> Dual {
834        1.0 / self.sinh()
835    }
836    fn sech(&self) -> Dual {
837        1.0 / self.cosh()
838    }
839    fn coth(&self) -> Dual {
840        1.0 / self.tanh()
841    }
842    fn asinh(&self) -> Dual {
843        <Dual as Float>::asinh(*self)
844    }
845    fn acosh(&self) -> Dual {
846        <Dual as Float>::acosh(*self)
847    }
848    fn atanh(&self) -> Dual {
849        <Dual as Float>::atanh(*self)
850    }
851    fn acsch(&self) -> Dual {
852        (Dual::new(1.0, 0.0) / *self).asinh()
853    }
854    fn asech(&self) -> Dual {
855        (Dual::new(1.0, 0.0) / *self).acosh()
856    }
857    fn acoth(&self) -> Dual {
858        (Dual::new(1.0, 0.0) / *self).atanh()
859    }
860}
861
862// --------
863// TESTING.
864// --------
865
866#[cfg(test)]
867mod tests {
868    use super::*;
869    use linalg_traits::Scalar;
870    use numtest::*;
871    use std::f64::consts::{E, FRAC_PI_4, FRAC_PI_6};
872
873    // Implementing the Compare trait exclusively for testing purposes.
874    impl Compare for Dual {
875        fn is_equal(&self, other: Self) -> bool {
876            let real_equal = self.get_real().is_equal(other.get_real());
877            let dual_equal = self.get_dual().is_equal(other.get_dual());
878            real_equal & dual_equal
879        }
880        fn is_equal_to_decimal(&self, other: Self, decimal: i32) -> (bool, i32) {
881            let (real_equal, real_decimal) = self
882                .get_real()
883                .is_equal_to_decimal(other.get_real(), decimal);
884            let (dual_equal, dual_decimal) = self
885                .get_dual()
886                .is_equal_to_decimal(other.get_dual(), decimal);
887            (real_equal & dual_equal, real_decimal.min(dual_decimal))
888        }
889        fn is_equal_to_atol(&self, other: Self, atol: Self) -> (bool, Self) {
890            let (real_equal, real_abs_diff) = self
891                .get_real()
892                .is_equal_to_atol(other.get_real(), atol.get_real());
893            let (dual_equal, dual_abs_diff) = self
894                .get_dual()
895                .is_equal_to_atol(other.get_dual(), atol.get_dual());
896            (
897                real_equal & dual_equal,
898                Dual::new(real_abs_diff, dual_abs_diff),
899            )
900        }
901        fn is_equal_to_rtol(&self, other: Self, rtol: Self) -> (bool, Self) {
902            let (real_equal, real_rel_diff) = self
903                .get_real()
904                .is_equal_to_rtol(other.get_real(), rtol.get_real());
905            let (dual_equal, dual_rel_diff) = self
906                .get_dual()
907                .is_equal_to_rtol(other.get_dual(), rtol.get_dual());
908            (
909                real_equal & dual_equal,
910                Dual::new(real_rel_diff, dual_rel_diff),
911            )
912        }
913    }
914
915    #[test]
916    fn test_new() {
917        let num1 = Dual::new(1.0, 2.0);
918        let num2 = Dual {
919            real: 1.0,
920            dual: 2.0,
921        };
922        assert_eq!(num1.real, num2.real);
923        assert_eq!(num1.dual, num2.dual);
924    }
925
926    #[test]
927    fn test_get_real() {
928        let num = Dual::new(1.0, 2.0);
929        assert_eq!(num.get_real(), 1.0);
930    }
931
932    #[test]
933    fn test_get_dual() {
934        let num = Dual::new(1.0, 2.0);
935        assert_eq!(num.get_dual(), 2.0);
936    }
937
938    #[test]
939    fn test_add_dual_dual() {
940        assert_eq!(
941            Dual::new(1.0, 2.0) + Dual::new(3.0, 4.0),
942            Dual::new(4.0, 6.0)
943        );
944    }
945
946    #[test]
947    fn test_sub_dual_dual() {
948        assert_eq!(
949            Dual::new(1.0, 2.0) - Dual::new(4.0, 3.0),
950            Dual::new(-3.0, -1.0)
951        );
952    }
953
954    #[test]
955    fn test_mul_dual_dual() {
956        assert_eq!(
957            Dual::new(1.0, 2.0) * Dual::new(3.0, -4.0),
958            Dual::new(3.0, 2.0)
959        );
960    }
961
962    #[test]
963    fn test_div_dual_dual() {
964        assert_eq!(
965            Dual::new(1.0, 2.0) / Dual::new(3.0, 4.0),
966            Dual::new(1.0 / 3.0, 2.0 / 9.0)
967        );
968    }
969
970    #[test]
971    fn test_rem_dual_dual() {
972        // Spot check.
973        assert_eq!(
974            Dual::new(5.0, 2.0) % Dual::new(3.0, 4.0),
975            Dual::new(2.0, -2.0)
976        );
977
978        // Check parity with the truncated definition of the remainder.
979        //  --> Reference: https://en.wikipedia.org/wiki/Modulo#In_programming_languages
980        let a = Dual::new(5.0, 2.0);
981        let n = Dual::new(3.0, 5.0);
982        assert_eq!(a % n, a - n * (a / n).trunc());
983    }
984
985    #[test]
986    fn test_zero() {
987        // Construction
988        assert_eq!(Dual::zero(), Dual::new(0.0, 0.0));
989
990        // Zero-check.
991        assert!(Dual::zero().is_zero());
992        assert!(Dual::new(0.0, 0.0).is_zero());
993
994        // Dual::zero() * Dual = Dual::zero().
995        assert_eq!(Dual::zero() * Dual::new(1.0, 2.0), Dual::zero());
996    }
997
998    #[test]
999    fn test_one() {
1000        // Construction.
1001        assert_eq!(Dual::one(), Dual::new(1.0, 0.0));
1002
1003        // Dual::one() * Dual = Dual.
1004        assert_eq!(Dual::one() * Dual::new(1.0, 2.0), Dual::new(1.0, 2.0));
1005
1006        // Dual::one() * scalar = Dual(scalar, scalar).
1007        assert_eq!(Dual::one() * 5.0, Dual::new(5.0, 0.0));
1008    }
1009
1010    #[test]
1011    fn test_from_str_radix() {
1012        assert_eq!(
1013            Dual::from_str_radix("2.125", 10).unwrap(),
1014            Dual::new(2.125, 0.0)
1015        );
1016    }
1017
1018    #[test]
1019    fn test_to_i64() {
1020        assert_eq!(Dual::new(1.0, 2.0).to_i64().unwrap(), 1_i64);
1021        assert_eq!(Dual::new(-1.0, 2.0).to_i64().unwrap(), -1_i64);
1022    }
1023
1024    #[test]
1025    fn test_to_u64() {
1026        assert_eq!(Dual::new(1.0, 2.0).to_u64().unwrap(), 1_u64);
1027        assert!(Dual::new(-1.0, 2.0).to_u64().is_none());
1028    }
1029
1030    #[test]
1031    fn test_to_f64() {
1032        assert_eq!(Dual::new(1.0, 2.0).to_f64().unwrap(), 1.0_f64);
1033        assert_eq!(Dual::new(-1.0, 2.0).to_f64().unwrap(), -1.0_f64);
1034    }
1035
1036    #[test]
1037    fn test_from_i64() {
1038        assert_eq!(<Dual as NumCast>::from(1_i64).unwrap(), Dual::new(1.0, 0.0));
1039        assert_eq!(
1040            <Dual as NumCast>::from(-1_i64).unwrap(),
1041            Dual::new(-1.0, 0.0)
1042        );
1043    }
1044
1045    #[test]
1046    fn test_from_u64() {
1047        assert_eq!(<Dual as NumCast>::from(1_u64).unwrap(), Dual::new(1.0, 0.0));
1048    }
1049
1050    #[test]
1051    fn test_from_f64() {
1052        assert_eq!(<Dual as NumCast>::from(1_f64).unwrap(), Dual::new(1.0, 0.0));
1053    }
1054
1055    #[test]
1056    fn test_partial_ord() {
1057        // Check <.
1058        assert!(Dual::new(1.0, 2.0) < Dual::new(3.0, 4.0));
1059        assert!(Dual::new(1.0, 4.0) < Dual::new(3.0, 2.0));
1060        assert!(Dual::new(-3.0, -4.0) < Dual::new(-1.0, -2.0));
1061        assert!(Dual::new(-3.0, -2.0) < Dual::new(-1.0, -4.0));
1062
1063        // Check >.
1064        assert!(Dual::new(3.0, 4.0) > Dual::new(1.0, 2.0));
1065        assert!(Dual::new(3.0, 2.0) > Dual::new(1.0, 4.0));
1066        assert!(Dual::new(-1.0, -2.0) > Dual::new(-3.0, -4.0));
1067        assert!(Dual::new(-1.0, -4.0) > Dual::new(-3.0, -2.0));
1068
1069        // Check <=.
1070        assert!(Dual::new(0.0, 2.0) <= Dual::new(1.0, 2.0));
1071        assert!(Dual::new(1.0, 2.0) <= Dual::new(1.0, 2.0));
1072
1073        // Check >=.
1074        assert!(Dual::new(2.0, 2.0) >= Dual::new(1.0, 2.0));
1075        assert!(Dual::new(1.0, 2.0) >= Dual::new(1.0, 2.0));
1076    }
1077
1078    #[test]
1079    fn test_neg() {
1080        assert_eq!(-Dual::new(1.0, 2.0), Dual::new(-1.0, -2.0));
1081        assert_eq!(-Dual::new(1.0, -2.0), Dual::new(-1.0, 2.0));
1082        assert_eq!(-Dual::new(-1.0, 2.0), Dual::new(1.0, -2.0));
1083        assert_eq!(-Dual::new(-1.0, -2.0), Dual::new(1.0, 2.0));
1084    }
1085
1086    #[test]
1087    fn test_nan() {
1088        let num = Dual::nan();
1089        assert!(num.get_real().is_nan());
1090        assert!(num.get_dual().is_nan());
1091    }
1092
1093    #[test]
1094    fn test_infinity() {
1095        let num = Dual::infinity();
1096        assert!(num.get_real().is_infinite() & (num.get_real() > 0.0));
1097        assert!(num.get_dual().is_infinite() & (num.get_dual() > 0.0));
1098    }
1099
1100    #[test]
1101    fn test_neg_infinity() {
1102        let num = Dual::neg_infinity();
1103        assert!(num.get_real().is_infinite() & (num.get_real() < 0.0));
1104        assert!(num.get_dual().is_infinite() & (num.get_dual() < 0.0));
1105    }
1106
1107    #[test]
1108    fn test_neg_zero() {
1109        let num = Dual::neg_zero();
1110        assert!(num.get_real().is_zero());
1111        assert!(num.get_dual().is_zero());
1112    }
1113
1114    #[test]
1115    fn test_min_value() {
1116        let num = Dual::min_value();
1117        assert!(num.get_real() == f64::MIN);
1118        assert!(num.get_dual() == f64::MIN);
1119    }
1120
1121    #[test]
1122    fn test_min_positive_value() {
1123        let num = Dual::min_positive_value();
1124        assert!(num.get_real() == f64::MIN_POSITIVE);
1125        assert!(num.get_dual() == f64::MIN_POSITIVE);
1126    }
1127
1128    #[test]
1129    fn test_max_value() {
1130        let num = Dual::max_value();
1131        assert!(num.get_real() == f64::MAX);
1132        assert!(num.get_dual() == f64::MAX);
1133    }
1134
1135    #[test]
1136    fn test_is_nan() {
1137        assert!(Dual::nan().is_nan());
1138        assert!(Dual::new(f64::NAN, 0.0).is_nan());
1139        assert!(!Dual::new(0.0, f64::NAN).is_nan());
1140        assert!(!Dual::new(0.0, 0.0).is_nan());
1141    }
1142
1143    #[test]
1144    fn test_is_infinite() {
1145        assert!(Dual::infinity().is_infinite());
1146        assert!(Dual::neg_infinity().is_infinite());
1147        assert!(Dual::new(f64::INFINITY, 0.0).is_infinite());
1148        assert!(Dual::new(f64::NEG_INFINITY, 0.0).is_infinite());
1149        assert!(!Dual::new(0.0, f64::INFINITY).is_infinite());
1150        assert!(!Dual::new(0.0, f64::NEG_INFINITY).is_infinite());
1151        assert!(!Dual::new(0.0, 0.0).is_infinite());
1152    }
1153
1154    #[test]
1155    fn test_is_finite() {
1156        assert!(!Dual::infinity().is_finite());
1157        assert!(!Dual::neg_infinity().is_finite());
1158        assert!(!Dual::new(f64::INFINITY, 0.0).is_finite());
1159        assert!(!Dual::new(f64::NEG_INFINITY, 0.0).is_finite());
1160        assert!(Dual::new(0.0, f64::INFINITY).is_finite());
1161        assert!(Dual::new(0.0, f64::NEG_INFINITY).is_finite());
1162        assert!(Dual::new(0.0, 0.0).is_finite());
1163    }
1164
1165    /// # References
1166    ///
1167    /// * https://docs.rs/num-traits/latest/num_traits/float/trait.Float.html#tymethod.is_normal
1168    ///
1169    /// # Note
1170    ///
1171    /// For each of these tests, we use a dual part of `f64::NAN` to ensure that `is_normal` is only
1172    /// checking the real part.
1173    #[test]
1174    fn test_is_normal() {
1175        // Normal (for these checks we use a not-normal dual part to ensure that only the real part
1176        // is being checked).
1177        assert!(Dual::new(1.0, f64::NAN).is_normal());
1178        assert!(Dual::new(f64::MIN_POSITIVE, f64::NAN).is_normal());
1179        assert!(Dual::new(f64::MAX, f64::NAN).is_normal());
1180
1181        // Not normal (for these checks we use a normal dual part to ensure that only the real part
1182        // is being checked).
1183        assert!(!Dual::new(0.0, 1.0).is_normal()); // Zero.
1184        assert!(!Dual::new(f64::NAN, 1.0).is_normal()); // NaN.
1185        assert!(!Dual::new(f64::INFINITY, 1.0).is_normal()); // Infinite.
1186        assert!(!Dual::new(f64::NEG_INFINITY, 1.0).is_normal()); // Infinite.
1187        assert!(!Dual::new(1.0e-308_f64, 1.0).is_normal()); // Subnormal (between 0 and f64::MIN).
1188    }
1189
1190    #[test]
1191    fn test_classify() {
1192        // Normal (for these checks we use a not-normal dual part to ensure that only the real part
1193        // is being checked).
1194        assert_eq!(Dual::new(1.0, f64::NAN).classify(), FpCategory::Normal);
1195        assert_eq!(
1196            Dual::new(f64::MIN_POSITIVE, f64::NAN).classify(),
1197            FpCategory::Normal
1198        );
1199        assert_eq!(Dual::new(f64::MAX, f64::NAN).classify(), FpCategory::Normal);
1200
1201        // Not normal (for these checks we use a normal dual part to ensure that only the real part
1202        // is being checked).
1203        assert_eq!(Dual::new(0.0, 1.0).classify(), FpCategory::Zero);
1204        assert_eq!(Dual::new(f64::NAN, 1.0).classify(), FpCategory::Nan);
1205        assert_eq!(
1206            Dual::new(f64::INFINITY, 1.0).classify(),
1207            FpCategory::Infinite
1208        );
1209        assert_eq!(
1210            Dual::new(f64::NEG_INFINITY, 1.0).classify(),
1211            FpCategory::Infinite
1212        );
1213        assert_eq!(
1214            Dual::new(1.0e-308_f64, 1.0).classify(),
1215            FpCategory::Subnormal
1216        );
1217    }
1218
1219    #[test]
1220    fn test_floor() {
1221        assert_eq!(Dual::new(2.7, 2.7).floor(), Dual::new(2.0, 0.0));
1222        assert_eq!(Dual::new(-2.7, -2.7).floor(), Dual::new(-3.0, 0.0));
1223    }
1224
1225    #[test]
1226    fn test_ceil() {
1227        assert_eq!(Dual::new(2.7, 2.7).ceil(), Dual::new(3.0, 0.0));
1228        assert_eq!(Dual::new(-2.7, -2.7).ceil(), Dual::new(-2.0, 0.0));
1229    }
1230
1231    #[test]
1232    fn test_round() {
1233        assert_eq!(Dual::new(2.3, 2.3).round(), Dual::new(2.0, 0.0));
1234        assert_eq!(Dual::new(2.7, 2.7).round(), Dual::new(3.0, 0.0));
1235        assert_eq!(Dual::new(-2.7, -2.7).round(), Dual::new(-3.0, 0.0));
1236        assert_eq!(Dual::new(-2.3, -2.3).round(), Dual::new(-2.0, 0.0));
1237    }
1238
1239    #[test]
1240    fn test_trunc() {
1241        assert_eq!(Dual::new(2.3, 2.3).trunc(), Dual::new(2.0, 0.0));
1242        assert_eq!(Dual::new(2.7, 2.7).trunc(), Dual::new(2.0, 0.0));
1243        assert_eq!(Dual::new(-2.7, -2.7).trunc(), Dual::new(-2.0, 0.0));
1244        assert_eq!(Dual::new(-2.3, -2.3).trunc(), Dual::new(-2.0, 0.0));
1245    }
1246
1247    #[test]
1248    fn test_fract() {
1249        assert_eq!(Dual::new(2.5, 2.5).fract(), Dual::new(0.5, 2.5));
1250        assert_eq!(Dual::new(-2.5, -2.5).fract(), Dual::new(-0.5, -2.5));
1251    }
1252
1253    #[test]
1254    fn test_abs() {
1255        assert_eq!(Dual::new(1.0, 2.0).abs(), Dual::new(1.0, 2.0));
1256        assert_eq!(Dual::new(-1.0, -2.0).abs(), Dual::new(1.0, 2.0));
1257        assert_eq!(Dual::new(-1.0, 2.0).abs(), Dual::new(1.0, -2.0));
1258    }
1259
1260    #[test]
1261    fn test_signum() {
1262        assert_eq!(Dual::new(2.0, 4.0).signum(), Dual::new(1.0, 0.0));
1263        assert_eq!(Dual::new(-2.0, -4.0).signum(), Dual::new(-1.0, 0.0));
1264    }
1265
1266    #[test]
1267    fn test_is_sign_positive() {
1268        assert!(Dual::new(2.0, -4.0).is_sign_positive());
1269        assert!(!Dual::new(-2.0, 4.0).is_sign_positive());
1270    }
1271
1272    #[test]
1273    fn test_is_sign_negative() {
1274        assert!(Dual::new(-2.0, 4.0).is_sign_negative());
1275        assert!(!Dual::new(2.0, -4.0).is_sign_negative());
1276    }
1277
1278    #[test]
1279    fn test_mul_add() {
1280        let a = Dual::new(1.0, 3.0);
1281        let b = Dual::new(-2.0, 5.0);
1282        let c = Dual::new(10.0, -4.0);
1283        assert_eq!(c.mul_add(a, b), (c * a) + b);
1284    }
1285
1286    #[test]
1287    fn test_recip() {
1288        assert_eq!(Dual::new(2.0, -5.0).recip(), Dual::new(0.5, 1.25));
1289    }
1290
1291    #[test]
1292    fn test_powi() {
1293        assert_eq!(Dual::new(2.0, -5.0).powi(0), Dual::new(1.0, 0.0));
1294        assert_eq!(Dual::new(2.0, -5.0).powi(1), Dual::new(2.0, -5.0));
1295        assert_eq!(Dual::new(2.0, -5.0).powi(2), Dual::new(4.0, -20.0));
1296        assert_eq!(Dual::new(2.0, -5.0).powi(3), Dual::new(8.0, -60.0));
1297    }
1298
1299    #[test]
1300    fn test_powf() {
1301        // Test against powi for integer powers.
1302        assert_eq!(
1303            Dual::new(2.0, -5.0).powf(Dual::new(0.0, 0.0)),
1304            Dual::new(2.0, -5.0).powi(0)
1305        );
1306        assert_eq!(
1307            Dual::new(2.0, -5.0).powf(Dual::new(1.0, 0.0)),
1308            Dual::new(2.0, -5.0).powi(1)
1309        );
1310        assert_eq!(
1311            Dual::new(2.0, -5.0).powf(Dual::new(2.0, 0.0)),
1312            Dual::new(2.0, -5.0).powi(2)
1313        );
1314        assert_equal_to_decimal!(
1315            Dual::new(2.0, -5.0).powf(Dual::new(3.0, 0.0)),
1316            Dual::new(2.0, -5.0).powi(3),
1317            14
1318        );
1319
1320        // Test against sqrt.
1321        assert_equal_to_decimal!(
1322            Dual::new(2.0, -5.0).powf(Dual::new(0.5, 0.0)),
1323            Dual::new(2.0, -5.0).sqrt(),
1324            15
1325        );
1326
1327        // Test against cbrt.
1328        assert_equal_to_decimal!(
1329            Dual::new(2.0, -5.0).powf(Dual::new(1.0 / 3.0, 0.0)),
1330            Dual::new(2.0, -5.0).cbrt(),
1331            15
1332        );
1333
1334        // Spot check.
1335        assert_eq!(
1336            Dual::new(2.0, -5.0).powf(Dual::new(-5.0, 4.0)),
1337            Dual::new(0.03125, 0.4772683975699932)
1338        );
1339    }
1340
1341    #[test]
1342    fn test_sqrt() {
1343        assert_eq!(Dual::new(4.0, 25.0).sqrt(), Dual::new(2.0, 6.25));
1344    }
1345
1346    #[test]
1347    fn test_exp() {
1348        assert_eq!(
1349            Dual::new(2.0, -3.0).exp(),
1350            Dual::new(2.0.exp(), -3.0 * 2.0.exp())
1351        );
1352    }
1353
1354    #[test]
1355    fn test_exp2() {
1356        assert_eq!(
1357            Dual::new(2.0, -3.0).exp2(),
1358            Dual::new(2.0.exp2(), -8.317766166719343)
1359        );
1360    }
1361
1362    #[test]
1363    fn test_ln() {
1364        assert_eq!(Dual::new(5.0, 8.0).ln(), Dual::new(5.0.ln(), 8.0 / 5.0));
1365    }
1366
1367    #[test]
1368    fn test_log() {
1369        assert_eq!(
1370            Dual::new(5.0, 8.0).log(Dual::new(4.5, 0.0)),
1371            Dual::new(5.0.log(4.5), 1.0637750447080176)
1372        );
1373    }
1374
1375    #[test]
1376    fn test_log2() {
1377        assert_eq!(
1378            Dual::new(5.0, 8.0).log2(),
1379            Dual::new(5.0.log2(), 2.3083120654223412)
1380        );
1381    }
1382
1383    #[test]
1384    fn test_log10() {
1385        assert_equal_to_decimal!(
1386            Dual::new(5.0, 8.0).log10(),
1387            Dual::new(5.0.log10(), 0.6948711710452028),
1388            16
1389        );
1390    }
1391
1392    #[test]
1393    fn test_max() {
1394        assert_eq!(
1395            Dual::new(1.0, 2.0).max(Dual::new(3.0, 4.0)),
1396            Dual::new(3.0, 4.0)
1397        );
1398        assert_eq!(
1399            Dual::new(3.0, 2.0).max(Dual::new(1.0, 4.0)),
1400            Dual::new(3.0, 2.0)
1401        );
1402        assert_eq!(
1403            Dual::new(3.0, 4.0).max(Dual::new(1.0, 2.0)),
1404            Dual::new(3.0, 4.0)
1405        );
1406        assert_eq!(
1407            Dual::new(-1.0, 2.0).max(Dual::new(-3.0, 4.0)),
1408            Dual::new(-1.0, 2.0)
1409        );
1410    }
1411
1412    #[test]
1413    fn test_min() {
1414        assert_eq!(
1415            Dual::new(1.0, 2.0).min(Dual::new(3.0, 4.0)),
1416            Dual::new(1.0, 2.0)
1417        );
1418        assert_eq!(
1419            Dual::new(3.0, 2.0).min(Dual::new(1.0, 4.0)),
1420            Dual::new(1.0, 4.0)
1421        );
1422        assert_eq!(
1423            Dual::new(3.0, 4.0).min(Dual::new(1.0, 2.0)),
1424            Dual::new(1.0, 2.0)
1425        );
1426        assert_eq!(
1427            Dual::new(-1.0, 2.0).min(Dual::new(-3.0, 4.0)),
1428            Dual::new(-3.0, 4.0)
1429        );
1430    }
1431
1432    #[test]
1433    fn test_abs_sub() {
1434        assert_eq!(
1435            Dual::new(4.0, 5.0).abs_sub(Dual::new(2.0, 8.0)),
1436            Dual::new(2.0, -3.0)
1437        );
1438    }
1439
1440    #[test]
1441    fn test_cbrt() {
1442        assert_eq!(Dual::new(8.0, 27.0).cbrt(), Dual::new(2.0, 2.25));
1443    }
1444
1445    #[test]
1446    fn test_hypot() {
1447        // Spot check.
1448        assert_eq!(
1449            Dual::new(1.0, 2.0).hypot(Dual::new(3.0, 4.0)),
1450            Dual::new(3.1622776601683795, 4.427188724235731)
1451        );
1452
1453        // Check parity with Euclidian norm.
1454        assert_eq!(
1455            Dual::new(1.0, 2.0).hypot(Dual::new(3.0, 4.0)),
1456            (Dual::new(1.0, 2.0).powi(2) + Dual::new(3.0, 4.0).powi(2)).sqrt()
1457        );
1458    }
1459
1460    #[test]
1461    fn test_sin() {
1462        assert_eq!(Dual::new(FRAC_PI_6, 2.0).sin(), Dual::new(0.5, 3.0.sqrt()));
1463    }
1464
1465    #[test]
1466    fn test_cos() {
1467        assert_eq!(
1468            Dual::new(FRAC_PI_6, 2.0).cos(),
1469            Dual::new(3.0.sqrt() / 2.0, -1.0)
1470        );
1471    }
1472
1473    #[test]
1474    fn test_tan() {
1475        assert_equal_to_decimal!(
1476            Dual::new(FRAC_PI_6, 2.0).tan(),
1477            Dual::new(3.0.sqrt() / 3.0, 8.0 / 3.0),
1478            15
1479        );
1480    }
1481
1482    #[test]
1483    fn test_asin() {
1484        assert_equal_to_decimal!(
1485            Dual::new(0.5, 3.0).asin(),
1486            Dual::new(FRAC_PI_6, 3.0 / 0.75.sqrt()),
1487            16
1488        );
1489    }
1490
1491    #[test]
1492    fn test_acos() {
1493        assert_equal_to_decimal!(
1494            Dual::new(3.0.sqrt() / 2.0, 3.0).acos(),
1495            Dual::new(FRAC_PI_6, -6.0),
1496            15
1497        );
1498    }
1499
1500    #[test]
1501    fn test_atan() {
1502        assert_eq!(Dual::new(1.0, 3.0).atan(), Dual::new(FRAC_PI_4, 1.5));
1503    }
1504
1505    #[test]
1506    fn test_atan2() {
1507        let x = Dual::new(3.0, 2.0);
1508        let y = Dual::new(-3.0, 5.0);
1509        let angle_expected = Dual::new(-FRAC_PI_4, 7.0 / 6.0);
1510        assert_eq!(y.atan2(x), angle_expected);
1511    }
1512
1513    #[test]
1514    fn test_sin_cos() {
1515        let (sin, cos) = Dual::new(FRAC_PI_6, 2.0).sin_cos();
1516        assert_eq!(sin, Dual::new(0.5, 3.0.sqrt()));
1517        assert_eq!(cos, Dual::new(3.0.sqrt() / 2.0, -1.0));
1518    }
1519
1520    #[test]
1521    fn test_exp_m1() {
1522        assert_eq!(
1523            Dual::new(3.0, 5.0).exp_m1(),
1524            Dual::new(3.0, 5.0).exp() - Dual::one()
1525        );
1526    }
1527
1528    #[test]
1529    fn test_ln_1p() {
1530        assert_eq!(
1531            Dual::new(3.0, 5.0).ln_1p(),
1532            (Dual::new(3.0, 5.0) + Dual::one()).ln()
1533        );
1534    }
1535
1536    #[test]
1537    fn test_sinh() {
1538        assert_equal_to_decimal!(
1539            Dual::new(1.0, 2.0).sinh(),
1540            Dual::new(((E * E) - 1.0) / (2.0 * E), ((E * E) + 1.0) / E),
1541            15
1542        );
1543    }
1544
1545    #[test]
1546    fn test_cosh() {
1547        assert_equal_to_decimal!(
1548            Dual::new(1.0, 2.0).cosh(),
1549            Dual::new(((E * E) + 1.0) / (2.0 * E), ((E * E) - 1.0) / E),
1550            15
1551        );
1552    }
1553
1554    #[test]
1555    fn test_tanh() {
1556        assert_equal_to_decimal!(
1557            Dual::new(1.0, 2.0).tanh(),
1558            Dual::new(
1559                (1.0 - E.powi(-2)) / (1.0 + E.powi(-2)),
1560                2.0 * ((2.0 * E) / (E.powi(2) + 1.0)).powi(2)
1561            ),
1562            15
1563        );
1564    }
1565
1566    #[test]
1567    fn test_asinh() {
1568        assert_eq!(Dual::new(1.0, 2.0).sinh().asinh(), Dual::new(1.0, 2.0));
1569    }
1570
1571    #[test]
1572    fn test_acosh() {
1573        assert_eq!(Dual::new(1.0, 2.0).cosh().acosh(), Dual::new(1.0, 2.0));
1574    }
1575
1576    #[test]
1577    fn test_atanh() {
1578        assert_equal_to_decimal!(Dual::new(1.0, 2.0).tanh().atanh(), Dual::new(1.0, 2.0), 16);
1579    }
1580
1581    #[test]
1582    fn test_integer_decode() {
1583        assert_eq!(
1584            Dual::new(1.2345e-5, 6.789e-7).integer_decode(),
1585            (1.2345e-5).integer_decode()
1586        );
1587    }
1588
1589    #[test]
1590    fn test_add_assign_dual_dual() {
1591        let mut a = Dual::new(1.0, 2.0);
1592        a += Dual::new(3.0, 4.0);
1593        assert_eq!(a, Dual::new(4.0, 6.0));
1594    }
1595
1596    #[test]
1597    fn test_sub_assign_dual_dual() {
1598        let mut a = Dual::new(1.0, 2.0);
1599        a -= Dual::new(4.0, 3.0);
1600        assert_eq!(a, Dual::new(-3.0, -1.0));
1601    }
1602
1603    #[test]
1604    fn test_mul_assign_dual_dual() {
1605        let mut a = Dual::new(1.0, 2.0);
1606        a *= Dual::new(3.0, -4.0);
1607        assert_eq!(a, Dual::new(3.0, 2.0));
1608    }
1609
1610    #[test]
1611    fn test_div_assign_dual_dual() {
1612        let mut a = Dual::new(1.0, 2.0);
1613        a /= Dual::new(3.0, 4.0);
1614        assert_eq!(a, Dual::new(1.0 / 3.0, 2.0 / 9.0));
1615    }
1616
1617    #[test]
1618    fn test_rem_assign_dual_dual() {
1619        let mut a = Dual::new(5.0, 2.0);
1620        let b = Dual::new(3.0, 4.0);
1621        a %= b;
1622        assert_eq!(a, Dual::new(2.0, -2.0));
1623    }
1624
1625    #[test]
1626    fn test_add_dual_f64() {
1627        assert_eq!(Dual::new(1.0, 2.0) + 3.0, Dual::new(4.0, 2.0));
1628    }
1629
1630    #[test]
1631    fn test_sub_dual_f64() {
1632        assert_eq!(Dual::new(1.0, 2.0) - 3.0, Dual::new(-2.0, 2.0));
1633    }
1634
1635    #[test]
1636    fn test_mul_dual_f64() {
1637        assert_eq!(Dual::new(1.0, -2.0) * 3.0, Dual::new(3.0, -6.0));
1638    }
1639
1640    #[test]
1641    fn test_div_dual_f64() {
1642        assert_eq!(Dual::new(1.0, 2.0) / 4.0, Dual::new(0.25, 0.5));
1643    }
1644
1645    #[test]
1646    fn test_rem_dual_f64() {
1647        // Spot check.
1648        assert_eq!(Dual::new(5.0, 1.0) % 3.0, Dual::new(2.0, 1.0));
1649
1650        // Check parity with the truncated definition of the remainder.
1651        //  --> Reference: https://en.wikipedia.org/wiki/Modulo#In_programming_languages
1652        let a = Dual::new(5.0, 1.0);
1653        let n = 3.0;
1654        assert_eq!(a % n, a - n * (a / n).trunc());
1655    }
1656
1657    #[test]
1658    fn test_add_assign_dual_f64() {
1659        let mut a = Dual::new(1.0, 2.0);
1660        a += 3.0;
1661        assert_eq!(a, Dual::new(4.0, 2.0));
1662    }
1663
1664    #[test]
1665    fn test_sub_assign_dual_f64() {
1666        let mut a = Dual::new(1.0, 2.0);
1667        a -= 3.0;
1668        assert_eq!(a, Dual::new(-2.0, 2.0));
1669    }
1670
1671    #[test]
1672    fn test_mul_assign_dual_f64() {
1673        let mut a = Dual::new(2.0, -3.0);
1674        a *= 5.0;
1675        assert_eq!(a, Dual::new(10.0, -15.0));
1676    }
1677
1678    #[test]
1679    fn test_div_assign_dual_f64() {
1680        let mut a = Dual::new(1.0, 2.0);
1681        a /= 4.0;
1682        assert_eq!(a, Dual::new(0.25, 0.5));
1683    }
1684
1685    #[test]
1686    fn test_rem_assign_dual_f64() {
1687        let mut a = Dual::new(5.0, 1.0);
1688        let n = 3.0;
1689        a %= n;
1690        assert_eq!(a, Dual::new(2.0, 1.0));
1691    }
1692
1693    #[test]
1694    fn test_add_f64_dual() {
1695        assert_eq!(1.0 + Dual::new(2.0, 3.0), Dual::new(3.0, 3.0));
1696    }
1697
1698    #[test]
1699    fn test_sub_f64_dual() {
1700        assert_eq!(1.0 - Dual::new(2.0, 3.0), Dual::new(-1.0, 3.0));
1701    }
1702
1703    #[test]
1704    fn test_mul_f64_dual() {
1705        assert_eq!(5.0 * Dual::new(2.0, -3.0), Dual::new(10.0, -15.0));
1706    }
1707
1708    #[test]
1709    fn test_div_f64_dual() {
1710        assert_eq!(5.0 / Dual::new(2.0, -3.0), Dual::new(2.5, 3.75));
1711    }
1712
1713    #[test]
1714    fn test_rem_f64_dual() {
1715        // Spot check.
1716        assert_eq!(5.0 % Dual::new(2.0, -3.0), Dual::new(1.0, 6.0));
1717
1718        // Check parity with "Dual % Dual" implementation.
1719        assert_eq!(
1720            5.0 % Dual::new(2.0, -3.0),
1721            Dual::new(5.0, 0.0) % Dual::new(2.0, -3.0)
1722        );
1723    }
1724
1725    // This just verifies that the scalar trait was fully implemented.
1726    #[test]
1727    fn test_scalar() {
1728        fn add_scalar<S: Scalar>(x: S, y: S) -> S {
1729            x + y
1730        }
1731        assert_eq!(
1732            add_scalar(Dual::new(5.0, 4.0), Dual::new(3.0, 2.0)),
1733            Dual::new(8.0, 6.0)
1734        );
1735    }
1736}
numdiff/automatic_differentiation/dual.rs

numdiff/automatic_differentiation/
dual.rs
