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#[derive(Copy, Clone, Debug)]
11#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
12pub struct HyperHyperDual<T> {
13 pub re: T,
15 pub eps1: T,
17 pub eps2: T,
19 pub eps3: T,
21 pub eps1eps2: T,
23 pub eps1eps3: T,
25 pub eps2eps3: T,
27 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 #[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 #[inline]
65 pub fn derivative1(mut self) -> Self {
66 self.eps1 = T::one();
67 self
68 }
69
70 #[inline]
72 pub fn derivative2(mut self) -> Self {
73 self.eps2 = T::one();
74 self
75 }
76
77 #[inline]
79 pub fn derivative3(mut self) -> Self {
80 self.eps3 = T::one();
81 self
82 }
83
84 #[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
166impl<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);