diff --git a/pixi.lock b/pixi.lock
index 0255d776..22b52fbe 100644
--- a/pixi.lock
+++ b/pixi.lock
@@ -3,8 +3,6 @@ environments:
   benchmark:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-7_kmp_llvm.conda
@@ -84,6 +82,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/matplotlib-base-3.10.8-py314h1194b4b_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mkl-2025.3.0-h0e700b2_463.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mkl-devel-2025.3.0-ha770c72_463.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mkl-include-2025.3.0-hf2ce2f3_463.conda
@@ -104,6 +103,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.14-8_cp314.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/qhull-2020.2-h434a139_5.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.17.0-py314hf07bd8e_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/seaborn-base-0.13.2-pyhd8ed1ab_3.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -202,6 +203,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/matplotlib-base-3.10.8-py314hd63e3f0_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/munkres-1.1.4-pyhd8ed1ab_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ncurses-6.5-h5e97a16_3.conda
@@ -219,6 +221,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.14-8_cp314.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/qhull-2020.2-h420ef59_5.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h46df422_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.17.0-py314hfc1f868_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/seaborn-base-0.13.2-pyhd8ed1ab_3.conda
@@ -298,6 +302,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/matplotlib-base-3.10.8-py314hfa45d96_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-2025.3.0-hac47afa_455.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-devel-2025.3.0-h57928b3_455.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-include-2025.3.0-h57928b3_455.conda
@@ -317,6 +322,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python-tzdata-2025.3-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.14-8_cp314.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/qhull-2020.2-hc790b64_5.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.17.0-py314h221f224_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/seaborn-base-0.13.2-pyhd8ed1ab_3.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -346,8 +353,6 @@ environments:
   default:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -472,6 +477,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/matplotlib-inline-0.2.1-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mypy-1.19.1-py314h5bd0f2a_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -509,6 +515,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/pyzmq-27.1.0-py312hfb55c3c_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/re2-2025.11.05-h5301d42_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.2-he8a4886_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.17.0-py314hf07bd8e_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -657,6 +665,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/matplotlib-inline-0.2.1-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mypy-1.19.1-py314hbdd0d06_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -694,6 +703,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/pyzmq-27.1.0-py312hd65ceae_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/re2-2025.11.05-h64b956e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h46df422_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.17.0-py314hfc1f868_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -806,6 +817,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/matplotlib-inline-0.2.1-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-2025.3.0-hac47afa_455.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mypy-1.19.1-py314h5a2d7ad_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
@@ -840,6 +852,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/pywin32-311-py314h8f8f202_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/pyzmq-27.1.0-py312hbb5da91_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/re2-2025.11.05-ha104f34_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.17.0-py314h221f224_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-9.2.2-pyhd8ed1ab_0.conda
@@ -871,8 +885,6 @@ environments:
   docs:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -1011,6 +1023,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/matplotlib-base-3.10.8-py314h1194b4b_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/matplotlib-inline-0.2.1-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mdit-py-plugins-0.5.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mdurl-0.1.2-pyhd8ed1ab_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mistune-3.2.0-pyhcf101f3_0.conda
@@ -1066,6 +1079,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/roman-numerals-4.1.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/roman-numerals-py-4.1.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/rpds-py-0.30.0-py314h2e6c369_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.17.0-py314hf07bd8e_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/seaborn-base-0.13.2-pyhd8ed1ab_3.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/send2trash-2.1.0-pyha191276_0.conda
@@ -1257,6 +1272,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/matplotlib-base-3.10.8-py314hd63e3f0_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/matplotlib-inline-0.2.1-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mdit-py-plugins-0.5.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mdurl-0.1.2-pyhd8ed1ab_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mistune-3.2.0-pyhcf101f3_0.conda
@@ -1314,6 +1330,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/roman-numerals-4.1.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/roman-numerals-py-4.1.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rpds-py-0.30.0-py314haad56a0_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.17.0-py314hfc1f868_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/seaborn-base-0.13.2-pyhd8ed1ab_3.conda
@@ -1481,6 +1499,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/matplotlib-base-3.10.8-py314hfa45d96_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/matplotlib-inline-0.2.1-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mdit-py-plugins-0.5.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mdurl-0.1.2-pyhd8ed1ab_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mistune-3.2.0-pyhcf101f3_0.conda
@@ -1535,6 +1554,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/roman-numerals-4.1.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/roman-numerals-py-4.1.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/rpds-py-0.30.0-py314h9f07db2_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.17.0-py314h221f224_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/seaborn-base-0.13.2-pyhd8ed1ab_3.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/send2trash-2.1.0-pyh6dadd2b_0.conda
@@ -1595,8 +1616,6 @@ environments:
   lint:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -1727,8 +1746,6 @@ environments:
   nightly:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -1840,6 +1857,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/make-4.4.1-hb9d3cd8_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mypy-1.19.1-py314h5bd0f2a_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -1869,6 +1887,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.14-8_cp314.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/re2-2025.11.05-h5301d42_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.2-he8a4886_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.17.0-py314hf07bd8e_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -1998,6 +2018,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/make-4.4.1-hc9fafa5_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mypy-1.19.1-py314hbdd0d06_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -2027,6 +2048,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.14-8_cp314.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/re2-2025.11.05-h64b956e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h46df422_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.17.0-py314hfc1f868_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -2121,6 +2144,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/make-4.4.1-h0e40799_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-2025.3.0-hac47afa_455.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mypy-1.19.1-py314h5a2d7ad_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
@@ -2148,6 +2172,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python-tzdata-2025.3-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.14-8_cp314.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/re2-2025.11.05-ha104f34_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.17.0-py314h221f224_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-9.2.2-pyhd8ed1ab_0.conda
@@ -2174,8 +2200,6 @@ environments:
   oldies:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -2289,6 +2313,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/make-4.4.1-hb9d3cd8_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/markupsafe-3.0.3-py310h3406613_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mypy-1.19.1-py310h7c4b9e2_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.0.1-pyhe01879c_0.conda
@@ -2319,6 +2344,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/pytz-2025.2-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/re2-2025.11.05-h5301d42_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.2-he8a4886_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.7.3-py310hea5193d_0.tar.bz2
       - conda: https://conda.anaconda.org/conda-forge/linux-64/setuptools-64.0.3-py310hff52083_0.tar.bz2
@@ -2447,6 +2474,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/make-4.4.1-hc9fafa5_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/markupsafe-3.0.3-py310hf4fd40f_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mypy-1.19.1-py310h230e4be_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.0.1-pyhe01879c_0.conda
@@ -2477,6 +2505,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/pytz-2025.2-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/re2-2025.11.05-h64b956e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h46df422_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.7.3-py310h6ecf4ae_0.tar.bz2
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/setuptools-64.0.3-py310hbe9552e_0.tar.bz2
@@ -2569,6 +2599,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/make-4.4.1-h013a479_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/markupsafe-3.0.3-py310hdb0e946_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-2024.2.2-h57928b3_16.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/msys2-conda-epoch-20160418-1.tar.bz2
       - conda: https://conda.anaconda.org/conda-forge/win-64/mypy-1.19.1-py310h29418f3_0.conda
@@ -2598,6 +2629,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.10-8_cp310.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/pytz-2025.2-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/re2-2025.11.05-ha104f34_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.7.3-py310h33db832_0.tar.bz2
       - conda: https://conda.anaconda.org/conda-forge/win-64/setuptools-64.0.3-py310h5588dad_0.tar.bz2
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-9.2.2-pyhd8ed1ab_0.conda
@@ -2624,8 +2657,6 @@ environments:
   py310:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -2738,6 +2769,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/make-4.4.1-hb9d3cd8_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/markupsafe-3.0.3-py310h3406613_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mypy-1.19.1-py310h7c4b9e2_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -2769,6 +2801,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/pytz-2025.2-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/re2-2025.11.05-h5301d42_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.2-he8a4886_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.15.2-py310h1d65ade_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -2897,6 +2931,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/make-4.4.1-hc9fafa5_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/markupsafe-3.0.3-py310hf4fd40f_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mypy-1.19.1-py310h230e4be_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -2928,6 +2963,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/pytz-2025.2-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/re2-2025.11.05-h64b956e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h46df422_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.15.2-py310h32ab4ed_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -3021,6 +3058,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/make-4.4.1-h0e40799_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/markupsafe-3.0.3-py310hdb0e946_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-2025.3.0-hac47afa_455.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mypy-1.19.1-py310h29418f3_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
@@ -3049,6 +3087,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.10-8_cp310.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/pytz-2025.2-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/re2-2025.11.05-ha104f34_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.15.2-py310h15c175c_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-9.2.2-pyhd8ed1ab_0.conda
@@ -3075,8 +3115,6 @@ environments:
   py311:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -3189,6 +3227,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/make-4.4.1-hb9d3cd8_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/markupsafe-3.0.3-py311h3778330_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mypy-1.19.1-py311h49ec1c0_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -3218,6 +3257,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.11-8_cp311.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/re2-2025.11.05-h5301d42_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.2-he8a4886_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.17.0-py311hbe70eeb_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -3346,6 +3387,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/make-4.4.1-hc9fafa5_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/markupsafe-3.0.3-py311ha9b3269_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mypy-1.19.1-py311h8b270aa_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -3375,6 +3417,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.11-8_cp311.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/re2-2025.11.05-h64b956e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h46df422_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.17.0-py311he9931d0_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -3468,6 +3512,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/make-4.4.1-h0e40799_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/markupsafe-3.0.3-py311h3f79411_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-2025.3.0-hac47afa_455.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mypy-1.19.1-py311h3485c13_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
@@ -3495,6 +3540,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python-tzdata-2025.3-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.11-8_cp311.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/re2-2025.11.05-ha104f34_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.17.0-py311h9c22a71_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-9.2.2-pyhd8ed1ab_0.conda
@@ -3521,8 +3568,6 @@ environments:
   py312:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -3635,6 +3680,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/make-4.4.1-hb9d3cd8_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/markupsafe-3.0.3-py312h8a5da7c_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mypy-1.19.1-py312h4c3975b_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -3664,6 +3710,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.12-8_cp312.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/re2-2025.11.05-h5301d42_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.2-he8a4886_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.17.0-py312h54fa4ab_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -3792,6 +3840,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/make-4.4.1-hc9fafa5_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/markupsafe-3.0.3-py312h5748b74_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mypy-1.19.1-py312hefc2c51_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -3821,6 +3870,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.12-8_cp312.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/re2-2025.11.05-h64b956e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h46df422_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.17.0-py312h0f234b1_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -3914,6 +3965,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/make-4.4.1-h0e40799_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/markupsafe-3.0.3-py312h05f76fc_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-2025.3.0-hac47afa_455.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mypy-1.19.1-py312he06e257_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
@@ -3941,6 +3993,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python-tzdata-2025.3-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.12-8_cp312.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/re2-2025.11.05-ha104f34_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.17.0-py312h9b3c559_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-9.2.2-pyhd8ed1ab_0.conda
@@ -3967,8 +4021,6 @@ environments:
   py313:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -4080,6 +4132,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/make-4.4.1-hb9d3cd8_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/markupsafe-3.0.3-py313h3dea7bd_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mypy-1.19.1-py313h07c4f96_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -4109,6 +4162,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.13-8_cp313.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/re2-2025.11.05-h5301d42_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.2-he8a4886_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.17.0-py313h4b8bb8b_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -4238,6 +4293,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/make-4.4.1-hc9fafa5_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/markupsafe-3.0.3-py313h7d74516_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mypy-1.19.1-py313hd3e6d80_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -4267,6 +4323,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.13-8_cp313.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/re2-2025.11.05-h64b956e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h46df422_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.17.0-py313hc753a45_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -4361,6 +4419,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/make-4.4.1-h0e40799_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/markupsafe-3.0.3-py313hd650c13_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-2025.3.0-hac47afa_455.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mypy-1.19.1-py313h5ea7bf4_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
@@ -4388,6 +4447,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python-tzdata-2025.3-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.13-8_cp313.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/re2-2025.11.05-ha104f34_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.17.0-py313he51e9a2_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-9.2.2-pyhd8ed1ab_0.conda
@@ -4414,8 +4475,6 @@ environments:
   py314:
     channels:
     - url: https://conda.anaconda.org/conda-forge/
-    options:
-      pypi-prerelease-mode: if-necessary-or-explicit
     packages:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
@@ -4527,6 +4586,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/make-4.4.1-hb9d3cd8_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/mypy-1.19.1-py314h5bd0f2a_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -4556,6 +4616,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.14-8_cp314.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/re2-2025.11.05-h5301d42_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.2-he8a4886_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/scipy-1.17.0-py314hf07bd8e_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -4685,6 +4747,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/make-4.4.1-hc9fafa5_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mypy-1.19.1-py314hbdd0d06_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/narwhals-2.15.0-pyhcf101f3_0.conda
@@ -4714,6 +4777,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.14-8_cp314.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/re2-2025.11.05-h64b956e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/readline-8.3-h46df422_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scipy-1.17.0-py314hfc1f868_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/sdkroot_env_osx-arm64-26.0-ha3f98da_6.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
@@ -4808,6 +4873,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/make-4.4.1-h0e40799_2.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mako-1.3.10-pyhcf101f3_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/markupsafe-3.0.3-pyh7db6752_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mkl-2025.3.0-hac47afa_455.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/mypy-1.19.1-py314h5a2d7ad_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/mypy_extensions-1.1.0-pyha770c72_0.conda
@@ -4835,6 +4901,8 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/noarch/python-tzdata-2025.3-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/python_abi-3.14-8_cp314.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/re2-2025.11.05-ha104f34_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/scipy-1.17.0-py314h221f224_1.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-80.10.2-pyh332efcf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/setuptools-scm-9.2.2-pyhd8ed1ab_0.conda
@@ -10727,6 +10795,51 @@ packages:
   license_family: BSD
   size: 15175
   timestamp: 1761214578417
+- conda: https://conda.anaconda.org/conda-forge/linux-64/maturin-1.11.5-py310hb4e1661_0.conda
+  noarch: python
+  sha256: 5123ae1d622a2d389f9af3a2212cc2cfe1a75b4e298e2bfaaeedec9da2e2a64e
+  md5: 3fc6cd19772b0033484de8d866dbf54d
+  depends:
+  - python
+  - tomli >=1.1.0
+  - libgcc >=14
+  - __glibc >=2.17,<3.0.a0
+  - openssl >=3.5.4,<4.0a0
+  constrains:
+  - __glibc >=2.17
+  license: MIT
+  license_family: MIT
+  size: 7582313
+  timestamp: 1767965835519
+- conda: https://conda.anaconda.org/conda-forge/osx-arm64/maturin-1.11.5-py310h1939404_0.conda
+  noarch: python
+  sha256: 74ffa13f5da5091962f2a23c81db8efcfc993b121c85f70006899091af66aab1
+  md5: aa879fb2ce3ff56518d54dc1e4f997a8
+  depends:
+  - python
+  - tomli >=1.1.0
+  - __osx >=11.0
+  - openssl >=3.5.4,<4.0a0
+  constrains:
+  - __osx >=11.0
+  license: MIT
+  license_family: MIT
+  size: 6851352
+  timestamp: 1767965954110
+- conda: https://conda.anaconda.org/conda-forge/win-64/maturin-1.11.5-py310h24db72e_0.conda
+  noarch: python
+  sha256: 65b55995de6f0420d6be61bd8f0585ee66050a710d53450ab5ea2d829b902d9c
+  md5: b58155ded427d3d4833503255d55b3b2
+  depends:
+  - python
+  - tomli >=1.1.0
+  - vc >=14.3,<15
+  - vc14_runtime >=14.44.35208
+  - ucrt >=10.0.20348.0
+  license: MIT
+  license_family: MIT
+  size: 6886894
+  timestamp: 1767965856695
 - conda: https://conda.anaconda.org/conda-forge/noarch/mdit-py-plugins-0.5.0-pyhd8ed1ab_0.conda
   sha256: 123cc004e2946879708cdb6a9eff24acbbb054990d6131bb94bca7a374ebebfc
   md5: 1997a083ef0b4c9331f9191564be275e
@@ -14762,6 +14875,71 @@ packages:
   license_family: MIT
   size: 9568276
   timestamp: 1769521017574
+- conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.92.0-h53717f1_0.conda
+  sha256: c82a58098e06e887e41c4a08591218ec38e11c0bb0890c9ad0bd28ab9f261810
+  md5: a78c3f096ec96b2b505a148fa3984101
+  depends:
+  - __glibc >=2.17,<3.0.a0
+  - gcc_impl_linux-64
+  - libgcc >=14
+  - libzlib >=1.3.1,<2.0a0
+  - rust-std-x86_64-unknown-linux-gnu 1.92.0 h2c6d0dc_0
+  - sysroot_linux-64 >=2.17
+  license: MIT
+  license_family: MIT
+  size: 232828649
+  timestamp: 1765821075605
+- conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.92.0-h4ff7c5d_0.conda
+  sha256: 7cc5407dc6d559ef90118931faa4063c282dfed0472be562eacb12bf09b096c9
+  md5: 0ea02a89903b4f23918ac8aa20500919
+  depends:
+  - rust-std-aarch64-apple-darwin 1.92.0 hf6ec828_0
+  license: MIT
+  license_family: MIT
+  size: 241496727
+  timestamp: 1765820634853
+- conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.92.0-hf8d6059_0.conda
+  sha256: 91546fe10a0dcfd0027542f195dcd172be11cba8eaecdeed73ba433fcd834fdd
+  md5: 612aa55a4027e6c29a73938453392240
+  depends:
+  - rust-std-x86_64-pc-windows-msvc 1.92.0 h17fc481_0
+  license: MIT
+  license_family: MIT
+  size: 259934274
+  timestamp: 1765822676115
+- conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.92.0-hf6ec828_0.conda
+  sha256: b86a91f07127469b5a8c490a8b791551f13bb67a5081958de033daa3d6ceb3d4
+  md5: fde071794782ac4359173f1ddd4ae8d2
+  depends:
+  - __unix
+  constrains:
+  - rust >=1.92.0,<1.92.1.0a0
+  license: MIT
+  license_family: MIT
+  size: 34887424
+  timestamp: 1765820242072
+- conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.92.0-h17fc481_0.conda
+  sha256: d6538f085668b6954eb0a2dc9e90865dc2614e823cffb41d55ea2d8d6dae355e
+  md5: 7875e6952005584eec69200211a90fe6
+  depends:
+  - __win
+  constrains:
+  - rust >=1.92.0,<1.92.1.0a0
+  license: MIT
+  license_family: MIT
+  size: 28703603
+  timestamp: 1765822434271
+- conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.92.0-h2c6d0dc_0.conda
+  sha256: 19570f26206e2635f78d987233ba8960c684576f8571298a6108eed4967e7c9a
+  md5: ee54789987e177271d9f95ef7fd7fa31
+  depends:
+  - __unix
+  constrains:
+  - rust >=1.92.0,<1.92.1.0a0
+  license: MIT
+  license_family: MIT
+  size: 38587633
+  timestamp: 1765820881154
 - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.2-he8a4886_1.conda
   sha256: dec76e9faa3173579d34d226dbc91892417a80784911daf8e3f0eb9bad19d7a6
   md5: bade189a194e66b93c03021bd36c337b
diff --git a/pixi.toml b/pixi.toml
index 0ebf0615..99ee8cdb 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -4,7 +4,16 @@ channels = ["conda-forge"]
 platforms = ["linux-64", "osx-arm64", "win-64"]
 
 [tasks]
-postinstall = "pip install --no-build-isolation --no-deps --disable-pip-version-check -e ."
+# Build the Rust extension using maturin (config in pyproject.toml [tool.maturin])
+build-rust = """
+    maturin build --release -m rust_ext/Cargo.toml \
+    && pip install --force-reinstall rust_ext/target/wheels/*.whl \
+    && cp $(python -c "import tabmat_rust_ext; print(tabmat_rust_ext.__file__.replace('__init__.py', 'tabmat_rust_ext.abi3.so'))") src/tabmat/tabmat_rust_ext/
+"""
+# Build only C++ extensions (skip Rust)
+postinstall-cpp = "pip install --no-build-isolation --no-deps --disable-pip-version-check -e ."
+# Full build: Rust extension first, then C++ extensions
+postinstall = { cmd = "pip install --no-build-isolation --no-deps --disable-pip-version-check -e .", depends-on = ["build-rust"] }
 
 [tasks.install-nightlies]
 cmd = """
@@ -88,6 +97,8 @@ mako = "*"
 setuptools = ">=62.0"
 setuptools-scm = ">=8.1"
 xsimd = "<11|>12.1"
+rust = ">=1.70"
+maturin = ">=1.0"
 [target.unix.dependencies]
 jemalloc-local = "*"
 
diff --git a/pyproject.toml b/pyproject.toml
index 56c85e21..2c70d39a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -92,3 +92,17 @@ before-all = [
   "cmake -DCMAKE_INSTALL_PREFIX=/usr/local ..",
   "make install"
 ]
+
+# =============================================================================
+# Maturin configuration for the Rust extension
+# =============================================================================
+# The Rust extension is built separately using maturin. This configuration
+# tells maturin where to find the Cargo.toml and how to build the extension.
+# Build with: maturin build --release -m rust_ext/Cargo.toml
+# =============================================================================
+
+[tool.maturin]
+# Path to the Cargo.toml for the Rust extension
+manifest-path = "rust_ext/Cargo.toml"
+# Use stable ABI for Python 3.10+
+features = ["pyo3/extension-module", "pyo3/abi3-py310"]
diff --git a/rust_ext/Cargo.lock b/rust_ext/Cargo.lock
new file mode 100644
index 00000000..bd7d640f
--- /dev/null
+++ b/rust_ext/Cargo.lock
@@ -0,0 +1,269 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+version = 4
+
+[[package]]
+name = "autocfg"
+version = "1.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8"
+
+[[package]]
+name = "cfg-if"
+version = "1.0.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801"
+
+[[package]]
+name = "heck"
+version = "0.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea"
+
+[[package]]
+name = "indoc"
+version = "2.0.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "79cf5c93f93228cf8efb3ba362535fb11199ac548a09ce117c9b1adc3030d706"
+dependencies = [
+ "rustversion",
+]
+
+[[package]]
+name = "libc"
+version = "0.2.180"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bcc35a38544a891a5f7c865aca548a982ccb3b8650a5b06d0fd33a10283c56fc"
+
+[[package]]
+name = "matrixmultiply"
+version = "0.3.10"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a06de3016e9fae57a36fd14dba131fccf49f74b40b7fbdb472f96e361ec71a08"
+dependencies = [
+ "autocfg",
+ "rawpointer",
+]
+
+[[package]]
+name = "memoffset"
+version = "0.9.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "488016bfae457b036d996092f6cb448677611ce4449e970ceaf42695203f218a"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
+name = "ndarray"
+version = "0.16.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "882ed72dce9365842bf196bdeedf5055305f11fc8c03dee7bb0194a6cad34841"
+dependencies = [
+ "matrixmultiply",
+ "num-complex",
+ "num-integer",
+ "num-traits",
+ "portable-atomic",
+ "portable-atomic-util",
+ "rawpointer",
+]
+
+[[package]]
+name = "num-complex"
+version = "0.4.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-integer"
+version = "0.1.46"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-traits"
+version = "0.2.19"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
+name = "numpy"
+version = "0.23.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b94caae805f998a07d33af06e6a3891e38556051b8045c615470a71590e13e78"
+dependencies = [
+ "libc",
+ "ndarray",
+ "num-complex",
+ "num-integer",
+ "num-traits",
+ "pyo3",
+ "rustc-hash",
+]
+
+[[package]]
+name = "once_cell"
+version = "1.21.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "42f5e15c9953c5e4ccceeb2e7382a716482c34515315f7b03532b8b4e8393d2d"
+
+[[package]]
+name = "portable-atomic"
+version = "1.13.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c33a9471896f1c69cecef8d20cbe2f7accd12527ce60845ff44c153bb2a21b49"
+
+[[package]]
+name = "portable-atomic-util"
+version = "0.2.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7a9db96d7fa8782dd8c15ce32ffe8680bbd1e978a43bf51a34d39483540495f5"
+dependencies = [
+ "portable-atomic",
+]
+
+[[package]]
+name = "proc-macro2"
+version = "1.0.106"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8fd00f0bb2e90d81d1044c2b32617f68fcb9fa3bb7640c23e9c748e53fb30934"
+dependencies = [
+ "unicode-ident",
+]
+
+[[package]]
+name = "pyo3"
+version = "0.23.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7778bffd85cf38175ac1f545509665d0b9b92a198ca7941f131f85f7a4f9a872"
+dependencies = [
+ "cfg-if",
+ "indoc",
+ "libc",
+ "memoffset",
+ "once_cell",
+ "portable-atomic",
+ "pyo3-build-config",
+ "pyo3-ffi",
+ "pyo3-macros",
+ "unindent",
+]
+
+[[package]]
+name = "pyo3-build-config"
+version = "0.23.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "94f6cbe86ef3bf18998d9df6e0f3fc1050a8c5efa409bf712e661a4366e010fb"
+dependencies = [
+ "once_cell",
+ "target-lexicon",
+]
+
+[[package]]
+name = "pyo3-ffi"
+version = "0.23.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e9f1b4c431c0bb1c8fb0a338709859eed0d030ff6daa34368d3b152a63dfdd8d"
+dependencies = [
+ "libc",
+ "pyo3-build-config",
+]
+
+[[package]]
+name = "pyo3-macros"
+version = "0.23.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fbc2201328f63c4710f68abdf653c89d8dbc2858b88c5d88b0ff38a75288a9da"
+dependencies = [
+ "proc-macro2",
+ "pyo3-macros-backend",
+ "quote",
+ "syn",
+]
+
+[[package]]
+name = "pyo3-macros-backend"
+version = "0.23.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fca6726ad0f3da9c9de093d6f116a93c1a38e417ed73bf138472cf4064f72028"
+dependencies = [
+ "heck",
+ "proc-macro2",
+ "pyo3-build-config",
+ "quote",
+ "syn",
+]
+
+[[package]]
+name = "quote"
+version = "1.0.44"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "21b2ebcf727b7760c461f091f9f0f539b77b8e87f2fd88131e7f1b433b3cece4"
+dependencies = [
+ "proc-macro2",
+]
+
+[[package]]
+name = "rawpointer"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
+
+[[package]]
+name = "rustc-hash"
+version = "2.1.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "357703d41365b4b27c590e3ed91eabb1b663f07c4c084095e60cbed4362dff0d"
+
+[[package]]
+name = "rustversion"
+version = "1.0.22"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b39cdef0fa800fc44525c84ccb54a029961a8215f9619753635a9c0d2538d46d"
+
+[[package]]
+name = "syn"
+version = "2.0.114"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d4d107df263a3013ef9b1879b0df87d706ff80f65a86ea879bd9c31f9b307c2a"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "unicode-ident",
+]
+
+[[package]]
+name = "tabmat_rust_ext"
+version = "0.1.0"
+dependencies = [
+ "numpy",
+ "pyo3",
+]
+
+[[package]]
+name = "target-lexicon"
+version = "0.12.16"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "61c41af27dd6d1e27b1b16b489db798443478cef1f06a660c96db617ba5de3b1"
+
+[[package]]
+name = "unicode-ident"
+version = "1.0.22"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9312f7c4f6ff9069b165498234ce8be658059c6728633667c526e27dc2cf1df5"
+
+[[package]]
+name = "unindent"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7264e107f553ccae879d21fbea1d6724ac785e8c3bfc762137959b5802826ef3"
diff --git a/rust_ext/Cargo.toml b/rust_ext/Cargo.toml
new file mode 100644
index 00000000..f63a7894
--- /dev/null
+++ b/rust_ext/Cargo.toml
@@ -0,0 +1,16 @@
+[package]
+name = "tabmat_rust_ext"
+version = "0.1.0"
+edition = "2021"
+
+[lib]
+name = "tabmat_rust_ext"
+crate-type = ["cdylib"]
+
+[dependencies]
+pyo3 = { version = "0.23", features = ["extension-module", "abi3-py310"] }
+numpy = "0.23"
+
+[profile.release]
+lto = true
+opt-level = 3
diff --git a/rust_ext/src/categorical.rs b/rust_ext/src/categorical.rs
new file mode 100644
index 00000000..eb910449
--- /dev/null
+++ b/rust_ext/src/categorical.rs
@@ -0,0 +1,164 @@
+//! Core categorical matrix operations implemented in Rust.
+//!
+//! These functions implement the categorical matrix operations using simple
+//! sequential iteration. The design prioritizes clarity and correctness.
+
+// =============================================================================
+// transpose_matvec
+// =============================================================================
+
+/// Transpose matrix-vector multiplication: X.T @ other
+///
+/// Computes X.T @ other where X is a categorical matrix represented by indices.
+/// For each row i, X[i, indices[i]] = 1 and all other elements are 0.
+/// So (X.T @ other)[j] = sum of other[i] for all i where indices[i] == j.
+///
+/// When drop_first is true, category 0 is skipped and indices are shifted.
+pub fn transpose_matvec(
+    indices: &[i32],
+    other: &[f64],
+    out_size: usize,
+    drop_first: bool,
+) -> Vec<f64> {
+    let mut result = vec![0.0; out_size];
+    let offset = if drop_first { 1 } else { 0 };
+
+    for (i, &idx) in indices.iter().enumerate() {
+        let col_idx = idx - offset;
+        if col_idx >= 0 {
+            result[col_idx as usize] += other[i];
+        }
+    }
+
+    result
+}
+
+// =============================================================================
+// matvec
+// =============================================================================
+
+/// Matrix-vector multiplication: out[i] += other[indices[i] - offset]
+///
+/// For categorical matrices, this is a simple gather operation.
+/// When drop_first is true, category 0 contributes nothing.
+pub fn matvec(indices: &[i32], other: &[f64], out: &mut [f64], drop_first: bool) {
+    let offset = if drop_first { 1 } else { 0 };
+
+    for (i, &idx) in indices.iter().enumerate() {
+        let col_idx = idx - offset;
+        if col_idx >= 0 {
+            out[i] += other[col_idx as usize];
+        }
+    }
+}
+
+// =============================================================================
+// sandwich_categorical
+// =============================================================================
+
+/// Sandwich product returning diagonal: X.T @ diag(d) @ X
+///
+/// For categorical matrices, the result is always diagonal.
+/// Result[j] = sum of d[k] for all rows k in `rows` where indices[k] == j.
+pub fn sandwich_diagonal(
+    indices: &[i32],
+    d: &[f64],
+    rows: &[i32],
+    n_cols: usize,
+    drop_first: bool,
+) -> Vec<f64> {
+    let mut result = vec![0.0; n_cols];
+    let offset = if drop_first { 1 } else { 0 };
+
+    for &row in rows {
+        let k = row as usize;
+        let col_idx = indices[k] - offset;
+        if col_idx >= 0 {
+            result[col_idx as usize] += d[k];
+        }
+    }
+
+    result
+}
+
+// =============================================================================
+// sandwich_cat_cat
+// =============================================================================
+
+/// Cat-cat sandwich: X1.T @ diag(d) @ X2
+///
+/// Result[i, j] = sum of d[k] for rows k where i_indices[k] == i and j_indices[k] == j.
+/// Returns a flat vector in row-major order.
+pub fn sandwich_cat_cat(
+    i_indices: &[i32],
+    j_indices: &[i32],
+    d: &[f64],
+    rows: &[i32],
+    i_ncol: usize,
+    j_ncol: usize,
+    i_drop_first: bool,
+    j_drop_first: bool,
+) -> Vec<f64> {
+    let mut result = vec![0.0; i_ncol * j_ncol];
+    let i_offset = if i_drop_first { 1 } else { 0 };
+    let j_offset = if j_drop_first { 1 } else { 0 };
+
+    for &row in rows {
+        let k = row as usize;
+        let i = i_indices[k] - i_offset;
+        let j = j_indices[k] - j_offset;
+
+        if i >= 0 && j >= 0 {
+            result[i as usize * j_ncol + j as usize] += d[k];
+        }
+    }
+
+    result
+}
+
+// =============================================================================
+// sandwich_cat_dense
+// =============================================================================
+
+/// Cat-dense sandwich: X_cat.T @ diag(d) @ X_dense
+///
+/// Result[i, j_idx] = sum over k of (d[k] * X_dense[k, j_cols[j_idx]]) where i_indices[k] == i.
+/// Supports both C-contiguous and Fortran-contiguous layouts.
+pub fn sandwich_cat_dense(
+    i_indices: &[i32],
+    d: &[f64],
+    mat_j: &[f64],
+    mat_j_shape: (usize, usize),
+    rows: &[i32],
+    j_cols: &[i32],
+    i_ncol: usize,
+    is_c_contiguous: bool,
+    drop_first: bool,
+) -> Vec<f64> {
+    let (mat_j_nrow, mat_j_ncol) = mat_j_shape;
+    let len_j_cols = j_cols.len();
+    let mut result = vec![0.0; i_ncol * len_j_cols];
+    let offset = if drop_first { 1 } else { 0 };
+
+    for &row in rows {
+        let k = row as usize;
+        let i = i_indices[k] - offset;
+
+        if i >= 0 {
+            let i = i as usize;
+            let d_k = d[k];
+
+            for (j_idx, &j_col) in j_cols.iter().enumerate() {
+                let j = j_col as usize;
+                let mat_val = if is_c_contiguous {
+                    mat_j[k * mat_j_ncol + j]
+                } else {
+                    mat_j[j * mat_j_nrow + k]
+                };
+                result[i * len_j_cols + j_idx] += d_k * mat_val;
+            }
+        }
+    }
+
+    result
+}
diff --git a/rust_ext/src/dense.rs b/rust_ext/src/dense.rs
new file mode 100644
index 00000000..8323d66d
--- /dev/null
+++ b/rust_ext/src/dense.rs
@@ -0,0 +1,160 @@
+//! Core dense matrix operations implemented in Rust.
+//!
+//! These functions implement dense matrix operations using simple sequential
+//! iteration.
+
+// =============================================================================
+// Helper: access dense matrix elements
+// =============================================================================
+
+/// Get element from dense matrix at (row, col).
+/// Handles both C-contiguous (row-major) and Fortran-contiguous (column-major) layouts.
+fn get_element(data: &[f64], row: usize, col: usize, n_rows: usize, n_cols: usize, is_c_contiguous: bool) -> f64 {
+    if is_c_contiguous {
+        data[row * n_cols + col]
+    } else {
+        data[col * n_rows + row]
+    }
+}
+
+// =============================================================================
+// dense_sandwich
+// =============================================================================
+
+/// Dense sandwich product: X.T @ diag(d) @ X
+///
+/// Computes the symmetric matrix X[rows, cols].T @ diag(d[rows]) @ X[rows, cols].
+/// Returns a flattened array of shape (len(cols), len(cols)) in row-major order.
+pub fn dense_sandwich(
+    x: &[f64],
+    x_shape: (usize, usize),
+    d: &[f64],
+    rows: &[i32],
+    cols: &[i32],
+    is_c_contiguous: bool,
+) -> Vec<f64> {
+    let (n_rows, n_cols) = x_shape;
+    let out_m = cols.len();
+    let mut out = vec![0.0; out_m * out_m];
+
+    if rows.is_empty() || out_m == 0 {
+        return out;
+    }
+
+    // Compute upper triangle (including diagonal)
+    for (ci, &col_i) in cols.iter().enumerate() {
+        let i = col_i as usize;
+
+        for (cj, &col_j) in cols.iter().enumerate().skip(ci) {
+            let j = col_j as usize;
+
+            // Sum over restricted rows: d[k] * X[k, i] * X[k, j]
+            out[ci * out_m + cj] = rows
+                .iter()
+                .map(|&row| {
+                    let k = row as usize;
+                    let x_ki = get_element(x, k, i, n_rows, n_cols, is_c_contiguous);
+                    let x_kj = get_element(x, k, j, n_rows, n_cols, is_c_contiguous);
+                    d[k] * x_ki * x_kj
+                })
+                .sum();
+        }
+    }
+
+    // Mirror upper triangle to lower triangle
+    for ci in 0..out_m {
+        for cj in (ci + 1)..out_m {
+            out[cj * out_m + ci] = out[ci * out_m + cj];
+        }
+    }
+
+    out
+}
+
+// =============================================================================
+// dense_rmatvec
+// =============================================================================
+
+/// Dense transpose matrix-vector multiplication: X[rows, cols].T @ v[rows]
+///
+/// Returns a vector of length cols.len().
+pub fn dense_rmatvec(
+    x: &[f64],
+    x_shape: (usize, usize),
+    v: &[f64],
+    rows: &[i32],
+    cols: &[i32],
+    is_c_contiguous: bool,
+) -> Vec<f64> {
+    let (n_rows, n_cols) = x_shape;
+
+    cols.iter()
+        .map(|&col| {
+            let j = col as usize;
+            rows.iter()
+                .map(|&row| {
+                    let i = row as usize;
+                    get_element(x, i, j, n_rows, n_cols, is_c_contiguous) * v[i]
+                })
+                .sum()
+        })
+        .collect()
+}
+
+// =============================================================================
+// dense_matvec
+// =============================================================================
+
+/// Dense matrix-vector multiplication: X[rows, cols] @ v[cols]
+///
+/// Returns a vector of length rows.len().
+pub fn dense_matvec(
+    x: &[f64],
+    x_shape: (usize, usize),
+    v: &[f64],
+    rows: &[i32],
+    cols: &[i32],
+    is_c_contiguous: bool,
+) -> Vec<f64> {
+    let (n_rows, n_cols) = x_shape;
+
+    rows.iter()
+        .map(|&row| {
+            let i = row as usize;
+            cols.iter()
+                .map(|&col| {
+                    let j = col as usize;
+                    get_element(x, i, j, n_rows, n_cols, is_c_contiguous) * v[j]
+                })
+                .sum()
+        })
+        .collect()
+}
+
+// =============================================================================
+// transpose_square_dot_weights
+// =============================================================================
+
+/// Compute weighted squared column norms with shift for dense matrix.
+///
+/// For each column j: out[j] = sum_i(weights[i] * (X[i, j] - shift[j])^2)
+pub fn dense_transpose_square_dot_weights(
+    x: &[f64],
+    x_shape: (usize, usize),
+    weights: &[f64],
+    shift: &[f64],
+    is_c_contiguous: bool,
+) -> Vec<f64> {
+    let (n_rows, n_cols) = x_shape;
+
+    (0..n_cols)
+        .map(|j| {
+            (0..n_rows)
+                .map(|i| {
+                    let x_ij = get_element(x, i, j, n_rows, n_cols, is_c_contiguous);
+                    weights[i] * (x_ij - shift[j]).powi(2)
+                })
+                .sum()
+        })
+        .collect()
+}
diff --git a/rust_ext/src/lib.rs b/rust_ext/src/lib.rs
new file mode 100644
index 00000000..c8143bd0
--- /dev/null
+++ b/rust_ext/src/lib.rs
@@ -0,0 +1,569 @@
+//! Rust backend for tabmat matrix operations.
+//!
+//! This module provides idiomatic Rust implementations of categorical and sparse
+//! matrix operations that are also available in the C++ backend. The implementations
+//! use simple sequential iteration for clarity and correctness.
+
+use numpy::{
+    PyArray1, PyArray2, PyArrayMethods, PyReadonlyArray1, PyReadonlyArray2, PyUntypedArrayMethods,
+};
+use pyo3::prelude::*;
+
+mod categorical;
+mod dense;
+mod sparse;
+
+/// Rust backend module for tabmat operations.
+#[pymodule]
+fn tabmat_rust_ext(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    // Categorical operations
+    m.add_function(wrap_pyfunction!(transpose_matvec, m)?)?;
+    m.add_function(wrap_pyfunction!(matvec, m)?)?;
+    m.add_function(wrap_pyfunction!(sandwich_categorical, m)?)?;
+    m.add_function(wrap_pyfunction!(sandwich_cat_cat, m)?)?;
+    m.add_function(wrap_pyfunction!(sandwich_cat_dense, m)?)?;
+
+    // Sparse operations
+    m.add_function(wrap_pyfunction!(csr_matvec_unrestricted, m)?)?;
+    m.add_function(wrap_pyfunction!(csr_matvec, m)?)?;
+    m.add_function(wrap_pyfunction!(csc_rmatvec_unrestricted, m)?)?;
+    m.add_function(wrap_pyfunction!(csc_rmatvec, m)?)?;
+    m.add_function(wrap_pyfunction!(sparse_sandwich, m)?)?;
+    m.add_function(wrap_pyfunction!(csr_dense_sandwich, m)?)?;
+    m.add_function(wrap_pyfunction!(transpose_square_dot_weights, m)?)?;
+
+    // Dense operations
+    m.add_function(wrap_pyfunction!(dense_sandwich, m)?)?;
+    m.add_function(wrap_pyfunction!(dense_rmatvec, m)?)?;
+    m.add_function(wrap_pyfunction!(dense_matvec, m)?)?;
+    m.add_function(wrap_pyfunction!(dense_transpose_square_dot_weights, m)?)?;
+    Ok(())
+}
+
+// =============================================================================
+// transpose_matvec
+// =============================================================================
+
+/// Transpose matrix-vector multiplication: out += X.T @ other
+///
+/// For categorical matrices, X[i, indices[i]] = 1 and all other elements are 0.
+/// So (X.T @ other)[j] = sum over all rows i where indices[i] == j of other[i].
+#[pyfunction]
+fn transpose_matvec<'py>(
+    _py: Python<'py>,
+    indices: PyReadonlyArray1<'py, i32>,
+    other: PyReadonlyArray1<'py, f64>,
+    out: &Bound<'py, PyArray1<f64>>,
+    drop_first: bool,
+) -> PyResult<()> {
+    let indices = indices.as_slice()?;
+    let other = other.as_slice()?;
+    let out_len = out.len();
+
+    let result = categorical::transpose_matvec(indices, other, out_len, drop_first);
+
+    // Add results to output array
+    {
+        let out_slice = unsafe { out.as_slice_mut()? };
+        for (o, r) in out_slice.iter_mut().zip(result.iter()) {
+            *o += *r;
+        }
+    }
+
+    Ok(())
+}
+
+// =============================================================================
+// matvec
+// =============================================================================
+
+/// Matrix-vector multiplication: out[i] += other[indices[i]]
+///
+/// For categorical matrices with one-hot encoding, this is just a lookup.
+#[pyfunction]
+fn matvec<'py>(
+    _py: Python<'py>,
+    indices: PyReadonlyArray1<'py, i32>,
+    other: PyReadonlyArray1<'py, f64>,
+    out: &Bound<'py, PyArray1<f64>>,
+    drop_first: bool,
+) -> PyResult<()> {
+    let indices = indices.as_slice()?;
+    let other = other.as_slice()?;
+
+    {
+        let out_slice = unsafe { out.as_slice_mut()? };
+        categorical::matvec(indices, other, out_slice, drop_first);
+    }
+
+    Ok(())
+}
+
+// =============================================================================
+// sandwich_categorical
+// =============================================================================
+
+/// Sandwich product for categorical matrix: X.T @ diag(d) @ X
+///
+/// For categorical matrices, the result is always diagonal.
+/// Result[j] = sum of d[i] for all rows i where indices[i] == j.
+#[pyfunction]
+fn sandwich_categorical<'py>(
+    py: Python<'py>,
+    indices: PyReadonlyArray1<'py, i32>,
+    d: PyReadonlyArray1<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    n_cols: usize,
+    drop_first: bool,
+) -> PyResult<Py<PyArray1<f64>>> {
+    let indices = indices.as_slice()?;
+    let d = d.as_slice()?;
+    let rows = rows.as_slice()?;
+
+    let result = categorical::sandwich_diagonal(indices, d, rows, n_cols, drop_first);
+
+    Ok(PyArray1::from_vec(py, result).into())
+}
+
+// =============================================================================
+// sandwich_cat_cat
+// =============================================================================
+
+/// Cat-cat sandwich: X1.T @ diag(d) @ X2
+///
+/// Result[i, j] = sum of d[k] for all rows k where i_indices[k] == i and j_indices[k] == j.
+#[pyfunction]
+fn sandwich_cat_cat<'py>(
+    py: Python<'py>,
+    i_indices: PyReadonlyArray1<'py, i32>,
+    j_indices: PyReadonlyArray1<'py, i32>,
+    d: PyReadonlyArray1<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    i_ncol: usize,
+    j_ncol: usize,
+    i_drop_first: bool,
+    j_drop_first: bool,
+) -> PyResult<Py<PyArray2<f64>>> {
+    let i_indices = i_indices.as_slice()?;
+    let j_indices = j_indices.as_slice()?;
+    let d = d.as_slice()?;
+    let rows = rows.as_slice()?;
+
+    let result = categorical::sandwich_cat_cat(
+        i_indices,
+        j_indices,
+        d,
+        rows,
+        i_ncol,
+        j_ncol,
+        i_drop_first,
+        j_drop_first,
+    );
+
+    // Convert flat Vec to 2D array
+    let array = PyArray2::from_vec2(
+        py,
+        &result
+            .chunks(j_ncol)
+            .map(|chunk| chunk.to_vec())
+            .collect::<Vec<_>>(),
+    )?;
+
+    Ok(array.into())
+}
+
+// =============================================================================
+// sandwich_cat_dense
+// =============================================================================
+
+/// Cat-dense sandwich: X_cat.T @ diag(d) @ X_dense
+///
+/// Result[i, j] = sum over k of (d[k] * X_dense[k, j]) where i_indices[k] == i.
+#[pyfunction]
+fn sandwich_cat_dense<'py>(
+    py: Python<'py>,
+    i_indices: PyReadonlyArray1<'py, i32>,
+    d: PyReadonlyArray1<'py, f64>,
+    mat_j: PyReadonlyArray2<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    j_cols: PyReadonlyArray1<'py, i32>,
+    i_ncol: usize,
+    is_c_contiguous: bool,
+    drop_first: bool,
+) -> PyResult<Py<PyArray2<f64>>> {
+    let i_indices = i_indices.as_slice()?;
+    let d = d.as_slice()?;
+    let rows = rows.as_slice()?;
+    let j_cols_slice = j_cols.as_slice()?;
+    let shape = mat_j.shape();
+    let mat_j_shape = (shape[0], shape[1]);
+
+    let mat_j_slice = mat_j.as_slice()?;
+
+    let n_j_cols = j_cols_slice.len();
+
+    // Handle edge case: empty j_cols produces empty result with correct shape
+    if n_j_cols == 0 || i_ncol == 0 {
+        let array = PyArray2::zeros(py, [i_ncol, n_j_cols], false);
+        return Ok(array.unbind());
+    }
+
+    let result = categorical::sandwich_cat_dense(
+        i_indices,
+        d,
+        mat_j_slice,
+        mat_j_shape,
+        rows,
+        j_cols_slice,
+        i_ncol,
+        is_c_contiguous,
+        drop_first,
+    );
+
+    let array = PyArray2::from_vec2(
+        py,
+        &result
+            .chunks(n_j_cols)
+            .map(|chunk| chunk.to_vec())
+            .collect::<Vec<_>>(),
+    )?;
+
+    Ok(array.into())
+}
+
+// =============================================================================
+// Sparse operations
+// =============================================================================
+
+/// CSR matrix-vector multiplication: out += X @ v (unrestricted)
+#[pyfunction]
+fn csr_matvec_unrestricted<'py>(
+    _py: Python<'py>,
+    data: PyReadonlyArray1<'py, f64>,
+    indices: PyReadonlyArray1<'py, i32>,
+    indptr: PyReadonlyArray1<'py, i32>,
+    v: PyReadonlyArray1<'py, f64>,
+    out: &Bound<'py, PyArray1<f64>>,
+) -> PyResult<()> {
+    let data = data.as_slice()?;
+    let indices = indices.as_slice()?;
+    let indptr = indptr.as_slice()?;
+    let v = v.as_slice()?;
+
+    {
+        let out_slice = unsafe { out.as_slice_mut()? };
+        sparse::csr_matvec_unrestricted(data, indices, indptr, v, out_slice);
+    }
+
+    Ok(())
+}
+
+/// CSR matrix-vector multiplication with row/column restrictions
+#[pyfunction]
+fn csr_matvec<'py>(
+    py: Python<'py>,
+    data: PyReadonlyArray1<'py, f64>,
+    indices: PyReadonlyArray1<'py, i32>,
+    indptr: PyReadonlyArray1<'py, i32>,
+    v: PyReadonlyArray1<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    cols: PyReadonlyArray1<'py, i32>,
+    n_cols_total: usize,
+) -> PyResult<Py<PyArray1<f64>>> {
+    let data = data.as_slice()?;
+    let indices = indices.as_slice()?;
+    let indptr = indptr.as_slice()?;
+    let v = v.as_slice()?;
+    let rows = rows.as_slice()?;
+    let cols = cols.as_slice()?;
+
+    let result = sparse::csr_matvec(data, indices, indptr, v, rows, cols, n_cols_total);
+
+    Ok(PyArray1::from_vec(py, result).into())
+}
+
+/// CSC transpose matrix-vector multiplication: out += XT.T @ v (unrestricted)
+#[pyfunction]
+fn csc_rmatvec_unrestricted<'py>(
+    _py: Python<'py>,
+    data: PyReadonlyArray1<'py, f64>,
+    indices: PyReadonlyArray1<'py, i32>,
+    indptr: PyReadonlyArray1<'py, i32>,
+    v: PyReadonlyArray1<'py, f64>,
+    out: &Bound<'py, PyArray1<f64>>,
+) -> PyResult<()> {
+    let data = data.as_slice()?;
+    let indices = indices.as_slice()?;
+    let indptr = indptr.as_slice()?;
+    let v = v.as_slice()?;
+
+    {
+        let out_slice = unsafe { out.as_slice_mut()? };
+        sparse::csc_rmatvec_unrestricted(data, indices, indptr, v, out_slice);
+    }
+
+    Ok(())
+}
+
+/// CSC transpose matrix-vector multiplication with restrictions
+#[pyfunction]
+fn csc_rmatvec<'py>(
+    py: Python<'py>,
+    data: PyReadonlyArray1<'py, f64>,
+    indices: PyReadonlyArray1<'py, i32>,
+    indptr: PyReadonlyArray1<'py, i32>,
+    v: PyReadonlyArray1<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    cols: PyReadonlyArray1<'py, i32>,
+    n_rows_total: usize,
+) -> PyResult<Py<PyArray1<f64>>> {
+    let data = data.as_slice()?;
+    let indices = indices.as_slice()?;
+    let indptr = indptr.as_slice()?;
+    let v = v.as_slice()?;
+    let rows = rows.as_slice()?;
+    let cols = cols.as_slice()?;
+
+    let result = sparse::csc_rmatvec(data, indices, indptr, v, rows, cols, n_rows_total);
+
+    Ok(PyArray1::from_vec(py, result).into())
+}
+
+/// Sparse sandwich product: AT @ diag(d) @ A
+#[pyfunction]
+fn sparse_sandwich<'py>(
+    py: Python<'py>,
+    a_data: PyReadonlyArray1<'py, f64>,
+    a_indices: PyReadonlyArray1<'py, i32>,
+    a_indptr: PyReadonlyArray1<'py, i32>,
+    at_data: PyReadonlyArray1<'py, f64>,
+    at_indices: PyReadonlyArray1<'py, i32>,
+    at_indptr: PyReadonlyArray1<'py, i32>,
+    d: PyReadonlyArray1<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    cols: PyReadonlyArray1<'py, i32>,
+    n_rows_total: usize,
+    n_cols_total: usize,
+) -> PyResult<Py<PyArray2<f64>>> {
+    let a_data = a_data.as_slice()?;
+    let a_indices = a_indices.as_slice()?;
+    let a_indptr = a_indptr.as_slice()?;
+    let at_data = at_data.as_slice()?;
+    let at_indices = at_indices.as_slice()?;
+    let at_indptr = at_indptr.as_slice()?;
+    let d = d.as_slice()?;
+    let rows = rows.as_slice()?;
+    let cols = cols.as_slice()?;
+
+    let m = cols.len();
+
+    // Handle edge case: empty cols produces empty result with correct shape
+    if m == 0 {
+        let array = PyArray2::zeros(py, [0, 0], false);
+        return Ok(array.unbind());
+    }
+
+    let result = sparse::sparse_sandwich(
+        a_data,
+        a_indices,
+        a_indptr,
+        at_data,
+        at_indices,
+        at_indptr,
+        d,
+        rows,
+        cols,
+        n_rows_total,
+        n_cols_total,
+    );
+
+    let array = PyArray2::from_vec2(
+        py,
+        &result
+            .chunks(m)
+            .map(|chunk| chunk.to_vec())
+            .collect::<Vec<_>>(),
+    )?;
+
+    Ok(array.into())
+}
+
+/// CSR-dense sandwich: A.T @ diag(d) @ B
+#[pyfunction]
+fn csr_dense_sandwich<'py>(
+    py: Python<'py>,
+    a_data: PyReadonlyArray1<'py, f64>,
+    a_indices: PyReadonlyArray1<'py, i32>,
+    a_indptr: PyReadonlyArray1<'py, i32>,
+    b: PyReadonlyArray2<'py, f64>,
+    d: PyReadonlyArray1<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    a_cols: PyReadonlyArray1<'py, i32>,
+    b_cols: PyReadonlyArray1<'py, i32>,
+    a_ncol: usize,
+    is_c_contiguous: bool,
+) -> PyResult<Py<PyArray2<f64>>> {
+    let a_data = a_data.as_slice()?;
+    let a_indices = a_indices.as_slice()?;
+    let a_indptr = a_indptr.as_slice()?;
+    let d = d.as_slice()?;
+    let rows = rows.as_slice()?;
+    let a_cols = a_cols.as_slice()?;
+    let b_cols = b_cols.as_slice()?;
+    let b_shape = (b.shape()[0], b.shape()[1]);
+    let b_slice = b.as_slice()?;
+
+    let n_a_cols = a_cols.len();
+    let n_b_cols = b_cols.len();
+
+    // Handle edge cases: empty dimensions produce empty result with correct shape
+    if n_a_cols == 0 || n_b_cols == 0 {
+        // Use zeros to create array with correct shape
+        let array = PyArray2::zeros(py, [n_a_cols, n_b_cols], false);
+        return Ok(array.unbind());
+    }
+
+    let result = sparse::csr_dense_sandwich(
+        a_data,
+        a_indices,
+        a_indptr,
+        b_slice,
+        b_shape,
+        d,
+        rows,
+        a_cols,
+        b_cols,
+        a_ncol,
+        is_c_contiguous,
+    );
+
+    let array = PyArray2::from_vec2(
+        py,
+        &result
+            .chunks(n_b_cols)
+            .map(|chunk| chunk.to_vec())
+            .collect::<Vec<_>>(),
+    )?;
+
+    Ok(array.into())
+}
+
+/// Compute weighted squared column norms for CSC matrix
+#[pyfunction]
+fn transpose_square_dot_weights<'py>(
+    py: Python<'py>,
+    data: PyReadonlyArray1<'py, f64>,
+    indices: PyReadonlyArray1<'py, i32>,
+    indptr: PyReadonlyArray1<'py, i32>,
+    weights: PyReadonlyArray1<'py, f64>,
+) -> PyResult<Py<PyArray1<f64>>> {
+    let data = data.as_slice()?;
+    let indices = indices.as_slice()?;
+    let indptr = indptr.as_slice()?;
+    let weights = weights.as_slice()?;
+
+    let result = sparse::transpose_square_dot_weights(data, indices, indptr, weights);
+
+    Ok(PyArray1::from_vec(py, result).into())
+}
+
+// =============================================================================
+// Dense operations
+// =============================================================================
+
+/// Dense sandwich product: X.T @ diag(d) @ X
+#[pyfunction]
+fn dense_sandwich<'py>(
+    py: Python<'py>,
+    x: PyReadonlyArray2<'py, f64>,
+    d: PyReadonlyArray1<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    cols: PyReadonlyArray1<'py, i32>,
+) -> PyResult<Py<PyArray2<f64>>> {
+    let x_slice = x.as_slice()?;
+    let x_shape = (x.shape()[0], x.shape()[1]);
+    let d = d.as_slice()?;
+    let rows = rows.as_slice()?;
+    let cols = cols.as_slice()?;
+    let is_c_contiguous = x.is_c_contiguous();
+
+    let out_m = cols.len();
+
+    if out_m == 0 {
+        let array = PyArray2::zeros(py, [0, 0], false);
+        return Ok(array.unbind());
+    }
+
+    let result = dense::dense_sandwich(x_slice, x_shape, d, rows, cols, is_c_contiguous);
+
+    let array = PyArray2::from_vec2(
+        py,
+        &result
+            .chunks(out_m)
+            .map(|chunk| chunk.to_vec())
+            .collect::<Vec<_>>(),
+    )?;
+
+    Ok(array.into())
+}
+
+/// Dense transpose matrix-vector multiplication: X.T @ v
+#[pyfunction]
+fn dense_rmatvec<'py>(
+    py: Python<'py>,
+    x: PyReadonlyArray2<'py, f64>,
+    v: PyReadonlyArray1<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    cols: PyReadonlyArray1<'py, i32>,
+) -> PyResult<Py<PyArray1<f64>>> {
+    let x_slice = x.as_slice()?;
+    let x_shape = (x.shape()[0], x.shape()[1]);
+    let v = v.as_slice()?;
+    let rows = rows.as_slice()?;
+    let cols = cols.as_slice()?;
+    let is_c_contiguous = x.is_c_contiguous();
+
+    let result = dense::dense_rmatvec(x_slice, x_shape, v, rows, cols, is_c_contiguous);
+
+    Ok(PyArray1::from_vec(py, result).into())
+}
+
+/// Dense matrix-vector multiplication: X @ v
+#[pyfunction]
+fn dense_matvec<'py>(
+    py: Python<'py>,
+    x: PyReadonlyArray2<'py, f64>,
+    v: PyReadonlyArray1<'py, f64>,
+    rows: PyReadonlyArray1<'py, i32>,
+    cols: PyReadonlyArray1<'py, i32>,
+) -> PyResult<Py<PyArray1<f64>>> {
+    let x_slice = x.as_slice()?;
+    let x_shape = (x.shape()[0], x.shape()[1]);
+    let v = v.as_slice()?;
+    let rows = rows.as_slice()?;
+    let cols = cols.as_slice()?;
+    let is_c_contiguous = x.is_c_contiguous();
+
+    let result = dense::dense_matvec(x_slice, x_shape, v, rows, cols, is_c_contiguous);
+
+    Ok(PyArray1::from_vec(py, result).into())
+}
+
+/// Compute weighted squared column norms with shift for dense matrix
+#[pyfunction]
+fn dense_transpose_square_dot_weights<'py>(
+    py: Python<'py>,
+    x: PyReadonlyArray2<'py, f64>,
+    weights: PyReadonlyArray1<'py, f64>,
+    shift: PyReadonlyArray1<'py, f64>,
+) -> PyResult<Py<PyArray1<f64>>> {
+    let x_slice = x.as_slice()?;
+    let x_shape = (x.shape()[0], x.shape()[1]);
+    let weights = weights.as_slice()?;
+    let shift = shift.as_slice()?;
+    let is_c_contiguous = x.is_c_contiguous();
+
+    let result = dense::dense_transpose_square_dot_weights(x_slice, x_shape, weights, shift, is_c_contiguous);
+
+    Ok(PyArray1::from_vec(py, result).into())
+}
diff --git a/rust_ext/src/sparse.rs b/rust_ext/src/sparse.rs
new file mode 100644
index 00000000..ce6ad750
--- /dev/null
+++ b/rust_ext/src/sparse.rs
@@ -0,0 +1,341 @@
+//! Core sparse matrix operations implemented in Rust.
+//!
+//! These functions implement sparse matrix operations (CSR/CSC format) using
+//! simple sequential iteration. The design prioritizes clarity and correctness.
+
+// =============================================================================
+// CSR matvec operations
+// =============================================================================
+
+/// CSR matrix-vector multiplication: out += X @ v (unrestricted)
+///
+/// For CSR format:
+/// - data[indptr[i]..indptr[i+1]] contains the non-zero values in row i
+/// - indices[indptr[i]..indptr[i+1]] contains the column indices for those values
+pub fn csr_matvec_unrestricted(
+    data: &[f64],
+    indices: &[i32],
+    indptr: &[i32],
+    v: &[f64],
+    out: &mut [f64],
+) {
+    let n_rows = indptr.len() - 1;
+
+    for i in 0..n_rows {
+        let start = indptr[i] as usize;
+        let end = indptr[i + 1] as usize;
+
+        for idx in start..end {
+            let j = indices[idx] as usize;
+            let x_val = data[idx];
+            out[i] += x_val * v[j];
+        }
+    }
+}
+
+/// CSR matrix-vector multiplication with row/column restrictions: out = X[rows, cols] @ v
+///
+/// Only considers specified rows and columns.
+/// Output has length rows.len().
+pub fn csr_matvec(
+    data: &[f64],
+    indices: &[i32],
+    indptr: &[i32],
+    v: &[f64],
+    rows: &[i32],
+    cols: &[i32],
+    n_cols_total: usize,
+) -> Vec<f64> {
+    let n = rows.len();
+    let mut out = vec![0.0; n];
+
+    // Build column inclusion mask
+    let mut col_included = vec![false; n_cols_total];
+    for &col in cols {
+        col_included[col as usize] = true;
+    }
+
+    for (ci, &row) in rows.iter().enumerate() {
+        let i = row as usize;
+        let start = indptr[i] as usize;
+        let end = indptr[i + 1] as usize;
+
+        for idx in start..end {
+            let j = indices[idx] as usize;
+            if !col_included[j] {
+                continue;
+            }
+            let x_val = data[idx];
+            out[ci] += x_val * v[j];
+        }
+    }
+
+    out
+}
+
+// =============================================================================
+// CSC rmatvec operations (transpose matvec)
+// =============================================================================
+
+/// CSC transpose matrix-vector multiplication: out += XT.T @ v (unrestricted)
+///
+/// XT is stored in CSC format (which is the same as X.T in CSR format).
+/// This computes X.T @ v where XT represents the transpose.
+///
+/// For CSC format:
+/// - data[indptr[j]..indptr[j+1]] contains the non-zero values in column j
+/// - indices[indptr[j]..indptr[j+1]] contains the row indices for those values
+pub fn csc_rmatvec_unrestricted(
+    data: &[f64],
+    indices: &[i32],
+    indptr: &[i32],
+    v: &[f64],
+    out: &mut [f64],
+) {
+    let n_cols = indptr.len() - 1;
+
+    for j in 0..n_cols {
+        let start = indptr[j] as usize;
+        let end = indptr[j + 1] as usize;
+
+        for idx in start..end {
+            let i = indices[idx] as usize;
+            let xt_val = data[idx];
+            out[j] += xt_val * v[i];
+        }
+    }
+}
+
+/// CSC transpose matrix-vector multiplication with restrictions: out = XT[rows, cols].T @ v
+///
+/// Only considers specified rows (of original matrix) and columns.
+/// Output has length cols.len().
+pub fn csc_rmatvec(
+    data: &[f64],
+    indices: &[i32],
+    indptr: &[i32],
+    v: &[f64],
+    rows: &[i32],
+    cols: &[i32],
+    n_rows_total: usize,
+) -> Vec<f64> {
+    let m = cols.len();
+    let mut out = vec![0.0; m];
+
+    // Build row inclusion mask
+    let mut row_included = vec![false; n_rows_total];
+    for &row in rows {
+        row_included[row as usize] = true;
+    }
+
+    for (cj, &col) in cols.iter().enumerate() {
+        let j = col as usize;
+        let start = indptr[j] as usize;
+        let end = indptr[j + 1] as usize;
+
+        for idx in start..end {
+            let i = indices[idx] as usize;
+            if !row_included[i] {
+                continue;
+            }
+            let xt_val = data[idx];
+            out[cj] += xt_val * v[i];
+        }
+    }
+
+    out
+}
+
+// =============================================================================
+// sparse_sandwich
+// =============================================================================
+
+/// Sparse sandwich product: AT @ diag(d) @ A
+///
+/// Both A and AT are in CSC format. AT is the transpose of A.
+/// This computes the symmetric matrix A.T @ diag(d) @ A.
+///
+/// Returns the upper triangle filled, then mirrors to lower triangle.
+pub fn sparse_sandwich(
+    a_data: &[f64],
+    a_indices: &[i32],
+    a_indptr: &[i32],
+    at_data: &[f64],
+    at_indices: &[i32],
+    at_indptr: &[i32],
+    d: &[f64],
+    rows: &[i32],
+    cols: &[i32],
+    n_rows_total: usize,
+    n_cols_total: usize,
+) -> Vec<f64> {
+    let m = cols.len();
+    let mut out = vec![0.0; m * m];
+
+    // Build row inclusion mask
+    let mut row_included = vec![false; n_rows_total];
+    for &row in rows {
+        row_included[row as usize] = true;
+    }
+
+    // Build column map: original col index -> output col index (or -1 if not included)
+    let mut col_map = vec![-1i32; n_cols_total];
+    for (cj, &col) in cols.iter().enumerate() {
+        col_map[col as usize] = cj as i32;
+    }
+
+    // For each output column Cj (corresponding to original column j)
+    for (cj, &col_j) in cols.iter().enumerate() {
+        let j = col_j as usize;
+
+        // Iterate over non-zeros in column j of A
+        let a_start = a_indptr[j] as usize;
+        let a_end = a_indptr[j + 1] as usize;
+
+        for a_idx in a_start..a_end {
+            let k = a_indices[a_idx] as usize;
+            if !row_included[k] {
+                continue;
+            }
+
+            let a_val = a_data[a_idx] * d[k];
+
+            // Now iterate over non-zeros in row k of AT (which is column k of AT in CSC)
+            let at_start = at_indptr[k] as usize;
+            let at_end = at_indptr[k + 1] as usize;
+
+            for at_idx in at_start..at_end {
+                let i = at_indices[at_idx] as usize;
+
+                // Only compute upper triangle (i <= j)
+                if i > j {
+                    break;
+                }
+
+                let ci = col_map[i];
+                if ci == -1 {
+                    continue;
+                }
+
+                let at_val = at_data[at_idx];
+                out[cj * m + ci as usize] += at_val * a_val;
+            }
+        }
+    }
+
+    // Mirror upper triangle to lower triangle
+    for i in 0..m {
+        for j in (i + 1)..m {
+            out[i * m + j] = out[j * m + i];
+        }
+    }
+
+    out
+}
+
+// =============================================================================
+// csr_dense_sandwich
+// =============================================================================
+
+/// CSR-dense sandwich: A.T @ diag(d) @ B
+///
+/// A is a sparse CSR matrix, B is a dense matrix.
+/// Returns a dense matrix of shape (len(A_cols), len(B_cols)).
+pub fn csr_dense_sandwich(
+    a_data: &[f64],
+    a_indices: &[i32],
+    a_indptr: &[i32],
+    b: &[f64],
+    b_shape: (usize, usize),
+    d: &[f64],
+    rows: &[i32],
+    a_cols: &[i32],
+    b_cols: &[i32],
+    a_ncol: usize,
+    is_c_contiguous: bool,
+) -> Vec<f64> {
+    let n_a_cols = a_cols.len();
+    let n_b_cols = b_cols.len();
+    let (_b_nrow, b_ncol) = b_shape;
+
+    let mut out = vec![0.0; n_a_cols * n_b_cols];
+
+    if rows.is_empty() || n_a_cols == 0 || n_b_cols == 0 {
+        return out;
+    }
+
+    // Build A column map
+    let mut a_col_map = vec![-1i32; a_ncol];
+    for (ci, &col) in a_cols.iter().enumerate() {
+        a_col_map[col as usize] = ci as i32;
+    }
+
+    // For each row k in the specified rows
+    for &row in rows {
+        let k = row as usize;
+        let d_k = d[k];
+
+        // Iterate over non-zeros in row k of A (CSR)
+        let start = a_indptr[k] as usize;
+        let end = a_indptr[k + 1] as usize;
+
+        for a_idx in start..end {
+            let i = a_indices[a_idx] as usize;
+            let ci = a_col_map[i];
+            if ci == -1 {
+                continue;
+            }
+            let ci = ci as usize;
+
+            let a_val = a_data[a_idx];
+            let q = a_val * d_k;
+
+            // Accumulate contribution for each B column
+            for (cj, &b_col) in b_cols.iter().enumerate() {
+                let j = b_col as usize;
+                let b_val = if is_c_contiguous {
+                    b[k * b_ncol + j]
+                } else {
+                    // Fortran contiguous
+                    b[j * b_shape.0 + k]
+                };
+                out[ci * n_b_cols + cj] += q * b_val;
+            }
+        }
+    }
+
+    out
+}
+
+// =============================================================================
+// transpose_square_dot_weights
+// =============================================================================
+
+/// Compute weighted squared column norms for CSC matrix.
+///
+/// For each column j: out[j] = sum over non-zeros in column j of (weights[i] * data[idx]^2)
+/// where i is the row index.
+///
+/// This is used for computing column-wise weighted squared norms.
+pub fn transpose_square_dot_weights(
+    data: &[f64],
+    indices: &[i32],
+    indptr: &[i32],
+    weights: &[f64],
+) -> Vec<f64> {
+    let n_cols = indptr.len() - 1;
+    let mut out = vec![0.0; n_cols];
+
+    for j in 0..n_cols {
+        let start = indptr[j] as usize;
+        let end = indptr[j + 1] as usize;
+
+        for idx in start..end {
+            let i = indices[idx] as usize;
+            let v = data[idx];
+            out[j] += weights[i] * v * v;
+        }
+    }
+
+    out
+}
diff --git a/setup.py b/setup.py
index 0c8af7f8..5a36f41f 100644
--- a/setup.py
+++ b/setup.py
@@ -156,6 +156,9 @@
     ],
     package_dir={"": "src"},
     packages=find_packages(where="src"),
+    package_data={
+        "tabmat.tabmat_rust_ext": ["*.so", "*.pyd", "*.dylib"],
+    },
     install_requires=["formulaic>=1.2", "narwhals", "numpy", "scipy"],
     python_requires=">=3.10",
     ext_modules=cythonize(
diff --git a/src/tabmat/categorical_matrix.py b/src/tabmat/categorical_matrix.py
index 3fbd5447..efcceef1 100644
--- a/src/tabmat/categorical_matrix.py
+++ b/src/tabmat/categorical_matrix.py
@@ -173,7 +173,7 @@ def matvec(mat, vec):
 from scipy import sparse as sps
 
 from .dense_matrix import DenseMatrix
-from .ext.categorical import (
+from .ext.categorical_dispatch import (
     matvec_complex,
     matvec_fast,
     multiply_complex,
@@ -183,7 +183,7 @@ def matvec(mat, vec):
     transpose_matvec_complex,
     transpose_matvec_fast,
 )
-from .ext.split import sandwich_cat_cat, sandwich_cat_dense
+from .ext.split_dispatch import sandwich_cat_cat, sandwich_cat_dense
 from .matrix_base import MatrixBase
 from .sparse_matrix import SparseMatrix
 from .util import (
diff --git a/src/tabmat/dense_matrix.py b/src/tabmat/dense_matrix.py
index 7daeeff6..4876359f 100644
--- a/src/tabmat/dense_matrix.py
+++ b/src/tabmat/dense_matrix.py
@@ -4,7 +4,7 @@
 
 import numpy as np
 
-from .ext.dense import (
+from .ext.dense_dispatch import (
     dense_matvec,
     dense_rmatvec,
     dense_sandwich,
diff --git a/src/tabmat/ext/backend.py b/src/tabmat/ext/backend.py
new file mode 100644
index 00000000..b506d136
--- /dev/null
+++ b/src/tabmat/ext/backend.py
@@ -0,0 +1,403 @@
+"""Backend selection for categorical and sparse matrix operations.
+
+This module provides a unified interface for categorical and sparse matrix
+operations that can use either the C++ (Cython) or Rust backend.
+
+Usage:
+    from tabmat.ext.backend import get_backend, set_backend
+
+    # Get current backend
+    backend = get_backend()
+
+    # Set backend to "rust" or "cpp"
+    set_backend("rust")
+
+    # Use functions from the backend
+    backend.transpose_matvec_fast(...)
+    backend.csr_matvec_unrestricted(...)
+"""
+
+import os
+from typing import Literal
+
+# Global backend setting
+_backend_env = os.environ.get("TABMAT_BACKEND", "cpp").lower()
+_BACKEND: Literal["cpp", "rust"] = "rust" if _backend_env == "rust" else "cpp"
+
+# Backend modules (lazy loaded)
+_cpp_backend = None
+_rust_backend = None
+
+
+def _load_cpp_backend():
+    """Load the C++ (Cython) backend."""
+    global _cpp_backend
+    if _cpp_backend is None:
+        from tabmat.ext import (  # type: ignore[attr-defined]
+            categorical as _cpp_categorical,
+        )
+        from tabmat.ext import sparse as _cpp_sparse  # type: ignore[attr-defined]
+        from tabmat.ext import split as _cpp_split  # type: ignore[attr-defined]
+
+        class CppBackend:
+            """C++ backend wrapper."""
+
+            # From categorical.pyx
+            transpose_matvec_fast = staticmethod(_cpp_categorical.transpose_matvec_fast)
+            transpose_matvec_complex = staticmethod(
+                _cpp_categorical.transpose_matvec_complex
+            )
+            matvec_fast = staticmethod(_cpp_categorical.matvec_fast)
+            matvec_complex = staticmethod(_cpp_categorical.matvec_complex)
+            sandwich_categorical_fast = staticmethod(
+                _cpp_categorical.sandwich_categorical_fast
+            )
+            sandwich_categorical_complex = staticmethod(
+                _cpp_categorical.sandwich_categorical_complex
+            )
+            get_col_included = staticmethod(_cpp_categorical.get_col_included)
+            multiply_complex = staticmethod(_cpp_categorical.multiply_complex)
+            subset_categorical_complex = staticmethod(
+                _cpp_categorical.subset_categorical_complex
+            )
+
+            # From split.pyx
+            sandwich_cat_dense = staticmethod(_cpp_split.sandwich_cat_dense)
+            sandwich_cat_cat = staticmethod(_cpp_split.sandwich_cat_cat)
+
+            # From sparse.pyx
+            csr_matvec_unrestricted = staticmethod(_cpp_sparse.csr_matvec_unrestricted)
+            csr_matvec = staticmethod(_cpp_sparse.csr_matvec)
+            csc_rmatvec_unrestricted = staticmethod(
+                _cpp_sparse.csc_rmatvec_unrestricted
+            )
+            csc_rmatvec = staticmethod(_cpp_sparse.csc_rmatvec)
+            sparse_sandwich = staticmethod(_cpp_sparse.sparse_sandwich)
+            csr_dense_sandwich = staticmethod(_cpp_sparse.csr_dense_sandwich)
+            transpose_square_dot_weights = staticmethod(
+                _cpp_sparse.transpose_square_dot_weights
+            )
+
+            name = "cpp"
+
+        _cpp_backend = CppBackend()
+    return _cpp_backend
+
+
+def _load_rust_backend():
+    """Load the Rust backend."""
+    global _rust_backend
+    if _rust_backend is None:
+        try:
+            from tabmat.tabmat_rust_ext import tabmat_rust_ext as _rust_ext
+        except ImportError as e:
+            raise ImportError(
+                "Rust backend not available. Build the Rust extension "
+                "or use set_backend('cpp')."
+            ) from e
+
+        class RustBackend:
+            """Rust backend wrapper.
+
+            The Rust backend uses unified functions with drop_first parameters.
+            These wrapper methods provide API compatibility with the C++ backend.
+
+            Note: Some functions still use the C++ backend because they're not
+            performance-critical or not yet implemented in Rust.
+            """
+
+            # Fall back to C++ for functions not yet in Rust
+            @staticmethod
+            def get_col_included(cols, n_cols):
+                return _load_cpp_backend().get_col_included(cols, n_cols)
+
+            @staticmethod
+            def multiply_complex(indices, d, ncols, dtype, drop_first):
+                return _load_cpp_backend().multiply_complex(
+                    indices, d, ncols, dtype, drop_first
+                )
+
+            @staticmethod
+            def subset_categorical_complex(indices, ncols, drop_first):
+                return _load_cpp_backend().subset_categorical_complex(
+                    indices, ncols, drop_first
+                )
+
+            # Wrappers for transpose_matvec that match C++ API
+            @staticmethod
+            def transpose_matvec_fast(
+                indices, other, n_cols, dtype, rows, other_col, out
+            ):
+                _rust_ext.transpose_matvec(indices, other, out, False)
+
+            @staticmethod
+            def transpose_matvec_complex(
+                indices, other, n_cols, dtype, rows, other_col, out, drop_first
+            ):
+                _rust_ext.transpose_matvec(indices, other, out, drop_first)
+
+            # Wrappers for matvec that match C++ API
+            @staticmethod
+            def matvec_fast(indices, other, n_rows, rows, n_cols, out):
+                _rust_ext.matvec(indices, other, out, False)
+
+            @staticmethod
+            def matvec_complex(indices, other, n_rows, rows, n_cols, out, drop_first):
+                _rust_ext.matvec(indices, other, out, drop_first)
+
+            # Wrappers for sandwich_categorical that match C++ API
+            @staticmethod
+            def sandwich_categorical_fast(indices, d, rows, dtype, n_cols):
+                return _rust_ext.sandwich_categorical(indices, d, rows, n_cols, False)
+
+            @staticmethod
+            def sandwich_categorical_complex(
+                indices, d, rows, dtype, n_cols, drop_first
+            ):
+                return _rust_ext.sandwich_categorical(
+                    indices, d, rows, n_cols, drop_first
+                )
+
+            # Wrapper for sandwich_cat_dense that matches C++ API
+            @staticmethod
+            def sandwich_cat_dense(
+                i_indices,
+                i_ncol,
+                d,
+                mat_j,
+                rows,
+                j_cols,
+                is_c_contiguous,
+                has_missings,
+                drop_first,
+            ):
+                return _rust_ext.sandwich_cat_dense(
+                    i_indices,
+                    d,
+                    mat_j,
+                    rows,
+                    j_cols,
+                    i_ncol,
+                    is_c_contiguous,
+                    drop_first,
+                )
+
+            # Wrapper for sandwich_cat_cat that matches C++ API
+            @staticmethod
+            def sandwich_cat_cat(
+                i_indices,
+                j_indices,
+                i_ncol,
+                j_ncol,
+                d,
+                rows,
+                dtype,
+                i_drop_first,
+                j_drop_first,
+                i_has_missings,
+                j_has_missings,
+            ):
+                return _rust_ext.sandwich_cat_cat(
+                    i_indices,
+                    j_indices,
+                    d,
+                    rows,
+                    i_ncol,
+                    j_ncol,
+                    i_drop_first,
+                    j_drop_first,
+                )
+
+            # Sparse operations
+            # Note: Rust backend only supports float64, so we convert if needed
+            @staticmethod
+            def csr_matvec_unrestricted(X, v, out, X_indices):
+                import numpy as np
+
+                orig_dtype = X.dtype
+                if orig_dtype != np.float64:
+                    data = X.data.astype(np.float64)
+                    v = v.astype(np.float64)
+                else:
+                    data = X.data
+
+                if out is None:
+                    out = np.zeros(X.shape[0], dtype=np.float64)
+                elif out.dtype != np.float64:
+                    out = out.astype(np.float64)
+                _rust_ext.csr_matvec_unrestricted(data, X.indices, X.indptr, v, out)
+                if orig_dtype != np.float64:
+                    out = out.astype(orig_dtype)
+                return out
+
+            @staticmethod
+            def csr_matvec(X, v, rows, cols):
+                import numpy as np
+
+                orig_dtype = X.dtype
+                if orig_dtype != np.float64:
+                    data = X.data.astype(np.float64)
+                    v = v.astype(np.float64)
+                else:
+                    data = X.data
+                result = _rust_ext.csr_matvec(
+                    data, X.indices, X.indptr, v, rows, cols, X.shape[1]
+                )
+                if orig_dtype != np.float64:
+                    result = result.astype(orig_dtype)
+                return result
+
+            @staticmethod
+            def csc_rmatvec_unrestricted(XT, v, out, XT_indices):
+                import numpy as np
+
+                orig_dtype = XT.dtype
+                if orig_dtype != np.float64:
+                    data = XT.data.astype(np.float64)
+                    v = v.astype(np.float64)
+                else:
+                    data = XT.data
+
+                if out is None:
+                    out = np.zeros(XT.shape[1], dtype=np.float64)
+                elif out.dtype != np.float64:
+                    out = out.astype(np.float64)
+                _rust_ext.csc_rmatvec_unrestricted(data, XT.indices, XT.indptr, v, out)
+                if orig_dtype != np.float64:
+                    out = out.astype(orig_dtype)
+                return out
+
+            @staticmethod
+            def csc_rmatvec(XT, v, rows, cols):
+                import numpy as np
+
+                orig_dtype = XT.dtype
+                if orig_dtype != np.float64:
+                    data = XT.data.astype(np.float64)
+                    v = v.astype(np.float64)
+                else:
+                    data = XT.data
+                result = _rust_ext.csc_rmatvec(
+                    data, XT.indices, XT.indptr, v, rows, cols, XT.shape[0]
+                )
+                if orig_dtype != np.float64:
+                    result = result.astype(orig_dtype)
+                return result
+
+            @staticmethod
+            def sparse_sandwich(A, AT, d, rows, cols):
+                import numpy as np
+
+                orig_dtype = A.dtype
+                if orig_dtype != np.float64:
+                    a_data = A.data.astype(np.float64)
+                    at_data = AT.data.astype(np.float64)
+                    d = d.astype(np.float64)
+                else:
+                    a_data = A.data
+                    at_data = AT.data
+                result = _rust_ext.sparse_sandwich(
+                    a_data,
+                    A.indices,
+                    A.indptr,
+                    at_data,
+                    AT.indices,
+                    AT.indptr,
+                    d,
+                    rows,
+                    cols,
+                    d.shape[0],
+                    A.shape[1],
+                )
+                if orig_dtype != np.float64:
+                    result = result.astype(orig_dtype)
+                return result
+
+            @staticmethod
+            def csr_dense_sandwich(A, B, d, rows, A_cols, B_cols):
+                import numpy as np
+
+                orig_dtype = A.dtype
+                is_c_contiguous = B.flags["C_CONTIGUOUS"]
+                if orig_dtype != np.float64:
+                    a_data = A.data.astype(np.float64)
+                    d = d.astype(np.float64)
+                    B = B.astype(np.float64)
+                else:
+                    a_data = A.data
+                result = _rust_ext.csr_dense_sandwich(
+                    a_data,
+                    A.indices,
+                    A.indptr,
+                    B,
+                    d,
+                    rows,
+                    A_cols,
+                    B_cols,
+                    A.shape[1],
+                    is_c_contiguous,
+                )
+                if orig_dtype != np.float64:
+                    result = result.astype(orig_dtype)
+                return result
+
+            @staticmethod
+            def transpose_square_dot_weights(data, indices, indptr, weights, dtype):
+                import numpy as np
+
+                orig_dtype = data.dtype
+                if orig_dtype != np.float64:
+                    data = data.astype(np.float64)
+                    weights = weights.astype(np.float64)
+                result = _rust_ext.transpose_square_dot_weights(
+                    data, indices, indptr, weights
+                )
+                if orig_dtype != np.float64:
+                    result = result.astype(orig_dtype)
+                return result
+
+            name = "rust"
+
+        _rust_backend = RustBackend()
+    return _rust_backend
+
+
+def get_backend():
+    """Get the currently active backend.
+
+    Returns:
+        Backend object with categorical matrix operations.
+    """
+    if _BACKEND == "rust":
+        return _load_rust_backend()
+    else:
+        return _load_cpp_backend()
+
+
+def set_backend(backend: Literal["cpp", "rust"]):
+    """Set the backend for categorical matrix operations.
+
+    Args:
+        backend: Either "cpp" for C++ (Cython) or "rust" for Rust backend.
+    """
+    global _BACKEND
+    if backend not in ("cpp", "rust"):
+        raise ValueError(f"Unknown backend: {backend}. Use 'cpp' or 'rust'.")
+    _BACKEND = backend
+
+
+def available_backends():
+    """Return list of available backends.
+
+    Returns:
+        List of backend names that are available.
+    """
+    available = ["cpp"]  # C++ is always available
+
+    try:
+        _load_rust_backend()
+        available.append("rust")
+    except ImportError:
+        pass
+
+    return available
diff --git a/src/tabmat/ext/categorical_dispatch.py b/src/tabmat/ext/categorical_dispatch.py
new file mode 100644
index 00000000..6e293cca
--- /dev/null
+++ b/src/tabmat/ext/categorical_dispatch.py
@@ -0,0 +1,195 @@
+"""Categorical operations with automatic backend dispatch.
+
+This module provides the same API as the C++ categorical module but automatically
+dispatches to either the C++ or Rust backend based on the current setting.
+
+Note: The Rust backend has simpler function signatures without row/column
+restrictions. When restrictions are present, we fall back to the C++ backend.
+"""
+
+import os
+from typing import Literal
+
+import numpy as np
+
+# Global backend setting
+_backend_env = os.environ.get("TABMAT_BACKEND", "cpp").lower()
+_BACKEND: Literal["cpp", "rust"] = "rust" if _backend_env == "rust" else "cpp"
+
+
+def set_backend(backend: Literal["cpp", "rust"]):
+    """Set the backend for categorical matrix operations."""
+    global _BACKEND
+    if backend not in ("cpp", "rust"):
+        raise ValueError(f"Unknown backend: {backend}. Use 'cpp' or 'rust'.")
+    _BACKEND = backend
+
+
+def get_backend_name() -> str:
+    """Get the name of the current backend."""
+    return _BACKEND
+
+
+# Lazy-loaded backend modules
+_cpp_module = None
+_rust_module = None
+
+
+def _get_cpp():
+    global _cpp_module
+    if _cpp_module is None:
+        from tabmat.ext import categorical as _cpp_module  # type: ignore[attr-defined]
+    return _cpp_module
+
+
+def _get_rust():
+    global _rust_module
+    if _rust_module is None:
+        from tabmat.tabmat_rust_ext import tabmat_rust_ext as _rust_module
+    return _rust_module
+
+
+def transpose_matvec_fast(indices, other, n_cols, dtype, rows, cols, out):
+    """Transpose matrix-vector multiplication (fast path, no drop_first)."""
+    # Rust doesn't support row/col restrictions, fall back to C++ if present
+    has_restrictions = rows is not None or cols is not None
+    if _BACKEND == "rust" and not has_restrictions:
+        rust = _get_rust()
+        orig_dtype = other.dtype
+        if orig_dtype != np.float64:
+            other = other.astype(np.float64)
+            out_f64 = out.astype(np.float64)
+        else:
+            out_f64 = out
+
+        rust.transpose_matvec(indices, other, out_f64, False)
+
+        if orig_dtype != np.float64:
+            out[:] = out_f64.astype(orig_dtype)
+        else:
+            out[:] = out_f64
+    else:
+        _get_cpp().transpose_matvec_fast(indices, other, n_cols, dtype, rows, cols, out)
+
+
+def transpose_matvec_complex(
+    indices, other, n_cols, dtype, rows, cols, out, drop_first
+):
+    """Transpose matrix-vector multiplication (complex path, with drop_first)."""
+    # Rust doesn't support row/col restrictions, fall back to C++ if present
+    has_restrictions = rows is not None or cols is not None
+    if _BACKEND == "rust" and not has_restrictions:
+        rust = _get_rust()
+        orig_dtype = other.dtype
+        if orig_dtype != np.float64:
+            other = other.astype(np.float64)
+            out_f64 = out.astype(np.float64)
+        else:
+            out_f64 = out
+
+        rust.transpose_matvec(indices, other, out_f64, drop_first)
+
+        if orig_dtype != np.float64:
+            out[:] = out_f64.astype(orig_dtype)
+        else:
+            out[:] = out_f64
+    else:
+        _get_cpp().transpose_matvec_complex(
+            indices, other, n_cols, dtype, rows, cols, out, drop_first
+        )
+
+
+def matvec_fast(indices, other, n_rows, cols, n_cols, out):
+    """Matrix-vector multiplication (fast path, no drop_first)."""
+    # Rust doesn't support col restrictions, fall back to C++ if present
+    if _BACKEND == "rust" and cols is None:
+        rust = _get_rust()
+        orig_dtype = other.dtype
+        if orig_dtype != np.float64:
+            other = other.astype(np.float64)
+            out_f64 = out.astype(np.float64)
+        else:
+            out_f64 = out
+
+        rust.matvec(indices, other, out_f64, False)
+
+        if orig_dtype != np.float64:
+            out[:] = out_f64.astype(orig_dtype)
+        else:
+            out[:] = out_f64
+    else:
+        _get_cpp().matvec_fast(indices, other, n_rows, cols, n_cols, out)
+
+
+def matvec_complex(indices, other, n_rows, cols, n_cols, out, drop_first):
+    """Matrix-vector multiplication (complex path, with drop_first)."""
+    # Rust doesn't support col restrictions, fall back to C++ if present
+    if _BACKEND == "rust" and cols is None:
+        rust = _get_rust()
+        orig_dtype = other.dtype
+        if orig_dtype != np.float64:
+            other = other.astype(np.float64)
+            out_f64 = out.astype(np.float64)
+        else:
+            out_f64 = out
+
+        rust.matvec(indices, other, out_f64, drop_first)
+
+        if orig_dtype != np.float64:
+            out[:] = out_f64.astype(orig_dtype)
+        else:
+            out[:] = out_f64
+    else:
+        _get_cpp().matvec_complex(indices, other, n_rows, cols, n_cols, out, drop_first)
+
+
+def sandwich_categorical_fast(indices, d, rows, dtype, n_cols):
+    """Sandwich product for categorical matrix (fast path, no drop_first)."""
+    if _BACKEND == "rust":
+        rust = _get_rust()
+        orig_dtype = d.dtype
+        if orig_dtype != np.float64:
+            d = d.astype(np.float64)
+
+        result = rust.sandwich_categorical(indices, d, rows, n_cols, False)
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().sandwich_categorical_fast(indices, d, rows, dtype, n_cols)
+
+
+def sandwich_categorical_complex(indices, d, rows, dtype, n_cols, drop_first):
+    """Sandwich product for categorical matrix (complex path, with drop_first)."""
+    if _BACKEND == "rust":
+        rust = _get_rust()
+        orig_dtype = d.dtype
+        if orig_dtype != np.float64:
+            d = d.astype(np.float64)
+
+        result = rust.sandwich_categorical(indices, d, rows, n_cols, drop_first)
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().sandwich_categorical_complex(
+            indices, d, rows, dtype, n_cols, drop_first
+        )
+
+
+def multiply_complex(indices, d, ncols, dtype, drop_first):
+    """Multiply a CategoricalMatrix by a vector d.
+
+    Note: This function is not implemented in Rust, always uses C++ backend.
+    """
+    return _get_cpp().multiply_complex(indices, d, ncols, dtype, drop_first)
+
+
+def subset_categorical_complex(indices, ncols, drop_first):
+    """Construct inputs to transform a CategoricalMatrix into a csr_matrix.
+
+    Note: This function is not implemented in Rust, always uses C++ backend.
+    """
+    return _get_cpp().subset_categorical_complex(indices, ncols, drop_first)
diff --git a/src/tabmat/ext/dense_dispatch.py b/src/tabmat/ext/dense_dispatch.py
new file mode 100644
index 00000000..24cf04cd
--- /dev/null
+++ b/src/tabmat/ext/dense_dispatch.py
@@ -0,0 +1,123 @@
+"""Dense operations with automatic backend dispatch.
+
+This module provides the same API as the C++ dense module but automatically
+dispatches to either the C++ or Rust backend based on the current setting.
+"""
+
+import os
+from typing import Literal
+
+import numpy as np
+
+# Global backend setting
+_backend_env = os.environ.get("TABMAT_BACKEND", "cpp").lower()
+_BACKEND: Literal["cpp", "rust"] = "rust" if _backend_env == "rust" else "cpp"
+
+
+def set_backend(backend: Literal["cpp", "rust"]):
+    """Set the backend for dense matrix operations."""
+    global _BACKEND
+    if backend not in ("cpp", "rust"):
+        raise ValueError(f"Unknown backend: {backend}. Use 'cpp' or 'rust'.")
+    _BACKEND = backend
+
+
+def get_backend_name() -> str:
+    """Get the name of the current backend."""
+    return _BACKEND
+
+
+# Lazy-loaded backend modules
+_cpp_module = None
+_rust_module = None
+
+
+def _get_cpp():
+    global _cpp_module
+    if _cpp_module is None:
+        from tabmat.ext import dense as _cpp_module  # type: ignore[attr-defined]
+    return _cpp_module
+
+
+def _get_rust():
+    global _rust_module
+    if _rust_module is None:
+        from tabmat.tabmat_rust_ext import tabmat_rust_ext as _rust_module
+    return _rust_module
+
+
+def dense_sandwich(X, d, rows, cols, thresh1d=32, kratio=16, innerblock=128):
+    """Dense sandwich product: X.T @ diag(d) @ X."""
+    if _BACKEND == "rust":
+        rust = _get_rust()
+        orig_dtype = X.dtype
+
+        if orig_dtype != np.float64:
+            X = X.astype(np.float64)
+            d = d.astype(np.float64)
+
+        result = rust.dense_sandwich(X, d, rows, cols)
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().dense_sandwich(X, d, rows, cols, thresh1d, kratio, innerblock)
+
+
+def dense_rmatvec(X, v, rows, cols):
+    """Dense transpose matrix-vector multiplication: X.T @ v."""
+    if _BACKEND == "rust":
+        rust = _get_rust()
+        orig_dtype = X.dtype
+
+        if orig_dtype != np.float64:
+            X = X.astype(np.float64)
+            v = v.astype(np.float64)
+
+        result = rust.dense_rmatvec(X, v, rows, cols)
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().dense_rmatvec(X, v, rows, cols)
+
+
+def dense_matvec(X, v, rows, cols):
+    """Dense matrix-vector multiplication: X @ v."""
+    if _BACKEND == "rust":
+        rust = _get_rust()
+        orig_dtype = X.dtype
+
+        if orig_dtype != np.float64:
+            X = X.astype(np.float64)
+            v = v.astype(np.float64)
+
+        result = rust.dense_matvec(X, v, rows, cols)
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().dense_matvec(X, v, rows, cols)
+
+
+def transpose_square_dot_weights(X, weights, shift):
+    """Compute weighted squared column norms with shift for dense matrix."""
+    if _BACKEND == "rust":
+        rust = _get_rust()
+        orig_dtype = X.dtype
+
+        if orig_dtype != np.float64:
+            X = X.astype(np.float64)
+            weights = weights.astype(np.float64)
+            shift = shift.astype(np.float64)
+
+        result = rust.dense_transpose_square_dot_weights(X, weights, shift)
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().transpose_square_dot_weights(X, weights, shift)
diff --git a/src/tabmat/ext/sparse_dispatch.py b/src/tabmat/ext/sparse_dispatch.py
new file mode 100644
index 00000000..6b4ae86d
--- /dev/null
+++ b/src/tabmat/ext/sparse_dispatch.py
@@ -0,0 +1,239 @@
+"""Sparse operations with automatic backend dispatch.
+
+This module provides the same API as the C++ sparse module but automatically
+dispatches to either the C++ or Rust backend based on the current setting.
+"""
+
+import os
+from typing import Literal
+
+# Global backend setting
+_backend_env = os.environ.get("TABMAT_BACKEND", "cpp").lower()
+_BACKEND: Literal["cpp", "rust"] = "rust" if _backend_env == "rust" else "cpp"
+
+
+def set_backend(backend: Literal["cpp", "rust"]):
+    """Set the backend for sparse matrix operations."""
+    global _BACKEND
+    if backend not in ("cpp", "rust"):
+        raise ValueError(f"Unknown backend: {backend}. Use 'cpp' or 'rust'.")
+    _BACKEND = backend
+
+
+def get_backend_name() -> str:
+    """Get the name of the current backend."""
+    return _BACKEND
+
+
+# Lazy-loaded backend modules
+_cpp_module = None
+_rust_module = None
+
+
+def _get_cpp():
+    global _cpp_module
+    if _cpp_module is None:
+        from tabmat.ext import sparse as _cpp_module  # type: ignore[attr-defined]
+    return _cpp_module
+
+
+def _get_rust():
+    global _rust_module
+    if _rust_module is None:
+        from tabmat.tabmat_rust_ext import tabmat_rust_ext as _rust_module
+    return _rust_module
+
+
+def csr_matvec_unrestricted(X, v, out, X_indices):
+    """CSR matrix-vector multiplication: out += X @ v."""
+    if _BACKEND == "rust":
+        import numpy as np
+
+        rust = _get_rust()
+        orig_dtype = X.dtype
+        if orig_dtype != np.float64:
+            data = X.data.astype(np.float64)
+            v = v.astype(np.float64)
+        else:
+            data = X.data
+
+        if out is None:
+            out = np.zeros(X.shape[0], dtype=np.float64)
+        elif out.dtype != np.float64:
+            out = out.astype(np.float64)
+
+        rust.csr_matvec_unrestricted(data, X.indices, X.indptr, v, out)
+
+        if orig_dtype != np.float64:
+            out = out.astype(orig_dtype)
+        return out
+    else:
+        return _get_cpp().csr_matvec_unrestricted(X, v, out, X_indices)
+
+
+def csr_matvec(X, v, rows, cols):
+    """CSR matrix-vector multiplication with row/column restrictions."""
+    if _BACKEND == "rust":
+        import numpy as np
+
+        rust = _get_rust()
+        orig_dtype = X.dtype
+        if orig_dtype != np.float64:
+            data = X.data.astype(np.float64)
+            v = v.astype(np.float64)
+        else:
+            data = X.data
+
+        result = rust.csr_matvec(data, X.indices, X.indptr, v, rows, cols, X.shape[1])
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().csr_matvec(X, v, rows, cols)
+
+
+def csc_rmatvec_unrestricted(XT, v, out, XT_indices):
+    """CSC transpose matrix-vector multiplication: out += XT.T @ v."""
+    if _BACKEND == "rust":
+        import numpy as np
+
+        rust = _get_rust()
+        orig_dtype = XT.dtype
+        if orig_dtype != np.float64:
+            data = XT.data.astype(np.float64)
+            v = v.astype(np.float64)
+        else:
+            data = XT.data
+
+        if out is None:
+            out = np.zeros(XT.shape[1], dtype=np.float64)
+        elif out.dtype != np.float64:
+            out = out.astype(np.float64)
+
+        rust.csc_rmatvec_unrestricted(data, XT.indices, XT.indptr, v, out)
+
+        if orig_dtype != np.float64:
+            out = out.astype(orig_dtype)
+        return out
+    else:
+        return _get_cpp().csc_rmatvec_unrestricted(XT, v, out, XT_indices)
+
+
+def csc_rmatvec(XT, v, rows, cols):
+    """CSC transpose matrix-vector multiplication with restrictions."""
+    if _BACKEND == "rust":
+        import numpy as np
+
+        rust = _get_rust()
+        orig_dtype = XT.dtype
+        if orig_dtype != np.float64:
+            data = XT.data.astype(np.float64)
+            v = v.astype(np.float64)
+        else:
+            data = XT.data
+
+        result = rust.csc_rmatvec(
+            data, XT.indices, XT.indptr, v, rows, cols, XT.shape[0]
+        )
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().csc_rmatvec(XT, v, rows, cols)
+
+
+def sparse_sandwich(A, AT, d, rows, cols):
+    """Sparse sandwich product: A.T @ diag(d) @ A."""
+    if _BACKEND == "rust":
+        import numpy as np
+
+        rust = _get_rust()
+        orig_dtype = A.dtype
+        if orig_dtype != np.float64:
+            a_data = A.data.astype(np.float64)
+            at_data = AT.data.astype(np.float64)
+            d = d.astype(np.float64)
+        else:
+            a_data = A.data
+            at_data = AT.data
+
+        result = rust.sparse_sandwich(
+            a_data,
+            A.indices,
+            A.indptr,
+            at_data,
+            AT.indices,
+            AT.indptr,
+            d,
+            rows,
+            cols,
+            d.shape[0],
+            A.shape[1],
+        )
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().sparse_sandwich(A, AT, d, rows, cols)
+
+
+def csr_dense_sandwich(A, B, d, rows, A_cols, B_cols):
+    """CSR-dense sandwich: A.T @ diag(d) @ B."""
+    if _BACKEND == "rust":
+        import numpy as np
+
+        rust = _get_rust()
+        orig_dtype = A.dtype
+        is_c_contiguous = B.flags["C_CONTIGUOUS"]
+
+        if orig_dtype != np.float64:
+            a_data = A.data.astype(np.float64)
+            d = d.astype(np.float64)
+            B = B.astype(np.float64)
+        else:
+            a_data = A.data
+
+        result = rust.csr_dense_sandwich(
+            a_data,
+            A.indices,
+            A.indptr,
+            B,
+            d,
+            rows,
+            A_cols,
+            B_cols,
+            A.shape[1],
+            is_c_contiguous,
+        )
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().csr_dense_sandwich(A, B, d, rows, A_cols, B_cols)
+
+
+def transpose_square_dot_weights(data, indices, indptr, weights, dtype):
+    """Compute weighted squared column norms for CSC matrix."""
+    if _BACKEND == "rust":
+        import numpy as np
+
+        rust = _get_rust()
+        orig_dtype = data.dtype
+
+        if orig_dtype != np.float64:
+            data = data.astype(np.float64)
+            weights = weights.astype(np.float64)
+
+        result = rust.transpose_square_dot_weights(data, indices, indptr, weights)
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().transpose_square_dot_weights(
+            data, indices, indptr, weights, dtype
+        )
diff --git a/src/tabmat/ext/split_dispatch.py b/src/tabmat/ext/split_dispatch.py
new file mode 100644
index 00000000..9cf2fc7d
--- /dev/null
+++ b/src/tabmat/ext/split_dispatch.py
@@ -0,0 +1,140 @@
+"""Split matrix operations with automatic backend dispatch.
+
+This module provides the same API as the C++ split module but automatically
+dispatches to either the C++ or Rust backend based on the current setting.
+"""
+
+import os
+from typing import Literal
+
+import numpy as np
+
+# Global backend setting
+_backend_env = os.environ.get("TABMAT_BACKEND", "cpp").lower()
+_BACKEND: Literal["cpp", "rust"] = "rust" if _backend_env == "rust" else "cpp"
+
+
+def set_backend(backend: Literal["cpp", "rust"]):
+    """Set the backend for split matrix operations."""
+    global _BACKEND
+    if backend not in ("cpp", "rust"):
+        raise ValueError(f"Unknown backend: {backend}. Use 'cpp' or 'rust'.")
+    _BACKEND = backend
+
+
+def get_backend_name() -> str:
+    """Get the name of the current backend."""
+    return _BACKEND
+
+
+# Lazy-loaded backend modules
+_cpp_module = None
+_rust_module = None
+
+
+def _get_cpp():
+    global _cpp_module
+    if _cpp_module is None:
+        from tabmat.ext import split as _cpp_module  # type: ignore[attr-defined]
+    return _cpp_module
+
+
+def _get_rust():
+    global _rust_module
+    if _rust_module is None:
+        from tabmat.tabmat_rust_ext import tabmat_rust_ext as _rust_module
+    return _rust_module
+
+
+def sandwich_cat_cat(
+    i_indices,
+    j_indices,
+    i_ncol,
+    j_ncol,
+    d,
+    rows,
+    dtype,
+    i_drop_first,
+    j_drop_first,
+    i_has_missings,
+    j_has_missings,
+):
+    """Cross-sandwich product between two categorical matrices."""
+    if _BACKEND == "rust":
+        rust = _get_rust()
+        orig_dtype = d.dtype
+        if orig_dtype != np.float64:
+            d = d.astype(np.float64)
+
+        result = rust.sandwich_cat_cat(
+            i_indices, j_indices, d, rows, i_ncol, j_ncol, i_drop_first, j_drop_first
+        )
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().sandwich_cat_cat(
+            i_indices,
+            j_indices,
+            i_ncol,
+            j_ncol,
+            d,
+            rows,
+            dtype,
+            i_drop_first,
+            j_drop_first,
+            i_has_missings,
+            j_has_missings,
+        )
+
+
+def sandwich_cat_dense(
+    i_indices,
+    i_ncol,
+    d,
+    mat_j,
+    rows,
+    j_cols,
+    is_c_contiguous,
+    has_missings,
+    drop_first,
+):
+    """Cross-sandwich product between categorical and dense matrices."""
+    if _BACKEND == "rust":
+        rust = _get_rust()
+        orig_dtype = d.dtype
+        if orig_dtype != np.float64:
+            d = d.astype(np.float64)
+            mat_j = mat_j.astype(np.float64)
+
+        result = rust.sandwich_cat_dense(
+            i_indices, d, mat_j, rows, j_cols, i_ncol, is_c_contiguous, drop_first
+        )
+
+        if orig_dtype != np.float64:
+            result = result.astype(orig_dtype)
+        return result
+    else:
+        return _get_cpp().sandwich_cat_dense(
+            i_indices,
+            i_ncol,
+            d,
+            mat_j,
+            rows,
+            j_cols,
+            is_c_contiguous,
+            has_missings,
+            drop_first,
+        )
+
+
+# Re-export functions that don't have Rust implementations
+def split_col_subsets(self, cols):
+    """Split column subsets - only available in C++ backend."""
+    return _get_cpp().split_col_subsets(self, cols)
+
+
+def is_sorted(a):
+    """Check if array is sorted - only available in C++ backend."""
+    return _get_cpp().is_sorted(a)
diff --git a/src/tabmat/sparse_matrix.py b/src/tabmat/sparse_matrix.py
index 324c2be4..9d31e99b 100644
--- a/src/tabmat/sparse_matrix.py
+++ b/src/tabmat/sparse_matrix.py
@@ -3,7 +3,7 @@
 import numpy as np
 from scipy import sparse as sps
 
-from .ext.sparse import (
+from .ext.sparse_dispatch import (
     csc_rmatvec,
     csc_rmatvec_unrestricted,
     csr_dense_sandwich,
diff --git a/src/tabmat/split_matrix.py b/src/tabmat/split_matrix.py
index 68a3d22e..223e0bb1 100644
--- a/src/tabmat/split_matrix.py
+++ b/src/tabmat/split_matrix.py
@@ -6,7 +6,7 @@
 from scipy import sparse as sps
 
 from .dense_matrix import DenseMatrix
-from .ext.split import is_sorted, split_col_subsets
+from .ext.split_dispatch import is_sorted, split_col_subsets
 from .matrix_base import MatrixBase
 from .sparse_matrix import SparseMatrix
 from .standardized_mat import StandardizedMatrix
diff --git a/src/tabmat/tabmat_rust_ext/__init__.py b/src/tabmat/tabmat_rust_ext/__init__.py
new file mode 100644
index 00000000..f458a7d7
--- /dev/null
+++ b/src/tabmat/tabmat_rust_ext/__init__.py
@@ -0,0 +1,35 @@
+from . import tabmat_rust_ext  # type: ignore[attr-defined]
+
+__all__ = [
+    "tabmat_rust_ext",
+    # Categorical operations
+    "transpose_matvec",
+    "matvec",
+    "sandwich_categorical",
+    "sandwich_cat_cat",
+    "sandwich_cat_dense",
+    # Sparse operations
+    "csr_matvec_unrestricted",
+    "csr_matvec",
+    "csc_rmatvec_unrestricted",
+    "csc_rmatvec",
+    "sparse_sandwich",
+    "csr_dense_sandwich",
+    "transpose_square_dot_weights",
+]
+
+# Re-export categorical functions
+transpose_matvec = tabmat_rust_ext.transpose_matvec
+matvec = tabmat_rust_ext.matvec
+sandwich_categorical = tabmat_rust_ext.sandwich_categorical
+sandwich_cat_cat = tabmat_rust_ext.sandwich_cat_cat
+sandwich_cat_dense = tabmat_rust_ext.sandwich_cat_dense
+
+# Re-export sparse functions
+csr_matvec_unrestricted = tabmat_rust_ext.csr_matvec_unrestricted
+csr_matvec = tabmat_rust_ext.csr_matvec
+csc_rmatvec_unrestricted = tabmat_rust_ext.csc_rmatvec_unrestricted
+csc_rmatvec = tabmat_rust_ext.csc_rmatvec
+sparse_sandwich = tabmat_rust_ext.sparse_sandwich
+csr_dense_sandwich = tabmat_rust_ext.csr_dense_sandwich
+transpose_square_dot_weights = tabmat_rust_ext.transpose_square_dot_weights