1#[macro_export]
2macro_rules! impl_derivatives {
3 ($deriv:ident, $nderiv:expr, $struct:ident, [$($im:ident),*]$(, [$($dim:tt),*]$(, [$($ddim:tt),*])*)?) => {
4 impl<T: DualNum<F>, F: DualNumFloat$($(, $dim: Dim)*)?> DualNum<F> for $struct<T, F$($(, $dim)*)?>
5 where
6 $($(DefaultAllocator: Allocator<$($ddim,)*>),*)?
7 {
8 const NDERIV: usize = T::NDERIV + $nderiv;
9
10 type InnerDual = T;
11 fn from_re(inner: Self::InnerDual) -> Self {
12 Self::from_re(inner)
13 }
14
15 #[inline]
16 fn recip(&self) -> Self {
17 let rec = self.re.recip();
18 let f0 = rec.clone();
19 first!($deriv, let f1 = -f0.clone() * &rec;);
20 second!($deriv, let f2 = f1.clone() * &rec * F::from(-2.0).unwrap(););
21 third!($deriv, let f3 = f2.clone() * rec * F::from(-3.0).unwrap(););
22 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
23 }
24
25 #[inline]
26 fn powi(&self, exp: i32) -> Self {
27 match exp {
28 0 => Self::one(),
29 1 => self.clone(),
30 2 => self * self,
31 _ => {
32 let pow3 = self.re.powi(exp - 3);
33 let f0 = pow3.clone() * &self.re * &self.re * &self.re;
34 first!($deriv, let f1 = pow3.clone() * &self.re * &self.re * F::from(exp).unwrap(););
35 second!($deriv, let f2 = pow3.clone() * &self.re * F::from(exp * (exp - 1)).unwrap(););
36 third!($deriv, let f3 = pow3 * F::from(exp * (exp - 1) * (exp - 2)).unwrap(););
37 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
38 }
39 }
40 }
41
42 #[inline]
43 fn powf(&self, n: F) -> Self {
44 if n.is_zero() {
45 Self::one()
46 } else if n.is_one() {
47 self.clone()
48 } else if (n - F::one() - F::one()).abs() < F::epsilon() {
49 self * self
50 } else {
51 let n1 = n - F::one();
52 let n2 = n1 - F::one();
53 let n3 = n2 - F::one();
54 let pow3 = self.re.powf(n3);
55 let f0 = pow3.clone() * &self.re * &self.re * &self.re;
56 first!($deriv, let f1 = pow3.clone() * &self.re * &self.re * n;);
57 second!($deriv, let f2 = pow3.clone() * &self.re * n * n1;);
58 third!($deriv, let f3 = pow3 * n * n1 * n2;);
59 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
60 }
61 }
62
63 #[inline]
64 fn sqrt(&self) -> Self {
65 first!($deriv, let rec = self.re.recip(););
66 first!($deriv, let half = F::from(0.5).unwrap(););
67 let f0 = self.re.sqrt();
68 first!($deriv, let f1 = f0.clone() * &rec * half;);
69 second!($deriv, let f2 = -f1.clone() * &rec * half;);
70 third!($deriv, let f3 = f2.clone() * rec * (-F::one() - half););
71 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
72 }
73
74 #[inline]
75 fn cbrt(&self) -> Self {
76 first!($deriv, let rec = self.re.recip(););
77 first!($deriv, let third = F::from(1.0 / 3.0).unwrap(););
78 let f0 = self.re.cbrt();
79 first!($deriv, let f1 = f0.clone() * &rec * third;);
80 second!($deriv, let f2 = f1.clone() * &rec * (third - F::one()););
81 third!($deriv, let f3 = f2.clone() * rec * (third - F::one() - F::one()););
82 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
83 }
84
85
86 #[inline]
87 fn exp(&self) -> Self {
88 let f = self.re.exp();
89 chain_rule!($deriv, Self::chain_rule(self, f.clone(), f.clone(), f.clone(), f))
90 }
91
92 #[inline]
93 fn exp2(&self) -> Self {
94 first!($deriv, let ln2 = F::from(2.0).unwrap().ln(););
95 let f0 = self.re.exp2();
96 first!($deriv, let f1 = f0.clone() * ln2;);
97 second!($deriv, let f2 = f1.clone() * ln2;);
98 third!($deriv, let f3 = f2.clone() * ln2;);
99 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
100 }
101
102 #[inline]
103 fn exp_m1(&self) -> Self {
104 let f0 = self.re.exp_m1();
105 first!($deriv, let f1 = self.re.exp(););
106 chain_rule!($deriv, Self::chain_rule(self, f0, f1.clone(), f1.clone(), f1))
107 }
108
109 #[inline]
110 fn ln(&self) -> Self {
111 first!($deriv, let rec = self.re.recip(););
112 let f0 = self.re.ln();
113 first!($deriv, let f1 = rec.clone(););
114 second!($deriv, let f2 = -f1.clone() * &rec;);
115 third!($deriv, let f3 = f2.clone() * rec * F::from(-2.0).unwrap(););
116 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
117 }
118
119 #[inline]
120 fn log(&self, base: F) -> Self {
121 first!($deriv, let rec = self.re.recip(););
122 let f0 = self.re.log(base);
123 first!($deriv, let f1 = rec.clone() / base.ln(););
124 second!($deriv, let f2 = -f1.clone() * &rec;);
125 third!($deriv, let f3 = f2.clone() * rec * F::from(-2.0).unwrap(););
126 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
127 }
128
129 #[inline]
130 fn log2(&self) -> Self {
131 first!($deriv, let rec = self.re.recip(););
132 let f0 = self.re.log2();
133 first!($deriv, let f1 = rec.clone() / (F::one() + F::one()).ln(););
134 second!($deriv, let f2 = -f1.clone() * &rec;);
135 third!($deriv, let f3 = f2.clone() * rec * F::from(-2.0).unwrap(););
136 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
137 }
138
139 #[inline]
140 fn log10(&self) -> Self {
141 first!($deriv, let rec = self.re.recip(););
142 let f0 = self.re.log10();
143 first!($deriv, let f1 = rec.clone() / F::from(10.0).unwrap().ln(););
144 second!($deriv, let f2 = -f1.clone() * &rec;);
145 third!($deriv, let f3 = f2.clone() * rec * F::from(-2.0).unwrap(););
146 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
147 }
148
149 #[inline]
150 fn ln_1p(&self) -> Self {
151 first!($deriv, let rec = (self.re.clone() + F::one()).recip(););
152 let f0 = self.re.ln_1p();
153 first!($deriv, let f1 = rec.clone(););
154 second!($deriv, let f2 = -f1.clone() * &rec;);
155 third!($deriv, let f3 = f2.clone() * rec * F::from(-2.0).unwrap(););
156 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
157 }
158
159 #[inline]
160 fn sin(&self) -> Self {
161 zeroth!($deriv, let s = self.re.sin(););
162 first!($deriv, let (s, c) = self.re.sin_cos(););
163 chain_rule!($deriv, Self::chain_rule(self, s.clone(), c.clone(), -s, -c))
164 }
165
166 #[inline]
167 fn cos(&self) -> Self {
168 zeroth!($deriv, let c = self.re.cos(););
169 first!($deriv, let (s, c) = self.re.sin_cos(););
170 chain_rule!($deriv, Self::chain_rule(self, c.clone(), -s.clone(), -c, s))
171 }
172
173 #[inline]
174 fn sin_cos(&self) -> (Self, Self) {
175 let (s, c) = self.re.sin_cos();
176 (
177 chain_rule!($deriv, Self::chain_rule(self, s.clone(), c.clone(), -s.clone(), -c.clone())),
178 chain_rule!($deriv, Self::chain_rule(self, c.clone(), -s.clone(), -c, s)))
179 }
180
181 #[inline]
182 fn tan(&self) -> Self {
183 let (sin, cos) = self.sin_cos();
184 sin / cos
185 }
186
187 #[inline]
188 fn asin(&self) -> Self {
189 first!($deriv, let rec = (T::one() - self.re.clone() * &self.re).recip(););
190 let f0 = self.re.asin();
191 first!($deriv, let f1 = rec.sqrt(););
192 second!($deriv, let f2 = self.re.clone() * &f1 * &rec;);
193 third!($deriv, let f3 = (self.re.clone() * &self.re * (F::one() + F::one()) + F::one()) * &f1 * &rec * rec;);
194 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
195 }
196
197 #[inline]
198 fn acos(&self) -> Self {
199 first!($deriv, let rec = (T::one() - self.re.clone() * &self.re).recip(););
200 let f0 = self.re.acos();
201 first!($deriv, let f1 = -rec.sqrt(););
202 second!($deriv, let f2 = self.re.clone() * &f1 * &rec;);
203 third!($deriv, let f3 = (self.re.clone() * &self.re * (F::one() + F::one()) + F::one()) * &f1 * &rec * rec;);
204 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
205 }
206
207 #[inline]
208 fn atan(&self) -> Self {
209 first!($deriv, let rec = (T::one() + self.re.clone() * &self.re).recip(););
210 let f0 = self.re.atan();
211 first!($deriv, let f1 = rec.clone(););
212 second!($deriv, let two = F::one() + F::one(););
213 second!($deriv, let f2 = -self.re.clone() * &f1 * &rec * two;);
214 third!($deriv, let f3 = (self.re.clone() * &self.re * F::from(6.0).unwrap() - two) * &f1 * &rec * rec;);
215 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
216 }
217
218 #[inline]
219 fn atan2(&self, other: Self) -> Self {
220 let mut res = (self / other.clone()).atan();
221 res.re = self.re.atan2(other.re);
222 res
223 }
224
225 #[inline]
226 fn sinh(&self) -> Self {
227 let s = self.re.sinh();
228 first!($deriv, let c = self.re.cosh(););
229 chain_rule!($deriv, Self::chain_rule(self, s.clone(), c.clone(), s, c))
230 }
231
232 #[inline]
233 fn cosh(&self) -> Self {
234 first!($deriv, let s = self.re.sinh(););
235 let c = self.re.cosh();
236 chain_rule!($deriv, Self::chain_rule(self, c.clone(), s.clone(), c, s))
237 }
238
239 #[inline]
240 fn tanh(&self) -> Self {
241 self.sinh() / self.cosh()
242 }
243
244 #[inline]
245 fn asinh(&self) -> Self {
246 first!($deriv, let rec = (T::one() + self.re.clone() * &self.re).recip(););
247 let f0 = self.re.asinh();
248 first!($deriv, let f1 = rec.sqrt(););
249 second!($deriv, let f2 = -self.re.clone() * &f1 * &rec;);
250 third!($deriv, let f3 = (self.re.clone() * &self.re * (F::one() + F::one()) - F::one()) * &f1 * &rec * rec;);
251 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
252 }
253
254 #[inline]
255 fn acosh(&self) -> Self {
256 first!($deriv, let rec = (self.re.clone() * &self.re - F::one()).recip(););
257 let f0 = self.re.acosh();
258 first!($deriv, let f1 = rec.sqrt(););
259 second!($deriv, let f2 = -self.re.clone() * &f1 * &rec;);
260 third!($deriv, let f3 = (self.re.clone() * &self.re * (F::one() + F::one()) + F::one()) * &f1 * &rec * rec;);
261 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
262 }
263
264 #[inline]
265 fn atanh(&self) -> Self {
266 first!($deriv, let rec = (T::one() - self.re.clone() * &self.re).recip(););
267 let f0 = self.re.atanh();
268 first!($deriv, let f1 = rec.clone(););
269 second!($deriv, let two = F::one() + F::one(););
270 second!($deriv, let f2 = self.re.clone() * &f1 * &rec * two;);
271 third!($deriv, let f3 = (self.re.clone() * &self.re * F::from(6.0).unwrap() + two) * &f1 * &rec * rec;);
272 chain_rule!($deriv, Self::chain_rule(self, f0, f1, f2, f3))
273 }
274
275 #[inline]
276 fn sph_j0(&self) -> Self {
277 if self.re() < F::epsilon() {
278 Self::one() - self * self / F::from(6.0).unwrap()
279 } else {
280 self.sin() / self
281 }
282 }
283
284 #[inline]
285 fn sph_j1(&self) -> Self {
286 if self.re() < F::epsilon() {
287 self.clone() / F::from(3.0).unwrap()
288 } else {
289 let (s, c) = self.sin_cos();
290 (s - self * c) / (self * self)
291 }
292 }
293
294 #[inline]
295 fn sph_j2(&self) -> Self {
296 if self.re() < F::epsilon() {
297 self * self / F::from(15.0).unwrap()
298 } else {
299 let (s, c) = self.sin_cos();
300 let s2 = self * self;
301 ((&s - self * c) * F::from(3.0).unwrap() - &s2 * s) / (s2 * self)
302 }
303 }
304 }
305 };
306}
307
308#[macro_export]
309macro_rules! zeroth {
310 (zeroth, $($code:tt)*) => {
311 $($code)*
312 };
313 (first, $($code:tt)*) => {};
314 (second, $($code:tt)*) => {};
315 (third, $($code:tt)*) => {};
316}
317
318#[macro_export]
319macro_rules! first {
320 (zeroth, $($code:tt)*) => {};
321 (first, $($code:tt)*) => {
322 $($code)*
323 };
324 (second, $($code:tt)*) => {
325 $($code)*
326 };
327 (third, $($code:tt)*) => {
328 $($code)*
329 };
330}
331
332#[macro_export]
333macro_rules! second {
334 (zeroth, $($code:tt)*) => {};
335 (first, $($code:tt)*) => {};
336 (second, $($code:tt)*) => {
337 $($code)*
338 };
339 (third, $($code:tt)*) => {
340 $($code)*
341 };
342}
343
344#[macro_export]
345macro_rules! third {
346 (zeroth, $($code:tt)*) => {};
347 (first, $($code:tt)*) => {};
348 (second, $($code:tt)*) => {};
349 (third, $($code:tt)*) => {
350 $($code)*
351 };
352}
353
354#[macro_export]
355macro_rules! chain_rule {
356 (zeroth, Self::chain_rule($self:ident, $f0:expr, $f1:expr, $f2:expr, $f3:expr)) => {
357 Self::chain_rule($self, $f0)
358 };
359 (first, Self::chain_rule($self:ident, $f0:expr, $f1:expr, $f2:expr, $f3:expr)) => {
360 Self::chain_rule($self, $f0, $f1)
361 };
362 (second, Self::chain_rule($self:ident, $f0:expr, $f1:expr, $f2:expr, $f3:expr)) => {
363 Self::chain_rule($self, $f0, $f1, $f2)
364 };
365 (third, Self::chain_rule($self:ident, $f0:expr, $f1:expr, $f2:expr, $f3:expr)) => {
366 Self::chain_rule($self, $f0, $f1, $f2, $f3)
367 };
368}
369
370#[macro_export]
371macro_rules! impl_zeroth_derivatives {
372 ($struct:ident, [$($im:ident),*]$(, [$($dim:tt),*]$(, [$($ddim:tt),*])*)?) => {
373 impl_derivatives!(zeroth, 0, $struct, [$($im),*]$(, [$($dim),*]$(, [$($ddim),*])*)?);
374 };
375}
376
377#[macro_export]
378macro_rules! impl_first_derivatives {
379 ($struct:ident, [$($im:ident),*]$(, [$($dim:tt),*]$(, [$($ddim:tt),*])*)?) => {
380 impl_derivatives!(first, 1, $struct, [$($im),*]$(, [$($dim),*]$(, [$($ddim),*])*)?);
381 };
382}
383
384#[macro_export]
385macro_rules! impl_second_derivatives {
386 ($struct:ident, [$($im:ident),*]$(, [$($dim:tt),*]$(, [$($ddim:tt),*])*)?) => {
387 impl_derivatives!(second, 2, $struct, [$($im),*]$(, [$($dim),*]$(, [$($ddim),*])*)?);
388 };
389}
390
391#[macro_export]
392macro_rules! impl_third_derivatives {
393 ($struct:ident, [$($im:ident),*]$(, [$($dim:tt),*]$(, [$($ddim:tt),*])*)?) => {
394 impl_derivatives!(third, 3, $struct, [$($im),*]$(, [$($dim),*]$(, [$($ddim),*])*)?);
395 };
396}