Repository: itt-ustutt/num-dual (language: python)

Issue:
Compatibility of 'DualNum' trait with `nalgebra`
By implementing [nalgebra::ComplexField](https://docs.rs/nalgebra/latest/nalgebra/trait.ComplexField.html) for dual numbers (analogous to the implementation for `f64`) the full linear algebra capabilities of nalgebra (in particular eigenvalues) can be used with arbitrary generalized (hyper-) dual numbers. Currently this is only possible with first order dual numbers.

In this context it is worth reconsidering whether the `Float` trait should also be part of the supertraits of `DualNum`.

If this is implemented, the `linalg` module can be discarded in favor of the well established functionalities in `nalgebra`.


Discussion / hints:
The required traits for using the eigenvalue function in nalgebra are so unfathomably restrictive (e.g., including `SimdValue`) that it is unlikely that they will ever be implemented for structs other than `f32` and `f64`. Therefore, it is unlikely that we can completely get rid of the eigenvalue algorithm in `num-dual`.

Fix patch:
diff --git a/Cargo.toml b/Cargo.toml
index 6ce203c..fab0df5 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -19,6 +19,8 @@ nalgebra = "0.32"
 pyo3 = { version = "0.18", optional = true, features = ["multiple-pymethods"] }
 ndarray = { version = "0.15", optional = true }
 numpy = { version = "0.18", optional = true }
+approx = "0.5"
+simba = "0.8"
 
 [profile.release]
 lto = true
@@ -30,7 +32,6 @@ linalg = ["ndarray"]
 
 [dev-dependencies]
 criterion = "0.4"
-approx = "0.5"
 
 [[bench]]
 name = "benchmark"
diff --git a/src/derivative.rs b/src/derivative.rs
index cb3e3b2..7824362 100644
--- a/src/derivative.rs
+++ b/src/derivative.rs
@@ -4,7 +4,9 @@ use nalgebra::constraint::{SameNumberOfRows, ShapeConstraint};
 use nalgebra::*;
 use std::fmt;
 use std::marker::PhantomData;
+use std::mem::MaybeUninit;
 use std::ops::{Add, AddAssign, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
+use num_traits::Zero;
 
 #[derive(PartialEq, Eq, Clone, Debug)]
 pub struct Derivative<T: DualNum<F>, F, R: Dim, C: Dim>(
@@ -35,6 +37,84 @@ where
         Self::new(None)
     }
 
+    pub(crate) fn map<T2, F2>(&self, mut f: impl FnMut(T) -> T2) -> Derivative<T2, F2, R, C>
+    where
+        T2: DualNum<F2>,
+        DefaultAllocator: Allocator<T2, R, C>,
+    {
+        let opt = self.0.as_ref().map(move |eps| eps.map(|e| f(e)));
+        Derivative::new(opt)
+    }
+
+    // A version of map that doesn't clone values before mapping. Useful for the SimdValue impl,
+    // which would be redundantly cloning all the lanes of each epsilon value before extracting
+    // just one of them.
+    //
+    // To implement, we inline a copy of Matrix::map, which implicitly clones values, and remove
+    // the cloning.
+    pub(crate) fn map_borrowed<T2, F2>(
+        &self,
+        mut f: impl FnMut(&T) -> T2,
+    ) -> Derivative<T2, F2, R, C>
+    where
+        T2: DualNum<F2>,
+        DefaultAllocator: Allocator<T2, R, C>,
+    {
+        let opt = self.0.as_ref().map(move |eps| {
+            let ref this = eps;
+            let mut f = |e| f(e);
+            let (nrows, ncols) = this.shape_generic();
+            let mut res: Matrix<MaybeUninit<T2>, R, C, _> = Matrix::uninit(nrows, ncols);
+
+            for j in 0..ncols.value() {
+                for i in 0..nrows.value() {
+                    // Safety: all indices are in range.
+                    unsafe {
+                        let a = this.data.get_unchecked(i, j);
+                        *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(a));
+                    }
+                }
+            }
+
+            // Safety: res is now fully initialized.
+            unsafe { res.assume_init() }
+        });
+        Derivative::new(opt)
+    }
+
+    /// Same but bails out if the closure returns None
+    pub(crate) fn try_map_borrowed<T2, F2>(
+        &self,
+        mut f: impl FnMut(&T) -> Option<T2>,
+    ) -> Option<Derivative<T2, F2, R, C>>
+    where
+        T2: DualNum<F2>,
+        DefaultAllocator: Allocator<T2, R, C>,
+    {
+        self.0
+            .as_ref()
+            .and_then(move |eps| {
+                let ref this = eps;
+                let mut f = |e| f(e);
+                let (nrows, ncols) = this.shape_generic();
+                let mut res: Matrix<MaybeUninit<T2>, R, C, _> = Matrix::uninit(nrows, ncols);
+
+                for j in 0..ncols.value() {
+                    for i in 0..nrows.value() {
+                        // Safety: all indices are in range.
+                        unsafe {
+                            let a = this.data.get_unchecked(i, j);
+                            *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(a)?);
+                        }
+                    }
+                }
+
+                // Safety: res is now fully initialized.
+                Some(unsafe { res.assume_init() })
+            })
+            .map(Derivative::some)
+    }
+
     pub fn derivative_generic(r: R, c: C, i: usize) -> Self {
         let mut m = OMatrix::zeros_generic(r, c);
         m[i] = T::one();
@@ -304,3 +384,167 @@ where
         }
     }
 }
+
+impl<T, R: Dim, C: Dim> nalgebra::SimdValue for Derivative<T, T::Element, R, C>
+where
+    DefaultAllocator: Allocator<T, R, C> + Allocator<T::Element, R, C>,
+    T: DualNum<T::Element> + SimdValue + Scalar,
+    T::Element: DualNum<T::Element> + Scalar + Zero,
+{
+    type Element = Derivative<T::Element, T::Element, R, C>;
+
+    type SimdBool = T::SimdBool;
+
+    #[inline]
+    fn lanes() -> usize {
+        T::lanes()
+    }
+
+    #[inline]
+    fn splat(val: Self::Element) -> Self {
+        val.map(|e| T::splat(e))
+    }
+
+    #[inline]
+    fn extract(&self, i: usize) -> Self::Element {
+        self.map_borrowed(|e| T::extract(e, i))
+    }
+
+    #[inline]
+    unsafe fn extract_unchecked(&self, i: usize) -> Self::Element {
+        let opt = self
+            .map_borrowed(|e| T::extract_unchecked(e, i))
+            .0
+            // Now check it's all zeros.
+            // Unfortunately there is no way to use the vectorized version of `is_zero`, which is
+            // only for matrices with statically known dimensions. Specialization would be
+            // required.
+            .filter(|x| Iterator::any(&mut x.iter(), |e| !e.is_zero()));
+        Derivative::new(opt)
+    }
+
+    // SIMD code will expect to be able to replace one lane with another Self::Element,
+    // even with a None Derivative, e.g.
+    //
+    // let single = Derivative::none();
+    // let mut x4 = Derivative::splat(single);
+    // let one = Derivative::some(...);
+    // x4.replace(1, one);
+    //
+    // So the implementation of `replace` will need to auto-upgrade to Some(zeros) in
+    // order to satisfy requests like that.
+    fn replace(&mut self, i: usize, val: Self::Element) {
+        match (&mut self.0, val.0) {
+            (Some(ours), Some(theirs)) => {
+                ours.zip_apply(&theirs, |e, replacement| e.replace(i, replacement));
+            }
+            (ours @ None, Some(theirs)) => {
+                let (r, c) = theirs.shape_generic();
+                let mut init: OMatrix<T, R, C> = OMatrix::zeros_generic(r, c);
+                init.zip_apply(&theirs, |e, replacement| e.replace(i, replacement));
+                *ours = Some(init);
+            }
+            (Some(ours), None) => {
+                ours.apply(|e| e.replace(i, T::Element::zero()));
+            }
+            _ => {}
+        }
+    }
+
+    unsafe fn replace_unchecked(&mut self, i: usize, val: Self::Element) {
+        match (&mut self.0, val.0) {
+            (Some(ours), Some(theirs)) => {
+                ours.zip_apply(&theirs, |e, replacement| {
+                    e.replace_unchecked(i, replacement)
+                });
+            }
+            (ours @ None, Some(theirs)) => {
+                let (r, c) = theirs.shape_generic();
+                let mut init: OMatrix<T, R, C> = OMatrix::zeros_generic(r, c);
+                init.zip_apply(&theirs, |e, replacement| {
+                    e.replace_unchecked(i, replacement)
+                });
+                *ours = Some(init);
+            }
+            (Some(ours), None) => {
+                ours.apply(|e| e.replace_unchecked(i, T::Element::zero()));
+            }
+            _ => {}
+        }
+    }
+
+    fn select(mut self, cond: Self::SimdBool, other: Self) -> Self {
+        // If cond is mixed, then we may need to generate big zero matrices to do the
+        // component-wise select on. So check if cond is all-true or all-first to avoid that.
+        if cond.all() {
+            self
+        } else if cond.none() {
+            other
+        } else {
+            match (&mut self.0, other.0) {
+                (Some(ours), Some(theirs)) => {
+                    ours.zip_apply(&theirs, |e, other_e| {
+                        // this will probably get optimized out
+                        let e_ = std::mem::replace(e, T::zero());
+                        *e = e_.select(cond, other_e)
+                    });
+                    self
+                }
+                (Some(ours), None) => {
+                    ours.apply(|e| {
+                        // this will probably get optimized out
+                        let e_ = std::mem::replace(e, T::zero());
+                        *e = e_.select(cond, T::zero());
+                    });
+                    self
+                }
+                (ours @ None, Some(mut theirs)) => {
+                    use std::ops::Not;
+                    let inverted: T::SimdBool = cond.not();
+                    theirs.apply(|e| {
+                        // this will probably get optimized out
+                        let e_ = std::mem::replace(e, T::zero());
+                        *e = e_.select(inverted, T::zero());
+                    });
+                    *ours = Some(theirs);
+                    self
+                }
+                _ => self,
+            }
+        }
+    }
+}
+
+use simba::scalar::{SubsetOf, SupersetOf};
+
+impl<TSuper, FSuper, T, F, R: Dim, C: Dim> SubsetOf<Derivative<TSuper, FSuper, R, C>>
+    for Derivative<T, F, R, C>
+where
+    TSuper: DualNum<FSuper> + SupersetOf<T>,
+    T: DualNum<F>,
+    DefaultAllocator: Allocator<T, R, C>,
+    DefaultAllocator: Allocator<TSuper, R, C>,
+    // DefaultAllocator: Allocator<TSuper, D>
+    //     + Allocator<TSuper, U1, D>
+    //     + Allocator<TSuper, D, U1>
+    //     + Allocator<TSuper, D, D>,
+{
+    #[inline(always)]
+    fn to_superset(&self) -> Derivative<TSuper, FSuper, R, C> {
+        self.map_borrowed(|elem| TSuper::from_subset(elem))
+    }
+    #[inline(always)]
+    fn from_superset(element: &Derivative<TSuper, FSuper, R, C>) -> Option<Self> {
+        element.try_map_borrowed(|elem| TSuper::to_subset(elem))
+    }
+    #[inline(always)]
+    fn from_superset_unchecked(element: &Derivative<TSuper, FSuper, R, C>) -> Self {
+        element.map_borrowed(|elem| TSuper::to_subset_unchecked(elem))
+    }
+    #[inline(always)]
+    fn is_in_subset(element: &Derivative<TSuper, FSuper, R, C>) -> bool {
+        element.0.as_ref().map_or(true, |matrix| {
+            matrix.iter().all(|elem| TSuper::is_in_subset(elem))
+        })
+    }
+}
diff --git a/src/dual.rs b/src/dual.rs
index 6ddd347..dcaf877 100644
--- a/src/dual.rs
+++ b/src/dual.rs
@@ -1,4 +1,5 @@
 use crate::{Derivative, DualNum, DualNumFloat};
+use approx::{AbsDiffEq, RelativeEq, UlpsEq};
 use nalgebra::allocator::Allocator;
 use nalgebra::*;
 use num_traits::{Float, FloatConst, FromPrimitive, Inv, Num, One, Signed, Zero};
@@ -11,7 +12,7 @@ use std::ops::{
 };
 
 /// A dual number for the calculations of gradients or Jacobians.
-#[derive(PartialEq, Eq, Clone, Debug)]
+#[derive(Clone, Debug)]
 pub struct DualVec<T: DualNum<F>, F, D: Dim>
 where
     DefaultAllocator: Allocator<T, D>,
@@ -283,3 +284,727 @@ where
 
 impl_first_derivatives!(DualVec, [eps], [D]);
 impl_dual!(DualVec, [eps], [D]);
+
+/**
+ * The SimdValue trait is for rearranging data into a form more suitable for Simd,
+ * and rearranging it back into a usable form. It is not documented particularly well.
+ *
+ * The primary job of this SimdValue impl is to allow people to use `simba::simd::f32x4` etc,
+ * instead of f32/f64. Those types implement nalgebra::SimdRealField/ComplexField, so they
+ * behave like scalars. When we use them, we would have `DualVec<f32x4, f32, N>` etc, with our
+ * F parameter set to <T as SimdValue>::Element. We will need to be able to split up that type
+ * into four of DualVec in order to get out of simd-land. That's what the SimdValue trait is for.
+ *
+ * Ultimately, someone will have to to implement SimdRealField on DualVec and call the
+ * simd_ functions of <T as SimdRealField>. That's future work for someone who finds
+ * num_dual is not fast enough.
+ *
+ * Unfortunately, doing anything with SIMD is blocked on
+ * <https://github.com/dimforge/simba/issues/44>.
+ *
+ */
+impl<T, D: Dim> nalgebra::SimdValue for DualVec<T, T::Element, D>
+where
+    DefaultAllocator: Allocator<T, D> + Allocator<T::Element, D>,
+    T: DualNum<T::Element> + SimdValue + Scalar,
+    T::Element: DualNum<T::Element> + Scalar,
+{
+    // Say T = simba::f32x4. T::Element is f32. T::SimdBool is AutoSimd<[bool; 4]>.
+    // AutoSimd<[f32; 4]> stores an actual [f32; 4], i.e. four floats in one slot.
+    // So our DualVec<AutoSimd<[f32; 4], f32, N> has 4 * (1+N) floats in it, stored in blocks of
+    // four. When we want to do any math on it but ignore its f32x4 storage mode, we need to break
+    // that type into FOUR of DualVec<f32, f32, N>; then we do math on it, then we bring it back
+    // together.
+    //
+    // Hence this definition of Element:
+    type Element = DualVec<T::Element, T::Element, D>;
+    type SimdBool = T::SimdBool;
+
+    #[inline]
+    fn lanes() -> usize {
+        T::lanes()
+    }
+
+    #[inline]
+    fn splat(val: Self::Element) -> Self {
+        // Need to make `lanes` copies of each of:
+        // - the real part
+        // - each of the N epsilon parts
+        let re = T::splat(val.re);
+        let eps = Derivative::splat(val.eps);
+        Self {
+            re,
+            eps,
+            f: PhantomData,
+        }
+    }
+
+    #[inline]
+    fn extract(&self, i: usize) -> Self::Element {
+        let re = self.re.extract(i);
+        let eps = self.eps.extract(i);
+        Self::Element {
+            re,
+            eps,
+            f: PhantomData,
+        }
+    }
+
+    #[inline]
+    unsafe fn extract_unchecked(&self, i: usize) -> Self::Element {
+        let re = self.re.extract_unchecked(i);
+        let eps = self.eps.extract_unchecked(i);
+        Self::Element {
+            re,
+            eps,
+            f: PhantomData,
+        }
+    }
+
+    #[inline]
+    fn replace(&mut self, i: usize, val: Self::Element) {
+        self.re.replace(i, val.re);
+        self.eps.replace(i, val.eps);
+    }
+
+    #[inline]
+    unsafe fn replace_unchecked(&mut self, i: usize, val: Self::Element) {
+        self.re.replace_unchecked(i, val.re);
+        self.eps.replace_unchecked(i, val.eps);
+    }
+
+    #[inline]
+    fn select(self, cond: Self::SimdBool, other: Self) -> Self {
+        let re = self.re.select(cond, other.re);
+        let eps = self.eps.select(cond, other.eps);
+        Self {
+            re,
+            eps,
+            f: PhantomData,
+        }
+    }
+}
+
+/// Comparisons are only made based on the real part. This allows the code to follow the
+/// same execution path as real-valued code would.
+impl<T: DualNum<F> + PartialEq, F: Float, D: Dim> PartialEq for DualVec<T, F, D>
+where
+    DefaultAllocator: Allocator<T, D>,
+{
+    #[inline]
+    fn eq(&self, other: &Self) -> bool {
+        self.re.eq(&other.re)
+    }
+}
+/// Like PartialEq, comparisons are only made based on the real part. This allows the code to follow the
+/// same execution path as real-valued code would.
+impl<T: DualNum<F> + PartialOrd, F: Float, D: Dim> PartialOrd for DualVec<T, F, D>
+where
+    DefaultAllocator: Allocator<T, D>,
+{
+    #[inline]
+    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
+        self.re.partial_cmp(&other.re)
+    }
+}
+/// Like PartialEq, comparisons are only made based on the real part. This allows the code to follow the
+/// same execution path as real-valued code would.
+impl<T: DualNum<F> + approx::AbsDiffEq<Epsilon = T>, F: Float, D: Dim> approx::AbsDiffEq
+    for DualVec<T, F, D>
+where
+    DefaultAllocator: Allocator<T, D>,
+{
+    type Epsilon = Self;
+    #[inline]
+    fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
+        self.re.abs_diff_eq(&other.re, epsilon.re)
+    }
+
+    #[inline]
+    fn default_epsilon() -> Self::Epsilon {
+        Self::from_re(T::default_epsilon())
+    }
+}
+/// Like PartialEq, comparisons are only made based on the real part. This allows the code to follow the
+/// same execution path as real-valued code would.
+impl<T: DualNum<F> + approx::RelativeEq<Epsilon = T>, F: Float, D: Dim> approx::RelativeEq
+    for DualVec<T, F, D>
+where
+    DefaultAllocator: Allocator<T, D>,
+{
+    #[inline]
+    fn default_max_relative() -> Self::Epsilon {
+        Self::from_re(T::default_max_relative())
+    }
+
+    #[inline]
+    fn relative_eq(
+        &self,
+        other: &Self,
+        epsilon: Self::Epsilon,
+        max_relative: Self::Epsilon,
+    ) -> bool {
+        self.re.relative_eq(&other.re, epsilon.re, max_relative.re)
+    }
+}
+impl<T: DualNum<F> + UlpsEq<Epsilon = T>, F: Float, D: Dim> UlpsEq for DualVec<T, F, D>
+where
+    DefaultAllocator: Allocator<T, D>,
+{
+    #[inline]
+    fn default_max_ulps() -> u32 {
+        T::default_max_ulps()
+    }
+
+    #[inline]
+    fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
+        T::ulps_eq(&self.re, &other.re, epsilon.re, max_ulps)
+    }
+}
+
+impl<T, D: Dim> nalgebra::Field for DualVec<T, T::Element, D>
+where
+    T: DualNum<T::Element> + SimdValue,
+    T::Element: DualNum<T::Element> + Scalar + Float,
+    DefaultAllocator:
+        Allocator<T, D> + Allocator<T, U1, D> + Allocator<T, D, U1> + Allocator<T, D, D>,
+    DefaultAllocator: Allocator<T::Element, D>
+        + Allocator<T::Element, U1, D>
+        + Allocator<T::Element, D, U1>
+        + Allocator<T::Element, D, D>,
+{
+}
+
+use simba::scalar::{SubsetOf, SupersetOf};
+
+impl<TSuper, FSuper, T, F, D: Dim> SubsetOf<DualVec<TSuper, FSuper, D>> for DualVec<T, F, D>
+where
+    TSuper: DualNum<FSuper> + SupersetOf<T>,
+    T: DualNum<F>,
+    DefaultAllocator:
+        Allocator<T, D> + Allocator<T, U1, D> + Allocator<T, D, U1> + Allocator<T, D, D>,
+    DefaultAllocator: Allocator<TSuper, D>
+        + Allocator<TSuper, U1, D>
+        + Allocator<TSuper, D, U1>
+        + Allocator<TSuper, D, D>,
+{
+    #[inline(always)]
+    fn to_superset(&self) -> DualVec<TSuper, FSuper, D> {
+        let re = TSuper::from_subset(&self.re);
+        let eps = Derivative::from_subset(&self.eps);
+        DualVec {
+            re,
+            eps,
+            f: PhantomData,
+        }
+    }
+    #[inline(always)]
+    fn from_superset(element: &DualVec<TSuper, FSuper, D>) -> Option<Self> {
+        let re = TSuper::to_subset(&element.re)?;
+        let eps = Derivative::to_subset(&element.eps)?;
+        Some(Self {
+            re,
+            eps,
+            f: PhantomData,
+        })
+    }
+    #[inline(always)]
+    fn from_superset_unchecked(element: &DualVec<TSuper, FSuper, D>) -> Self {
+        let re = TSuper::to_subset_unchecked(&element.re);
+        let eps = Derivative::to_subset_unchecked(&element.eps);
+        Self {
+            re,
+            eps,
+            f: PhantomData,
+        }
+    }
+    #[inline(always)]
+    fn is_in_subset(element: &DualVec<TSuper, FSuper, D>) -> bool {
+        TSuper::is_in_subset(&element.re)
+            && <Derivative<_, _, _, _> as SupersetOf<Derivative<_, _, _, _>>>::is_in_subset(
+                &element.eps,
+            )
+    }
+}
+
+impl<TSuper, FSuper, D: Dim> SupersetOf<f32> for DualVec<TSuper, FSuper, D>
+where
+    TSuper: DualNum<FSuper> + SupersetOf<f32>,
+    DefaultAllocator: Allocator<TSuper, D>
+        + Allocator<TSuper, U1, D>
+        + Allocator<TSuper, D, U1>
+        + Allocator<TSuper, D, D>,
+{
+    #[inline(always)]
+    fn is_in_subset(&self) -> bool {
+        self.re.is_in_subset()
+    }
+
+    #[inline(always)]
+    fn to_subset_unchecked(&self) -> f32 {
+        self.re.to_subset_unchecked()
+    }
+
+    #[inline(always)]
+    fn from_subset(element: &f32) -> Self {
+        // Interpret as a purely real number
+        let re = TSuper::from_subset(element);
+        let eps = Derivative::none();
+        Self {
+            re,
+            eps,
+            f: PhantomData,
+        }
+    }
+}
+
+impl<TSuper, FSuper, D: Dim> SupersetOf<f64> for DualVec<TSuper, FSuper, D>
+where
+    TSuper: DualNum<FSuper> + SupersetOf<f64>,
+    DefaultAllocator: Allocator<TSuper, D>
+        + Allocator<TSuper, U1, D>
+        + Allocator<TSuper, D, U1>
+        + Allocator<TSuper, D, D>,
+{
+    #[inline(always)]
+    fn is_in_subset(&self) -> bool {
+        self.re.is_in_subset()
+    }
+
+    #[inline(always)]
+    fn to_subset_unchecked(&self) -> f64 {
+        self.re.to_subset_unchecked()
+    }
+
+    #[inline(always)]
+    fn from_subset(element: &f64) -> Self {
+        // Interpret as a purely real number
+        let re = TSuper::from_subset(element);
+        let eps = Derivative::none();
+        Self {
+            re,
+            eps,
+            f: PhantomData,
+        }
+    }
+}
+
+// We can't do a simd implementation until simba lets us implement SimdPartialOrd
+// using _T_'s SimdBool. The blanket impl gets in the way. So we must constrain
+// T to SimdValue<Element = T, SimdBool = bool>, which is basically the same as
+// saying f32 or f64 only.
+//
+// Limitation of simba. See https://github.com/dimforge/simba/issues/44
+
+use nalgebra::{ComplexField, RealField};
+// This impl is modelled on `impl ComplexField for f32`. The imaginary part is nothing.
+impl<T, D: Dim> ComplexField for DualVec<T, T::Element, D>
+where
+    T: DualNum<T::Element> + SupersetOf<T> + AbsDiffEq<Epsilon = T> + Sync + Send,
+    T::Element: DualNum<T::Element> + Scalar + DualNumFloat + Sync + Send,
+    T: SupersetOf<T::Element>,
+    T: SupersetOf<f64>,
+    T: SimdPartialOrd + PartialOrd,
+    T: SimdValue<Element = T, SimdBool = bool>,
+    T: RelativeEq + UlpsEq + AbsDiffEq,
+    DefaultAllocator:
+        Allocator<T, D> + Allocator<T, U1, D> + Allocator<T, D, U1> + Allocator<T, D, D>,
+    <DefaultAllocator as Allocator<T, D>>::Buffer: Sync + Send,
+{
+    type RealField = Self;
+
+    #[inline]
+    fn from_real(re: Self::RealField) -> Self {
+        re
+    }
+
+    #[inline]
+    fn real(self) -> Self::RealField {
+        self
+    }
+
+    #[inline]
+    fn imaginary(self) -> Self::RealField {
+        Self::zero()
+    }
+
+    #[inline]
+    fn modulus(self) -> Self::RealField {
+        self.abs()
+    }
+
+    #[inline]
+    fn modulus_squared(self) -> Self::RealField {
+        &self * &self
+    }
+
+    #[inline]
+    fn argument(self) -> Self::RealField {
+        Self::zero()
+    }
+
+    #[inline]
+    fn norm1(self) -> Self::RealField {
+        self.abs()
+    }
+
+    #[inline]
+    fn scale(self, factor: Self::RealField) -> Self {
+        self * factor
+    }
+
+    #[inline]
+    fn unscale(self, factor: Self::RealField) -> Self {
+        self / factor
+    }
+
+    #[inline]
+    fn floor(self) -> Self {
+        panic!("called floor() on a dual number")
+    }
+
+    #[inline]
+    fn ceil(self) -> Self {
+        panic!("called ceil() on a dual number")
+    }
+
+    #[inline]
+    fn round(self) -> Self {
+        panic!("called round() on a dual number")
+    }
+
+    #[inline]
+    fn trunc(self) -> Self {
+        panic!("called trunc() on a dual number")
+    }
+
+    #[inline]
+    fn fract(self) -> Self {
+        panic!("called fract() on a dual number")
+    }
+
+    #[inline]
+    fn mul_add(self, a: Self, b: Self) -> Self {
+        DualNum::mul_add(&self, a, b)
+    }
+
+    #[inline]
+    fn abs(self) -> Self::RealField {
+        Signed::abs(&self)
+    }
+
+    #[inline]
+    fn hypot(self, other: Self) -> Self::RealField {
+        let sum_sq = self.powi(2) + other.powi(2);
+        DualNum::sqrt(&sum_sq)
+    }
+
+    #[inline]
+    fn recip(self) -> Self {
+        DualNum::recip(&self)
+    }
+
+    #[inline]
+    fn conjugate(self) -> Self {
+        self
+    }
+
+    #[inline]
+    fn sin(self) -> Self {
+        DualNum::sin(&self)
+    }
+
+    #[inline]
+    fn cos(self) -> Self {
+        DualNum::cos(&self)
+    }
+
+    #[inline]
+    fn sin_cos(self) -> (Self, Self) {
+        DualNum::sin_cos(&self)
+    }
+
+    #[inline]
+    fn tan(self) -> Self {
+        DualNum::tan(&self)
+    }
+
+    #[inline]
+    fn asin(self) -> Self {
+        DualNum::asin(&self)
+    }
+
+    #[inline]
+    fn acos(self) -> Self {
+        DualNum::acos(&self)
+    }
+
+    #[inline]
+    fn atan(self) -> Self {
+        DualNum::atan(&self)
+    }
+
+    #[inline]
+    fn sinh(self) -> Self {
+        DualNum::sinh(&self)
+    }
+
+    #[inline]
+    fn cosh(self) -> Self {
+        DualNum::cosh(&self)
+    }
+
+    #[inline]
+    fn tanh(self) -> Self {
+        DualNum::tanh(&self)
+    }
+
+    #[inline]
+    fn asinh(self) -> Self {
+        DualNum::asinh(&self)
+    }
+
+    #[inline]
+    fn acosh(self) -> Self {
+        DualNum::acosh(&self)
+    }
+
+    #[inline]
+    fn atanh(self) -> Self {
+        DualNum::atanh(&self)
+    }
+
+    #[inline]
+    fn log(self, base: Self::RealField) -> Self {
+        DualNum::ln(&self) / DualNum::ln(&base)
+    }
+
+    #[inline]
+    fn log2(self) -> Self {
+        DualNum::log2(&self)
+    }
+
+    #[inline]
+    fn log10(self) -> Self {
+        DualNum::log10(&self)
+    }
+
+    #[inline]
+    fn ln(self) -> Self {
+        DualNum::ln(&self)
+    }
+
+    #[inline]
+    fn ln_1p(self) -> Self {
+        DualNum::ln_1p(&self)
+    }
+
+    #[inline]
+    fn sqrt(self) -> Self {
+        DualNum::sqrt(&self)
+    }
+
+    #[inline]
+    fn exp(self) -> Self {
+        DualNum::exp(&self)
+    }
+
+    #[inline]
+    fn exp2(self) -> Self {
+        DualNum::exp2(&self)
+    }
+
+    #[inline]
+    fn exp_m1(self) -> Self {
+        DualNum::exp_m1(&self)
+    }
+
+    #[inline]
+    fn powi(self, n: i32) -> Self {
+        DualNum::powi(&self, n)
+    }
+
+    #[inline]
+    fn powf(self, n: Self::RealField) -> Self {
+        // n could be a dual.
+        DualNum::powd(&self, n)
+    }
+
+    #[inline]
+    fn powc(self, n: Self) -> Self {
+        // same as powf, Self isn't complex
+        self.powf(n)
+    }
+
+    #[inline]
+    fn cbrt(self) -> Self {
+        DualNum::cbrt(&self)
+    }
+
+    #[inline]
+    fn is_finite(&self) -> bool {
+        self.re.is_finite()
+    }
+
+    #[inline]
+    fn try_sqrt(self) -> Option<Self> {
+        if self > Self::zero() {
+            Some(DualNum::sqrt(&self))
+        } else {
+            None
+        }
+    }
+}
+
+impl<T, D: Dim> RealField for DualVec<T, T::Element, D>
+where
+    T: DualNum<T::Element> + SupersetOf<T> + Sync + Send,
+    T::Element: DualNum<T::Element> + Scalar + DualNumFloat,
+    T: SupersetOf<T::Element>,
+    T: SupersetOf<f64>,
+    T: SimdPartialOrd + PartialOrd,
+    T: RelativeEq + AbsDiffEq<Epsilon = T>,
+    T: SimdValue<Element = T, SimdBool = bool>,
+    T: UlpsEq,
+    T: AbsDiffEq,
+    DefaultAllocator:
+        Allocator<T, D> + Allocator<T, U1, D> + Allocator<T, D, U1> + Allocator<T, D, D>,
+    <DefaultAllocator as Allocator<T, D>>::Buffer: Sync + Send,
+{
+    #[inline]
+    fn copysign(self, _sign: Self) -> Self {
+        todo!("copysign not yet implemented on dual numbers")
+    }
+
+    #[inline]
+    fn atan2(self, _other: Self) -> Self {
+        todo!("atan2 not yet implemented on dual numbers")
+    }
+
+    #[inline]
+    fn pi() -> Self {
+        Self::from_re(<T as FloatConst>::PI())
+    }
+
+    #[inline]
+    fn two_pi() -> Self {
+        Self::from_re(<T as FloatConst>::TAU())
+    }
+
+    #[inline]
+    fn frac_pi_2() -> Self {
+        Self::from_re(<T as FloatConst>::FRAC_PI_4())
+    }
+
+    #[inline]
+    fn frac_pi_3() -> Self {
+        Self::from_re(<T as FloatConst>::FRAC_PI_3())
+    }
+
+    #[inline]
+    fn frac_pi_4() -> Self {
+        Self::from_re(<T as FloatConst>::FRAC_PI_4())
+    }
+
+    #[inline]
+    fn frac_pi_6() -> Self {
+        Self::from_re(<T as FloatConst>::FRAC_PI_6())
+    }
+
+    #[inline]
+    fn frac_pi_8() -> Self {
+        Self::from_re(<T as FloatConst>::FRAC_PI_8())
+    }
+
+    #[inline]
+    fn frac_1_pi() -> Self {
+        Self::from_re(<T as FloatConst>::FRAC_1_PI())
+    }
+
+    #[inline]
+    fn frac_2_pi() -> Self {
+        Self::from_re(<T as FloatConst>::FRAC_2_PI())
+    }
+
+    #[inline]
+    fn frac_2_sqrt_pi() -> Self {
+        Self::from_re(<T as FloatConst>::FRAC_2_SQRT_PI())
+    }
+
+    #[inline]
+    fn e() -> Self {
+        Self::from_re(<T as FloatConst>::E())
+    }
+
+    #[inline]
+    fn log2_e() -> Self {
+        Self::from_re(<T as FloatConst>::LOG2_E())
+    }
+
+    #[inline]
+    fn log10_e() -> Self {
+        Self::from_re(<T as FloatConst>::LOG10_E())
+    }
+
+    #[inline]
+    fn ln_2() -> Self {
+        Self::from_re(<T as FloatConst>::LN_2())
+    }
+
+    #[inline]
+    fn ln_10() -> Self {
+        Self::from_re(<T as FloatConst>::LN_10())
+    }
+
+    #[inline]
+    fn is_sign_positive(&self) -> bool {
+        self.re.is_sign_positive()
+    }
+
+    #[inline]
+    fn is_sign_negative(&self) -> bool {
+        self.re.is_sign_negative()
+    }
+
+    /// Got to be careful using this, because it throws away the derivatives of the one not chosen
+    #[inline]
+    fn max(self, other: Self) -> Self {
+        if other > self {
+            other
+        } else {
+            self
+        }
+    }
+
+    /// Got to be careful using this, because it throws away the derivatives of the one not chosen
+    #[inline]
+    fn min(self, other: Self) -> Self {
+        if other < self {
+            other
+        } else {
+            self
+        }
+    }
+
+    /// If the min/max values are constants and the clamping has an effect, you lose your gradients.
+    #[inline]
+    fn clamp(self, min: Self, max: Self) -> Self {
+        if self < min {
+            min
+        } else if self > max {
+            max
+        } else {
+            self
+        }
+    }
+
+    #[inline]
+    fn min_value() -> Option<Self> {
+        Some(Self::from_re(T::min_value()))
+    }
+
+    #[inline]
+    fn max_value() -> Option<Self> {
+        Some(Self::from_re(T::max_value()))
+    }
+}
diff --git a/src/lib.rs b/src/lib.rs
index ccbf725..fc73411 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -42,7 +42,7 @@
 #![warn(clippy::all)]
 #![allow(clippy::needless_range_loop)]
 
-use num_traits::{Float, FromPrimitive, Inv, NumAssignOps, NumOps, Signed};
+use num_traits::{Float, FloatConst, FromPrimitive, Inv, NumAssignOps, NumOps, Signed};
 use std::fmt;
 use std::iter::{Product, Sum};
 
@@ -212,11 +212,11 @@ pub trait DualNum<F>:
 
 /// The underlying data type of individual derivatives. Usually f32 or f64.
 pub trait DualNumFloat:
-    Float + FromPrimitive + Signed + fmt::Display + fmt::Debug + Sync + Send + 'static
+    Float + FloatConst + FromPrimitive + Signed + fmt::Display + fmt::Debug + Sync + Send + 'static
 {
 }
 impl<T> DualNumFloat for T where
-    T: Float + FromPrimitive + Signed + fmt::Display + fmt::Debug + Sync + Send + 'static
+    T: Float + FloatConst + FromPrimitive + Signed + fmt::Display + fmt::Debug + Sync + Send + 'static
 {
 }
 
diff --git a/src/macros.rs b/src/macros.rs
index 1f73583..d6baf87 100644
--- a/src/macros.rs
+++ b/src/macros.rs
@@ -557,66 +557,82 @@ macro_rules! impl_float_const {
             $($(DefaultAllocator: Allocator<T, $dim> + Allocator<T, U1, $dim> + Allocator<T, $dim, $dim>,)*
             DefaultAllocator: Allocator<T$(, $dim)*>)?
         {
+            #[allow(non_snake_case)]
             fn E() -> Self {
                 Self::from(F::E())
             }
 
+            #[allow(non_snake_case)]
             fn FRAC_1_PI() -> Self {
                 Self::from(F::FRAC_1_PI())
             }
 
+            #[allow(non_snake_case)]
             fn FRAC_1_SQRT_2() -> Self {
                 Self::from(F::FRAC_1_SQRT_2())
             }
 
+            #[allow(non_snake_case)]
             fn FRAC_2_PI() -> Self {
                 Self::from(F::FRAC_2_PI())
             }
 
+            #[allow(non_snake_case)]
             fn FRAC_2_SQRT_PI() -> Self {
                 Self::from(F::FRAC_2_SQRT_PI())
             }
 
+            #[allow(non_snake_case)]
             fn FRAC_PI_2() -> Self {
                 Self::from(F::FRAC_PI_2())
             }
 
+            #[allow(non_snake_case)]
             fn FRAC_PI_3() -> Self {
                 Self::from(F::FRAC_PI_3())
             }
 
+            #[allow(non_snake_case)]
             fn FRAC_PI_4() -> Self {
                 Self::from(F::FRAC_PI_4())
             }
 
+            #[allow(non_snake_case)]
             fn FRAC_PI_6() -> Self {
                 Self::from(F::FRAC_PI_6())
             }
 
+            #[allow(non_snake_case)]
             fn FRAC_PI_8() -> Self {
                 Self::from(F::FRAC_PI_8())
             }
 
+            #[allow(non_snake_case)]
             fn LN_10() -> Self {
                 Self::from(F::LN_10())
             }
 
+            #[allow(non_snake_case)]
             fn LN_2() -> Self {
                 Self::from(F::LN_2())
             }
 
+            #[allow(non_snake_case)]
             fn LOG10_E() -> Self {
                 Self::from(F::LOG10_E())
             }
 
+            #[allow(non_snake_case)]
             fn LOG2_E() -> Self {
                 Self::from(F::LOG2_E())
             }
 
+            #[allow(non_snake_case)]
             fn PI() -> Self {
                 Self::from(F::PI())
             }
 
+            #[allow(non_snake_case)]
             fn SQRT_2() -> Self {
                 Self::from(F::SQRT_2())
             }


Test patch:
diff --git a/tests/test_dual.rs b/tests/test_dual.rs
index 49c1ef0..9c6eb8b 100644
--- a/tests/test_dual.rs
+++ b/tests/test_dual.rs
@@ -370,3 +370,53 @@ fn test_dual_bessel_j2_4() {
     assert!((res.re - -0.279979741339189).abs() < 1e-12);
     assert!((res.eps.unwrap() - -0.132099570594364).abs() < 1e-12);
 }
+
+mod nalgebra_api {
+    use nalgebra::{Point2, Point3, UnitQuaternion, Vector2, Vector3};
+    use num_dual::*;
+    use num_traits::Zero;
+    use std::f32::consts::*;
+
+    fn unit_circle(t: Dual32) -> Point2<Dual32> {
+        let x_dir = |x: Dual32| Vector2::new(x, Dual32::zero());
+        let y_dir = |y: Dual32| Vector2::new(Dual32::zero(), y);
+        let theta = t * TAU;
+        Point2::from(x_dir(theta.cos()) + y_dir(theta.sin()))
+    }
+
+    // This is testing that you can type-check code that whacks DualVec in
+    // nalgebra structures and tries to use them.
+    #[test]
+    fn use_nalgebra_2d() {
+        // 1 radian around the circle
+        let t = Dual32::from_re(0.25).derivative();
+        let point = unit_circle(t);
+        let real = point.map(|x| x.re);
+        let grad = point.map(|x| x.eps.unwrap());
+        println!("point: {}", point.coords);
+        approx::assert_relative_eq!(real, Point2::new(0., 1.), epsilon = 1e-3);
+        approx::assert_relative_eq!(grad, Point2::new(-TAU, 0.), epsilon = 1e-3);
+    }
+
+    #[test]
+    fn use_nalgebra_3d() {
+        // First one does nothing to the gradient. Still got no y or z direction.
+        let rot1 = UnitQuaternion::from_axis_angle(&Vector3::x_axis(), FRAC_PI_8);
+        // Second one goes pi/8 (22.5deg) about the y axis, which should shift some of the steepness
+        // in the x direction into the z direction, but not much.
+        let rot2 = UnitQuaternion::from_axis_angle(&Vector3::y_axis(), FRAC_PI_8);
+        let rotation = (rot2 * rot1).cast::<Dual32>();
+        let lifted_3d_circle = |t: Dual32| {
+            let xy = unit_circle(t); // [0, 1] with derivatives [-1, 0]
+            Point3::new(xy.x, xy.y, Dual32::zero())
+        };
+        let function = |t: Dual32| rotation * lifted_3d_circle(t);
+        let point = function(Dual32::from_re(0.25).derivative());
+        let real = point.map(|x| x.re);
+        let grad = point.map(|x| x.eps.unwrap());
+        println!("rotated point: {}", point.coords);
+        approx::assert_relative_eq!(real.coords, real.coords.normalize(), epsilon = 1e-3);
+        approx::assert_relative_eq!(real, Point3::new(0.146, 0.924, 0.354), epsilon = 1e-3);
+        approx::assert_relative_eq!(grad, Point3::new(-5.805, 0., 2.404), epsilon = 1e-3);
+    }
+}