frechet/
lib.rs

1//! `frechet` is a library providing simple dual number types with interfaces
2//! similar to that of the standard floating point types. Additionally, provides
3//! a small interface to abstract the computation of derivatives.
4//!
5//! # Refresher on dual numbers
6//!
7//! A dual number $z$ is represented as the sum of a "real part" $x$ and an "imaginary part" $y$,
8//! and written $z = x + yj$, where $j$ is purely symbolic and follows the rule $j^2 = 0$.
9//! As an example, evaluating the polynomial $P(X)=3X^2-X+1$ at the dual number $X + j$ yields
10//! $$P(X+j)=3X^2-X+1 + (6X-1)j = P(X) + P^\prime(X)j.$$
11//! This motivates the following extension of any (differentiable) real function $f$ to dual numbers:
12//! $$f(x+yj) = f(x) + yf^\prime(x)j.$$
13//!
14//! # Example
15//! ```
16//! use frechet::*;
17//!
18//! fn p(x: dual32) -> dual32 { x.powf(2.5).atanh() + 1.0  }
19//! fn p_derivative(x: f32) -> f32 { -2.5 * x.powf(1.5)/(x.powi(5) - 1.0) }
20//!
21//! // using the `derivative` function
22//! let z1 = derivative(p, 2.0);
23//!
24//! // manually
25//! let z2 = p(2.0.as_dual_variable()).d;
26//!
27//! // exact derivative
28//! let z3 = p_derivative(2.0);
29//!
30//! assert!((z1 - z3).abs() < f32::EPSILON);
31//! assert!((z2 - z3).abs() < f32::EPSILON);
32//! ```
33
34/// Type representing a simple 32-bit precision dual number.
35#[allow(non_camel_case_types)]
36#[derive(Clone, Copy, Debug)]
37pub struct dual32 {
38    pub x: f32,
39    pub d: f32,
40}
41
42impl dual32 {
43    pub const ZERO: dual32 = dual32 { x: 0.0, d: 0.0 };
44    pub const ONE: dual32 = dual32 { x: 1.0, d: 0.0 };
45
46    /// The unit "imaginary" dual number. Satisfies  `J*J == ZERO`.
47    pub const J: dual32 = dual32 { x: 0.0, d: 1.0 };
48
49    #[inline]
50    pub const fn new(x: f32, d: f32) -> dual32 {
51        Self { x, d }
52    }
53}
54
55impl core::ops::Add<dual32> for dual32 {
56    type Output = dual32;
57
58    #[inline]
59    fn add(self, rhs: dual32) -> Self::Output {
60        Self::Output::new(self.x + rhs.x, self.d + rhs.d)
61    }
62}
63
64impl core::ops::Add<f32> for dual32 {
65    type Output = dual32;
66
67    #[inline]
68    fn add(self, rhs: f32) -> Self::Output {
69        Self::Output::new(self.x + rhs, self.d)
70    }
71}
72
73impl core::ops::Add<dual32> for f32 {
74    type Output = dual32;
75
76    #[inline]
77    fn add(self, rhs: dual32) -> Self::Output {
78        Self::Output::new(self + rhs.x, rhs.d)
79    }
80}
81
82impl core::ops::Sub<dual32> for dual32 {
83    type Output = dual32;
84
85    #[inline]
86    fn sub(self, rhs: dual32) -> Self::Output {
87        Self::Output::new(self.x - rhs.x, self.d - rhs.d)
88    }
89}
90
91impl core::ops::Sub<f32> for dual32 {
92    type Output = dual32;
93
94    #[inline]
95    fn sub(self, rhs: f32) -> Self::Output {
96        Self::Output::new(self.x - rhs, self.d)
97    }
98}
99
100impl core::ops::Sub<dual32> for f32 {
101    type Output = dual32;
102
103    #[inline]
104    fn sub(self, rhs: dual32) -> Self::Output {
105        Self::Output::new(self - rhs.x, rhs.d)
106    }
107}
108
109impl core::ops::Mul<dual32> for dual32 {
110    type Output = dual32;
111
112    #[inline]
113    fn mul(self, rhs: dual32) -> Self::Output {
114        Self::Output::new(self.x * rhs.x, self.x * rhs.d + self.d * rhs.x)
115    }
116}
117
118impl core::ops::Mul<f32> for dual32 {
119    type Output = dual32;
120
121    #[inline]
122    fn mul(self, rhs: f32) -> Self::Output {
123        Self::Output::new(self.x * rhs, self.x * rhs)
124    }
125}
126
127impl core::ops::Mul<dual32> for f32 {
128    type Output = dual32;
129
130    #[inline]
131    fn mul(self, rhs: dual32) -> Self::Output {
132        Self::Output::new(self * rhs.x, self * rhs.d)
133    }
134}
135
136impl core::ops::Div<dual32> for dual32 {
137    type Output = dual32;
138
139    #[inline]
140    fn div(self, rhs: dual32) -> Self::Output {
141        Self::Output::new(
142            self.x / rhs.x,
143            (self.d * rhs.x - self.x * rhs.d) / (rhs.x * rhs.x),
144        )
145    }
146}
147
148impl core::ops::Div<f32> for dual32 {
149    type Output = dual32;
150
151    #[inline]
152    fn div(self, rhs: f32) -> Self::Output {
153        Self::Output::new(self.x / rhs, self.x / rhs)
154    }
155}
156
157impl core::ops::Div<dual32> for f32 {
158    type Output = dual32;
159
160    #[inline]
161    fn div(self, rhs: dual32) -> Self::Output {
162        Self::Output::new(self / rhs.x, (-self * rhs.d) / (rhs.x * rhs.x))
163    }
164}
165
166// Power functions
167impl dual32 {
168    #[inline]
169    pub fn recip(self) -> dual32 {
170        Self::new(self.x.recip(), -self.d * self.x.recip().powi(2))
171    }
172
173    #[inline]
174    pub fn powf(self, n: f32) -> dual32 {
175        Self::new(self.x.powf(n), self.d * n * self.x.powf(n - 1.0))
176    }
177
178    #[inline]
179    pub fn powi(self, n: i32) -> dual32 {
180        Self::new(self.x.powi(n), self.d * n as f32 * self.x.powi(n - 1))
181    }
182
183    #[inline]
184    pub fn sqrt(self) -> dual32 {
185        Self::new(self.x.sqrt(), self.d * 0.5 / self.x.sqrt())
186    }
187
188    #[inline]
189    pub fn cbrt(self) -> dual32 {
190        Self::new(self.x.cbrt(), self.d * (1.0 / 3.0) / self.x.cbrt().powi(2))
191    }
192}
193
194// Exponentials and logarithms
195impl dual32 {
196    #[inline]
197    pub fn exp(self) -> dual32 {
198        Self::new(self.x.exp(), self.d * self.x.exp())
199    }
200
201    #[inline]
202    pub fn exp2(self) -> dual32 {
203        Self::new(
204            self.x.exp2(),
205            self.d * core::f32::consts::LN_2 * self.x.exp2(),
206        )
207    }
208
209    #[inline]
210    pub fn exp_m1(self) -> dual32 {
211        Self::new(self.x.exp_m1(), self.d * self.x.exp())
212    }
213
214    #[inline]
215    pub fn ln(self) -> dual32 {
216        Self::new(self.x.ln(), self.d / self.x)
217    }
218
219    #[inline]
220    pub fn ln_1p(self) -> dual32 {
221        Self::new(self.x.ln_1p(), self.d / (1.0 + self.x))
222    }
223
224    #[inline]
225    pub fn log(self, base: f32) -> dual32 {
226        Self::new(self.x.log(base), self.d / (base.ln() * self.x))
227    }
228
229    #[inline]
230    pub fn log10(self) -> dual32 {
231        Self::new(self.x.log10(), self.d / (core::f32::consts::LN_10 * self.x))
232    }
233
234    #[inline]
235    pub fn log2(self) -> dual32 {
236        Self::new(self.x.log2(), self.d / (core::f32::consts::LN_2 * self.x))
237    }
238}
239
240// Trigonometric functions
241// atan2 ?
242// sin_cos ?
243impl dual32 {
244    #[inline]
245    pub fn cos(self) -> dual32 {
246        Self::new(self.x.cos(), -self.d * self.x.sin())
247    }
248
249    #[inline]
250    pub fn sin(self) -> dual32 {
251        Self::new(self.x.sin(), self.d * self.x.cos())
252    }
253
254    #[inline]
255    pub fn tan(self) -> dual32 {
256        Self::new(self.x.tan(), self.d / self.x.cos().powi(2))
257    }
258
259    #[inline]
260    pub fn acos(self) -> dual32 {
261        Self::new(self.x.acos(), -self.d / (1.0 - self.x * self.x).sqrt())
262    }
263
264    #[inline]
265    pub fn asin(self) -> dual32 {
266        Self::new(self.x.asin(), self.d / (1.0 - self.x * self.x).sqrt())
267    }
268
269    #[inline]
270    pub fn atan(self) -> dual32 {
271        Self::new(self.x.tan(), self.d / (1.0 + self.x * self.x))
272    }
273}
274
275// Hyperbolic functions
276impl dual32 {
277    #[inline]
278    pub fn cosh(self) -> dual32 {
279        Self::new(self.x.cosh(), self.d * self.x.sinh())
280    }
281
282    #[inline]
283    pub fn sinh(self) -> dual32 {
284        Self::new(self.x.sinh(), self.d * self.x.cosh())
285    }
286
287    #[inline]
288    pub fn tanh(self) -> dual32 {
289        Self::new(self.x.tanh(), self.d / self.x.cosh().powi(2))
290    }
291
292    #[inline]
293    pub fn acosh(self) -> dual32 {
294        Self::new(self.x.acosh(), self.d / (self.x * self.x - 1.0).sqrt())
295    }
296
297    #[inline]
298    pub fn asinh(self) -> dual32 {
299        Self::new(self.x.asinh(), self.d / (self.x * self.x + 1.0).sqrt())
300    }
301
302    #[inline]
303    pub fn atanh(self) -> dual32 {
304        Self::new(self.x.atanh(), self.d / (1.0 - self.x * self.x))
305    }
306}
307
308// Continuous differentiable almost-everywhere functions
309impl dual32 {
310    #[inline]
311    pub fn abs(self) -> dual32 {
312        Self::new(self.x.abs(), self.d * self.x.signum())
313    }
314
315    #[inline]
316    pub fn min(self, other: f32) -> dual32 {
317        Self::new(
318            self.x.min(other),
319            if self.x < other { self.d * 1.0 } else { 0.0 },
320        )
321    }
322
323    #[inline]
324    pub fn max(self, other: f32) -> dual32 {
325        Self::new(
326            self.x.max(other),
327            if other < self.x { self.d * 1.0 } else { 0.0 },
328        )
329    }
330
331    #[inline]
332    pub fn clamp(self, min: f32, max: f32) -> dual32 {
333        Self::new(
334            self.x.clamp(min, max),
335            if min < self.x && self.x < max {
336                self.d * 1.0
337            } else {
338                0.0
339            },
340        )
341    }
342}
343
344/// Piece-wise constant functions
345impl dual32 {
346    #[inline]
347    pub fn ceil(self) -> dual32 {
348        Self::new(self.x.ceil(), 0.0)
349    }
350
351    #[inline]
352    pub fn floor(self) -> dual32 {
353        Self::new(self.x.floor(), 0.0)
354    }
355
356    // not constant but x.fract() == x - x.floor()
357    #[inline]
358    pub fn fract(self) -> dual32 {
359        Self::new(self.x.fract(), self.d)
360    }
361
362    #[inline]
363    pub fn round(self) -> dual32 {
364        Self::new(self.x.round(), 0.0)
365    }
366
367    #[inline]
368    pub fn signum(self) -> dual32 {
369        Self::new(self.x.signum(), 0.0)
370    }
371
372    #[inline]
373    pub fn trunc(self) -> dual32 {
374        Self::new(self.x.trunc(), 0.0)
375    }
376}
377
378/// Sporadic functions
379impl dual32 {
380    #[inline]
381    pub fn mul_add(self, a: f32, b: f32) -> dual32 {
382        Self::new(self.x.mul_add(a, b), self.d * a * self.x)
383    }
384}
385
386/// Type representing a simple 64-bit precision dual number.
387#[allow(non_camel_case_types)]
388#[derive(Clone, Copy, Debug)]
389pub struct dual64 {
390    pub x: f64,
391    pub d: f64,
392}
393
394impl dual64 {
395    pub const ZERO: dual64 = dual64 { x: 0.0, d: 0.0 };
396    pub const ONE: dual64 = dual64 { x: 1.0, d: 0.0 };
397
398    /// The unit "imaginary" dual number. Satisfies  `J*J == ZERO`.
399    pub const J: dual64 = dual64 { x: 0.0, d: 1.0 };
400
401    #[inline]
402    pub const fn new(x: f64, d: f64) -> dual64 {
403        Self { x, d }
404    }
405}
406
407impl core::ops::Add<dual64> for dual64 {
408    type Output = dual64;
409
410    #[inline]
411    fn add(self, rhs: dual64) -> Self::Output {
412        Self::Output::new(self.x + rhs.x, self.d + rhs.d)
413    }
414}
415
416impl core::ops::Add<f64> for dual64 {
417    type Output = dual64;
418
419    #[inline]
420    fn add(self, rhs: f64) -> Self::Output {
421        Self::Output::new(self.x + rhs, self.d)
422    }
423}
424
425impl core::ops::Add<dual64> for f64 {
426    type Output = dual64;
427
428    #[inline]
429    fn add(self, rhs: dual64) -> Self::Output {
430        Self::Output::new(self + rhs.x, rhs.d)
431    }
432}
433
434impl core::ops::Sub<dual64> for dual64 {
435    type Output = dual64;
436
437    #[inline]
438    fn sub(self, rhs: dual64) -> Self::Output {
439        Self::Output::new(self.x - rhs.x, self.d - rhs.d)
440    }
441}
442
443impl core::ops::Sub<f64> for dual64 {
444    type Output = dual64;
445
446    #[inline]
447    fn sub(self, rhs: f64) -> Self::Output {
448        Self::Output::new(self.x - rhs, self.d)
449    }
450}
451
452impl core::ops::Sub<dual64> for f64 {
453    type Output = dual64;
454
455    #[inline]
456    fn sub(self, rhs: dual64) -> Self::Output {
457        Self::Output::new(self - rhs.x, rhs.d)
458    }
459}
460
461impl core::ops::Mul<dual64> for dual64 {
462    type Output = dual64;
463
464    #[inline]
465    fn mul(self, rhs: dual64) -> Self::Output {
466        Self::Output::new(self.x * rhs.x, self.x * rhs.d + self.d * rhs.x)
467    }
468}
469
470impl core::ops::Mul<f64> for dual64 {
471    type Output = dual64;
472
473    #[inline]
474    fn mul(self, rhs: f64) -> Self::Output {
475        Self::Output::new(self.x * rhs, self.x * rhs)
476    }
477}
478
479impl core::ops::Mul<dual64> for f64 {
480    type Output = dual64;
481
482    #[inline]
483    fn mul(self, rhs: dual64) -> Self::Output {
484        Self::Output::new(self * rhs.x, self * rhs.d)
485    }
486}
487
488impl core::ops::Div<dual64> for dual64 {
489    type Output = dual64;
490
491    #[inline]
492    fn div(self, rhs: dual64) -> Self::Output {
493        Self::Output::new(
494            self.x / rhs.x,
495            (self.d * rhs.x - self.x * rhs.d) / (rhs.x * rhs.x),
496        )
497    }
498}
499
500impl core::ops::Div<f64> for dual64 {
501    type Output = dual64;
502
503    #[inline]
504    fn div(self, rhs: f64) -> Self::Output {
505        Self::Output::new(self.x / rhs, self.x / rhs)
506    }
507}
508
509impl core::ops::Div<dual64> for f64 {
510    type Output = dual64;
511
512    #[inline]
513    fn div(self, rhs: dual64) -> Self::Output {
514        Self::Output::new(self / rhs.x, (-self * rhs.d) / (rhs.x * rhs.x))
515    }
516}
517
518// Power functions
519impl dual64 {
520    #[inline]
521    pub fn recip(self) -> dual64 {
522        Self::new(self.x.recip(), -self.d * self.x.recip().powi(2))
523    }
524
525    #[inline]
526    pub fn powf(self, n: f64) -> dual64 {
527        Self::new(self.x.powf(n), self.d * n * self.x.powf(n - 1.0))
528    }
529
530    #[inline]
531    pub fn powi(self, n: i32) -> dual64 {
532        Self::new(self.x.powi(n), self.d * n as f64 * self.x.powi(n - 1))
533    }
534
535    #[inline]
536    pub fn sqrt(self) -> dual64 {
537        Self::new(self.x.sqrt(), self.d * 0.5 / self.x.sqrt())
538    }
539
540    #[inline]
541    pub fn cbrt(self) -> dual64 {
542        Self::new(self.x.cbrt(), self.d * (1.0 / 3.0) / self.x.cbrt().powi(2))
543    }
544}
545
546// Exponentials and logarithms
547impl dual64 {
548    #[inline]
549    pub fn exp(self) -> dual64 {
550        Self::new(self.x.exp(), self.d * self.x.exp())
551    }
552
553    #[inline]
554    pub fn exp2(self) -> dual64 {
555        Self::new(
556            self.x.exp2(),
557            self.d * core::f64::consts::LN_2 * self.x.exp2(),
558        )
559    }
560
561    #[inline]
562    pub fn exp_m1(self) -> dual64 {
563        Self::new(self.x.exp_m1(), self.d * self.x.exp())
564    }
565
566    #[inline]
567    pub fn ln(self) -> dual64 {
568        Self::new(self.x.ln(), self.d / self.x)
569    }
570
571    #[inline]
572    pub fn ln_1p(self) -> dual64 {
573        Self::new(self.x.ln_1p(), self.d / (1.0 + self.x))
574    }
575
576    #[inline]
577    pub fn log(self, base: f64) -> dual64 {
578        Self::new(self.x.log(base), self.d / (base.ln() * self.x))
579    }
580
581    #[inline]
582    pub fn log10(self) -> dual64 {
583        Self::new(self.x.log10(), self.d / (core::f64::consts::LN_10 * self.x))
584    }
585
586    #[inline]
587    pub fn log2(self) -> dual64 {
588        Self::new(self.x.log2(), self.d / (core::f64::consts::LN_2 * self.x))
589    }
590}
591
592// Trigonometric functions
593// atan2 ?
594// sin_cos ?
595impl dual64 {
596    #[inline]
597    pub fn cos(self) -> dual64 {
598        Self::new(self.x.cos(), -self.d * self.x.sin())
599    }
600
601    #[inline]
602    pub fn sin(self) -> dual64 {
603        Self::new(self.x.sin(), self.d * self.x.cos())
604    }
605
606    #[inline]
607    pub fn tan(self) -> dual64 {
608        Self::new(self.x.tan(), self.d / self.x.cos().powi(2))
609    }
610
611    #[inline]
612    pub fn acos(self) -> dual64 {
613        Self::new(self.x.acos(), -self.d / (1.0 - self.x * self.x).sqrt())
614    }
615
616    #[inline]
617    pub fn asin(self) -> dual64 {
618        Self::new(self.x.asin(), self.d / (1.0 - self.x * self.x).sqrt())
619    }
620
621    #[inline]
622    pub fn atan(self) -> dual64 {
623        Self::new(self.x.tan(), self.d / (1.0 + self.x * self.x))
624    }
625}
626
627// Hyperbolic functions
628impl dual64 {
629    #[inline]
630    pub fn cosh(self) -> dual64 {
631        Self::new(self.x.cosh(), self.d * self.x.sinh())
632    }
633
634    #[inline]
635    pub fn sinh(self) -> dual64 {
636        Self::new(self.x.sinh(), self.d * self.x.cosh())
637    }
638
639    #[inline]
640    pub fn tanh(self) -> dual64 {
641        Self::new(self.x.tanh(), self.d / self.x.cosh().powi(2))
642    }
643
644    #[inline]
645    pub fn acosh(self) -> dual64 {
646        Self::new(self.x.acosh(), self.d / (self.x * self.x - 1.0).sqrt())
647    }
648
649    #[inline]
650    pub fn asinh(self) -> dual64 {
651        Self::new(self.x.asinh(), self.d / (self.x * self.x + 1.0).sqrt())
652    }
653
654    #[inline]
655    pub fn atanh(self) -> dual64 {
656        Self::new(self.x.atanh(), self.d / (1.0 - self.x * self.x))
657    }
658}
659
660// Continuous differentiable almost-everywhere functions
661impl dual64 {
662    #[inline]
663    pub fn abs(self) -> dual64 {
664        Self::new(self.x.abs(), self.d * self.x.signum())
665    }
666
667    #[inline]
668    pub fn min(self, other: f64) -> dual64 {
669        Self::new(
670            self.x.min(other),
671            if self.x < other { self.d * 1.0 } else { 0.0 },
672        )
673    }
674
675    #[inline]
676    pub fn max(self, other: f64) -> dual64 {
677        Self::new(
678            self.x.max(other),
679            if other < self.x { self.d * 1.0 } else { 0.0 },
680        )
681    }
682
683    #[inline]
684    pub fn clamp(self, min: f64, max: f64) -> dual64 {
685        Self::new(
686            self.x.clamp(min, max),
687            if min < self.x && self.x < max {
688                self.d * 1.0
689            } else {
690                0.0
691            },
692        )
693    }
694}
695
696/// Piece-wise constant functions
697impl dual64 {
698    #[inline]
699    pub fn ceil(self) -> dual64 {
700        Self::new(self.x.ceil(), 0.0)
701    }
702
703    #[inline]
704    pub fn floor(self) -> dual64 {
705        Self::new(self.x.floor(), 0.0)
706    }
707
708    // not constant but x.fract() == x - x.floor()
709    #[inline]
710    pub fn fract(self) -> dual64 {
711        Self::new(self.x.fract(), self.d)
712    }
713
714    #[inline]
715    pub fn round(self) -> dual64 {
716        Self::new(self.x.round(), 0.0)
717    }
718
719    #[inline]
720    pub fn signum(self) -> dual64 {
721        Self::new(self.x.signum(), 0.0)
722    }
723
724    #[inline]
725    pub fn trunc(self) -> dual64 {
726        Self::new(self.x.trunc(), 0.0)
727    }
728}
729
730/// Sporadic functions
731impl dual64 {
732    #[inline]
733    pub fn mul_add(self, a: f64, b: f64) -> dual64 {
734        Self::new(self.x.mul_add(a, b), self.d * a * self.x)
735    }
736}
737
738#[cfg(test)]
739mod tests {
740    use super::*;
741
742    #[test]
743    fn identity32() {
744        let dual32 { x, d } = dual32::J.powi(2);
745        assert!(x.abs() < f32::EPSILON);
746        assert!(d.abs() < f32::EPSILON);
747    }
748
749    #[test]
750    fn identity64() {
751        let dual64 { x, d } = dual64::J.powi(2);
752        assert!(x.abs() < f64::EPSILON);
753        assert!(d.abs() < f64::EPSILON);
754    }
755
756    #[test]
757    fn polynomial32() {
758        fn p(x: dual32) -> dual32 {
759            4.0 * x * x - 3.0 * x + 3.0
760        }
761        let x = 3.0;
762        let dual32 { x: px, d: pdx } = p(x + dual32::J);
763
764        assert!((px - 30.0).abs() <= f32::EPSILON);
765        assert!((pdx - 21.0).abs() <= f32::EPSILON);
766    }
767
768    #[test]
769    fn polynomial64() {
770        fn p(x: dual64) -> dual64 {
771            4.0 * x * x - 3.0 * x + 3.0
772        }
773        let x = 3.0;
774        let dual64 { x: px, d: pdx } = p(x + dual64::J);
775
776        assert!((px - 30.0).abs() <= f64::EPSILON);
777        assert!((pdx - 21.0).abs() <= f64::EPSILON);
778    }
779}
780
781// --------------------
782
783/// Trait providing generic getters for dual32 and dual64
784pub trait Dual<F> {
785    fn real(self) -> F;
786    fn imag(self) -> F;
787}
788
789impl Dual<f32> for dual32 {
790    fn real(self) -> f32 {
791        self.x
792    }
793
794    fn imag(self) -> f32 {
795        self.d
796    }
797}
798
799impl Dual<f64> for dual64 {
800    fn real(self) -> f64 {
801        self.x
802    }
803
804    fn imag(self) -> f64 {
805        self.d
806    }
807}
808
809/// Trait used to convert a non-dual numeric type into a dual "variable",
810/// i.e. a dual number with unit imaginary part.
811pub trait AsDualVariable<Dual> {
812    fn as_dual_variable(self) -> Dual;
813    fn from_dual(z: Dual) -> Self;
814}
815
816impl AsDualVariable<dual32> for f32 {
817    fn as_dual_variable(self) -> dual32 {
818        dual32::new(self, 1.0)
819    }
820
821    fn from_dual(z: dual32) -> Self {
822        z.x
823    }
824}
825
826impl AsDualVariable<dual32> for dual32 {
827    fn as_dual_variable(self) -> dual32 {
828        self
829    }
830
831    fn from_dual(z: dual32) -> Self {
832        z
833    }
834}
835
836impl AsDualVariable<dual64> for f64 {
837    fn as_dual_variable(self) -> dual64 {
838        dual64::new(self, 1.0)
839    }
840
841    fn from_dual(z: dual64) -> Self {
842        z.x
843    }
844}
845
846impl AsDualVariable<dual64> for dual64 {
847    fn as_dual_variable(self) -> dual64 {
848        self
849    }
850
851    fn from_dual(z: dual64) -> Self {
852        z
853    }
854}
855
856/// Computes the derivative of a function `f` at a point `x`
857pub fn derivative<F, D>(f: impl Fn(D) -> D, x: F) -> F
858where
859    D: Dual<F>,
860    F: AsDualVariable<D>,
861{
862    f(x.as_dual_variable()).imag()
863}
frechet/lib.rs

frechet/
lib.rs
