Merge release/0.2.0

Scalar is now num-dual's trait, blanket-implemented, rather than a
hand-rolled one declaring the same methods over a fixed list of types.

Breaking: the `autodiff` feature is gone, since num-dual is no longer
optional; `Scalar::re`, `exp`, `ln`, `log10`, `sqrt`, `sin`, `cos`,
`powi`, `powf`, `abs` and `powt` are gone, replaced by DualNum,
DualStruct and Signed, with `powt` now `powd`; and a call built only
from float literals must name its scalar, as gtd7d::<f64>(..).

Not breaking, but the reason it is worth doing: a caller's own dual type
is now covered automatically, and the dependency tree drops from about
twenty crates to two.

Cargo.lock

@@ -325,11 +325,10 @@ dependencies = [
 
 [[package]]
 name = "nrlmsise00"
-version = "0.1.0"
+version = "0.2.0"
 dependencies = [
  "approx",
  "criterion",
- "nalgebra",
  "num-dual",
  "num-traits",
  "serde_json",

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "nrlmsise00"
-version = "0.1.0"
+version = "0.2.0"
 edition = "2021"
 authors = ["Connor Johnstone"]
 description = "NRLMSISE-00 atmosphere model, generic over its scalar type so density can be differentiated"
@@ -25,9 +25,7 @@ all-features = true
 members = ["nrlmsis2"]
 
 [features]
-default = ["autodiff", "drivers"]
+default = ["drivers"]
-# Dual-number scalars, so density can be differentiated.
-autodiff = ["dep:num-dual", "dep:nalgebra"]
 # The compiled-in solar and geomagnetic record, and date lookups into it.
 drivers = []
 # Trades the reference's exact bits for speed: the angle-sum identity for the expansions'
@@ -36,13 +34,14 @@ drivers = []
 fast = []
 
 [dependencies]
-nalgebra = { version = "0.35", optional = true }
+num-dual = { version = "0.15", default-features = false }
-num-dual = { version = "0.15", optional = true }
 num-traits = { version = "0.2", default-features = false }
 
 [dev-dependencies]
 approx = "0.5"
 criterion = "0.7"
+# With default features, for the benchmark: it measures the vector dual types, which are
+# the ones behind num-dual's nalgebra feature.
 num-dual = "0.15"
 # arbitrary_precision keeps each number's own digits, so a fixture value is parsed by
 # str::parse rather than by serde_json's own float parser, which reads some f64s back one
@@ -53,10 +52,10 @@ serde_json = { version = "1.0", features = ["arbitrary_precision"] }
 [[bench]]
 name = "density"
 harness = false
-# Measures the dual-number scalars and the driver lookup, so it needs both.
+# Measures the driver lookup as well as the model.
-required-features = ["autodiff", "drivers"]
+required-features = ["drivers"]
 
 [[example]]
 name = "density"
-# Reads the compiled-in driver record and seeds a dual number, so it needs both.
+# Reads the compiled-in driver record.
-required-features = ["drivers", "autodiff"]
+required-features = ["drivers"]

readme.md

@@ -12,7 +12,8 @@ implementations exist and both are `f64` only.
 use nrlmsise00::{Input, MagneticActivity, Model};
 
 let model = Model::new();
-let output = model.gtd7d(&Input::consistent(
+// The scalar is named, because a call built only from literals has nothing else to fix it.
+let output = model.gtd7d::<f64>(&Input::consistent(
     172,                                   // day of year
     29000.0,                               // seconds of universal time
     400.0,                                 // altitude, km
@@ -32,8 +33,6 @@ assert!((output.temperature - 1098.2).abs() < 0.1);
 And the derivative, through the model rather than around it:
 
 ```rust
-# #[cfg(feature = "autodiff")]
-# fn main() {
 use num_dual::{Dual64, DualNum};
 use nrlmsise00::{Input, MagneticActivity, Model, Scalar};
 
@@ -53,15 +52,49 @@ let output = model.gtd7d(&Input::consistent(
 // d(rho)/dz, in kg/m^3 per km, and the scale height it implies.
 let scale_height = -output.total_mass.re / output.total_mass.eps;
 assert!((scale_height - 55.8).abs() < 0.1);
-# }
-# #[cfg(not(feature = "autodiff"))]
-# fn main() {}
 ```
 
 `cargo run --example density` prints a vertical profile, the same altitude across two
 solar cycles, and the altitude derivative a drag calculation needs, checked against a
 central difference.
 
+## The scalar type is `num-dual`'s
+
+`Scalar` is a marker over [`num_dual::DualNum`] with one method of its own, `lift`, and it
+is blanket-implemented:
+
+```rust,ignore
+pub trait Scalar: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy {
+    fn lift(x: f64) -> Self { Self::from(x) }
+}
+
+impl<T: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy> Scalar for T {}
+```
+
+So every type that implements `num-dual`'s trait at `f64` precision is a scalar here —
+`f64`, each of its dual numbers, nestings of them, and anything a caller writes. There is
+nothing to implement and nothing to enumerate, which is the point: an earlier version
+declared `exp`, `ln`, `powf` and the rest itself and implemented them for a fixed list of
+types, so a caller whose scalar was a dual type the list had missed was stopped by the
+orphan rule with no way out, and a caller generic over another crate's scalar trait of the
+same shape was adapting between two names for one thing.
+
+Two things follow for a caller:
+
+- **Generic code imports only `Scalar`.** The other methods are its supertraits', and a
+  type parameter bounded by `Scalar` carries them. `DualNum`, `DualStruct` and `Signed` are
+  re-exported from the `scalar` module for code that calls them on a concrete dual number.
+- **A call built only from float literals has to name its scalar**, as
+  `model.gtd7d::<f64>(..)`. Inference will not settle a literal on `f64` when the blanket
+  bound admits every dual number as well. Code that gets its position from somewhere typed
+  never meets this.
+
+`num-dual` is therefore not optional, and this crate takes it with `default-features =
+false`: that drops `nalgebra`, which only the vector dual types need and which this crate
+never names. A caller who wants `DualSVec64` enables it from their own manifest, and the
+blanket implementation covers those types without this crate knowing about them.
+
+
 ## `gtd7` or `gtd7d`
 
 They differ in one number. `gtd7`'s total mass density sums helium, atomic oxygen,
@@ -94,7 +127,9 @@ about the model rather than anyone's code.
 below says what it costs.
 
 At a dual type the value agrees to about 1e-14 rather than exactly, because `num-dual`
-computes a division as a reciprocal and a multiply. See the `scalar` module docs.
+computes a division as a reciprocal and a multiply. At `f64` every method is a thin forward
+to the `std` function of the same name, which is why that path is the exact one. See the
+`scalar` module docs.
 
 ## Where it deliberately differs from the reference
 
@@ -252,7 +287,9 @@ use nrlmsise00::{drivers, Model};
 
 let model = Model::new();
 let when = drivers::mjd(2015, 3, 15, 8, 3, 20.0);   // or hifitime's to_mjd_utc_days()
-let output = model.gtd7d_at(when, 400.0, 60.0, -70.0).expect("inside the record");
+let output = model
+    .gtd7d_at::<f64>(when, 400.0, 60.0, -70.0)
+    .expect("inside the record");
 
 assert!((output.total_mass - 2.3e-12).abs() < 0.5e-12);
 # }

src/expansion.rs

@@ -548,18 +548,17 @@ pub(crate) fn globe7<T: Scalar>(
                 * (plg[0][1] * p[119] * swc[5] * cd14 + 1.0)
                 * ((plg[0][1] * p[68] + plg[0][3] * p[69] + plg[0][5] * p[70])
                     * shifted(work.universal, phase[P71], SR, input.seconds_of_day, p[71]));
-            t[11] = t[11]
+            t[11] += (plg[2][3] * p[76] + plg[2][5] * p[77] + plg[2][7] * p[78])
-                + (plg[2][3] * p[76] + plg[2][5] * p[77] + plg[2][7] * p[78])
+                * swc[11]
-                    * swc[11]
+                * mixed(
-                    * mixed(
+                    work.universal,
-                        work.universal,
+                    work.east2,
-                        work.east2,
+                    phase[P79],
-                        phase[P79],
+                    input.seconds_of_day,
-                        input.seconds_of_day,
+                    p[79],
-                        p[79],
+                    input.longitude_deg,
-                        input.longitude_deg,
+                )
-                    )
+                * (dfa * p[137] * swc[1] + 1.0);
-                    * (dfa * p[137] * swc[1] + 1.0);
         }
 
         // Mixed magnetic activity, universal time and longitude.
@@ -621,7 +620,7 @@ pub(crate) fn globe7<T: Scalar>(
     let mut tinf = T::lift(p[30]);
     for (index, term) in t.iter().enumerate() {
         let switch = if index == 8 { sw9 } else { sw[index + 1] };
-        tinf = tinf + *term * switch.abs();
+        tinf += *term * switch.abs();
     }
     tinf
 }
@@ -735,7 +734,7 @@ pub(crate) fn glob7s<T: Scalar>(
     let mut tt = T::lift(0.0);
     for (index, term) in t.iter().enumerate() {
         let switch = if index == 8 { sw9 } else { sw[index + 1] };
-        tt = tt + *term * switch.abs();
+        tt += *term * switch.abs();
     }
     tt
 }

src/helpers.rs

@@ -225,15 +225,14 @@ pub(crate) fn splini<T: Scalar>(xa: &[T], ya: &[T], y2a: &[T], x: T) -> T {
         let b = (xx - xa[klo]) / h;
         let a2 = a * a;
         let b2 = b * b;
-        yi = yi
+        yi += ((T::lift(1.0) - a2) * ya[klo] / 2.0
-            + ((T::lift(1.0) - a2) * ya[klo] / 2.0
+            + b2 * ya[khi] / 2.0
-                + b2 * ya[khi] / 2.0
+            + ((-(a2 * a2 + 1.0) / 4.0 + a2 / 2.0) * y2a[klo]
-                + ((-(a2 * a2 + 1.0) / 4.0 + a2 / 2.0) * y2a[klo]
+                + (b2 * b2 / 4.0 - b2 / 2.0) * y2a[khi])
-                    + (b2 * b2 / 4.0 - b2 / 2.0) * y2a[khi])
+                * h
-                    * h
+                * h
-                    * h
+                / 6.0)
-                    / 6.0)
+            * h;
-                * h;
         klo += 1;
         khi += 1;
     }
@@ -271,8 +270,8 @@ mod tests {
     /// coefficient puts there.
     #[test]
     fn test_gravity_varies_with_latitude_as_the_coefficient_says() {
-        let (equator, _) = glatf(0.0);
+        let (equator, _): (f64, f64) = glatf(0.0);
-        let (pole, _) = glatf(90.0);
+        let (pole, _): (f64, f64) = glatf(90.0);
         assert!(pole > equator);
         let fraction = (pole - equator) / equator;
         // 1 - 0.0026373 * cos(0) against 1 - 0.0026373 * cos(pi), so twice the coefficient.

src/model.rs

@@ -43,7 +43,7 @@ const NUCLEON_MASS: f64 = 1.66e-24;
 /// use nrlmsise00::{Input, MagneticActivity, Model};
 ///
 /// let model = Model::new();
-/// let output = model.gtd7d(&Input::consistent(
+/// let output = model.gtd7d::<f64>(&Input::consistent(
 ///     172,                                     // day of year
 ///     29000.0,                                 // seconds of universal time
 ///     400.0,                                   // km
@@ -138,7 +138,7 @@ impl Model {
         let mut output = self.evaluate(input, &mut work);
         output.total_mass = mass_density(&output, true);
         if self.switches.units == Units::Si {
-            output.total_mass = output.total_mass / 1000.0;
+            output.total_mass /= 1000.0;
         }
         output
     }
@@ -363,7 +363,7 @@ impl Model {
         // the output units, so only the thousand is applied here.
         output.total_mass = mass_density(&output, false);
         if metric {
-            output.total_mass = output.total_mass / 1000.0;
+            output.total_mass /= 1000.0;
         }
 
         // The temperature at the altitude, from the same profile with no mass.
@@ -530,7 +530,7 @@ impl Model {
             let rl = (b28 * PDM[0][1] / b04).ln();
             let zc04 = PDM[0][4] * PDL[1][0];
             let hc04 = PDM[0][5] * PDL[1][1];
-            output.helium = output.helium * ccor(z, rl, hc04, zc04);
+            output.helium *= ccor(z, rl, hc04, zc04);
         }
 
         // Atomic oxygen.
@@ -563,12 +563,12 @@ impl Model {
             let hc16 = PDM[1][5] * PDL[1][3];
             let zc16 = PDM[1][4] * PDL[1][2];
             let hc216 = PDM[1][5] * PDL[1][4];
-            output.atomic_oxygen = output.atomic_oxygen * ccor2(z, rl, hc16, zc16, hc216);
+            output.atomic_oxygen *= ccor2(z, rl, hc16, zc16, hc216);
             // Chemistry correction.
             let hcc16 = PDM[1][7] * PDL[1][13];
             let zcc16 = PDM[1][6] * PDL[1][12];
             let rc16 = PDM[1][3] * PDL[1][14];
-            output.atomic_oxygen = output.atomic_oxygen * ccor(z, T::lift(rc16), hcc16, zcc16);
+            output.atomic_oxygen *= ccor(z, T::lift(rc16), hcc16, zcc16);
         }
 
         // Molecular oxygen.
@@ -600,7 +600,7 @@ impl Model {
                 let rl = (b28 * PDM[3][1] / b32).ln();
                 let hc32 = PDM[3][5] * PDL[1][7];
                 let zc32 = PDM[3][4] * PDL[1][6];
-                output.molecular_oxygen = output.molecular_oxygen * ccor(z, rl, hc32, zc32);
+                output.molecular_oxygen *= ccor(z, rl, hc32, zc32);
             }
             // Correction for a general departure from diffusive equilibrium above the
             // lower boundary, which applies at every altitude.
@@ -609,8 +609,7 @@ impl Model {
             let zcc32 = PDM[3][6] * PDL[1][21];
             let rc32 = ((input.f107_average - 150.0) * (sw[1] * PDL[0][23]) + 1.0)
                 * (PDM[3][3] * PDL[1][23]);
-            output.molecular_oxygen =
+            output.molecular_oxygen *= ccor2(z, rc32, hcc32, zcc32, hcc232);
-                output.molecular_oxygen * ccor2(z, rc32, hcc32, zcc32, hcc232);
         }
 
         // Argon.
@@ -641,7 +640,7 @@ impl Model {
             let rl = (b28 * PDM[4][1] / b40).ln();
             let hc40 = PDM[4][5] * PDL[1][9];
             let zc40 = PDM[4][4] * PDL[1][8];
-            output.argon = output.argon * ccor(z, rl, hc40, zc40);
+            output.argon *= ccor(z, rl, hc40, zc40);
         }
 
         // Hydrogen.
@@ -673,12 +672,12 @@ impl Model {
             let rl = (b28 * PDM[5][1] * PDL[1][17].abs() / b01).ln();
             let hc01 = PDM[5][5] * PDL[1][11];
             let zc01 = PDM[5][4] * PDL[1][10];
-            output.hydrogen = output.hydrogen * ccor(z, rl, hc01, zc01);
+            output.hydrogen *= ccor(z, rl, hc01, zc01);
             // Chemistry correction.
             let hcc01 = PDM[5][7] * PDL[1][19];
             let zcc01 = PDM[5][6] * PDL[1][18];
             let rc01 = PDM[5][3] * PDL[1][20];
-            output.hydrogen = output.hydrogen * ccor(z, T::lift(rc01), hcc01, zcc01);
+            output.hydrogen *= ccor(z, T::lift(rc01), hcc01, zcc01);
         }
 
         // Atomic nitrogen.
@@ -709,12 +708,12 @@ impl Model {
             let rl = (b28 * PDM[6][1] * PDL[0][2].abs() / b14).ln();
             let hc14 = PDM[6][5] * PDL[0][1];
             let zc14 = PDM[6][4] * PDL[0][0];
-            output.atomic_nitrogen = output.atomic_nitrogen * ccor(z, rl, hc14, zc14);
+            output.atomic_nitrogen *= ccor(z, rl, hc14, zc14);
             // Chemistry correction.
             let hcc14 = PDM[6][7] * PDL[0][4];
             let zcc14 = PDM[6][6] * PDL[0][3];
             let rc14 = PDM[6][3] * PDL[0][5];
-            output.atomic_nitrogen = output.atomic_nitrogen * ccor(z, T::lift(rc14), hcc14, zcc14);
+            output.atomic_nitrogen *= ccor(z, T::lift(rc14), hcc14, zcc14);
         }
 
         // Anomalous oxygen: hot atomic oxygen, on its own barometric profile rather than
@@ -793,14 +792,14 @@ fn mass_density<T: Scalar>(output: &Output<T>, anomalous: bool) -> T {
 /// mass density is divided by a thousand. Both operations, in that order: g/cm^3 scaled by
 /// 1e6 and then by 1e-3 is kg/m^3, and doing it in one step would not round the same way.
 fn to_si<T: Scalar>(output: &mut Output<T>) {
-    output.helium = output.helium * 1.0e6;
+    output.helium *= 1.0e6;
-    output.atomic_oxygen = output.atomic_oxygen * 1.0e6;
+    output.atomic_oxygen *= 1.0e6;
-    output.molecular_nitrogen = output.molecular_nitrogen * 1.0e6;
+    output.molecular_nitrogen *= 1.0e6;
-    output.molecular_oxygen = output.molecular_oxygen * 1.0e6;
+    output.molecular_oxygen *= 1.0e6;
-    output.argon = output.argon * 1.0e6;
+    output.argon *= 1.0e6;
-    output.total_mass = output.total_mass * 1.0e6;
+    output.total_mass *= 1.0e6;
-    output.hydrogen = output.hydrogen * 1.0e6;
+    output.hydrogen *= 1.0e6;
-    output.atomic_nitrogen = output.atomic_nitrogen * 1.0e6;
+    output.atomic_nitrogen *= 1.0e6;
-    output.anomalous_oxygen = output.anomalous_oxygen * 1.0e6;
+    output.anomalous_oxygen *= 1.0e6;
-    output.total_mass = output.total_mass / 1000.0;
+    output.total_mass /= 1000.0;
 }

src/profile.rs

@@ -106,7 +106,7 @@ pub(crate) fn densu<T: Scalar>(
     if tt.re() <= 0.0 {
         expl = T::lift(50.0);
     }
-    densu_temp = dlb * (tlb / tt).powt(gamma + (1.0 + alpha)) * expl;
+    densu_temp = dlb * (tlb / tt).powd(gamma + (1.0 + alpha)) * expl;
     if alt.re() >= za {
         return densu_temp;
     }

src/scalar.rs

@@ -5,27 +5,90 @@
 //! exists: density falls by a factor of e about every 50 km, so a drag partial that drops
 //! its dependence on position is badly wrong rather than slightly wrong.
 //!
-//! [`Scalar`] is implemented for `f64` always, and for the dual numbers of the
+//! # It is `num-dual`'s trait, not one of ours
-//! [`num_dual`] crate under the `autodiff` feature, which is on by default.
+//!
+//! [`Scalar`] declares one method of its own and inherits the rest from
+//! [`num_dual::DualNum`], which `num-dual` implements for `f64` and `f32` as well as for
+//! its dual types. It is then blanket-implemented, so **every** type satisfying that trait
+//! at `f64` precision is a scalar here, including nestings and types this crate has never
+//! heard of.
+//!
+//! That is deliberate. An earlier version declared `exp`, `ln`, `powf` and the rest itself
+//! and implemented them for a list of types. Two things were wrong with that. The list
+//! could not be complete, and nobody else could extend it: a caller whose scalar was a
+//! dual type the list had missed was stopped by the orphan rule, with no way to implement
+//! this crate's trait for it. And declaring the same methods a second time made this trait
+//! incompatible with every other crate's trait over the same types — `ordinary-diffeq`
+//! declares one of the same shape — so a caller generic over both was adapting between
+//! two names for one thing. Taking the bound from `num-dual` leaves one name.
+//!
+//! A caller generic over another crate's scalar trait still writes both bounds, since
+//! neither trait implies the other, but the second one is now a formality that any dual
+//! number satisfies rather than a wall that some do not.
+//!
+//! # What has to be in scope
+//!
+//! Nothing, to write generic code. The methods belong to [`Scalar`]'s supertraits, and a
+//! type parameter bounded by `Scalar` carries them without the supertraits being imported:
+//!
+//! ```
+//! use nrlmsise00::Scalar;
+//!
+//! fn scale_height<T: Scalar>(density: T, slope: T) -> T {
+//!     // ln from DualNum, abs from Signed, re from DualStruct, lift from Scalar, and only
+//!     // Scalar is imported.
+//!     let _ = (density.ln(), density.abs(), density.re());
+//!     -density / slope + T::lift(0.0)
+//! }
+//!
+//! assert!((scale_height(1.0_f64, -0.02) - 50.0).abs() < 1e-12);
+//! ```
+//!
+//! Calling the same methods on a *concrete* dual number does need the trait that declares
+//! them, and the three are re-exported here for that: [`DualNum`] for the arithmetic,
+//! [`DualStruct`] for [`DualStruct::re`], [`Signed`] for [`Signed::abs`].
+//!
+//! # A literal-only call has to name its scalar
+//!
+//! The one wrinkle of the bound being a blanket one. Inference will not settle a float
+//! literal on `f64` when the only thing constraining it is `T: Scalar`, because the
+//! blanket implementation admits every dual number too, so a call built entirely from
+//! literals is ambiguous and says so at the first method call on its result. Name the
+//! scalar once and the whole call follows:
+//!
+//! ```
+//! # use nrlmsise00::{Input, MagneticActivity, Model};
+//! let model = Model::new();
+//! let input = Input::consistent(
+//!     172, 29000.0, 400.0, 60.0, -70.0, 150.0, 150.0,
+//!     MagneticActivity::Daily { ap: 4.0 },
+//! );
+//! let output = model.gtd7d::<f64>(&input);
+//! assert!(output.total_mass > 0.0);
+//! ```
+//!
+//! Code that gets its position and time from somewhere with a type never meets this.
 //!
 //! # Writing against it
 //!
 //! Three rules, each of which the model's own code obeys throughout:
 //!
-//! - **Magnitude comparisons branch on [`Scalar::re`].** A dual number carries a
+//! - **Magnitude comparisons branch on [`DualStruct::re`].** `num-dual` does order its
-//!   derivative that an ordering cannot account for, and every comparison in the model —
+//!   dual numbers, and it does so on the real part alone, so `x > y` would agree with
-//!   the altitude bands, the exponent clamps, the guards in `ccor` and `dnet` — is asking
+//!   `x.re() > y.re()`. The model writes the second anyway: every comparison in it — the
-//!   about a real quantity, where the real part is the right thing to test.
+//!   altitude bands, the exponent clamps, the guards in `ccor` and `dnet` — is a question
+//!   about a real quantity, and saying so at each site is worth the few characters.
 //! - **Coefficients stay `f64`.** `T op f64` is available, so lifting three thousand
 //!   coefficients to `T` would only cost memory and cache.
 //! - **`f64` on the left of an operator does not exist.** `T * 2.0` works; `2.0 * T`
 //!   cannot, because of the orphan rule. Write the scalar on the left, or
-//!   [`Scalar::lift`] the constant.
+//!   [`Scalar::lift`] the constant. Where the reference puts a constant on the left of a
+//!   *non-commutative* operator, lifting is the only way to keep its arithmetic.
 //!
 //! # A dual evaluation's value is not bit for bit the `f64` one
 //!
-//! It agrees to about 1e-14 relative. The difference is `num-dual`'s arithmetic rather
+//! It agrees to about 1e-14. The difference is `num-dual`'s arithmetic rather than
-//! than anything this crate does, and two places account for it:
+//! anything this crate does, and two places account for it:
 //!
 //! - **Division goes through the reciprocal.** `Dual / Dual` is computed as
 //!   `self.re * other.re.recip()`, two roundings where a division is one, so a dual
@@ -39,245 +102,87 @@
 //! Both were measured rather than assumed. Correcting them here was tried and dropped: a
 //! real part can be fixed exactly, by adding back a difference that is itself exactly
 //! representable, but the same cannot be done for division without replacing every `/` in
-//! the model, so the result was a slower model that was still not exact. The `f64` path is
+//! the model, so the result was a slower model that was still not exact. At `f64` every
-//! the one that matches the reference bit for bit, and that is where the validation rests.
+//! method is a thin forward to the `std` function of the same name, so that path is the
-//!
+//! one that matches the reference bit for bit, and it is where the validation rests.
-//! # A consumer generic over its own scalar trait
-//!
-//! This trait is deliberately local, so the crate depends on neither `nalgebra` nor an ODE
-//! solver. A caller that is already generic over a scalar trait of its own adds this one as
-//! a second bound:
-//!
-//! ```ignore
-//! fn drag<T: my_crate::Scalar + nrlmsise00::Scalar>(...) { ... }
-//! ```
-//!
-//! Both are satisfied by `f64` and by `num-dual`'s types, so nothing else changes.
 
-use core::fmt::Debug;
+pub use num_dual::{DualNum, DualStruct};
-use core::ops::{Add, Div, Mul, Neg, Sub};
+pub use num_traits::Signed;
 
 /// A scalar the model can be evaluated in.
 ///
-/// The methods are the ones NRLMSISE-00 actually computes with, and no more.
+/// Satisfied by `f64`, by `num-dual`'s dual numbers over it, and by anything else that
-pub trait Scalar:
+/// implements [`DualNum`] at `f64` precision and is [`Copy`]. There is nothing to
-    Copy
+/// implement: the blanket implementation below covers every such type.
-    + Debug
+///
-    + Add<Output = Self>
+/// The methods the model computes with are [`DualNum`]'s, [`DualStruct::re`] and
-    + Sub<Output = Self>
+/// [`Signed::abs`]; see the module documentation for what to import.
-    + Mul<Output = Self>
+pub trait Scalar: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy {
-    + Div<Output = Self>
-    + Neg<Output = Self>
-    + Add<f64, Output = Self>
-    + Sub<f64, Output = Self>
-    + Mul<f64, Output = Self>
-    + Div<f64, Output = Self>
-{
-    /// This value with all derivative information discarded.
-    ///
-    /// Every magnitude comparison in the model goes through here.
-    fn re(&self) -> f64;
-
     /// A real constant as a value of this type, with zero derivative.
-    fn lift(x: f64) -> Self;
-
-    /// `e^self`.
-    fn exp(self) -> Self;
-
-    /// The natural logarithm.
-    fn ln(self) -> Self;
-
-    /// The base-10 logarithm.
-    fn log10(self) -> Self;
-
-    /// The square root.
-    fn sqrt(self) -> Self;
-
-    /// The sine, in radians.
-    fn sin(self) -> Self;
-
-    /// The cosine, in radians.
-    fn cos(self) -> Self;
-
-    /// The magnitude.
     ///
-    /// The reference writes `sqrt(x*x)` for this. The two agree bit for bit at `f64`, and
+    /// `num-dual` spells this `Self::from`, which is a less obvious thing to read in the
-    /// this one also has a derivative at zero, which `sqrt(x*x)` does not: the model takes
+    /// middle of an expansion than `T::lift`, and which the model writes a hundred and
-    /// the magnitude of the altitude and of the latitude, and both are sampled at zero.
+    /// twenty times.
-    fn abs(self) -> Self;
-
-    /// `self` raised to an integer power.
-    fn powi(self, n: i32) -> Self;
-
-    /// `self` raised to a real power.
-    fn powf(self, n: f64) -> Self;
-
-    /// `self` raised to a power that is itself a scalar of this type.
-    ///
-    /// At `f64` this is the platform `pow`, which is what the reference calls.
-    fn powt(self, n: Self) -> Self;
-}
-
-impl Scalar for f64 {
-    #[inline(always)]
-    fn re(&self) -> f64 {
-        *self
-    }
-
-    #[inline(always)]
     fn lift(x: f64) -> Self {
-        x
+        Self::from(x)
-    }
-
-    #[inline(always)]
-    fn exp(self) -> Self {
-        f64::exp(self)
-    }
-
-    #[inline(always)]
-    fn ln(self) -> Self {
-        f64::ln(self)
-    }
-
-    #[inline(always)]
-    fn log10(self) -> Self {
-        f64::log10(self)
-    }
-
-    #[inline(always)]
-    fn sqrt(self) -> Self {
-        f64::sqrt(self)
-    }
-
-    #[inline(always)]
-    fn sin(self) -> Self {
-        f64::sin(self)
-    }
-
-    #[inline(always)]
-    fn cos(self) -> Self {
-        f64::cos(self)
-    }
-
-    #[inline(always)]
-    fn abs(self) -> Self {
-        f64::abs(self)
-    }
-
-    #[inline(always)]
-    fn powi(self, n: i32) -> Self {
-        f64::powi(self, n)
-    }
-
-    #[inline(always)]
-    fn powf(self, n: f64) -> Self {
-        f64::powf(self, n)
-    }
-
-    #[inline(always)]
-    fn powt(self, n: Self) -> Self {
-        f64::powf(self, n)
     }
 }
 
-#[cfg(feature = "autodiff")]
+impl<T: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy> Scalar for T {}
-mod dual {
-    use super::Scalar;
-    use num_dual::{DualNum, DualStruct};
 
-    /// Implements `Scalar` for one of `num-dual`'s dual number types.
+#[cfg(test)]
-    ///
+mod tests {
-    /// The bounds are stated on `Self` rather than on the type parameters so that one arm
+    use super::*;
-    /// covers every nesting depth those types admit.
+
-    macro_rules! impl_scalar {
+    /// The blanket implementation reaches the primitive, the dual types, and a nesting of
-        ($ty:ident) => {
+    /// them that this crate never names — which is the point of it being blanket.
-            impl<T> Scalar for num_dual::$ty<T>
+    #[test]
-            where
+    fn test_the_blanket_implementation_reaches_what_it_should() {
-                Self: Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
+        fn assert_scalar<T: Scalar>(x: T) -> f64 {
-            {
+            // One method from each of the four traits, to show they are all available.
-                impl_scalar_body!();
+            (x.exp().ln().abs() + T::lift(1.0)).re()
-            }
+        }
-        };
+
-        ($ty:ident, $($n:ident),+) => {
+        assert_eq!(assert_scalar(2.0_f64), 3.0);
-            impl<T: nalgebra::Scalar $(, const $n: usize)+> Scalar
+        assert_eq!(assert_scalar(num_dual::Dual64::from(2.0)), 3.0);
-                for num_dual::$ty<T $(, nalgebra::Const<$n>)+>
+        assert_eq!(assert_scalar(num_dual::Dual2_64::from(2.0)), 3.0);
-            where
+        assert_eq!(assert_scalar(num_dual::HyperDual64::from(2.0)), 3.0);
-                Self: Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
+        // A dual of a dual, which the old list of implementations did not include.
-            {
+        assert_eq!(
-                impl_scalar_body!();
+            assert_scalar(num_dual::Dual::<num_dual::Dual64>::from(2.0)),
-            }
+            3.0
-        };
+        );
     }
 
-    /// The method bodies, which are the same for every one of those types.
+    /// `lift` is a constant: the value with no derivative, whatever the order.
-    macro_rules! impl_scalar_body {
+    #[test]
-        () => {
+    fn test_lift_carries_no_derivative() {
-            #[inline(always)]
+        use num_dual::Dual64;
-            fn re(&self) -> f64 {
-                DualStruct::re(self)
-            }
 
-            #[inline(always)]
+        let constant = Dual64::lift(3.5);
-            fn lift(x: f64) -> Self {
+        assert_eq!(constant.re, 3.5);
-                Self::from(x)
+        assert_eq!(constant.eps, 0.0);
-            }
 
-            #[inline(always)]
+        let seeded = Dual64::from(3.5).derivative();
-            fn exp(self) -> Self {
+        assert_eq!(seeded.re(), 3.5);
-                DualNum::exp(&self)
+        assert_eq!(seeded.eps, 1.0);
-            }
-
-            #[inline(always)]
-            fn ln(self) -> Self {
-                DualNum::ln(&self)
-            }
-
-            #[inline(always)]
-            fn log10(self) -> Self {
-                DualNum::log10(&self)
-            }
-
-            #[inline(always)]
-            fn sqrt(self) -> Self {
-                DualNum::sqrt(&self)
-            }
-
-            #[inline(always)]
-            fn sin(self) -> Self {
-                DualNum::sin(&self)
-            }
-
-            #[inline(always)]
-            fn cos(self) -> Self {
-                DualNum::cos(&self)
-            }
-
-            #[inline(always)]
-            fn abs(self) -> Self {
-                num_traits::Signed::abs(&self)
-            }
-
-            #[inline(always)]
-            fn powi(self, n: i32) -> Self {
-                DualNum::powi(&self, n)
-            }
-
-            #[inline(always)]
-            fn powf(self, n: f64) -> Self {
-                DualNum::powf(&self, n)
-            }
-
-            #[inline(always)]
-            fn powt(self, n: Self) -> Self {
-                DualNum::powd(&self, n)
-            }
-        };
     }
 
-    impl_scalar!(Dual);
+    /// Ordering a dual number compares its real part, which is what the model's branches
-    impl_scalar!(Dual2);
+    /// want, so the two forms agree.
-    impl_scalar!(Dual3);
+    #[test]
-    impl_scalar!(HyperDual);
+    fn test_ordering_a_dual_compares_the_real_part() {
-    impl_scalar!(DualVec, N);
+        use num_dual::Dual64;
-    impl_scalar!(Dual2Vec, N);
+
+        let small = Dual64::from(1.0).derivative();
+        let large = Dual64::from(2.0);
+        assert!(small < large);
+        assert!(small.re() < large.re());
+
+        // Equal real parts compare equal even though the derivatives differ, where a
+        // lexicographic ordering would have made one of them the smaller.
+        let other = Dual64::from(1.0);
+        assert_ne!(small.eps, other.eps);
+        assert_eq!(small.partial_cmp(&other), Some(core::cmp::Ordering::Equal));
+    }
 }

tests/derivatives.rs

@@ -1,5 +1,3 @@
-#![cfg(feature = "autodiff")]
-
 //! That the model differentiates.
 //!
 //! This is the one thing a port buys that depending on an existing `f64` implementation
@@ -116,7 +114,9 @@ fn test_the_model_jumps_at_the_species_threshold_altitudes() {
     // join between the Bates profile above and the spline below, and an eighth at
     // 428.475 km, too small to pin here, is where atomic nitrogen's turbopause blend
     // saturates and `dnet` switches to the mixed density with a step of about exp(-10).
-    let thresholds = [
+    // Annotated, because `Scalar` is blanket-implemented over `num-dual`'s trait and a
+    // bare float literal flowing into the generic model no longer falls back to `f64`.
+    let thresholds: [(f64, f64); 7] = [
         (72.5, 1.0e-3),
         (123.435, 1.5e-3),
         (160.0, 7.1e-5),
@@ -169,7 +169,7 @@ fn test_the_model_jumps_at_the_species_threshold_altitudes() {
     }
 
     // Above 500 km there are none at all, which is the range most of a mission sits in.
-    let mut z = 500.0;
+    let mut z: f64 = 500.0;
     while z < 1000.0 {
         let step = 1.0e-9;
         let below = model.gtd7d(&input(z - step)).total_mass;