Skip to main content

num_dual/datatypes/
hyperhyperdual.rs

1use crate::{DualNum, DualNumFloat, DualStruct};
2use num_traits::{FloatConst, FromPrimitive, Inv, Num, One, Signed, Zero};
3#[cfg(feature = "serde")]
4use serde::{Deserialize, Serialize};
5use std::fmt;
6use std::iter::{Product, Sum};
7use std::ops::*;
8
9/// A scalar hyper-hyper-dual number for the calculation of third partial derivatives.
10#[derive(Copy, Clone, Debug)]
11#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
12pub struct HyperHyperDual<T> {
13    /// Real part of the hyper-hyper-dual number
14    pub re: T,
15    /// First partial derivative part of the hyper-hyper-dual number
16    pub eps1: T,
17    /// First partial derivative part of the hyper-hyper-dual number
18    pub eps2: T,
19    /// First partial derivative part of the hyper-hyper-dual number
20    pub eps3: T,
21    /// Second partial derivative part of the hyper-hyper-dual number
22    pub eps1eps2: T,
23    /// Second partial derivative part of the hyper-hyper-dual number
24    pub eps1eps3: T,
25    /// Second partial derivative part of the hyper-hyper-dual number
26    pub eps2eps3: T,
27    /// Third partial derivative part of the hyper-hyper-dual number
28    pub eps1eps2eps3: T,
29}
30
31#[cfg(feature = "ndarray")]
32impl<T: DualNum> ndarray::ScalarOperand for HyperHyperDual<T> {}
33
34pub type HyperHyperDual32 = HyperHyperDual<f32>;
35pub type HyperHyperDual64 = HyperHyperDual<f64>;
36
37impl<T: DualNum> HyperHyperDual<T> {
38    /// Create a new hyper-hyper-dual number from its fields.
39    #[inline]
40    #[expect(clippy::too_many_arguments)]
41    pub fn new(
42        re: T,
43        eps1: T,
44        eps2: T,
45        eps3: T,
46        eps1eps2: T,
47        eps1eps3: T,
48        eps2eps3: T,
49        eps1eps2eps3: T,
50    ) -> Self {
51        Self {
52            re,
53            eps1,
54            eps2,
55            eps3,
56            eps1eps2,
57            eps1eps3,
58            eps2eps3,
59            eps1eps2eps3,
60        }
61    }
62
63    /// Set the partial derivative part w.r.t. the 1st variable to 1.
64    #[inline]
65    pub fn derivative1(mut self) -> Self {
66        self.eps1 = T::one();
67        self
68    }
69
70    /// Set the partial derivative part w.r.t. the 2nd variable to 1.
71    #[inline]
72    pub fn derivative2(mut self) -> Self {
73        self.eps2 = T::one();
74        self
75    }
76
77    /// Set the partial derivative part w.r.t. the 3rd variable to 1.
78    #[inline]
79    pub fn derivative3(mut self) -> Self {
80        self.eps3 = T::one();
81        self
82    }
83
84    /// Create a new hyper-hyper-dual number from the real part.
85    #[inline]
86    pub fn from_re(re: T) -> Self {
87        Self::new(
88            re,
89            T::zero(),
90            T::zero(),
91            T::zero(),
92            T::zero(),
93            T::zero(),
94            T::zero(),
95            T::zero(),
96        )
97    }
98}
99
100impl<T: DualNum> HyperHyperDual<T> {
101    #[inline]
102    fn chain_rule(&self, f0: T, f1: T, f2: T, f3: T) -> Self {
103        Self::new(
104            f0,
105            f1.clone() * &self.eps1,
106            f1.clone() * &self.eps2,
107            f1.clone() * &self.eps3,
108            f1.clone() * &self.eps1eps2 + f2.clone() * &self.eps1 * &self.eps2,
109            f1.clone() * &self.eps1eps3 + f2.clone() * &self.eps1 * &self.eps3,
110            f1.clone() * &self.eps2eps3 + f2.clone() * &self.eps2 * &self.eps3,
111            f1 * &self.eps1eps2eps3
112                + f2 * (self.eps1.clone() * &self.eps2eps3
113                    + self.eps2.clone() * &self.eps1eps3
114                    + self.eps3.clone() * &self.eps1eps2)
115                + f3 * self.eps1.clone() * &self.eps2 * &self.eps3,
116        )
117    }
118}
119
120impl<T: DualNum> Mul<&HyperHyperDual<T>> for &HyperHyperDual<T> {
121    type Output = HyperHyperDual<T>;
122    #[inline]
123    fn mul(self, rhs: &HyperHyperDual<T>) -> HyperHyperDual<T> {
124        HyperHyperDual::new(
125            self.re.clone() * &rhs.re,
126            self.eps1.clone() * &rhs.re + self.re.clone() * &rhs.eps1,
127            self.eps2.clone() * &rhs.re + self.re.clone() * &rhs.eps2,
128            self.eps3.clone() * &rhs.re + self.re.clone() * &rhs.eps3,
129            self.eps1eps2.clone() * &rhs.re
130                + self.eps1.clone() * &rhs.eps2
131                + self.eps2.clone() * &rhs.eps1
132                + self.re.clone() * &rhs.eps1eps2,
133            self.eps1eps3.clone() * &rhs.re
134                + self.eps1.clone() * &rhs.eps3
135                + self.eps3.clone() * &rhs.eps1
136                + self.re.clone() * &rhs.eps1eps3,
137            self.eps2eps3.clone() * &rhs.re
138                + self.eps2.clone() * &rhs.eps3
139                + self.eps3.clone() * &rhs.eps2
140                + self.re.clone() * &rhs.eps2eps3,
141            self.eps1eps2eps3.clone() * &rhs.re
142                + self.eps1.clone() * &rhs.eps2eps3
143                + self.eps2.clone() * &rhs.eps1eps3
144                + self.eps3.clone() * &rhs.eps1eps2
145                + self.eps2eps3.clone() * &rhs.eps1
146                + self.eps1eps3.clone() * &rhs.eps2
147                + self.eps1eps2.clone() * &rhs.eps3
148                + self.re.clone() * &rhs.eps1eps2eps3,
149        )
150    }
151}
152
153impl<T: DualNum> Div<&HyperHyperDual<T>> for &HyperHyperDual<T> {
154    type Output = HyperHyperDual<T>;
155    #[inline]
156    fn div(self, rhs: &HyperHyperDual<T>) -> HyperHyperDual<T> {
157        let rec = T::one() / &rhs.re;
158        let f0 = rec.clone();
159        let f1 = -f0.clone() * &rec;
160        let f2 = -f1.clone() * &rec * T::Primitive::TWO;
161        let f3 = -f2.clone() * rec * T::Primitive::THREE;
162        self * rhs.chain_rule(f0, f1, f2, f3)
163    }
164}
165
166/* string conversions */
167impl<T: fmt::Display> fmt::Display for HyperHyperDual<T> {
168    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
169        write!(
170            f,
171            "{} + {}ε1 + {}ε2 + {}ε3 + {}ε1ε2 + {}ε1ε3 + {}ε2ε3 + {}ε1ε2ε3",
172            self.re,
173            self.eps1,
174            self.eps2,
175            self.eps3,
176            self.eps1eps2,
177            self.eps1eps3,
178            self.eps2eps3,
179            self.eps1eps2eps3
180        )
181    }
182}
183
184impl_third_derivatives!(
185    HyperHyperDual,
186    [eps1, eps2, eps3, eps1eps2, eps1eps3, eps2eps3, eps1eps2eps3]
187);
188impl_dual!(
189    HyperHyperDual,
190    [eps1, eps2, eps3, eps1eps2, eps1eps3, eps2eps3, eps1eps2eps3]
191);
192#[cfg(feature = "nalgebra")]
193impl_nalgebra!(
194    HyperHyperDual,
195    [eps1, eps2, eps3, eps1eps2, eps1eps3, eps2eps3, eps1eps2eps3]
196);