diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml
index d3cae4f..5f0d33c 100644
--- a/.github/workflows/rust.yml
+++ b/.github/workflows/rust.yml
@@ -32,3 +32,15 @@ jobs:
       - uses: actions/checkout@v4
       - name: Run tests (ndarray)
         run: cargo test --release --features "ndarray approx"
+  test_nalgebra:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - name: Run tests (nalgebra)
+        run: cargo test --release --features "nalgebra approx"
+  test_num_dual:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - name: Run tests (num-dual)
+        run: cargo test --release --features "num-dual approx"
diff --git a/.github/workflows/test_python.yml b/.github/workflows/test_python.yml
index 18be9ca..ca6ced5 100644
--- a/.github/workflows/test_python.yml
+++ b/.github/workflows/test_python.yml
@@ -35,4 +35,4 @@ jobs:
           maturin build --release --out dist
           pip install si-units --no-index --find-links dist --force-reinstall
       - name: Test with pytest
-        run: pytest example/extend_quantity/test.py
+        run: pytest --doctest-modules example/extend_quantity/test.py
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c657e2a..8e85bc2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,8 +6,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.11.0] - 2025-09-28
+### Added
+- Added core functionalities for quantities with underlying data structures from `nalgebra`. [#92](https://github.com/itt-ustutt/quantity/pull/92)
+- Added automatic differentiation capabilities from the `num-dual` crate. [#92](https://github.com/itt-ustutt/quantity/pull/92)
+
 ## [0.10.6] - 2025-06-24
-## Changed
+### Changed
 - Updated the optional num-dual dependency to 0.11.2. [#91](https://github.com/itt-ustutt/quantity/pull/91)
 
 ## [0.10.5] - 2025-06-03
diff --git a/Cargo.toml b/Cargo.toml
index 53710b6..8aa2dd7 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,12 +1,12 @@
 [package]
 name = "quantity"
-version = "0.10.6"
+version = "0.11.0"
 authors = [
     "Philipp Rehner <prehner@ethz.ch>",
     "Gernot Bauer <bauer@itt.uni-stuttgart.de>",
 ]
-rust-version = "1.81"
-edition = "2021"
+rust-version = "1.87"
+edition = "2024"
 license = "MIT OR Apache-2.0"
 description = "Representation of quantites, i.e. of unit valued scalars and arrays."
 homepage = "https://github.com/itt-ustutt/quantity"
@@ -29,17 +29,20 @@ num-traits = "0.2"
 document-features = "0.2"
 ## Use N-dimensional arrays from the [ndarray] crate as value of a quantity.
 ndarray = { version = "0.16", optional = true }
+## Use dynamic or static arrays from the [nalgebra] crate as value of a quantity.
+nalgebra = { version = "0.34", optional = true }
 approx = { version = "0.5", optional = true }
-pyo3 = { version = "0.23", optional = true }
-numpy = { version = "0.23", optional = true }
-## Use generalized (hyper-)dual numbers from the [num-dual] crate as value of a quantity.
-num-dual = { version = "0.11.2", optional = true }
+pyo3 = { version = "0.26", optional = true }
+numpy = { version = "0.26", optional = true }
+num-dual = { version = "0.12", optional = true }
 
 [features]
 default = []
+## Use generalized (hyper-)dual numbers from the [num-dual] crate as value of a quantity.
+num-dual = ["dep:num-dual", "nalgebra"]
 ## Directly use (scalar) quantities in Python interfaces through [pyo3] and the [si-units](https://pypi.org/project/si-units/) package.
 python = ["pyo3"]
 ## Use scalar and array quantities in Python interfaces through [pyo3], [numpy], and the [si-units](https://pypi.org/project/si-units/) package.
-python_numpy = ["python", "numpy", "ndarray"]
+python_numpy = ["python", "numpy/nalgebra", "ndarray", "nalgebra"]
 ## Enable approximate comparisons through the [approx] crate.
 approx = ["dep:approx", "ndarray?/approx"]
diff --git a/README.md b/README.md
index 78aeac1..631260e 100644
--- a/README.md
+++ b/README.md
@@ -15,7 +15,7 @@ Add this to your `Cargo.toml`:
 
 ```
 [dependencies]
-quantity = "0.10"
+quantity = "0.11"
 ```
 
 ## Examples
@@ -24,7 +24,7 @@ Calculate pressure of an ideal gas.
 
 ```rust
 let temperature = 25.0 * CELSIUS;
-let volume = 1.5 * METER.powi(3);
+let volume = 1.5 * METER.powi::<P3>();
 let moles = 75.0 * MOL;
 let pressure = moles * RGAS * temperature / volume;
 println!("{:.5}", pressure);            // 123.94785 kPa
@@ -36,19 +36,19 @@ Calculate the gravitational pull of the moon on the earth.
 let mass_earth = 5.9724e24 * KILOGRAM;
 let mass_moon = 7.346e22 * KILOGRAM;
 let distance = 383.398 * KILO * METER;
-let force = G * mass_earth * mass_moon / distance.powi(2);
+let force = G * mass_earth * mass_moon / distance.powi::<P2>();
 println!("{:.5e}", force);              // 1.99208e26 N
 ```
 
 Calculate the pressure distribution in the atmosphere using the barometric formula.
 
 ```rust
-let z = SIArray1::linspace(1.0 * METER, 70.0 * KILO * METER, 10)?;
-let g = 9.81 * METER / SECOND.powi(2);
+let z = Quantity::linspace(1.0 * METER, 70.0 * KILO * METER, 10);
+let g = 9.81 * METER / SECOND.powi::<P2>();
 let m = 28.949 * GRAM / MOL;
 let t = 10.0 * CELSIUS;
 let p0 = BAR;
-let pressure = p0 * (-&z * m * g).to_reduced(RGAS * t)?.mapv(f64::exp);
+let pressure = ((-z.clone() * m * g) / (RGAS * t)).mapv(f64::exp) * p0;
 for i in 0..10 {
     println!("z = {:8.5}   p = {:9.5}", z.get(i), pressure.get(i));
 }
diff --git a/example/extend_quantity/Cargo.toml b/example/extend_quantity/Cargo.toml
index 886503e..31fd5de 100644
--- a/example/extend_quantity/Cargo.toml
+++ b/example/extend_quantity/Cargo.toml
@@ -1,13 +1,14 @@
 [package]
 name = "extend_quantity"
 version = "0.1.0"
-edition = "2021"
+edition = "2024"
 
 [lib]
 name = "extend_quantity"
 crate-type = ["cdylib"]
 
 [dependencies]
-pyo3 = { version = "0.23", features = ["extension-module", "abi3-py39"] }
+pyo3 = { version = "0.26", features = ["extension-module", "abi3-py39"] }
 quantity = { version = "*", path = "../../", features = ["python_numpy"] }
 ndarray = "0.16"
+nalgebra = "0.34"
\ No newline at end of file
diff --git a/example/extend_quantity/src/lib.rs b/example/extend_quantity/src/lib.rs
index e8e9384..bd6b677 100644
--- a/example/extend_quantity/src/lib.rs
+++ b/example/extend_quantity/src/lib.rs
@@ -2,6 +2,7 @@ use pyo3::pymodule;
 
 #[pymodule]
 mod extend_quantity {
+    use nalgebra::{DMatrix, DVector};
     use ndarray::Array1;
     use pyo3::pyfunction;
     use quantity::*;
@@ -34,4 +35,12 @@ mod extend_quantity {
     fn law_of_cosines2(a: Length, b: Length, c: Length) -> Angle {
         Angle::acos((a * a + b * b - c * c).convert_into(2.0 * a * b))
     }
+
+    #[pyfunction]
+    fn test_nalgebra(
+        pressure: Pressure<DMatrix<f64>>,
+        volume: Volume<DVector<f64>>,
+    ) -> Energy<DVector<f64>> {
+        pressure * volume
+    }
 }
diff --git a/example/extend_quantity/test.py b/example/extend_quantity/test.py
index 72909f6..63cf441 100644
--- a/example/extend_quantity/test.py
+++ b/example/extend_quantity/test.py
@@ -18,3 +18,14 @@ def test_array():
     >>> ideal_gas_array(T, V, N)
     array([ 826.31900987, 1058.09851279]) Pa
     """
+
+def test_nalgebra():
+    """
+    >>> from extend_quantity import test_nalgebra
+    >>> import si_units as si
+    >>> import numpy as np
+    >>> p = np.array([[1, 2], [3, 4]]) * si.BAR
+    >>> V = np.array([1, 2]) * si.LITER
+    >>> test_nalgebra(p, V)
+    array([ 500., 1100.]) J
+    """
diff --git a/si-units/Cargo.toml b/si-units/Cargo.toml
index bc45d70..ca2dae0 100644
--- a/si-units/Cargo.toml
+++ b/si-units/Cargo.toml
@@ -5,8 +5,8 @@ authors = [
     "Philipp Rehner <prehner@ethz.ch>",
     "Gernot Bauer <bauer@itt.uni-stuttgart.de>",
 ]
-rust-version = "1.81"
-edition = "2021"
+rust-version = "1.87"
+edition = "2024"
 license = "MIT OR Apache-2.0"
 description = "Representation of SI unit valued scalars and arrays."
 homepage = "https://github.com/itt-ustutt/quantity/tree/master/si-units"
@@ -22,7 +22,7 @@ crate-type = ["cdylib"]
 
 [dependencies]
 ndarray = "0.16"
-numpy = "0.23"
-pyo3 = { version = "0.23", features = ["extension-module", "abi3-py39"] }
+numpy = "0.26"
+pyo3 = { version = "0.26", features = ["extension-module", "abi3-py39"] }
 regex = "1.11"
 thiserror = "2.0"
diff --git a/si-units/src/lib.rs b/si-units/src/lib.rs
index 109ee35..39ed6c6 100644
--- a/si-units/src/lib.rs
+++ b/si-units/src/lib.rs
@@ -28,12 +28,12 @@ pub enum QuantityError {
 
 #[pyclass(name = "SIObject", module = "si_units._core", frozen)]
 pub struct PySIObject {
-    value: PyObject,
+    value: Py<PyAny>,
     unit: SIUnit,
 }
 
 impl PySIObject {
-    fn new(value: PyObject, unit: SIUnit) -> Self {
+    fn new(value: Py<PyAny>, unit: SIUnit) -> Self {
         Self { value, unit }
     }
 
@@ -93,7 +93,7 @@ impl PySIObject {
             .ok()
     }
 
-    fn __richcmp__(&self, py: Python, other: &Self, op: CompareOp) -> PyResult<PyObject> {
+    fn __richcmp__(&self, py: Python, other: &Self, op: CompareOp) -> PyResult<Py<PyAny>> {
         self.check_units(other).and_then(|_| match op {
             CompareOp::Eq => self.value.call_method1(py, "__eq__", (&other.value,)),
             CompareOp::Ne => self.value.call_method1(py, "__ne__", (&other.value,)),
@@ -240,7 +240,7 @@ impl PySIObject {
     }
 
     #[getter]
-    fn get_shape(&self, py: Python) -> PyResult<PyObject> {
+    fn get_shape(&self, py: Python) -> PyResult<Py<PyAny>> {
         self.value.getattr(py, "shape")
     }
 
diff --git a/src/ad.rs b/src/ad.rs
new file mode 100644
index 0000000..367bfb2
--- /dev/null
+++ b/src/ad.rs
@@ -0,0 +1,454 @@
+use super::Quantity;
+use nalgebra::{DefaultAllocator, Dim, OMatrix, OVector, U1, allocator::Allocator};
+use num_dual::{
+    Dual, Dual2, Dual2Vec, Dual3, DualNum, DualStruct, DualVec, Gradients, HyperDual, HyperDualVec,
+    HyperHyperDual, Real,
+};
+use std::ops::Sub;
+use typenum::Diff;
+
+impl<D, F, T: DualStruct<D, F>, U> DualStruct<D, F> for Quantity<T, U> {
+    type Real = Quantity<T::Real, U>;
+    type Inner = Quantity<T::Inner, U>;
+
+    fn re(&self) -> Self::Real {
+        Quantity::new(self.0.re())
+    }
+
+    fn from_inner(inner: &Self::Inner) -> Self {
+        Quantity::new(T::from_inner(&inner.0))
+    }
+}
+
+pub fn zeroth_derivative<G, T: DualNum<f64>, UX, UY>(g: G, x: Quantity<T, UX>) -> Quantity<T, UY>
+where
+    G: Fn(Quantity<Real<T, f64>, UX>) -> Quantity<Real<T, f64>, UY>,
+    UY: Sub<UX>,
+{
+    let r = num_dual::zeroth_derivative(|x| g(Quantity::new(x)).0, x.0);
+    Quantity::new(r)
+}
+
+pub fn first_derivative<G, T: DualNum<f64>, UX, UY>(
+    g: G,
+    x: Quantity<T, UX>,
+) -> (Quantity<T, UY>, Quantity<T, Diff<UY, UX>>)
+where
+    G: Fn(Quantity<Dual<T, f64>, UX>) -> Quantity<Dual<T, f64>, UY>,
+    UY: Sub<UX>,
+{
+    let r = num_dual::first_derivative(|x| g(Quantity::new(x)).0, x.0);
+    (Quantity::new(r.0), Quantity::new(r.1))
+}
+
+#[expect(clippy::type_complexity)]
+pub fn gradient<G, T: DualNum<f64>, UX, UY, N: Dim>(
+    g: G,
+    x: &Quantity<OVector<T, N>, UX>,
+) -> (Quantity<T, UY>, Quantity<OVector<T, N>, Diff<UY, UX>>)
+where
+    DefaultAllocator: Allocator<N>,
+    G: Fn(Quantity<OVector<DualVec<T, f64, N>, N>, UX>) -> Quantity<DualVec<T, f64, N>, UY>,
+    UY: Sub<UX>,
+{
+    let r = num_dual::gradient(|x| g(Quantity::new(x)).0, &x.0);
+    (Quantity::new(r.0), Quantity::new(r.1))
+}
+
+#[expect(clippy::type_complexity)]
+pub fn gradient_copy<G, T: DualNum<f64> + Copy, UX, UY, N: Gradients>(
+    g: G,
+    x: &Quantity<OVector<T, N>, UX>,
+) -> (Quantity<T, UY>, Quantity<OVector<T, N>, Diff<UY, UX>>)
+where
+    DefaultAllocator: Allocator<N>,
+    G: Fn(Quantity<OVector<N::Dual<T, f64>, N>, UX>) -> Quantity<N::Dual<T, f64>, UY>,
+    UY: Sub<UX>,
+{
+    let r = N::gradient(|x, _: &()| g(Quantity::new(x)).0, &x.0, &());
+    (Quantity::new(r.0), Quantity::new(r.1))
+}
+
+#[expect(clippy::type_complexity)]
+pub fn jacobian<G, T: DualNum<f64>, UX, UY, M: Dim, N: Dim>(
+    g: G,
+    x: &Quantity<OVector<T, N>, UX>,
+) -> (
+    Quantity<OVector<T, M>, UY>,
+    Quantity<OMatrix<T, M, N>, Diff<UY, UX>>,
+)
+where
+    DefaultAllocator: Allocator<M> + Allocator<N> + Allocator<M, N> + Allocator<U1, N>,
+    G: Fn(
+        Quantity<OVector<DualVec<T, f64, N>, N>, UX>,
+    ) -> Quantity<OVector<DualVec<T, f64, N>, M>, UY>,
+    UY: Sub<UX>,
+{
+    let r = num_dual::jacobian(|x| g(Quantity::new(x)).0, &x.0);
+    (Quantity::new(r.0), Quantity::new(r.1))
+}
+
+#[expect(clippy::type_complexity)]
+pub fn jacobian_copy<G, T: DualNum<f64> + Copy, UX, UY, N: Gradients>(
+    g: G,
+    x: &Quantity<OVector<T, N>, UX>,
+) -> (
+    Quantity<OVector<T, N>, UY>,
+    Quantity<OMatrix<T, N, N>, Diff<UY, UX>>,
+)
+where
+    DefaultAllocator: Allocator<N> + Allocator<N> + Allocator<N, N> + Allocator<U1, N>,
+    G: Fn(Quantity<OVector<N::Dual<T, f64>, N>, UX>) -> Quantity<OVector<N::Dual<T, f64>, N>, UY>,
+    UY: Sub<UX>,
+{
+    let r = N::jacobian(|x, _: &()| g(Quantity::new(x)).0, &x.0, &());
+    (Quantity::new(r.0), Quantity::new(r.1))
+}
+
+#[expect(clippy::type_complexity)]
+pub fn second_derivative<G, T: DualNum<f64>, UX, UY>(
+    g: G,
+    x: Quantity<T, UX>,
+) -> (
+    Quantity<T, UY>,
+    Quantity<T, Diff<UY, UX>>,
+    Quantity<T, Diff<Diff<UY, UX>, UX>>,
+)
+where
+    G: Fn(Quantity<Dual2<T, f64>, UX>) -> Quantity<Dual2<T, f64>, UY>,
+    UY: Sub<UX>,
+    Diff<UY, UX>: Sub<UX>,
+{
+    let r = num_dual::second_derivative(|x| g(Quantity::new(x)).0, x.0);
+    (Quantity::new(r.0), Quantity::new(r.1), Quantity::new(r.2))
+}
+
+#[expect(clippy::type_complexity)]
+pub fn hessian<G, T: DualNum<f64>, UX, UY, N: Dim>(
+    g: G,
+    x: &Quantity<OVector<T, N>, UX>,
+) -> (
+    Quantity<T, UY>,
+    Quantity<OVector<T, N>, Diff<UY, UX>>,
+    Quantity<OMatrix<T, N, N>, Diff<Diff<UY, UX>, UX>>,
+)
+where
+    DefaultAllocator: Allocator<N> + Allocator<U1, N> + Allocator<N, N>,
+    G: Fn(Quantity<OVector<Dual2Vec<T, f64, N>, N>, UX>) -> Quantity<Dual2Vec<T, f64, N>, UY>,
+    UY: Sub<UX>,
+    Diff<UY, UX>: Sub<UX>,
+{
+    let r = num_dual::hessian(|x| g(Quantity::new(x)).0, &x.0);
+    (Quantity::new(r.0), Quantity::new(r.1), Quantity::new(r.2))
+}
+
+#[expect(clippy::type_complexity)]
+pub fn hessian_copy<G, T: DualNum<f64> + Copy, UX, UY, N: Gradients>(
+    g: G,
+    x: &Quantity<OVector<T, N>, UX>,
+) -> (
+    Quantity<T, UY>,
+    Quantity<OVector<T, N>, Diff<UY, UX>>,
+    Quantity<OMatrix<T, N, N>, Diff<Diff<UY, UX>, UX>>,
+)
+where
+    DefaultAllocator: Allocator<N> + Allocator<N, N>,
+    G: Fn(Quantity<OVector<N::Dual2<T, f64>, N>, UX>) -> Quantity<N::Dual2<T, f64>, UY>,
+    UY: Sub<UX>,
+    Diff<UY, UX>: Sub<UX>,
+{
+    let r = N::hessian(|x, _: &()| g(Quantity::new(x)).0, &x.0, &());
+    (Quantity::new(r.0), Quantity::new(r.1), Quantity::new(r.2))
+}
+
+#[expect(clippy::type_complexity)]
+pub fn second_partial_derivative<G, T: DualNum<f64>, UX, UY, UZ>(
+    g: G,
+    (x, y): (Quantity<T, UX>, Quantity<T, UY>),
+) -> (
+    Quantity<T, UZ>,
+    Quantity<T, Diff<UZ, UX>>,
+    Quantity<T, Diff<UZ, UY>>,
+    Quantity<T, Diff<Diff<UZ, UX>, UY>>,
+)
+where
+    G: Fn(
+        (
+            Quantity<HyperDual<T, f64>, UX>,
+            Quantity<HyperDual<T, f64>, UY>,
+        ),
+    ) -> Quantity<HyperDual<T, f64>, UZ>,
+    UZ: Sub<UX>,
+    UZ: Sub<UY>,
+    Diff<UZ, UX>: Sub<UY>,
+{
+    let r = num_dual::second_partial_derivative(
+        |(x, y)| g((Quantity::new(x), Quantity::new(y))).0,
+        (x.0, y.0),
+    );
+    (
+        Quantity::new(r.0),
+        Quantity::new(r.1),
+        Quantity::new(r.2),
+        Quantity::new(r.3),
+    )
+}
+
+#[expect(clippy::type_complexity)]
+pub fn partial_hessian<G, T: DualNum<f64>, UX, UY, UZ, M: Dim, N: Dim>(
+    g: G,
+    (x, y): (&Quantity<OVector<T, M>, UX>, &Quantity<OVector<T, N>, UY>),
+) -> (
+    Quantity<T, UZ>,
+    Quantity<OVector<T, M>, Diff<UZ, UX>>,
+    Quantity<OVector<T, N>, Diff<UZ, UY>>,
+    Quantity<OMatrix<T, M, N>, Diff<Diff<UZ, UX>, UY>>,
+)
+where
+    G: Fn(
+        (
+            Quantity<OVector<HyperDualVec<T, f64, M, N>, M>, UX>,
+            Quantity<OVector<HyperDualVec<T, f64, M, N>, N>, UY>,
+        ),
+    ) -> Quantity<HyperDualVec<T, f64, M, N>, UZ>,
+    DefaultAllocator: Allocator<N> + Allocator<M> + Allocator<M, N> + Allocator<U1, N>,
+    UZ: Sub<UX>,
+    UZ: Sub<UY>,
+    Diff<UZ, UX>: Sub<UY>,
+{
+    let r = num_dual::partial_hessian(
+        |(x, y)| g((Quantity::new(x), Quantity::new(y))).0,
+        (&x.0, &y.0),
+    );
+    (
+        Quantity::new(r.0),
+        Quantity::new(r.1),
+        Quantity::new(r.2),
+        Quantity::new(r.3),
+    )
+}
+
+#[expect(clippy::type_complexity)]
+pub fn partial_hessian_copy<G, T: DualNum<f64> + Copy, UX, UY, UZ, N: Gradients>(
+    g: G,
+    (x, y): (&Quantity<OVector<T, N>, UX>, Quantity<T, UY>),
+) -> (
+    Quantity<T, UZ>,
+    Quantity<OVector<T, N>, Diff<UZ, UX>>,
+    Quantity<T, Diff<UZ, UY>>,
+    Quantity<OVector<T, N>, Diff<Diff<UZ, UX>, UY>>,
+)
+where
+    G: Fn(
+        (
+            Quantity<OVector<N::HyperDual<T, f64>, N>, UX>,
+            Quantity<N::HyperDual<T, f64>, UY>,
+        ),
+    ) -> Quantity<N::HyperDual<T, f64>, UZ>,
+    DefaultAllocator: Allocator<N>,
+    UZ: Sub<UX>,
+    UZ: Sub<UY>,
+    Diff<UZ, UX>: Sub<UY>,
+{
+    let r = N::partial_hessian(
+        |x, y, _: &()| g((Quantity::new(x), Quantity::new(y))).0,
+        &x.0,
+        y.0,
+        &(),
+    );
+    (
+        Quantity::new(r.0),
+        Quantity::new(r.1),
+        Quantity::new(r.2),
+        Quantity::new(r.3),
+    )
+}
+
+#[expect(clippy::type_complexity)]
+pub fn third_derivative<G, T: DualNum<f64>, UX, UY>(
+    g: G,
+    x: Quantity<T, UX>,
+) -> (
+    Quantity<T, UY>,
+    Quantity<T, Diff<UY, UX>>,
+    Quantity<T, Diff<Diff<UY, UX>, UX>>,
+    Quantity<T, Diff<Diff<Diff<UY, UX>, UX>, UX>>,
+)
+where
+    G: Fn(Quantity<Dual3<T, f64>, UX>) -> Quantity<Dual3<T, f64>, UY>,
+    UY: Sub<UX>,
+    Diff<UY, UX>: Sub<UX>,
+    Diff<Diff<UY, UX>, UX>: Sub<UX>,
+{
+    let r = num_dual::third_derivative(|x| g(Quantity::new(x)).0, x.0);
+    (
+        Quantity::new(r.0),
+        Quantity::new(r.1),
+        Quantity::new(r.2),
+        Quantity::new(r.3),
+    )
+}
+
+#[expect(clippy::type_complexity)]
+pub fn third_partial_derivative<G, T: DualNum<f64>, UX, UY, UZ, U>(
+    g: G,
+    (x, y, z): (Quantity<T, UX>, Quantity<T, UY>, Quantity<T, UZ>),
+) -> (
+    Quantity<T, U>,
+    Quantity<T, Diff<U, UX>>,
+    Quantity<T, Diff<U, UY>>,
+    Quantity<T, Diff<U, UZ>>,
+    Quantity<T, Diff<Diff<U, UX>, UY>>,
+    Quantity<T, Diff<Diff<U, UX>, UY>>,
+    Quantity<T, Diff<Diff<U, UX>, UY>>,
+    Quantity<T, Diff<Diff<Diff<U, UX>, UY>, UZ>>,
+)
+where
+    G: Fn(
+        (
+            Quantity<HyperHyperDual<T, f64>, UX>,
+            Quantity<HyperHyperDual<T, f64>, UY>,
+            Quantity<HyperHyperDual<T, f64>, UZ>,
+        ),
+    ) -> Quantity<HyperHyperDual<T, f64>, UY>,
+    U: Sub<UX>,
+    U: Sub<UY>,
+    U: Sub<UZ>,
+    Diff<U, UX>: Sub<UY>,
+    Diff<U, UY>: Sub<UX>,
+    Diff<U, UX>: Sub<UX>,
+    Diff<Diff<U, UX>, UY>: Sub<UZ>,
+{
+    let r = num_dual::third_partial_derivative(
+        |(x, y, z)| g((Quantity::new(x), Quantity::new(y), Quantity::new(z))).0,
+        (x.0, y.0, z.0),
+    );
+    (
+        Quantity::new(r.0),
+        Quantity::new(r.1),
+        Quantity::new(r.2),
+        Quantity::new(r.3),
+        Quantity::new(r.4),
+        Quantity::new(r.5),
+        Quantity::new(r.6),
+        Quantity::new(r.7),
+    )
+}
+
+#[cfg(test)]
+mod test_num_dual {
+    use super::*;
+    use crate::{Area, Length, METER, Temperature, Volume};
+    use approx::assert_relative_eq;
+    use nalgebra::{SMatrix, SVector, vector};
+    use num_dual::{Dual64, ImplicitDerivative, ImplicitFunction};
+    use typenum::{P2, P3};
+
+    struct MyArgs<D> {
+        temperature: Temperature<D>,
+    }
+
+    impl<D: DualNum<f64> + Copy> DualStruct<D, f64> for MyArgs<D> {
+        type Real = MyArgs<f64>;
+        type Inner = MyArgs<D::Inner>;
+
+        fn re(&self) -> Self::Real {
+            MyArgs {
+                temperature: self.temperature.re(),
+            }
+        }
+
+        fn from_inner(inner: &Self::Inner) -> Self {
+            MyArgs {
+                temperature: Temperature::from_inner(&inner.temperature),
+            }
+        }
+    }
+
+    struct AreaImplicit;
+    impl ImplicitFunction<f64> for AreaImplicit {
+        type Parameters<D> = Area<D>;
+        type Variable<D> = D;
+        fn residual<D: DualNum<f64> + Copy>(x: D, &area: &Area<D>) -> D {
+            let l = Length::new(x);
+            (area - l * l).convert_into(area)
+        }
+    }
+
+    #[test]
+    fn test_implicit() {
+        let a = Area::new(Dual64::from(25.0).derivative());
+        let implicit = ImplicitDerivative::new(AreaImplicit, a);
+        let x = implicit.implicit_derivative(5.0);
+        println!("{x}");
+        assert_eq!(x, a.0.sqrt())
+    }
+
+    fn volume<D: DualNum<f64> + Copy>(x: Length<D>) -> Volume<D> {
+        x * x * x
+    }
+
+    fn distance<D: DualNum<f64>, N: Dim>(x: Length<OVector<D, N>>) -> Length<D>
+    where
+        DefaultAllocator: Allocator<N>,
+    {
+        x.dot(&x).sqrt()
+    }
+
+    fn volume2<D: DualNum<f64> + Copy>((x, h): (Length<D>, Length<D>)) -> Volume<D> {
+        x * x * h
+    }
+
+    #[test]
+    fn test_derivative() {
+        let (v, dv) = first_derivative(volume, 5.0 * METER);
+        println!("{v}\t{dv:3}");
+        assert_eq!(v, 125.0 * METER.powi::<P3>());
+        assert_eq!(dv, 75.0 * METER.powi::<P2>());
+
+        let (v, dv, d2v) = second_derivative(volume, 5.0 * METER);
+        println!("{v}\t{dv:3}\t\t{d2v}");
+        assert_eq!(v, 125.0 * METER.powi::<P3>(),);
+        assert_eq!(dv, 75.0 * METER.powi::<P2>(),);
+        assert_eq!(d2v, 30.0 * METER);
+
+        let (v, dv_dx, dv_dh, d2v) =
+            second_partial_derivative(volume2, (5.0 * METER, 20.0 * METER));
+        println!("{v}\t{dv_dx:3}\t{dv_dh:3}\t{d2v}");
+        assert_eq!(v, 500.0 * METER.powi::<P3>(),);
+        assert_eq!(dv_dx, 200.0 * METER.powi::<P2>(),);
+        assert_eq!(dv_dh, 25.0 * METER.powi::<P2>(),);
+        assert_eq!(d2v, 10.0 * METER);
+
+        let (v, dv, d2v, d3v) = third_derivative(volume, 5.0 * METER);
+        println!("{v}\t{dv:3}\t\t{d2v}\t{d3v}");
+        assert_eq!(v, 125.0 * METER.powi::<P3>(),);
+        assert_eq!(dv, 75.0 * METER.powi::<P2>(),);
+        assert_eq!(d2v, 30.0 * METER);
+        assert_eq!(d3v.into_value(), 6.0);
+    }
+
+    #[test]
+    fn test_gradient_and_hessian() {
+        let x = vector![1.0, 5.0, 5.0, 7.0];
+        let (d, grad) = gradient(distance, &Length::new(x));
+        println!("{d}\t{grad:.1}");
+        assert_eq!(d, 10.0 * METER);
+        assert_relative_eq!(grad.into_value(), SVector::from([0.1, 0.5, 0.5, 0.7]));
+
+        let x = vector![1.0, 5.0, 5.0, 7.0];
+        let (d, grad, hess) = hessian(distance, &Length::new(x));
+        println!("{d}\t{grad:.1}\t{hess:.3?}");
+        assert_eq!(d, 10.0 * METER);
+        assert_relative_eq!(grad.into_value(), SVector::from([0.1, 0.5, 0.5, 0.7]));
+        assert_relative_eq!(
+            (hess * METER).into_value(),
+            SMatrix::from([
+                [0.099, -0.005, -0.005, -0.007],
+                [-0.005, 0.075, -0.025, -0.035],
+                [-0.005, -0.025, 0.075, -0.035],
+                [-0.007, -0.035, -0.035, 0.051]
+            ])
+        );
+    }
+}
diff --git a/src/array.rs b/src/array.rs
index 9e1571d..966202b 100644
--- a/src/array.rs
+++ b/src/array.rs
@@ -142,7 +142,7 @@ impl<T, S: Data<Elem = T>, U, D: Dimension> Quantity<ArrayBase<S, D>, U> {
     }
 
     /// Return a producer and iterable that traverses over all 1D lanes pointing in the direction of axis.
-    pub fn lanes_mut(&mut self, axis: Axis) -> LanesMut<T, D::Smaller>
+    pub fn lanes_mut(&mut self, axis: Axis) -> LanesMut<'_, T, D::Smaller>
     where
         S: DataMut,
     {
diff --git a/src/fmt.rs b/src/fmt.rs
index 3e35ef5..6b96098 100644
--- a/src/fmt.rs
+++ b/src/fmt.rs
@@ -4,20 +4,20 @@ use ndarray::{Array, Dimension};
 use std::collections::HashMap;
 use std::fmt;
 use std::sync::LazyLock;
-use typenum::{Quot, N1, N2, N3, P2, P3, P4};
+use typenum::{N1, N2, N3, P2, P3, P4, Quot};
 
 const UNIT_SYMBOLS: [&str; 7] = ["s", "m", "kg", "A", "K", "mol", "cd"];
 
 impl<
-        Inner: fmt::Debug,
-        T: Integer,
-        L: Integer,
-        M: Integer,
-        I: Integer,
-        THETA: Integer,
-        N: Integer,
-        J: Integer,
-    > fmt::Debug for Quantity<Inner, SIUnit<T, L, M, I, THETA, N, J>>
+    Inner: fmt::Debug,
+    T: Integer,
+    L: Integer,
+    M: Integer,
+    I: Integer,
+    THETA: Integer,
+    N: Integer,
+    J: Integer,
+> fmt::Debug for Quantity<Inner, SIUnit<T, L, M, I, THETA, N, J>>
 {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         self.0.fmt(f)?;
diff --git a/src/lib.rs b/src/lib.rs
index b3cd56f..bb9f58f 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -146,11 +146,15 @@
 use ndarray::{Array, ArrayBase, Data, Dimension};
 use std::marker::PhantomData;
 use std::ops::{Deref, Div, Mul};
-use typenum::{ATerm, Diff, Integer, Negate, Quot, Sum, TArr, N1, N2, P1, P3, Z0};
+use typenum::{ATerm, Diff, Integer, N1, N2, Negate, P1, P3, Quot, Sum, TArr, Z0};
 
+#[cfg(feature = "num-dual")]
+pub mod ad;
 #[cfg(feature = "ndarray")]
 mod array;
 mod fmt;
+#[cfg(feature = "nalgebra")]
+mod nalgebra;
 mod ops;
 #[cfg(feature = "python")]
 mod python;
@@ -519,31 +523,6 @@ impl<T> Deref for Dimensionless<T> {
     }
 }
 
-#[cfg(feature = "num-dual")]
-mod num_dual {
-    use super::Quantity;
-    use num_dual::{DualNum, DualStruct};
-
-    impl<D: DualNum<f64>, U> Quantity<D, U> {
-        pub fn re(&self) -> Quantity<f64, U> {
-            Quantity::new(self.0.re())
-        }
-    }
-
-    impl<D, F, T: DualStruct<D, F>, U> DualStruct<D, F> for Quantity<T, U> {
-        type Real = Quantity<T::Real, U>;
-        type Lifted<D2: DualNum<F, Inner = D>> = Quantity<T::Lifted<D2>, U>;
-
-        fn real(&self) -> Self::Real {
-            Quantity::new(self.0.real())
-        }
-
-        fn lift<D2: DualNum<F, Inner = D>>(&self) -> Self::Lifted<D2> {
-            Quantity::new(self.0.lift())
-        }
-    }
-}
-
 #[cfg(test)]
 mod test {
     use super::*;
@@ -555,51 +534,3 @@ mod test {
         assert_eq!(x, 1.0135e5_f64.ln())
     }
 }
-
-#[cfg(test)]
-#[cfg(feature = "num-dual")]
-mod hmm {
-    use super::*;
-    use ::num_dual::{Dual64, DualNum, DualStruct, ImplicitDerivative, ImplicitFunction};
-
-    struct MyArgs<D> {
-        temperature: Temperature<D>,
-    }
-
-    impl<D: DualNum<f64> + Copy> DualStruct<D, f64> for MyArgs<D> {
-        type Real = MyArgs<f64>;
-
-        type Lifted<D2: DualNum<f64, Inner = D>> = MyArgs<D::Lifted<D2>>;
-
-        fn real(&self) -> Self::Real {
-            MyArgs {
-                temperature: self.temperature.real(),
-            }
-        }
-
-        fn lift<D2: DualNum<f64, Inner = D>>(&self) -> Self::Lifted<D2> {
-            MyArgs {
-                temperature: self.temperature.lift(),
-            }
-        }
-    }
-
-    struct AreaImplicit;
-    impl ImplicitFunction<f64> for AreaImplicit {
-        type Parameters<D: DualNum<f64>> = Area<D>;
-        type Variable<D> = D;
-        fn residual<D: DualNum<f64> + Copy>(&area: &Area<D>, x: D) -> D {
-            let l = Length::new(x);
-            (area - l * l).convert_into(area)
-        }
-    }
-
-    #[test]
-    fn test_implicit() {
-        let a = Area::new(Dual64::from(25.0).derivative());
-        let implicit = ImplicitDerivative::new(AreaImplicit, a);
-        let x = implicit.implicit_derivative(5.0);
-        println!("{x}");
-        assert_eq!(x, a.0.sqrt())
-    }
-}
diff --git a/src/nalgebra.rs b/src/nalgebra.rs
new file mode 100644
index 0000000..aa6014d
--- /dev/null
+++ b/src/nalgebra.rs
@@ -0,0 +1,131 @@
+use super::Quantity;
+use nalgebra::allocator::Allocator;
+use nalgebra::constraint::{DimEq, ShapeConstraint};
+use nalgebra::{ClosedAddAssign, ClosedMulAssign, DMatrix, DefaultAllocator, Dim, OMatrix, Scalar};
+use num_traits::Zero;
+use std::marker::PhantomData;
+use std::ops::Add;
+use typenum::Sum;
+
+impl<R: Dim, C: Dim, U, T: Scalar> Quantity<OMatrix<T, R, C>, U>
+where
+    DefaultAllocator: Allocator<R, C>,
+{
+    /// Return the total number of elements in the matrix.
+    pub fn len(&self) -> usize {
+        self.0.len()
+    }
+
+    /// Return whether the matrix has any elements
+    pub fn is_empty(&self) -> bool {
+        self.0.is_empty()
+    }
+
+    /// Return the sum of all elements in the matrix.
+    ///
+    /// # Example
+    /// ```
+    /// # use quantity::BAR;
+    /// # use nalgebra::dvector;
+    /// # use approx::assert_relative_eq;
+    /// let x = dvector![1.5, 2.5] * BAR;
+    /// assert_relative_eq!(x.sum(), &(4.0 * BAR));
+    /// ```
+    pub fn sum(&self) -> Quantity<T, U>
+    where
+        T: Zero + ClosedAddAssign,
+    {
+        Quantity(self.0.sum(), PhantomData)
+    }
+
+    pub fn get(&self, index: usize) -> Quantity<T, U>
+    where
+        T: Copy,
+    {
+        Quantity(self.0[index], PhantomData)
+    }
+
+    pub fn set(&mut self, index: usize, value: Quantity<T, U>) {
+        self.0[index] = value.0;
+    }
+
+    pub fn get2(&self, i: usize, j: usize) -> Quantity<T, U>
+    where
+        T: Copy,
+    {
+        Quantity(self.0[(i, j)], PhantomData)
+    }
+
+    pub fn set2(&mut self, i: usize, j: usize, value: Quantity<T, U>) {
+        self.0[(i, j)] = value.0;
+    }
+
+    pub fn shape_generic(&self) -> (R, C) {
+        self.0.shape_generic()
+    }
+}
+
+impl<T: Scalar, R: Dim, C: Dim, U> Quantity<OMatrix<T, R, C>, U>
+where
+    DefaultAllocator: Allocator<R, C>,
+{
+    pub fn add_scalar(&self, rhs: Quantity<T, U>) -> Self
+    where
+        T: ClosedAddAssign,
+    {
+        Self(self.0.add_scalar(rhs.0), PhantomData)
+    }
+
+    pub fn component_mul<U2>(
+        &self,
+        rhs: &Quantity<OMatrix<T, R, C>, U2>,
+    ) -> Quantity<OMatrix<T, R, C>, Sum<U, U2>>
+    where
+        T: ClosedMulAssign,
+        U: Add<U2>,
+    {
+        Quantity(self.0.component_mul(&rhs.0), PhantomData)
+    }
+
+    pub fn dot<U2, R2: Dim, C2: Dim>(
+        &self,
+        rhs: &Quantity<OMatrix<T, R2, C2>, U2>,
+    ) -> Quantity<T, Sum<U, U2>>
+    where
+        DefaultAllocator: Allocator<R2, C2>,
+        T: Zero + ClosedAddAssign + ClosedMulAssign,
+        U: Add<U2>,
+        ShapeConstraint: DimEq<R, R2> + DimEq<C, C2>,
+    {
+        Quantity(self.0.dot(&rhs.0), PhantomData)
+    }
+
+    pub fn from_fn_generic<F>(nrows: R, ncols: C, mut f: F) -> Self
+    where
+        F: FnMut(usize, usize) -> Quantity<T, U>,
+    {
+        Self(
+            OMatrix::from_fn_generic(nrows, ncols, |i, j| f(i, j).0),
+            PhantomData,
+        )
+    }
+
+    pub fn from_element_generic(nrows: R, ncols: C, elem: Quantity<T, U>) -> Self {
+        Self(
+            OMatrix::from_element_generic(nrows, ncols, elem.0),
+            PhantomData,
+        )
+    }
+}
+
+impl<T: Scalar, U> Quantity<DMatrix<T>, U> {
+    pub fn from_fn<F>(nrows: usize, ncols: usize, mut f: F) -> Self
+    where
+        F: FnMut(usize, usize) -> Quantity<T, U>,
+    {
+        Self(
+            DMatrix::from_fn(nrows, ncols, |i, j| f(i, j).0),
+            PhantomData,
+        )
+    }
+}
diff --git a/src/ops.rs b/src/ops.rs
index 4ef63b0..8aaa9b1 100644
--- a/src/ops.rs
+++ b/src/ops.rs
@@ -1,12 +1,18 @@
 use super::Quantity;
 #[cfg(feature = "approx")]
 use approx::{AbsDiffEq, RelativeEq};
+#[cfg(feature = "nalgebra")]
+use nalgebra::allocator::Allocator;
+#[cfg(feature = "nalgebra")]
+use nalgebra::{DefaultAllocator, Dim, OMatrix};
 #[cfg(feature = "ndarray")]
 use ndarray::{Array, ArrayBase, Data, DataMut, DataOwned, Dimension};
+#[cfg(feature = "num-dual")]
+use num_dual::DualNum;
 use num_traits::{Inv, Signed};
 use std::marker::PhantomData;
 use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
-use typenum::{Diff, Integer, Negate, Prod, Quot, Sum, P2, P3};
+use typenum::{Diff, Integer, Negate, P2, P3, Prod, Quot, Sum};
 
 // Multiplication
 impl<T1, T2, U1, U2> Mul<Quantity<T2, U2>> for Quantity<T1, U1>
@@ -98,6 +104,28 @@ impl<U, S: DataOwned<Elem = f64> + DataMut, D: Dimension> Mul<Quantity<f64, U>>
     }
 }
 
+#[cfg(feature = "nalgebra")]
+impl<U, R: Dim, C: Dim> Mul<Quantity<f64, U>> for &OMatrix<f64, R, C>
+where
+    DefaultAllocator: Allocator<R, C>,
+{
+    type Output = Quantity<OMatrix<f64, R, C>, U>;
+    fn mul(self, other: Quantity<f64, U>) -> Self::Output {
+        Quantity(self * other.0, PhantomData)
+    }
+}
+
+#[cfg(feature = "nalgebra")]
+impl<U, R: Dim, C: Dim> Mul<Quantity<f64, U>> for OMatrix<f64, R, C>
+where
+    DefaultAllocator: Allocator<R, C>,
+{
+    type Output = Quantity<OMatrix<f64, R, C>, U>;
+    fn mul(self, other: Quantity<f64, U>) -> Self::Output {
+        Quantity(self * other.0, PhantomData)
+    }
+}
+
 impl<U, T1, T2> MulAssign<T2> for Quantity<T1, U>
 where
     T1: MulAssign<T2>,
@@ -336,6 +364,7 @@ where
     }
 }
 
+#[cfg(not(feature = "num-dual"))]
 impl<U> Quantity<f64, U> {
     /// Calculate the integer power of self.
     ///
@@ -404,6 +433,75 @@ impl<U> Quantity<f64, U> {
     }
 }
 
+#[cfg(feature = "num-dual")]
+impl<D: DualNum<f64>, U> Quantity<D, U> {
+    /// Calculate the integer power of self.
+    ///
+    /// # Example
+    /// ```
+    /// # use quantity::METER;
+    /// # use approx::assert_relative_eq;
+    /// # use typenum::P2;
+    /// let x = 3.0 * METER;
+    /// assert_relative_eq!(x.powi::<P2>(), 9.0 * METER * METER);
+    /// ```
+    pub fn powi<E: Integer>(self) -> Quantity<D, Prod<U, E>>
+    where
+        U: Mul<E>,
+    {
+        Quantity(self.0.powi(E::I32), PhantomData)
+    }
+
+    /// Calculate the square root of self.
+    ///
+    /// # Example
+    /// ```
+    /// # use quantity::METER;
+    /// # use approx::assert_relative_eq;
+    /// let x = 9.0 * METER * METER;
+    /// assert_relative_eq!(x.sqrt(), 3.0 * METER);
+    /// ```
+    pub fn sqrt(self) -> Quantity<D, Quot<U, P2>>
+    where
+        U: Div<P2>,
+    {
+        Quantity(self.0.sqrt(), PhantomData)
+    }
+
+    /// Calculate the cubic root of self.
+    ///
+    /// # Example
+    /// ```
+    /// # use quantity::METER;
+    /// # use approx::assert_relative_eq;
+    /// let x = 27.0 * METER * METER * METER;
+    /// assert_relative_eq!(x.cbrt(), 3.0 * METER);
+    /// ```
+    pub fn cbrt(self) -> Quantity<D, Quot<U, P3>>
+    where
+        U: Div<P3>,
+    {
+        Quantity(self.0.cbrt(), PhantomData)
+    }
+
+    /// Calculate the integer root of self.
+    ///
+    /// # Example
+    /// ```
+    /// # use quantity::METER;
+    /// # use approx::assert_relative_eq;
+    /// # use typenum::P4;
+    /// let x = 81.0 * METER * METER * METER * METER;
+    /// assert_relative_eq!(x.root::<P4>(), 3.0 * METER);
+    /// ```
+    pub fn root<R: Integer>(self) -> Quantity<D, Quot<U, R>>
+    where
+        U: Div<R>,
+    {
+        Quantity(self.0.powf(1.0 / R::I32 as f64), PhantomData)
+    }
+}
+
 impl<T, U> Quantity<T, U> {
     /// Return the absolute value of `self`.
     ///
diff --git a/src/python.rs b/src/python.rs
index 226ed3b..c1d145c 100644
--- a/src/python.rs
+++ b/src/python.rs
@@ -1,15 +1,21 @@
 use super::{Angle, Quantity, SIUnit};
 use crate::fmt::PrintUnit;
+#[cfg(feature = "nalgebra")]
+use nalgebra::{DMatrix, DVector, Dyn};
 #[cfg(feature = "ndarray")]
 use ndarray::{Array, Dimension};
 #[cfg(feature = "ndarray")]
-use numpy::{IntoPyArray, PyReadonlyArray};
+use numpy::IntoPyArray;
+#[cfg(any(feature = "nalgebra", feature = "ndarray"))]
+use numpy::PyReadonlyArray;
+#[cfg(feature = "nalgebra")]
+use numpy::{PyReadonlyArray1, PyReadonlyArray2, ToPyArray};
 use pyo3::{exceptions::PyValueError, prelude::*};
 use std::{marker::PhantomData, sync::LazyLock};
 use typenum::Integer;
 
-static SIOBJECT: LazyLock<PyObject> = LazyLock::new(|| {
-    Python::with_gil(|py| {
+static SIOBJECT: LazyLock<Py<PyAny>> = LazyLock::new(|| {
+    Python::attach(|py| {
         PyModule::import(py, "si_units")
             .unwrap()
             .getattr("SIObject")
@@ -18,16 +24,8 @@ static SIOBJECT: LazyLock<PyObject> = LazyLock::new(|| {
     })
 });
 
-impl<
-        'py,
-        T: Integer,
-        L: Integer,
-        M: Integer,
-        I: Integer,
-        THETA: Integer,
-        N: Integer,
-        J: Integer,
-    > IntoPyObject<'py> for Quantity<f64, SIUnit<T, L, M, I, THETA, N, J>>
+impl<'py, T: Integer, L: Integer, M: Integer, I: Integer, THETA: Integer, N: Integer, J: Integer>
+    IntoPyObject<'py> for Quantity<f64, SIUnit<T, L, M, I, THETA, N, J>>
 {
     type Target = PyAny;
     type Output = Bound<'py, PyAny>;
@@ -41,16 +39,16 @@ impl<
 
 #[cfg(feature = "ndarray")]
 impl<
-        'py,
-        T: Integer,
-        L: Integer,
-        M: Integer,
-        I: Integer,
-        THETA: Integer,
-        N: Integer,
-        J: Integer,
-        D: Dimension,
-    > IntoPyObject<'py> for Quantity<Array<f64, D>, SIUnit<T, L, M, I, THETA, N, J>>
+    'py,
+    T: Integer,
+    L: Integer,
+    M: Integer,
+    I: Integer,
+    THETA: Integer,
+    N: Integer,
+    J: Integer,
+    D: Dimension,
+> IntoPyObject<'py> for Quantity<Array<f64, D>, SIUnit<T, L, M, I, THETA, N, J>>
 {
     type Target = PyAny;
     type Output = Bound<'py, PyAny>;
@@ -63,16 +61,38 @@ impl<
     }
 }
 
-impl<
-        'py,
-        T: Integer,
-        L: Integer,
-        M: Integer,
-        I: Integer,
-        THETA: Integer,
-        N: Integer,
-        J: Integer,
-    > FromPyObject<'py> for Quantity<f64, SIUnit<T, L, M, I, THETA, N, J>>
+#[cfg(feature = "nalgebra")]
+impl<'py, T: Integer, L: Integer, M: Integer, I: Integer, THETA: Integer, N: Integer, J: Integer>
+    IntoPyObject<'py> for Quantity<DMatrix<f64>, SIUnit<T, L, M, I, THETA, N, J>>
+{
+    type Target = PyAny;
+    type Output = Bound<'py, PyAny>;
+    type Error = PyErr;
+
+    fn into_pyobject(self, py: Python<'py>) -> PyResult<Bound<'py, PyAny>> {
+        let unit = [L::I8, M::I8, T::I8, I::I8, N::I8, THETA::I8, J::I8];
+        let value = self.0.to_pyarray(py).into_any();
+        SIOBJECT.bind(py).call1((value, unit))
+    }
+}
+
+#[cfg(feature = "nalgebra")]
+impl<'py, T: Integer, L: Integer, M: Integer, I: Integer, THETA: Integer, N: Integer, J: Integer>
+    IntoPyObject<'py> for Quantity<DVector<f64>, SIUnit<T, L, M, I, THETA, N, J>>
+{
+    type Target = PyAny;
+    type Output = Bound<'py, PyAny>;
+    type Error = PyErr;
+
+    fn into_pyobject(self, py: Python<'py>) -> PyResult<Bound<'py, PyAny>> {
+        let unit = [L::I8, M::I8, T::I8, I::I8, N::I8, THETA::I8, J::I8];
+        let value = numpy::PyArray1::from_slice(py, self.0.data.as_vec()).into_any();
+        SIOBJECT.bind(py).call1((value, unit))
+    }
+}
+
+impl<'py, T: Integer, L: Integer, M: Integer, I: Integer, THETA: Integer, N: Integer, J: Integer>
+    FromPyObject<'py> for Quantity<f64, SIUnit<T, L, M, I, THETA, N, J>>
 where
     Self: PrintUnit,
 {
@@ -102,16 +122,16 @@ where
 
 #[cfg(feature = "ndarray")]
 impl<
-        'py,
-        T: Integer,
-        L: Integer,
-        M: Integer,
-        I: Integer,
-        THETA: Integer,
-        N: Integer,
-        J: Integer,
-        D: Dimension,
-    > FromPyObject<'py> for Quantity<Array<f64, D>, SIUnit<T, L, M, I, THETA, N, J>>
+    'py,
+    T: Integer,
+    L: Integer,
+    M: Integer,
+    I: Integer,
+    THETA: Integer,
+    N: Integer,
+    J: Integer,
+    D: Dimension,
+> FromPyObject<'py> for Quantity<Array<f64, D>, SIUnit<T, L, M, I, THETA, N, J>>
 where
     Self: PrintUnit,
 {
@@ -140,8 +160,74 @@ where
     }
 }
 
-static ANGLE: LazyLock<PyObject> = LazyLock::new(|| {
-    Python::with_gil(|py| {
+#[cfg(feature = "nalgebra")]
+impl<'py, T: Integer, L: Integer, M: Integer, I: Integer, THETA: Integer, N: Integer, J: Integer>
+    FromPyObject<'py> for Quantity<DVector<f64>, SIUnit<T, L, M, I, THETA, N, J>>
+where
+    Self: PrintUnit,
+{
+    fn extract_bound(ob: &Bound<'py, PyAny>) -> PyResult<Self> {
+        let Ok((value, unit_from)) = ob.call_method0("__getnewargs__").and_then(|raw| {
+            raw.extract::<(PyReadonlyArray1<f64>, [i8; 7])>()
+                .map(|(m, u)| {
+                    let m: nalgebra::DVectorView<f64, Dyn, Dyn> = m.try_as_matrix().unwrap();
+                    (m.clone_owned(), u)
+                })
+        }) else {
+            return Err(PyErr::new::<PyValueError, _>(format!(
+                "Missing units! Expected {}, got {}.",
+                Self::UNIT,
+                ob.call_method0("__repr__")?
+            )));
+        };
+        let unit_into = [L::I8, M::I8, T::I8, I::I8, N::I8, THETA::I8, J::I8];
+        if unit_into == unit_from {
+            Ok(Quantity(value, PhantomData))
+        } else {
+            Err(PyErr::new::<PyValueError, _>(format!(
+                "Wrong units! Expected {}, got {}.",
+                Self::UNIT,
+                ob.call_method0("__repr__")?
+            )))
+        }
+    }
+}
+
+#[cfg(feature = "nalgebra")]
+impl<'py, T: Integer, L: Integer, M: Integer, I: Integer, THETA: Integer, N: Integer, J: Integer>
+    FromPyObject<'py> for Quantity<DMatrix<f64>, SIUnit<T, L, M, I, THETA, N, J>>
+where
+    Self: PrintUnit,
+{
+    fn extract_bound(ob: &Bound<'py, PyAny>) -> PyResult<Self> {
+        let Ok((value, unit_from)) = ob.call_method0("__getnewargs__").and_then(|raw| {
+            raw.extract::<(PyReadonlyArray2<f64>, [i8; 7])>()
+                .map(|(m, u)| {
+                    let m: nalgebra::DMatrixView<f64, Dyn, Dyn> = m.try_as_matrix().unwrap();
+                    (m.clone_owned(), u)
+                })
+        }) else {
+            return Err(PyErr::new::<PyValueError, _>(format!(
+                "Missing units! Expected {}, got {}.",
+                Self::UNIT,
+                ob.call_method0("__repr__")?
+            )));
+        };
+        let unit_into = [L::I8, M::I8, T::I8, I::I8, N::I8, THETA::I8, J::I8];
+        if unit_into == unit_from {
+            Ok(Quantity(value, PhantomData))
+        } else {
+            Err(PyErr::new::<PyValueError, _>(format!(
+                "Wrong units! Expected {}, got {}.",
+                Self::UNIT,
+                ob.call_method0("__repr__")?
+            )))
+        }
+    }
+}
+
+static ANGLE: LazyLock<Py<PyAny>> = LazyLock::new(|| {
+    Python::attach(|py| {
         PyModule::import(py, "si_units")
             .unwrap()
             .getattr("Angle")