1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
//! Scalar types that every numerical kernel is generic over.
//!
//! The numerical kernels (the Mie solution and the special functions) are generic over
//! [`Real`], so one implementation gives:
//!
//! * values, with `f64`;
//! * extended-precision values, with [`DoubleDouble`] (about 32 significant digits);
//! * values with N first derivatives, with [`Dual<f64, N>`](Dual);
//! * extended-precision derivatives, with `Dual<DoubleDouble, N>`;
//! * second derivatives, with nested duals `Dual<Dual<T, M>, N>`.
//!
//! ```
//! use scatrs::Real;
//! use scatrs::scalar::{DoubleDouble, Dual};
//!
//! // A generic function, evaluated in double-double and with a derivative.
//! fn f<T: Real>(x: T) -> T {
//! (x * 2.0).sin() * x.exp()
//! }
//! let value = f(DoubleDouble::from_f64(0.5));
//! assert!((value.approx_f64() - (1.0_f64).sin() * (0.5_f64).exp()).abs() < 1e-15);
//!
//! let [x] = Dual::<f64, 1>::variables([0.5]);
//! let y = f(x); // y.eps()[0] = d/dx [sin 2x · eˣ] = (2 cos 2x + sin 2x) eˣ
//! let exact = (2.0 * (1.0_f64).cos() + (1.0_f64).sin()) * (0.5_f64).exp();
//! assert!((y.eps()[0] - exact).abs() < 1e-14);
//! ```
//!
//! [`Real`] lists what generic code may rely on, and [`DoubleDouble`] and [`Dual`] state their
//! accuracy and conventions.
use Debug;
use ;
use ;
pub use ;
pub use Dual;
/// A real scalar type for generic kernels.
///
/// Implemented for `f64`, `DoubleDouble` (about 32 significant digits) and `Dual<T, N>`
/// (forward-mode derivatives). The `num_traits::Float` supertrait makes
/// `num_complex::Complex<T>` usable for every `T: Real`.
/// Mixed operations with `f64` take the `f64` on the right (`x * 2.0`, not `2.0 * x`);
/// for `Complex<T>` use `z.scale(T::from_f64(2.0))`.
///
/// # Notes for generic code
///
/// `DoubleDouble` follows `f64` closely, but generic code must not assume more than `Float`
/// promises (see [`DoubleDouble`]'s "Use as `num_traits::Float`"):
///
/// * `mul_add` is fused for `f64` but not for `DoubleDouble`, which rounds the product before
/// adding. Error-free transformations must not rely on it: Kahan's idiom
/// `e = b.mul_add(c, -(b * c))` gives exactly 0 for `DoubleDouble`, not the rounding error
/// of `b * c`.
/// * `epsilon()` is a unit (2⁻⁵² for `f64`, 2⁻¹⁰⁴ for `DoubleDouble`), not a bound on the
/// error of one operation: a `DoubleDouble` division errs by up to about 2.5·epsilon.
/// * `NumCast::from` another `Real` keeps only its leading word (`to_f64`); use the value
/// itself, or [`from_f64_pair`](Real::from_f64_pair), to keep full precision.
/// * For [`Dual`], `==` and `<` compare values only, while `is_zero`, `is_one`, `is_nan` and
/// `is_finite` also inspect the derivatives (see its "Comparisons and zero tests"): branch on
/// comparisons of values, as `x == T::zero()`, where the derivatives should not matter.
/// * `Dual<T, N>` takes the values of `sin` and `cos` from `T::sin_cos`. They equal `T::sin` and
/// `T::cos` where `sin_cos` does: always for `DoubleDouble`. For `f64` on Apple targets, LLVM
/// fuses `sin` and `cos` of one argument into one call whose sine can be an ulp off a lone
/// `sin`, in `Dual` as in any `f64` code computing both.