diff --git a/include/xsf/.clang-format b/.clang-format
similarity index 100%
rename from include/xsf/.clang-format
rename to .clang-format
diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml
index c0d7aac52c..c18a5d5d13 100644
--- a/.github/workflows/linux.yml
+++ b/.github/workflows/linux.yml
@@ -1,10 +1,65 @@
 name: Linux Tests
 
-on:
-  push:
-    branches:
-      - main
-      - catch-test-setup
-    pull_request:
-      - main
-      - catch-test-setup
+on: [ push, pull_request ]
+
+permissions:
+   contents: read  # to fetch code (actions/checkout)
+
+env:
+  CCACHE_DIR: "${{ github.workspace }}/.ccache"
+  INSTALLDIR: "build-install"
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.head_ref || github.run_id }}
+  cancel-in-progress: true
+
+jobs:
+  gcc_x86_64_linux:
+    name: GCC X86_64 Linux
+    runs-on: ubuntu-22.04
+    env:
+      XSREF_TABLES_PATH: "${{ github.workspace }}/xsref/tables"
+
+    steps:
+      - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2
+      - name: Checkout xsref
+        uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683
+        with:
+          repository: 'scipy/xsref'
+          path: 'xsref'
+          ref: 'main'
+
+      - name: Checkout Catch2
+        uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683
+        with:
+          repository: 'catchorg/Catch2'
+          path: 'Catch2'
+          ref: 'v3.8.0'
+
+      - name: Install Ubuntu dependencies
+        run: |
+          sudo apt-get update
+          wget https://apache.jfrog.io/artifactory/arrow/$(lsb_release --id --short | tr 'A-Z' 'a-z')/apache-arrow-apt-source-latest-$(lsb_release --codename --short).deb
+          sudo apt-get install -y -V ./apache-arrow-apt-source-latest-$(lsb_release --codename --short).deb
+          sudo apt-get update
+          sudo apt-get install -y -V libarrow-dev # For C++
+          sudo apt-get install -y -V libarrow-glib-dev # For GLib (C)
+          sudo apt-get install -y -V libarrow-dataset-dev # For Apache Arrow Dataset C++
+          sudo apt-get install -y -V libparquet-dev # For Apache Parquet C++
+          sudo apt-get install -y -V libparquet-glib-dev # For Apache Parquet GLib (C)
+          sudo apt-get install -y libzstd-dev ccache
+
+      - name: Install Catch2
+        run: |
+          cmake -Bbuild -H. -DBUILD_TESTING=OFF
+          sudo cmake --build build/ --target install
+        working-directory: ./Catch2
+
+      - name: Configure and build
+        run: |
+          cmake .
+          make
+
+      - name: run tests
+        run: |
+          ctest --output-on-failure
diff --git a/.gitignore b/.gitignore
index 259148fa18..7ceaf2f66a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -30,3 +30,22 @@
 *.exe
 *.out
 *.app
+
+# CMake stuff
+CMakeLists.txt.user
+CMakeCache.txt
+CMakeFiles
+CMakeScripts
+Testing
+Makefile
+cmake_install.cmake
+install_manifest.txt
+compile_commands.json
+CTestTestfile.cmake
+_deps
+CMakeUserPresets.json
+*.cmake
+DartConfiguration.tcl
+
+# test executable extension
+*.test
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000000..ca1ae8c746
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 3.14)
+project(xsf)
+
+# Scipy is restricted to C++17
+# https://docs.scipy.org/doc/scipy/dev/toolchain.html#c-language-standards
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+# Tests
+enable_testing()
+add_subdirectory(tests)
diff --git a/include/xsf/cephes/cbrt.h b/include/xsf/cephes/cbrt.h
index 3e9fbd4eab..ee7c6284bc 100644
--- a/include/xsf/cephes/cbrt.h
+++ b/include/xsf/cephes/cbrt.h
@@ -52,80 +52,81 @@ namespace cephes {
 
     namespace detail {
 
-        constexpr double CBRT2 = 1.2599210498948731647672;
-        constexpr double CBRT4 = 1.5874010519681994747517;
-        constexpr double CBRT2I = 0.79370052598409973737585;
-        constexpr double CBRT4I = 0.62996052494743658238361;
+        inline constexpr double CBRT2 = 1.2599210498948731647672;
+        inline constexpr double CBRT4 = 1.5874010519681994747517;
+        inline constexpr double CBRT2I = 0.79370052598409973737585;
+        inline constexpr double CBRT4I = 0.62996052494743658238361;
 
-        XSF_HOST_DEVICE inline double cbrt(double x) {
-            int e, rem, sign;
-            double z;
+    }
 
-            if (!std::isfinite(x)) {
-                return x;
-            }
-            if (x == 0) {
-                return (x);
-            }
-            if (x > 0) {
-                sign = 1;
-            } else {
-                sign = -1;
-                x = -x;
-            }
+    XSF_HOST_DEVICE inline double cbrt(double x) {
+	int e, rem, sign;
+	double z;
 
-            z = x;
-            /* extract power of 2, leaving
-             * mantissa between 0.5 and 1
-             */
-            x = std::frexp(x, &e);
+	if (!std::isfinite(x)) {
+	    return x;
+	}
+	if (x == 0) {
+	    return (x);
+	}
+	if (x > 0) {
+	    sign = 1;
+	} else {
+	    sign = -1;
+	    x = -x;
+	}
+	
+	z = x;
+	/* extract power of 2, leaving
+	 * mantissa between 0.5 and 1
+	 */
+	x = std::frexp(x, &e);
+	
+	/* Approximate cube root of number between .5 and 1,
+	 * peak relative error = 9.2e-6
+	 */
+	x = (((-1.3466110473359520655053e-1 * x + 5.4664601366395524503440e-1) * x - 9.5438224771509446525043e-1) *
+	     x +
+	     1.1399983354717293273738e0) *
+	    x +
+	    4.0238979564544752126924e-1;
+	
+	/* exponent divided by 3 */
+	if (e >= 0) {
+	    rem = e;
+	    e /= 3;
+	    rem -= 3 * e;
+	    if (rem == 1) {
+		x *= detail::CBRT2;
+	    } else if (rem == 2) {
+		x *= detail::CBRT4;
+	    }
+	}
+	/* argument less than 1 */
+	else {
+	    e = -e;
+	    rem = e;
+	    e /= 3;
+	    rem -= 3 * e;
+	    if (rem == 1) {
+		x *= detail::CBRT2I;
+	    } else if (rem == 2) {
+		x *= detail::CBRT4I;
+	    }
+	    e = -e;
+	}
 
-            /* Approximate cube root of number between .5 and 1,
-             * peak relative error = 9.2e-6
-             */
-            x = (((-1.3466110473359520655053e-1 * x + 5.4664601366395524503440e-1) * x - 9.5438224771509446525043e-1) *
-                     x +
-                 1.1399983354717293273738e0) *
-                    x +
-                4.0238979564544752126924e-1;
+	/* multiply by power of 2 */
+	x = std::ldexp(x, e);
+	
+	/* Newton iteration */
+	x -= (x - (z / (x * x))) * 0.33333333333333333333;
+	x -= (x - (z / (x * x))) * 0.33333333333333333333;
 
-            /* exponent divided by 3 */
-            if (e >= 0) {
-                rem = e;
-                e /= 3;
-                rem -= 3 * e;
-                if (rem == 1) {
-                    x *= CBRT2;
-                } else if (rem == 2) {
-                    x *= CBRT4;
-                }
-            }
-            /* argument less than 1 */
-            else {
-                e = -e;
-                rem = e;
-                e /= 3;
-                rem -= 3 * e;
-                if (rem == 1) {
-                    x *= CBRT2I;
-                } else if (rem == 2) {
-                    x *= CBRT4I;
-                }
-                e = -e;
-            }
-
-            /* multiply by power of 2 */
-            x = std::ldexp(x, e);
-
-            /* Newton iteration */
-            x -= (x - (z / (x * x))) * 0.33333333333333333333;
-            x -= (x - (z / (x * x))) * 0.33333333333333333333;
-
-            if (sign < 0)
-                x = -x;
-            return (x);
-        }
-    } // namespace detail
+	if (sign < 0)
+	    x = -x;
+	return (x);
+    }
 
 } // namespace cephes
 } // namespace xsf
diff --git a/include/xsf/cephes/exp10.h b/include/xsf/cephes/exp10.h
index 56e8e62866..2ae728c53f 100644
--- a/include/xsf/cephes/exp10.h
+++ b/include/xsf/cephes/exp10.h
@@ -57,6 +57,8 @@
 #include "../config.h"
 #include "../error.h"
 
+#include "polevl.h"
+
 namespace xsf {
 namespace cephes {
 
diff --git a/include/xsf/cephes/gdtr.h b/include/xsf/cephes/gdtr.h
index 5123e16904..dd70d951dd 100644
--- a/include/xsf/cephes/gdtr.h
+++ b/include/xsf/cephes/gdtr.h
@@ -111,7 +111,7 @@ namespace cephes {
     XSF_HOST_DEVICE inline double gdtr(double a, double b, double x) {
 
         if (x < 0.0) {
-            sf_error("gdtr", SF_ERROR_DOMAIN, NULL);
+            set_error("gdtr", SF_ERROR_DOMAIN, NULL);
             return (std::numeric_limits<double>::quiet_NaN());
         }
         return (igam(b, a * x));
@@ -129,7 +129,7 @@ namespace cephes {
     XSF_HOST_DEVICE inline double gdtri(double a, double b, double y) {
 
         if ((y < 0.0) || (y > 1.0) || (a <= 0.0) || (b < 0.0)) {
-            sf_error("gdtri", SF_ERROR_DOMAIN, NULL);
+            set_error("gdtri", SF_ERROR_DOMAIN, NULL);
             return (std::numeric_limits<double>::quiet_NaN());
         }
 
diff --git a/include/xsf/cephes/ndtr.h b/include/xsf/cephes/ndtr.h
index a3611d26ba..21bfc12210 100644
--- a/include/xsf/cephes/ndtr.h
+++ b/include/xsf/cephes/ndtr.h
@@ -135,6 +135,7 @@
 #pragma once
 
 #include "../config.h"
+#include "../error.h"
 
 #include "const.h"
 #include "polevl.h"
diff --git a/include/xsf/cephes/shichi.h b/include/xsf/cephes/shichi.h
index fcdd2d7986..f249ca4b04 100644
--- a/include/xsf/cephes/shichi.h
+++ b/include/xsf/cephes/shichi.h
@@ -151,7 +151,7 @@ namespace cephes {
     } // namespace detail
 
     /* Sine and cosine integrals */
-    XSF_HOST_DEVICE inline int shichi(double x, double *si, double *ci) {
+    XSF_HOST_DEVICE inline int shichi(double x, double &si, double &ci) {
         double k, z, c, s, a, b;
         short sign;
 
@@ -163,8 +163,8 @@ namespace cephes {
         }
 
         if (x == 0.0) {
-            *si = 0.0;
-            *ci = -std::numeric_limits<double>::infinity();
+            si = 0.0;
+            ci = -std::numeric_limits<double>::infinity();
             return (0);
         }
 
@@ -216,8 +216,8 @@ namespace cephes {
 
     asymp:
         if (x > 1000) {
-            *si = std::numeric_limits<double>::infinity();
-            *ci = std::numeric_limits<double>::infinity();
+            si = std::numeric_limits<double>::infinity();
+            ci = std::numeric_limits<double>::infinity();
         } else {
             /* Asymptotic expansions
              * http://functions.wolfram.com/GammaBetaErf/CoshIntegral/06/02/
@@ -225,11 +225,11 @@ namespace cephes {
              */
             a = detail::hyp3f0(0.5, 1, 1, 4.0 / (x * x));
             b = detail::hyp3f0(1, 1, 1.5, 4.0 / (x * x));
-            *si = std::cosh(x) / x * a + std::sinh(x) / (x * x) * b;
-            *ci = std::sinh(x) / x * a + std::cosh(x) / (x * x) * b;
+            si = std::cosh(x) / x * a + std::sinh(x) / (x * x) * b;
+            ci = std::sinh(x) / x * a + std::cosh(x) / (x * x) * b;
         }
         if (sign) {
-            *si = -*si;
+            si = -si;
         }
         return 0;
 
@@ -238,9 +238,9 @@ namespace cephes {
             s = -s;
         }
 
-        *si = s;
+        si = s;
 
-        *ci = detail::SCIPY_EULER + std::log(x) + c;
+        ci = detail::SCIPY_EULER + std::log(x) + c;
         return (0);
     }
 
diff --git a/include/xsf/cephes/sici.h b/include/xsf/cephes/sici.h
index c22612ccc9..9031d9fef6 100644
--- a/include/xsf/cephes/sici.h
+++ b/include/xsf/cephes/sici.h
@@ -138,7 +138,7 @@ namespace cephes {
 
     } // namespace detail
 
-    XSF_HOST_DEVICE inline int sici(double x, double *si, double *ci) {
+    XSF_HOST_DEVICE inline int sici(double x, double &si, double &ci) {
         double z, c, s, f, g;
         short sign;
 
@@ -150,24 +150,24 @@ namespace cephes {
         }
 
         if (x == 0.0) {
-            *si = 0.0;
-            *ci = -std::numeric_limits<double>::infinity();
+            si = 0.0;
+            ci = -std::numeric_limits<double>::infinity();
             return (0);
         }
 
         if (x > 1.0e9) {
             if (std::isinf(x)) {
                 if (sign == -1) {
-                    *si = -M_PI_2;
-                    *ci = std::numeric_limits<double>::quiet_NaN();
+                    si = -M_PI_2;
+                    ci = std::numeric_limits<double>::quiet_NaN();
                 } else {
-                    *si = M_PI_2;
-                    *ci = 0;
+                    si = M_PI_2;
+                    ci = 0;
                 }
                 return 0;
             }
-            *si = M_PI_2 - std::cos(x) / x;
-            *ci = std::sin(x) / x;
+            si = M_PI_2 - std::cos(x) / x;
+            ci = std::sin(x) / x;
         }
 
         if (x > 4.0) {
@@ -181,8 +181,8 @@ namespace cephes {
         if (sign) {
             s = -s;
         }
-        *si = s;
-        *ci = detail::SCIPY_EULER + std::log(x) + c; /* real part if x < 0 */
+        si = s;
+        ci = detail::SCIPY_EULER + std::log(x) + c; /* real part if x < 0 */
         return (0);
 
         /* The auxiliary functions are:
@@ -211,11 +211,11 @@ namespace cephes {
             f = polevl(z, detail::sici_FN8, 8) / (x * p1evl(z, detail::sici_FD8, 8));
             g = z * polevl(z, detail::sici_GN8, 8) / p1evl(z, detail::sici_GD8, 9);
         }
-        *si = M_PI_2 - f * c - g * s;
+        si = M_PI_2 - f * c - g * s;
         if (sign) {
-            *si = -(*si);
+            si = -(si);
         }
-        *ci = f * s - g * c;
+        ci = f * s - g * c;
 
         return (0);
     }
diff --git a/include/xsf/exp.h b/include/xsf/exp.h
index 0d85511090..d2ef26d9a9 100644
--- a/include/xsf/exp.h
+++ b/include/xsf/exp.h
@@ -1,7 +1,11 @@
 #pragma once
 
-#include "xsf/cephes/exp10.h"
-#include "xsf/cephes/exp2.h"
+#include "config.h"
+#include "trig.h"
+
+#include "cephes/exp10.h"
+#include "cephes/exp2.h"
+#include "cephes/unity.h"
 
 namespace xsf {
 
@@ -45,12 +49,12 @@ inline std::complex<float> expm1(std::complex<float> z) {
     return static_cast<std::complex<float>>(expm1(static_cast<std::complex<double>>(z)));
 }
 
-double exp2(double x) { return cephes::exp2(x); }
+inline double exp2(double x) { return cephes::exp2(x); }
 
-float exp2(float x) { return exp2(static_cast<double>(x)); }
+inline float exp2(float x) { return exp2(static_cast<double>(x)); }
 
-double exp10(double x) { return cephes::exp10(x); }
+inline double exp10(double x) { return cephes::exp10(x); }
 
-float exp10(float x) { return exp10(static_cast<double>(x)); }
+inline float exp10(float x) { return exp10(static_cast<double>(x)); }
 
 } // namespace xsf
diff --git a/include/xsf/mathieu.h b/include/xsf/mathieu.h
index 7fe82a5ec2..a74b6eea19 100644
--- a/include/xsf/mathieu.h
+++ b/include/xsf/mathieu.h
@@ -1,5 +1,6 @@
 #pragma once
 
+#include "error.h"
 #include "specfun/specfun.h"
 
 namespace xsf {
diff --git a/include/xsf/par_cyl.h b/include/xsf/par_cyl.h
index 6092ab2dee..0c6734d47d 100644
--- a/include/xsf/par_cyl.h
+++ b/include/xsf/par_cyl.h
@@ -1,5 +1,6 @@
 #pragma once
 
+#include "error.h"
 #include "specfun/specfun.h"
 
 namespace xsf {
diff --git a/include/xsf/sici.h b/include/xsf/sici.h
index 4d26b64e02..94812d1787 100644
--- a/include/xsf/sici.h
+++ b/include/xsf/sici.h
@@ -26,23 +26,23 @@ namespace xsf {
 namespace detail {
     
     XSF_HOST_DEVICE inline void sici_power_series(int sgn, std::complex<double> z,
-						       std::complex<double> *s, std::complex<double> *c) {
+						       std::complex<double> &s, std::complex<double> &c) {
 	/* DLMF 6.6.5 and 6.6.6. If sgn = -1 computes si/ci, and if sgn = 1
 	 * computes shi/chi.
 	 */        
 	std::complex<double> fac = z;
-	*s = fac;
-	*c = 0;
+	s = fac;
+	c = 0;
 	std::complex<double> term1, term2;
 	for (int n = 1; n < 100; n++) {
 	    fac *= static_cast<double>(sgn)*z/(2.0*n);
 	    term2 = fac/(2.0*n);
-	    *c += term2;
+	    c += term2;
 	    fac *= z/(2.0*n + 1.0);
 	    term1 = fac/(2.0*n + 1.0);
-	    *s += term1;
+	    s += term1;
 	    constexpr double tol = std::numeric_limits<double>::epsilon();
-	    if (std::abs(term1) < tol*std::abs(*s) && std::abs(term2) < tol*std::abs(*c)) {
+	    if (std::abs(term1) < tol*std::abs(s) && std::abs(term2) < tol*std::abs(c)) {
 		break;
 	    }
 	}
@@ -52,7 +52,7 @@ namespace detail {
 
     
 XSF_HOST_DEVICE inline int sici(std::complex<double> z,
-				    std::complex<double> *si, std::complex<double> *ci) {
+				    std::complex<double> &si, std::complex<double> &ci) {
     /* Compute sin/cos integrals at complex arguments. The algorithm
      * largely follows that of [1].
      */
@@ -60,13 +60,13 @@ XSF_HOST_DEVICE inline int sici(std::complex<double> z,
     constexpr double EULER = xsf::cephes::detail::SCIPY_EULER;
 
     if (z == std::numeric_limits<double>::infinity()) {
-        *si = M_PI_2;
-        *ci = 0;
+        si = M_PI_2;
+        ci = 0;
         return 0;
     }
     if (z == -std::numeric_limits<double>::infinity()) {
-	*si = -M_PI_2;
-        *ci = {0.0, M_PI};
+	si = -M_PI_2;
+        ci = {0.0, M_PI};
         return 0;
     }
 
@@ -76,9 +76,9 @@ XSF_HOST_DEVICE inline int sici(std::complex<double> z,
 
         if (z == 0.0) {
             set_error("sici", SF_ERROR_DOMAIN, NULL);
-            *ci = {-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::quiet_NaN()};
+            ci = {-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::quiet_NaN()};
         } else {
-            *ci += EULER + std::log(z);
+            ci += EULER + std::log(z);
 	}
         return 0;
     }
@@ -87,51 +87,51 @@ XSF_HOST_DEVICE inline int sici(std::complex<double> z,
     std::complex<double> jz = std::complex<double>(0.0, 1.0) * z;
     std::complex<double> term1 = expi(jz);
     std::complex<double> term2 = expi(-jz);
-    *si = std::complex<double>(0.0, -0.5)*(term1 - term2);
-    *ci = 0.5*(term1 + term2);
+    si = std::complex<double>(0.0, -0.5)*(term1 - term2);
+    ci = 0.5*(term1 + term2);
     if (z.real() == 0) {
         if (z.imag() > 0) {
-            *ci += std::complex<double>(0.0, M_PI_2);
+            ci += std::complex<double>(0.0, M_PI_2);
 	} else if (z.imag() < 0) {
-            *ci -= std::complex<double>(0.0, M_PI_2);
+            ci -= std::complex<double>(0.0, M_PI_2);
 	}
     } else if (z.real() > 0) {
-        *si -= M_PI_2;
+        si -= M_PI_2;
     } else {
-        *si += M_PI_2;
+        si += M_PI_2;
         if (z.imag() >= 0) {
-            *ci += std::complex<double>(0.0, M_PI);
+            ci += std::complex<double>(0.0, M_PI);
         } else {
-            *ci -= std::complex<double>(0.0, M_PI);
+            ci -= std::complex<double>(0.0, M_PI);
 	}
     }
     return 0;
 }
 
 XSF_HOST_DEVICE inline int sici(std::complex<float> z,
-				    std::complex<float> *si_f, std::complex<float> *ci_f) {
+				    std::complex<float> &si_f, std::complex<float> &ci_f) {
     std::complex<double> si;
     std::complex<double> ci;
-    int res = sici(z, &si, &ci);
-    *si_f = si;
-    *ci_f = ci;
+    int res = sici(z, si, ci);
+    si_f = si;
+    ci_f = ci;
     return res;
 }
 
 XSF_HOST_DEVICE inline int shichi(std::complex<double> z,
-				       std::complex<double> *shi, std::complex<double> *chi) {
+				       std::complex<double> &shi, std::complex<double> &chi) {
     /* Compute sinh/cosh integrals at complex arguments. The algorithm
      * largely follows that of [1].
      */
     constexpr double EULER = xsf::cephes::detail::SCIPY_EULER;
     if (z == std::numeric_limits<double>::infinity()) {
-        *shi = std::numeric_limits<double>::infinity();
-        *chi = std::numeric_limits<double>::infinity();
+        shi = std::numeric_limits<double>::infinity();
+        chi = std::numeric_limits<double>::infinity();
         return 0;
     }
     if (z == -std::numeric_limits<double>::infinity()) {
-        *shi = -std::numeric_limits<double>::infinity();
-        *chi = std::numeric_limits<double>::infinity();
+        shi = -std::numeric_limits<double>::infinity();
+        chi = std::numeric_limits<double>::infinity();
         return 0;
     }
     if (std::abs(z) < 0.8) {
@@ -139,62 +139,62 @@ XSF_HOST_DEVICE inline int shichi(std::complex<double> z,
 	detail::sici_power_series(1, z, shi, chi);
         if (z == 0.0) {
             set_error("shichi", SF_ERROR_DOMAIN, NULL);
-            *chi = {-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::quiet_NaN()};
+            chi = {-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::quiet_NaN()};
         } else {
-            *chi += EULER + std::log(z);
+            chi += EULER + std::log(z);
 	}
 	return 0;
     }
 
     std::complex<double> term1 = expi(z);
     std::complex<double> term2 = expi(-z);
-    *shi = 0.5*(term1 - term2);
-    *chi = 0.5*(term1 + term2);
+    shi = 0.5*(term1 - term2);
+    chi = 0.5*(term1 + term2);
     if (z.imag() > 0) {
-        *shi -= std::complex<double>(0.0, 0.5*M_PI);
-        *chi += std::complex<double>(0.0, 0.5*M_PI);
+        shi -= std::complex<double>(0.0, 0.5*M_PI);
+        chi += std::complex<double>(0.0, 0.5*M_PI);
     } else if (z.imag() < 0) {
-        *shi += std::complex<double>(0.0, 0.5*M_PI);
-        *chi -= std::complex<double>(0.0, 0.5*M_PI);
+        shi += std::complex<double>(0.0, 0.5*M_PI);
+        chi -= std::complex<double>(0.0, 0.5*M_PI);
     } else if (z.real() < 0) {
-        *chi += std::complex<double>(0.0, M_PI);
+        chi += std::complex<double>(0.0, M_PI);
     }
     return 0;
 }
 
 XSF_HOST_DEVICE inline int shichi(std::complex<float> z,
-				    std::complex<float> *shi_f, std::complex<float> *chi_f) {
+				    std::complex<float> &shi_f, std::complex<float> &chi_f) {
     std::complex<double> shi;
     std::complex<double> chi;
-    int res = shichi(z, &shi, &chi);
-    *shi_f = shi;
-    *chi_f = chi;
+    int res = shichi(z, shi, chi);
+    shi_f = shi;
+    chi_f = chi;
     return res;
 }
 
-XSF_HOST_DEVICE inline int sici(double x, double *si, double *ci) {
+XSF_HOST_DEVICE inline int sici(double x, double &si, double &ci) {
     return cephes::sici(x, si, ci);
 }
 
-XSF_HOST_DEVICE inline int shichi(double x, double *shi, double *chi) {
+XSF_HOST_DEVICE inline int shichi(double x, double &shi, double &chi) {
     return cephes::shichi(x, shi, chi);
 }
 
-XSF_HOST_DEVICE inline int sici(float x, float *si_f, float *ci_f) {
+XSF_HOST_DEVICE inline int sici(float x, float &si_f, float &ci_f) {
     double si;
     double ci;
-    int res = cephes::sici(x, &si, &ci);
-    *si_f = si;
-    *ci_f = ci;
+    int res = cephes::sici(x, si, ci);
+    si_f = si;
+    ci_f = ci;
     return res;
 }
 
-XSF_HOST_DEVICE inline int shichi(float x, float *shi_f, float *chi_f) {
+XSF_HOST_DEVICE inline int shichi(float x, float &shi_f, float &chi_f) {
     double shi;
     double chi;
-    int res = cephes::shichi(x, &shi, &chi);
-    *shi_f = shi;
-    *chi_f = chi;
+    int res = cephes::shichi(x, shi, chi);
+    shi_f = shi;
+    chi_f = chi;
     return res;
 }
 }
diff --git a/include/xsf/specfun.h b/include/xsf/specfun.h
index 305416f3e9..278625a709 100644
--- a/include/xsf/specfun.h
+++ b/include/xsf/specfun.h
@@ -52,7 +52,7 @@ inline std::complex<double> chyp2f1(double a, double b, double c, std::complex<d
     return outz;
 }
 
-inline std::complex<double> chyp1f1(double a, double b, std::complex<double> z) {
+inline std::complex<double> hyp1f1(double a, double b, std::complex<double> z) {
     std::complex<double> outz = specfun::cchg(a, b, z);
     if (outz.real() == 1e300) {
         set_error("chyp1f1", SF_ERROR_OVERFLOW, NULL);
diff --git a/include/xsf/stats.h b/include/xsf/stats.h
index 267bf3c6d0..c639797743 100644
--- a/include/xsf/stats.h
+++ b/include/xsf/stats.h
@@ -127,7 +127,7 @@ inline std::complex<double> log_ndtr(std::complex<double> z) {
     std::complex<double> val1 = std::complex<double>(mRe_z2, im);
 
     std::complex<double> val2 = log(xsf::wofz(complex<double>(-y, x)));
-    std::complex<double> result = val1 + val2 - NPY_LOGE2;
+    std::complex<double> result = val1 + val2 - M_LN2;
 
     /* Again, select the principal branch: log(z) = log|z| + i arg(z), thus
      * the imaginary part of the result should belong to [-pi, pi].
@@ -155,6 +155,8 @@ inline double nbdtri(int k, int n, double p) { return cephes::nbdtri(k, n, p); }
 
 inline double ndtri(double x) { return cephes::ndtri(x); }
 
+inline float ndtri(float x) { return static_cast<float>(cephes::ndtri(x)); }
+
 inline double owens_t(double h, double a) { return cephes::owens_t(h, a); }
 
 inline double pdtr(double k, double m) { return cephes::pdtr(k, m); }
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
new file mode 100644
index 0000000000..2397699a4d
--- /dev/null
+++ b/tests/CMakeLists.txt
@@ -0,0 +1,23 @@
+if(NOT DEFINED ENV{XSREF_TABLES_PATH})
+    message(FATAL_ERROR "XSREF_TABLES_PATH environment variable is not set. Specify the path to xsref/tables in your environment.")
+endif()
+
+set(XSREF_TABLES_PATH $ENV{XSREF_TABLES_PATH})
+
+find_package(Catch2 3 REQUIRED)
+find_package(Arrow REQUIRED)
+find_package(Parquet REQUIRED)
+
+add_library(xsf INTERFACE)
+target_include_directories(xsf INTERFACE ${CMAKE_SOURCE_DIR}/include)
+
+file(GLOB TEST_SOURCES "*/test_*.cpp")
+foreach(test_file ${TEST_SOURCES})
+  get_filename_component(test_name ${test_file} NAME_WE)
+  add_executable(${test_name}.test ${test_file})
+  target_link_libraries(${test_name}.test PRIVATE Catch2::Catch2WithMain Arrow::arrow_shared Parquet::parquet_shared xsf)
+  target_compile_definitions(${test_name}.test PRIVATE XSREF_TABLES_PATH="${XSREF_TABLES_PATH}")
+  include(CTest)
+  include(Catch)
+  catch_discover_tests(${test_name}.test)
+endforeach()
diff --git a/tests/scipy_special_tests/test_airy.cpp b/tests/scipy_special_tests/test_airy.cpp
new file mode 100644
index 0000000000..68254fe15a
--- /dev/null
+++ b/tests/scipy_special_tests/test_airy.cpp
@@ -0,0 +1,51 @@
+#include "../testing_utils.h"
+
+#include <xsf/airy.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "airy"};
+
+TEST_CASE("airy D->DDDD scipy_special_tests", "[airy][D->DDDD][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::complex<double>,
+            std::tuple<std::complex<double>, std::complex<double>, std::complex<double>, std::complex<double>, bool>,
+            std::tuple<double, double, double, double>>(
+            tables_path / "In_cd-cd_cd_cd_cd.parquet", tables_path / "Out_cd-cd_cd_cd_cd.parquet",
+            tables_path / ("Err_cd-cd_cd_cd_cd_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto z = input;
+    auto [desired0, desired1, desired2, desired3, fallback] = output;
+
+    std::complex<double> out0;
+    std::complex<double> out1;
+    std::complex<double> out2;
+    std::complex<double> out3;
+
+    xsf::airy(z, out0, out1, out2, out3);
+    auto [tol0, tol1, tol2, tol3] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(z, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(z, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+
+    auto error2 = xsf::extended_relative_error(out2, desired2);
+    tol2 = adjust_tolerance(tol2);
+    CAPTURE(z, out2, desired2, error2, tol2, fallback);
+    REQUIRE(error2 <= tol2);
+
+    auto error3 = xsf::extended_relative_error(out3, desired3);
+    tol3 = adjust_tolerance(tol3);
+    CAPTURE(z, out3, desired3, error3, tol3, fallback);
+    REQUIRE(error3 <= tol3);
+}
diff --git a/tests/scipy_special_tests/test_airye.cpp b/tests/scipy_special_tests/test_airye.cpp
new file mode 100644
index 0000000000..6b45832fc9
--- /dev/null
+++ b/tests/scipy_special_tests/test_airye.cpp
@@ -0,0 +1,51 @@
+#include "../testing_utils.h"
+
+#include <xsf/airy.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "airye"};
+
+TEST_CASE("airye D->DDDD scipy_special_tests", "[airye][D->DDDD][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::complex<double>,
+            std::tuple<std::complex<double>, std::complex<double>, std::complex<double>, std::complex<double>, bool>,
+            std::tuple<double, double, double, double>>(
+            tables_path / "In_cd-cd_cd_cd_cd.parquet", tables_path / "Out_cd-cd_cd_cd_cd.parquet",
+            tables_path / ("Err_cd-cd_cd_cd_cd_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto z = input;
+    auto [desired0, desired1, desired2, desired3, fallback] = output;
+
+    std::complex<double> out0;
+    std::complex<double> out1;
+    std::complex<double> out2;
+    std::complex<double> out3;
+
+    xsf::airye(z, out0, out1, out2, out3);
+    auto [tol0, tol1, tol2, tol3] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(z, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(z, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+
+    auto error2 = xsf::extended_relative_error(out2, desired2);
+    tol2 = adjust_tolerance(tol2);
+    CAPTURE(z, out2, desired2, error2, tol2, fallback);
+    REQUIRE(error2 <= tol2);
+
+    auto error3 = xsf::extended_relative_error(out3, desired3);
+    tol3 = adjust_tolerance(tol3);
+    CAPTURE(z, out3, desired3, error3, tol3, fallback);
+    REQUIRE(error3 <= tol3);
+}
diff --git a/tests/scipy_special_tests/test_bdtr.cpp b/tests/scipy_special_tests/test_bdtr.cpp
new file mode 100644
index 0000000000..24ecfd0ab2
--- /dev/null
+++ b/tests/scipy_special_tests/test_bdtr.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "bdtr"};
+
+TEST_CASE("bdtr dpd->d scipy_special_tests", "[bdtr][dpd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, std::ptrdiff_t, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_p_d-d.parquet", tables_path / "Out_d_p_d-d.parquet",
+            tables_path / ("Err_d_p_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [k, n, p] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::bdtr(k, n, p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(k, n, p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_bdtrc.cpp b/tests/scipy_special_tests/test_bdtrc.cpp
new file mode 100644
index 0000000000..81a23d798c
--- /dev/null
+++ b/tests/scipy_special_tests/test_bdtrc.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "bdtrc"};
+
+TEST_CASE("bdtrc dpd->d scipy_special_tests", "[bdtrc][dpd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, std::ptrdiff_t, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_p_d-d.parquet", tables_path / "Out_d_p_d-d.parquet",
+            tables_path / ("Err_d_p_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [k, n, p] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::bdtrc(k, n, p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(k, n, p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_bdtri.cpp b/tests/scipy_special_tests/test_bdtri.cpp
new file mode 100644
index 0000000000..a2078ae143
--- /dev/null
+++ b/tests/scipy_special_tests/test_bdtri.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "bdtri"};
+
+TEST_CASE("bdtri dpd->d scipy_special_tests", "[bdtri][dpd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, std::ptrdiff_t, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_p_d-d.parquet", tables_path / "Out_d_p_d-d.parquet",
+            tables_path / ("Err_d_p_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [k, n, y] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::bdtri(k, n, y);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(k, n, y, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_bei.cpp b/tests/scipy_special_tests/test_bei.cpp
new file mode 100644
index 0000000000..09de66e3f4
--- /dev/null
+++ b/tests/scipy_special_tests/test_bei.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/kelvin.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "bei"};
+
+TEST_CASE("bei d->d scipy_special_tests", "[bei][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::bei(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_beip.cpp b/tests/scipy_special_tests/test_beip.cpp
new file mode 100644
index 0000000000..36b68845f6
--- /dev/null
+++ b/tests/scipy_special_tests/test_beip.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/kelvin.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "beip"};
+
+TEST_CASE("beip d->d scipy_special_tests", "[beip][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::beip(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_ber.cpp b/tests/scipy_special_tests/test_ber.cpp
new file mode 100644
index 0000000000..154c16fc97
--- /dev/null
+++ b/tests/scipy_special_tests/test_ber.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/kelvin.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ber"};
+
+TEST_CASE("ber d->d scipy_special_tests", "[ber][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::ber(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_berp.cpp b/tests/scipy_special_tests/test_berp.cpp
new file mode 100644
index 0000000000..72c428922f
--- /dev/null
+++ b/tests/scipy_special_tests/test_berp.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/kelvin.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "berp"};
+
+TEST_CASE("berp d->d scipy_special_tests", "[berp][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::berp(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_besselpoly.cpp b/tests/scipy_special_tests/test_besselpoly.cpp
new file mode 100644
index 0000000000..dd6ae93ca5
--- /dev/null
+++ b/tests/scipy_special_tests/test_besselpoly.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "besselpoly"};
+
+TEST_CASE("besselpoly ddd->d scipy_special_tests", "[besselpoly][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [a, lmb, nu] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::besselpoly(a, lmb, nu);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, lmb, nu, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_beta.cpp b/tests/scipy_special_tests/test_beta.cpp
new file mode 100644
index 0000000000..a40dbb7b37
--- /dev/null
+++ b/tests/scipy_special_tests/test_beta.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/beta.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "beta"};
+
+TEST_CASE("beta dd->d scipy_special_tests", "[beta][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [a, b] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::beta(a, b);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, b, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_betaln.cpp b/tests/scipy_special_tests/test_betaln.cpp
new file mode 100644
index 0000000000..173f2fcd34
--- /dev/null
+++ b/tests/scipy_special_tests/test_betaln.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/beta.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "betaln"};
+
+TEST_CASE("betaln dd->d scipy_special_tests", "[betaln][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [a, b] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::betaln(a, b);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, b, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_binom.cpp b/tests/scipy_special_tests/test_binom.cpp
new file mode 100644
index 0000000000..6941e830c8
--- /dev/null
+++ b/tests/scipy_special_tests/test_binom.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/binom.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "binom"};
+
+TEST_CASE("binom dd->d scipy_special_tests", "[binom][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [n, k] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::binom(n, k);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(n, k, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cbrt.cpp b/tests/scipy_special_tests/test_cbrt.cpp
new file mode 100644
index 0000000000..a091ef0728
--- /dev/null
+++ b/tests/scipy_special_tests/test_cbrt.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/cephes/cbrt.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cbrt"};
+
+TEST_CASE("cephes::cbrt d->d scipy_special_tests", "[cephes::cbrt][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cephes::cbrt(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cem.cpp b/tests/scipy_special_tests/test_cem.cpp
new file mode 100644
index 0000000000..80fd4e01a4
--- /dev/null
+++ b/tests/scipy_special_tests/test_cem.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/mathieu.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cem"};
+
+TEST_CASE("cem ddd->dd scipy_special_tests", "[cem][ddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 std::tuple<double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d-d_d.parquet", tables_path / "Out_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [m, q, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::cem(m, q, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, q, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, q, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_cem_cva.cpp b/tests/scipy_special_tests/test_cem_cva.cpp
new file mode 100644
index 0000000000..9efe785752
--- /dev/null
+++ b/tests/scipy_special_tests/test_cem_cva.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/mathieu.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cem_cva"};
+
+TEST_CASE("cem_cva dd->d scipy_special_tests", "[cem_cva][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [m, q] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cem_cva(m, q);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(m, q, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_chdtr.cpp b/tests/scipy_special_tests/test_chdtr.cpp
new file mode 100644
index 0000000000..b3771abd3c
--- /dev/null
+++ b/tests/scipy_special_tests/test_chdtr.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "chdtr"};
+
+TEST_CASE("chdtr dd->d scipy_special_tests", "[chdtr][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::chdtr(v, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_chdtrc.cpp b/tests/scipy_special_tests/test_chdtrc.cpp
new file mode 100644
index 0000000000..50aa2bffc3
--- /dev/null
+++ b/tests/scipy_special_tests/test_chdtrc.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "chdtrc"};
+
+TEST_CASE("chdtrc dd->d scipy_special_tests", "[chdtrc][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::chdtrc(v, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_chdtri.cpp b/tests/scipy_special_tests/test_chdtri.cpp
new file mode 100644
index 0000000000..7b8d2f119f
--- /dev/null
+++ b/tests/scipy_special_tests/test_chdtri.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "chdtri"};
+
+TEST_CASE("chdtri dd->d scipy_special_tests", "[chdtri][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, p] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::chdtri(v, p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cosdg.cpp b/tests/scipy_special_tests/test_cosdg.cpp
new file mode 100644
index 0000000000..1e806b2ed0
--- /dev/null
+++ b/tests/scipy_special_tests/test_cosdg.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/trig.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cosdg"};
+
+TEST_CASE("cosdg d->d scipy_special_tests", "[cosdg][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cosdg(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cosm1.cpp b/tests/scipy_special_tests/test_cosm1.cpp
new file mode 100644
index 0000000000..26f64a7a50
--- /dev/null
+++ b/tests/scipy_special_tests/test_cosm1.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/trig.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cosm1"};
+
+TEST_CASE("cosm1 d->d scipy_special_tests", "[cosm1][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cosm1(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cospi.cpp b/tests/scipy_special_tests/test_cospi.cpp
new file mode 100644
index 0000000000..dcee3392f3
--- /dev/null
+++ b/tests/scipy_special_tests/test_cospi.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/trig.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cospi"};
+
+TEST_CASE("cospi D->D scipy_special_tests", "[cospi][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cospi(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cotdg.cpp b/tests/scipy_special_tests/test_cotdg.cpp
new file mode 100644
index 0000000000..2e364d405a
--- /dev/null
+++ b/tests/scipy_special_tests/test_cotdg.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/trig.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cotdg"};
+
+TEST_CASE("cotdg d->d scipy_special_tests", "[cotdg][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cotdg(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_i.cpp b/tests/scipy_special_tests/test_cyl_bessel_i.cpp
new file mode 100644
index 0000000000..de52287de6
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_i.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_i"};
+
+TEST_CASE("cyl_bessel_i dd->d scipy_special_tests", "[cyl_bessel_i][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_i(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_i0.cpp b/tests/scipy_special_tests/test_cyl_bessel_i0.cpp
new file mode 100644
index 0000000000..7884dcb5cc
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_i0.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_i0"};
+
+TEST_CASE("cyl_bessel_i0 f->f scipy_special_tests", "[cyl_bessel_i0][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_i0(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_i0e.cpp b/tests/scipy_special_tests/test_cyl_bessel_i0e.cpp
new file mode 100644
index 0000000000..47b4e7b4da
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_i0e.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_i0e"};
+
+TEST_CASE("cyl_bessel_i0e f->f scipy_special_tests", "[cyl_bessel_i0e][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_i0e(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_i1.cpp b/tests/scipy_special_tests/test_cyl_bessel_i1.cpp
new file mode 100644
index 0000000000..1524dacd0a
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_i1.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_i1"};
+
+TEST_CASE("cyl_bessel_i1 f->f scipy_special_tests", "[cyl_bessel_i1][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_i1(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_i1e.cpp b/tests/scipy_special_tests/test_cyl_bessel_i1e.cpp
new file mode 100644
index 0000000000..74e8b49db0
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_i1e.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_i1e"};
+
+TEST_CASE("cyl_bessel_i1e f->f scipy_special_tests", "[cyl_bessel_i1e][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_i1e(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_ie.cpp b/tests/scipy_special_tests/test_cyl_bessel_ie.cpp
new file mode 100644
index 0000000000..5c364c49af
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_ie.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_ie"};
+
+TEST_CASE("cyl_bessel_ie dd->d scipy_special_tests", "[cyl_bessel_ie][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_ie(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_j.cpp b/tests/scipy_special_tests/test_cyl_bessel_j.cpp
new file mode 100644
index 0000000000..6c391978ab
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_j.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_j"};
+
+TEST_CASE("cyl_bessel_j dd->d scipy_special_tests", "[cyl_bessel_j][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_j(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_j0.cpp b/tests/scipy_special_tests/test_cyl_bessel_j0.cpp
new file mode 100644
index 0000000000..1c704c72b7
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_j0.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_j0"};
+
+TEST_CASE("cyl_bessel_j0 d->d scipy_special_tests", "[cyl_bessel_j0][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_j0(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_j1.cpp b/tests/scipy_special_tests/test_cyl_bessel_j1.cpp
new file mode 100644
index 0000000000..43a11b2e08
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_j1.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_j1"};
+
+TEST_CASE("cyl_bessel_j1 d->d scipy_special_tests", "[cyl_bessel_j1][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_j1(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_je.cpp b/tests/scipy_special_tests/test_cyl_bessel_je.cpp
new file mode 100644
index 0000000000..9753338ed6
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_je.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_je"};
+
+TEST_CASE("cyl_bessel_je dd->d scipy_special_tests", "[cyl_bessel_je][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_je(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_k.cpp b/tests/scipy_special_tests/test_cyl_bessel_k.cpp
new file mode 100644
index 0000000000..fbe1a474ad
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_k.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_k"};
+
+TEST_CASE("cyl_bessel_k dd->d scipy_special_tests", "[cyl_bessel_k][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_k(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_k0.cpp b/tests/scipy_special_tests/test_cyl_bessel_k0.cpp
new file mode 100644
index 0000000000..5d6a9ca8f2
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_k0.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_k0"};
+
+TEST_CASE("cyl_bessel_k0 d->d scipy_special_tests", "[cyl_bessel_k0][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_k0(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_k0e.cpp b/tests/scipy_special_tests/test_cyl_bessel_k0e.cpp
new file mode 100644
index 0000000000..878cfac673
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_k0e.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_k0e"};
+
+TEST_CASE("cyl_bessel_k0e d->d scipy_special_tests", "[cyl_bessel_k0e][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_k0e(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_k1.cpp b/tests/scipy_special_tests/test_cyl_bessel_k1.cpp
new file mode 100644
index 0000000000..c1db7516c0
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_k1.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_k1"};
+
+TEST_CASE("cyl_bessel_k1 d->d scipy_special_tests", "[cyl_bessel_k1][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_k1(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_k1e.cpp b/tests/scipy_special_tests/test_cyl_bessel_k1e.cpp
new file mode 100644
index 0000000000..61f8ad54f8
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_k1e.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_k1e"};
+
+TEST_CASE("cyl_bessel_k1e d->d scipy_special_tests", "[cyl_bessel_k1e][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_k1e(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_ke.cpp b/tests/scipy_special_tests/test_cyl_bessel_ke.cpp
new file mode 100644
index 0000000000..af40ef8196
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_ke.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_ke"};
+
+TEST_CASE("cyl_bessel_ke dd->d scipy_special_tests", "[cyl_bessel_ke][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_ke(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_y.cpp b/tests/scipy_special_tests/test_cyl_bessel_y.cpp
new file mode 100644
index 0000000000..afb8f9330c
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_y.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_y"};
+
+TEST_CASE("cyl_bessel_y dd->d scipy_special_tests", "[cyl_bessel_y][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_y(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_y0.cpp b/tests/scipy_special_tests/test_cyl_bessel_y0.cpp
new file mode 100644
index 0000000000..0161c5cc7f
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_y0.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_y0"};
+
+TEST_CASE("cyl_bessel_y0 d->d scipy_special_tests", "[cyl_bessel_y0][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_y0(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_y1.cpp b/tests/scipy_special_tests/test_cyl_bessel_y1.cpp
new file mode 100644
index 0000000000..64399bcc55
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_y1.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_y1"};
+
+TEST_CASE("cyl_bessel_y1 d->d scipy_special_tests", "[cyl_bessel_y1][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_y1(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_bessel_ye.cpp b/tests/scipy_special_tests/test_cyl_bessel_ye.cpp
new file mode 100644
index 0000000000..36d8c2eecb
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_bessel_ye.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_bessel_ye"};
+
+TEST_CASE("cyl_bessel_ye dd->d scipy_special_tests", "[cyl_bessel_ye][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_bessel_ye(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_hankel_1.cpp b/tests/scipy_special_tests/test_cyl_hankel_1.cpp
new file mode 100644
index 0000000000..eaf926bea1
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_hankel_1.cpp
@@ -0,0 +1,25 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_hankel_1"};
+
+TEST_CASE("cyl_hankel_1 dD->D scipy_special_tests", "[cyl_hankel_1][dD->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<std::tuple<double, std::complex<double>>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_d_cd-cd.parquet", tables_path / "Out_d_cd-cd.parquet",
+            tables_path / ("Err_d_cd-cd_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_hankel_1(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_hankel_1e.cpp b/tests/scipy_special_tests/test_cyl_hankel_1e.cpp
new file mode 100644
index 0000000000..34a2e9f607
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_hankel_1e.cpp
@@ -0,0 +1,25 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_hankel_1e"};
+
+TEST_CASE("cyl_hankel_1e dD->D scipy_special_tests", "[cyl_hankel_1e][dD->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<std::tuple<double, std::complex<double>>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_d_cd-cd.parquet", tables_path / "Out_d_cd-cd.parquet",
+            tables_path / ("Err_d_cd-cd_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_hankel_1e(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_hankel_2.cpp b/tests/scipy_special_tests/test_cyl_hankel_2.cpp
new file mode 100644
index 0000000000..649b836786
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_hankel_2.cpp
@@ -0,0 +1,25 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_hankel_2"};
+
+TEST_CASE("cyl_hankel_2 dD->D scipy_special_tests", "[cyl_hankel_2][dD->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<std::tuple<double, std::complex<double>>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_d_cd-cd.parquet", tables_path / "Out_d_cd-cd.parquet",
+            tables_path / ("Err_d_cd-cd_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_hankel_2(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_cyl_hankel_2e.cpp b/tests/scipy_special_tests/test_cyl_hankel_2e.cpp
new file mode 100644
index 0000000000..49d91b14f0
--- /dev/null
+++ b/tests/scipy_special_tests/test_cyl_hankel_2e.cpp
@@ -0,0 +1,25 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "cyl_hankel_2e"};
+
+TEST_CASE("cyl_hankel_2e dD->D scipy_special_tests", "[cyl_hankel_2e][dD->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<std::tuple<double, std::complex<double>>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_d_cd-cd.parquet", tables_path / "Out_d_cd-cd.parquet",
+            tables_path / ("Err_d_cd-cd_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [v, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cyl_hankel_2e(v, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_dawsn.cpp b/tests/scipy_special_tests/test_dawsn.cpp
new file mode 100644
index 0000000000..8a713b73d1
--- /dev/null
+++ b/tests/scipy_special_tests/test_dawsn.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/erf.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "dawsn"};
+
+TEST_CASE("dawsn D->D scipy_special_tests", "[dawsn][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::dawsn(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_digamma.cpp b/tests/scipy_special_tests/test_digamma.cpp
new file mode 100644
index 0000000000..cffcb62c01
--- /dev/null
+++ b/tests/scipy_special_tests/test_digamma.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/digamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "digamma"};
+
+TEST_CASE("digamma D->D scipy_special_tests", "[digamma][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::digamma(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_ellipe.cpp b/tests/scipy_special_tests/test_ellipe.cpp
new file mode 100644
index 0000000000..39577aa14b
--- /dev/null
+++ b/tests/scipy_special_tests/test_ellipe.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/ellip.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ellipe"};
+
+TEST_CASE("ellipe d->d scipy_special_tests", "[ellipe][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto m = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::ellipe(m);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(m, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_ellipeinc.cpp b/tests/scipy_special_tests/test_ellipeinc.cpp
new file mode 100644
index 0000000000..301104c975
--- /dev/null
+++ b/tests/scipy_special_tests/test_ellipeinc.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/ellip.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ellipeinc"};
+
+TEST_CASE("ellipeinc dd->d scipy_special_tests", "[ellipeinc][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [phi, m] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::ellipeinc(phi, m);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(phi, m, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_ellipj.cpp b/tests/scipy_special_tests/test_ellipj.cpp
new file mode 100644
index 0000000000..cb524521c6
--- /dev/null
+++ b/tests/scipy_special_tests/test_ellipj.cpp
@@ -0,0 +1,48 @@
+#include "../testing_utils.h"
+
+#include <xsf/ellip.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ellipj"};
+
+TEST_CASE("ellipj dd->dddd scipy_special_tests", "[ellipj][dd->dddd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<
+                                         std::tuple<double, double>, std::tuple<double, double, double, double, bool>,
+                                         std::tuple<double, double, double, double>>(
+        tables_path / "In_d_d-d_d_d_d.parquet", tables_path / "Out_d_d-d_d_d_d.parquet",
+        tables_path / ("Err_d_d-d_d_d_d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [u, m] = input;
+    auto [desired0, desired1, desired2, desired3, fallback] = output;
+
+    double out0;
+    double out1;
+    double out2;
+    double out3;
+
+    xsf::ellipj(u, m, out0, out1, out2, out3);
+    auto [tol0, tol1, tol2, tol3] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(u, m, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(u, m, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+
+    auto error2 = xsf::extended_relative_error(out2, desired2);
+    tol2 = adjust_tolerance(tol2);
+    CAPTURE(u, m, out2, desired2, error2, tol2, fallback);
+    REQUIRE(error2 <= tol2);
+
+    auto error3 = xsf::extended_relative_error(out3, desired3);
+    tol3 = adjust_tolerance(tol3);
+    CAPTURE(u, m, out3, desired3, error3, tol3, fallback);
+    REQUIRE(error3 <= tol3);
+}
diff --git a/tests/scipy_special_tests/test_ellipk.cpp b/tests/scipy_special_tests/test_ellipk.cpp
new file mode 100644
index 0000000000..fe6c4e9705
--- /dev/null
+++ b/tests/scipy_special_tests/test_ellipk.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/ellip.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ellipk"};
+
+TEST_CASE("ellipk d->d scipy_special_tests", "[ellipk][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto m = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::ellipk(m);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(m, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_ellipkinc.cpp b/tests/scipy_special_tests/test_ellipkinc.cpp
new file mode 100644
index 0000000000..a2c436d07d
--- /dev/null
+++ b/tests/scipy_special_tests/test_ellipkinc.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/ellip.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ellipkinc"};
+
+TEST_CASE("ellipkinc dd->d scipy_special_tests", "[ellipkinc][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [phi, m] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::ellipkinc(phi, m);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(phi, m, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_ellipkm1.cpp b/tests/scipy_special_tests/test_ellipkm1.cpp
new file mode 100644
index 0000000000..1f52934888
--- /dev/null
+++ b/tests/scipy_special_tests/test_ellipkm1.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/ellip.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ellipkm1"};
+
+TEST_CASE("ellipkm1 d->d scipy_special_tests", "[ellipkm1][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto p = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::ellipkm1(p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_erf.cpp b/tests/scipy_special_tests/test_erf.cpp
new file mode 100644
index 0000000000..10384012db
--- /dev/null
+++ b/tests/scipy_special_tests/test_erf.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/erf.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "erf"};
+
+TEST_CASE("erf f->f scipy_special_tests", "[erf][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::erf(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_erfc.cpp b/tests/scipy_special_tests/test_erfc.cpp
new file mode 100644
index 0000000000..10d108a4ad
--- /dev/null
+++ b/tests/scipy_special_tests/test_erfc.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/erf.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "erfc"};
+
+TEST_CASE("erfc f->f scipy_special_tests", "[erfc][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::erfc(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_erfcinv.cpp b/tests/scipy_special_tests/test_erfcinv.cpp
new file mode 100644
index 0000000000..c74e60a22c
--- /dev/null
+++ b/tests/scipy_special_tests/test_erfcinv.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/cephes/erfinv.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "erfcinv"};
+
+TEST_CASE("cephes::erfcinv d->d scipy_special_tests", "[cephes::erfcinv][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cephes::erfcinv(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_erfcx.cpp b/tests/scipy_special_tests/test_erfcx.cpp
new file mode 100644
index 0000000000..68854fc648
--- /dev/null
+++ b/tests/scipy_special_tests/test_erfcx.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/erf.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "erfcx"};
+
+TEST_CASE("erfcx D->D scipy_special_tests", "[erfcx][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::erfcx(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_erfi.cpp b/tests/scipy_special_tests/test_erfi.cpp
new file mode 100644
index 0000000000..076ef8d537
--- /dev/null
+++ b/tests/scipy_special_tests/test_erfi.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/erf.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "erfi"};
+
+TEST_CASE("erfi D->D scipy_special_tests", "[erfi][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::erfi(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_exp1.cpp b/tests/scipy_special_tests/test_exp1.cpp
new file mode 100644
index 0000000000..a5cf9cd0a6
--- /dev/null
+++ b/tests/scipy_special_tests/test_exp1.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/expint.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "exp1"};
+
+TEST_CASE("exp1 D->D scipy_special_tests", "[exp1][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::exp1(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_exp10.cpp b/tests/scipy_special_tests/test_exp10.cpp
new file mode 100644
index 0000000000..04b1f84e44
--- /dev/null
+++ b/tests/scipy_special_tests/test_exp10.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/exp.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "exp10"};
+
+TEST_CASE("exp10 d->d scipy_special_tests", "[exp10][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::exp10(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_exp2.cpp b/tests/scipy_special_tests/test_exp2.cpp
new file mode 100644
index 0000000000..38847b554d
--- /dev/null
+++ b/tests/scipy_special_tests/test_exp2.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/exp.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "exp2"};
+
+TEST_CASE("exp2 d->d scipy_special_tests", "[exp2][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::exp2(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_expi.cpp b/tests/scipy_special_tests/test_expi.cpp
new file mode 100644
index 0000000000..4251d20062
--- /dev/null
+++ b/tests/scipy_special_tests/test_expi.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/expint.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "expi"};
+
+TEST_CASE("expi D->D scipy_special_tests", "[expi][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::expi(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_expit.cpp b/tests/scipy_special_tests/test_expit.cpp
new file mode 100644
index 0000000000..a1a86ca608
--- /dev/null
+++ b/tests/scipy_special_tests/test_expit.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/log_exp.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "expit"};
+
+TEST_CASE("expit f->f scipy_special_tests", "[expit][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::expit(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_expm1.cpp b/tests/scipy_special_tests/test_expm1.cpp
new file mode 100644
index 0000000000..50ab2e5007
--- /dev/null
+++ b/tests/scipy_special_tests/test_expm1.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/exp.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "expm1"};
+
+TEST_CASE("expm1 D->D scipy_special_tests", "[expm1][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::expm1(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_expn.cpp b/tests/scipy_special_tests/test_expn.cpp
new file mode 100644
index 0000000000..840c970a8d
--- /dev/null
+++ b/tests/scipy_special_tests/test_expn.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/cephes/expn.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "expn"};
+
+TEST_CASE("cephes::expn dd->d scipy_special_tests", "[cephes::expn][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [n, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cephes::expn(n, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(n, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_exprel.cpp b/tests/scipy_special_tests/test_exprel.cpp
new file mode 100644
index 0000000000..733f339f85
--- /dev/null
+++ b/tests/scipy_special_tests/test_exprel.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/log_exp.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "exprel"};
+
+TEST_CASE("exprel d->d scipy_special_tests", "[exprel][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::exprel(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_fdtr.cpp b/tests/scipy_special_tests/test_fdtr.cpp
new file mode 100644
index 0000000000..c02afd357a
--- /dev/null
+++ b/tests/scipy_special_tests/test_fdtr.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "fdtr"};
+
+TEST_CASE("fdtr ddd->d scipy_special_tests", "[fdtr][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [dfn, dfd, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::fdtr(dfn, dfd, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(dfn, dfd, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_fdtrc.cpp b/tests/scipy_special_tests/test_fdtrc.cpp
new file mode 100644
index 0000000000..5d77d48107
--- /dev/null
+++ b/tests/scipy_special_tests/test_fdtrc.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "fdtrc"};
+
+TEST_CASE("fdtrc ddd->d scipy_special_tests", "[fdtrc][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [dfn, dfd, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::fdtrc(dfn, dfd, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(dfn, dfd, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_fdtri.cpp b/tests/scipy_special_tests/test_fdtri.cpp
new file mode 100644
index 0000000000..64a5883e58
--- /dev/null
+++ b/tests/scipy_special_tests/test_fdtri.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "fdtri"};
+
+TEST_CASE("fdtri ddd->d scipy_special_tests", "[fdtri][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [dfn, dfd, p] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::fdtri(dfn, dfd, p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(dfn, dfd, p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_fresnel.cpp b/tests/scipy_special_tests/test_fresnel.cpp
new file mode 100644
index 0000000000..80a2e60d7b
--- /dev/null
+++ b/tests/scipy_special_tests/test_fresnel.cpp
@@ -0,0 +1,35 @@
+#include "../testing_utils.h"
+
+#include <xsf/fresnel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "fresnel"};
+
+TEST_CASE("fresnel d->dd scipy_special_tests", "[fresnel][d->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<double, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d-d_d.parquet", tables_path / "Out_d-d_d.parquet",
+            tables_path / ("Err_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::fresnel(x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_gamma.cpp b/tests/scipy_special_tests/test_gamma.cpp
new file mode 100644
index 0000000000..8c71ba196d
--- /dev/null
+++ b/tests/scipy_special_tests/test_gamma.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/gamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gamma"};
+
+TEST_CASE("gamma D->D scipy_special_tests", "[gamma][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gamma(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_gammainc.cpp b/tests/scipy_special_tests/test_gammainc.cpp
new file mode 100644
index 0000000000..4164340a9c
--- /dev/null
+++ b/tests/scipy_special_tests/test_gammainc.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/gamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gammainc"};
+
+TEST_CASE("gammainc dd->d scipy_special_tests", "[gammainc][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [a, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gammainc(a, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_gammaincc.cpp b/tests/scipy_special_tests/test_gammaincc.cpp
new file mode 100644
index 0000000000..b8a0f237b0
--- /dev/null
+++ b/tests/scipy_special_tests/test_gammaincc.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/gamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gammaincc"};
+
+TEST_CASE("gammaincc dd->d scipy_special_tests", "[gammaincc][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [a, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gammaincc(a, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_gammainccinv.cpp b/tests/scipy_special_tests/test_gammainccinv.cpp
new file mode 100644
index 0000000000..5eacf597ae
--- /dev/null
+++ b/tests/scipy_special_tests/test_gammainccinv.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/gamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gammainccinv"};
+
+TEST_CASE("gammainccinv dd->d scipy_special_tests", "[gammainccinv][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [a, y] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gammainccinv(a, y);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, y, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_gammaincinv.cpp b/tests/scipy_special_tests/test_gammaincinv.cpp
new file mode 100644
index 0000000000..ce06c9fc06
--- /dev/null
+++ b/tests/scipy_special_tests/test_gammaincinv.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/gamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gammaincinv"};
+
+TEST_CASE("gammaincinv dd->d scipy_special_tests", "[gammaincinv][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [a, y] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gammaincinv(a, y);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, y, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_gammaln.cpp b/tests/scipy_special_tests/test_gammaln.cpp
new file mode 100644
index 0000000000..e2234188d7
--- /dev/null
+++ b/tests/scipy_special_tests/test_gammaln.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/gamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gammaln"};
+
+TEST_CASE("gammaln f->f scipy_special_tests", "[gammaln][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gammaln(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_gammasgn.cpp b/tests/scipy_special_tests/test_gammasgn.cpp
new file mode 100644
index 0000000000..f3339507f8
--- /dev/null
+++ b/tests/scipy_special_tests/test_gammasgn.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/gamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gammasgn"};
+
+TEST_CASE("gammasgn d->d scipy_special_tests", "[gammasgn][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gammasgn(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_gdtr.cpp b/tests/scipy_special_tests/test_gdtr.cpp
new file mode 100644
index 0000000000..72e9cc1678
--- /dev/null
+++ b/tests/scipy_special_tests/test_gdtr.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gdtr"};
+
+TEST_CASE("gdtr ddd->d scipy_special_tests", "[gdtr][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [a, b, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gdtr(a, b, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, b, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_gdtrc.cpp b/tests/scipy_special_tests/test_gdtrc.cpp
new file mode 100644
index 0000000000..c7d0fae62b
--- /dev/null
+++ b/tests/scipy_special_tests/test_gdtrc.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gdtrc"};
+
+TEST_CASE("gdtrc ddd->d scipy_special_tests", "[gdtrc][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [a, b, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gdtrc(a, b, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, b, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_gdtrib.cpp b/tests/scipy_special_tests/test_gdtrib.cpp
new file mode 100644
index 0000000000..9d4b60d518
--- /dev/null
+++ b/tests/scipy_special_tests/test_gdtrib.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/cdflib.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "gdtrib"};
+
+TEST_CASE("gdtrib ddd->d scipy_special_tests", "[gdtrib][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [a, p, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::gdtrib(a, p, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, p, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_hyp1f1.cpp b/tests/scipy_special_tests/test_hyp1f1.cpp
new file mode 100644
index 0000000000..57374e7907
--- /dev/null
+++ b/tests/scipy_special_tests/test_hyp1f1.cpp
@@ -0,0 +1,25 @@
+#include "../testing_utils.h"
+
+#include <xsf/specfun.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "hyp1f1"};
+
+TEST_CASE("hyp1f1 ddD->D scipy_special_tests", "[hyp1f1][ddD->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 std::tuple<double, double, std::complex<double>>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_d_d_cd-cd.parquet", tables_path / "Out_d_d_cd-cd.parquet",
+            tables_path / ("Err_d_d_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [a, b, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::hyp1f1(a, b, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, b, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_hyp2f1.cpp b/tests/scipy_special_tests/test_hyp2f1.cpp
new file mode 100644
index 0000000000..deeff0c235
--- /dev/null
+++ b/tests/scipy_special_tests/test_hyp2f1.cpp
@@ -0,0 +1,26 @@
+#include "../testing_utils.h"
+
+#include <xsf/hyp2f1.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "hyp2f1"};
+
+TEST_CASE("hyp2f1 dddD->D scipy_special_tests", "[hyp2f1][dddD->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<double, double, double, std::complex<double>>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_d_d_d_cd-cd.parquet", tables_path / "Out_d_d_d_cd-cd.parquet",
+            tables_path / ("Err_d_d_d_cd-cd_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [a, b, c, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::hyp2f1(a, b, c, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, b, c, z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_it1i0k0.cpp b/tests/scipy_special_tests/test_it1i0k0.cpp
new file mode 100644
index 0000000000..b69f5bdfb4
--- /dev/null
+++ b/tests/scipy_special_tests/test_it1i0k0.cpp
@@ -0,0 +1,35 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "it1i0k0"};
+
+TEST_CASE("it1i0k0 d->dd scipy_special_tests", "[it1i0k0][d->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<double, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d-d_d.parquet", tables_path / "Out_d-d_d.parquet",
+            tables_path / ("Err_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::it1i0k0(x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_it1j0y0.cpp b/tests/scipy_special_tests/test_it1j0y0.cpp
new file mode 100644
index 0000000000..692690c2ca
--- /dev/null
+++ b/tests/scipy_special_tests/test_it1j0y0.cpp
@@ -0,0 +1,35 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "it1j0y0"};
+
+TEST_CASE("it1j0y0 d->dd scipy_special_tests", "[it1j0y0][d->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<double, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d-d_d.parquet", tables_path / "Out_d-d_d.parquet",
+            tables_path / ("Err_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::it1j0y0(x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_it2i0k0.cpp b/tests/scipy_special_tests/test_it2i0k0.cpp
new file mode 100644
index 0000000000..7374ebfe90
--- /dev/null
+++ b/tests/scipy_special_tests/test_it2i0k0.cpp
@@ -0,0 +1,35 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "it2i0k0"};
+
+TEST_CASE("it2i0k0 d->dd scipy_special_tests", "[it2i0k0][d->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<double, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d-d_d.parquet", tables_path / "Out_d-d_d.parquet",
+            tables_path / ("Err_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::it2i0k0(x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_it2j0y0.cpp b/tests/scipy_special_tests/test_it2j0y0.cpp
new file mode 100644
index 0000000000..9332e18daf
--- /dev/null
+++ b/tests/scipy_special_tests/test_it2j0y0.cpp
@@ -0,0 +1,35 @@
+#include "../testing_utils.h"
+
+#include <xsf/bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "it2j0y0"};
+
+TEST_CASE("it2j0y0 d->dd scipy_special_tests", "[it2j0y0][d->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<double, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d-d_d.parquet", tables_path / "Out_d-d_d.parquet",
+            tables_path / ("Err_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::it2j0y0(x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_it2struve0.cpp b/tests/scipy_special_tests/test_it2struve0.cpp
new file mode 100644
index 0000000000..0786d9ef2b
--- /dev/null
+++ b/tests/scipy_special_tests/test_it2struve0.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/struve.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "it2struve0"};
+
+TEST_CASE("it2struve0 d->d scipy_special_tests", "[it2struve0][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::it2struve0(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_itairy.cpp b/tests/scipy_special_tests/test_itairy.cpp
new file mode 100644
index 0000000000..290565947f
--- /dev/null
+++ b/tests/scipy_special_tests/test_itairy.cpp
@@ -0,0 +1,48 @@
+#include "../testing_utils.h"
+
+#include <xsf/airy.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "itairy"};
+
+TEST_CASE("itairy d->dddd scipy_special_tests", "[itairy][d->dddd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 double, std::tuple<double, double, double, double, bool>, std::tuple<double, double, double, double>>(
+            tables_path / "In_d-d_d_d_d.parquet", tables_path / "Out_d-d_d_d_d.parquet",
+            tables_path / ("Err_d-d_d_d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, desired2, desired3, fallback] = output;
+
+    double out0;
+    double out1;
+    double out2;
+    double out3;
+
+    xsf::itairy(x, out0, out1, out2, out3);
+    auto [tol0, tol1, tol2, tol3] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+
+    auto error2 = xsf::extended_relative_error(out2, desired2);
+    tol2 = adjust_tolerance(tol2);
+    CAPTURE(x, out2, desired2, error2, tol2, fallback);
+    REQUIRE(error2 <= tol2);
+
+    auto error3 = xsf::extended_relative_error(out3, desired3);
+    tol3 = adjust_tolerance(tol3);
+    CAPTURE(x, out3, desired3, error3, tol3, fallback);
+    REQUIRE(error3 <= tol3);
+}
diff --git a/tests/scipy_special_tests/test_itmodstruve0.cpp b/tests/scipy_special_tests/test_itmodstruve0.cpp
new file mode 100644
index 0000000000..e231b76cdf
--- /dev/null
+++ b/tests/scipy_special_tests/test_itmodstruve0.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/struve.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "itmodstruve0"};
+
+TEST_CASE("itmodstruve0 d->d scipy_special_tests", "[itmodstruve0][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::itmodstruve0(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_itstruve0.cpp b/tests/scipy_special_tests/test_itstruve0.cpp
new file mode 100644
index 0000000000..88c6b8c548
--- /dev/null
+++ b/tests/scipy_special_tests/test_itstruve0.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/struve.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "itstruve0"};
+
+TEST_CASE("itstruve0 d->d scipy_special_tests", "[itstruve0][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::itstruve0(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_iv_ratio.cpp b/tests/scipy_special_tests/test_iv_ratio.cpp
new file mode 100644
index 0000000000..721bafd79d
--- /dev/null
+++ b/tests/scipy_special_tests/test_iv_ratio.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/iv_ratio.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "iv_ratio"};
+
+TEST_CASE("iv_ratio dd->d scipy_special_tests", "[iv_ratio][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::iv_ratio(v, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_iv_ratio_c.cpp b/tests/scipy_special_tests/test_iv_ratio_c.cpp
new file mode 100644
index 0000000000..103d0661f1
--- /dev/null
+++ b/tests/scipy_special_tests/test_iv_ratio_c.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/iv_ratio.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "iv_ratio_c"};
+
+TEST_CASE("iv_ratio_c dd->d scipy_special_tests", "[iv_ratio_c][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::iv_ratio_c(v, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_kei.cpp b/tests/scipy_special_tests/test_kei.cpp
new file mode 100644
index 0000000000..b883d4a036
--- /dev/null
+++ b/tests/scipy_special_tests/test_kei.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/kelvin.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "kei"};
+
+TEST_CASE("kei d->d scipy_special_tests", "[kei][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::kei(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_keip.cpp b/tests/scipy_special_tests/test_keip.cpp
new file mode 100644
index 0000000000..6f12def2d5
--- /dev/null
+++ b/tests/scipy_special_tests/test_keip.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/kelvin.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "keip"};
+
+TEST_CASE("keip d->d scipy_special_tests", "[keip][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::keip(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_kelvin.cpp b/tests/scipy_special_tests/test_kelvin.cpp
new file mode 100644
index 0000000000..d1afa114cb
--- /dev/null
+++ b/tests/scipy_special_tests/test_kelvin.cpp
@@ -0,0 +1,51 @@
+#include "../testing_utils.h"
+
+#include <xsf/kelvin.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "kelvin"};
+
+TEST_CASE("kelvin d->DDDD scipy_special_tests", "[kelvin][d->DDDD][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            double,
+            std::tuple<std::complex<double>, std::complex<double>, std::complex<double>, std::complex<double>, bool>,
+            std::tuple<double, double, double, double>>(
+            tables_path / "In_d-cd_cd_cd_cd.parquet", tables_path / "Out_d-cd_cd_cd_cd.parquet",
+            tables_path / ("Err_d-cd_cd_cd_cd_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto x = input;
+    auto [desired0, desired1, desired2, desired3, fallback] = output;
+
+    std::complex<double> out0;
+    std::complex<double> out1;
+    std::complex<double> out2;
+    std::complex<double> out3;
+
+    xsf::kelvin(x, out0, out1, out2, out3);
+    auto [tol0, tol1, tol2, tol3] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+
+    auto error2 = xsf::extended_relative_error(out2, desired2);
+    tol2 = adjust_tolerance(tol2);
+    CAPTURE(x, out2, desired2, error2, tol2, fallback);
+    REQUIRE(error2 <= tol2);
+
+    auto error3 = xsf::extended_relative_error(out3, desired3);
+    tol3 = adjust_tolerance(tol3);
+    CAPTURE(x, out3, desired3, error3, tol3, fallback);
+    REQUIRE(error3 <= tol3);
+}
diff --git a/tests/scipy_special_tests/test_ker.cpp b/tests/scipy_special_tests/test_ker.cpp
new file mode 100644
index 0000000000..d59055c902
--- /dev/null
+++ b/tests/scipy_special_tests/test_ker.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/kelvin.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ker"};
+
+TEST_CASE("ker d->d scipy_special_tests", "[ker][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::ker(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_kerp.cpp b/tests/scipy_special_tests/test_kerp.cpp
new file mode 100644
index 0000000000..c36fcc0331
--- /dev/null
+++ b/tests/scipy_special_tests/test_kerp.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/kelvin.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "kerp"};
+
+TEST_CASE("kerp d->d scipy_special_tests", "[kerp][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::kerp(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_kolmogc.cpp b/tests/scipy_special_tests/test_kolmogc.cpp
new file mode 100644
index 0000000000..aa311ff394
--- /dev/null
+++ b/tests/scipy_special_tests/test_kolmogc.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "kolmogc"};
+
+TEST_CASE("kolmogc d->d scipy_special_tests", "[kolmogc][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::kolmogc(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_kolmogci.cpp b/tests/scipy_special_tests/test_kolmogci.cpp
new file mode 100644
index 0000000000..29718f80ef
--- /dev/null
+++ b/tests/scipy_special_tests/test_kolmogci.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "kolmogci"};
+
+TEST_CASE("kolmogci d->d scipy_special_tests", "[kolmogci][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::kolmogci(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_kolmogi.cpp b/tests/scipy_special_tests/test_kolmogi.cpp
new file mode 100644
index 0000000000..d171d39a5b
--- /dev/null
+++ b/tests/scipy_special_tests/test_kolmogi.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "kolmogi"};
+
+TEST_CASE("kolmogi d->d scipy_special_tests", "[kolmogi][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::kolmogi(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_kolmogorov.cpp b/tests/scipy_special_tests/test_kolmogorov.cpp
new file mode 100644
index 0000000000..65310cac6d
--- /dev/null
+++ b/tests/scipy_special_tests/test_kolmogorov.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "kolmogorov"};
+
+TEST_CASE("kolmogorov d->d scipy_special_tests", "[kolmogorov][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::kolmogorov(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_kolmogp.cpp b/tests/scipy_special_tests/test_kolmogp.cpp
new file mode 100644
index 0000000000..1d07a6916d
--- /dev/null
+++ b/tests/scipy_special_tests/test_kolmogp.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "kolmogp"};
+
+TEST_CASE("kolmogp d->d scipy_special_tests", "[kolmogp][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::kolmogp(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_lambertw.cpp b/tests/scipy_special_tests/test_lambertw.cpp
new file mode 100644
index 0000000000..51c5b4c0b4
--- /dev/null
+++ b/tests/scipy_special_tests/test_lambertw.cpp
@@ -0,0 +1,26 @@
+#include "../testing_utils.h"
+
+#include <xsf/lambertw.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "lambertw"};
+
+TEST_CASE("lambertw Dpd->D scipy_special_tests", "[lambertw][Dpd->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<std::complex<double>, std::ptrdiff_t, double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd_p_d-cd.parquet", tables_path / "Out_cd_p_d-cd.parquet",
+            tables_path / ("Err_cd_p_d-cd_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [z, k, eps] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::lambertw(z, k, eps);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, k, eps, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_lanczos_sum_expg_scaled.cpp b/tests/scipy_special_tests/test_lanczos_sum_expg_scaled.cpp
new file mode 100644
index 0000000000..0ce2f3715e
--- /dev/null
+++ b/tests/scipy_special_tests/test_lanczos_sum_expg_scaled.cpp
@@ -0,0 +1,26 @@
+#include "../testing_utils.h"
+
+#include <xsf/cephes/lanczos.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "lanczos_sum_expg_scaled"};
+
+TEST_CASE(
+    "cephes::lanczos_sum_expg_scaled d->d scipy_special_tests",
+    "[cephes::lanczos_sum_expg_scaled][d->d][scipy_special_tests]"
+) {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cephes::lanczos_sum_expg_scaled(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_lgam1p.cpp b/tests/scipy_special_tests/test_lgam1p.cpp
new file mode 100644
index 0000000000..1c6ef2084b
--- /dev/null
+++ b/tests/scipy_special_tests/test_lgam1p.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/cephes/unity.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "lgam1p"};
+
+TEST_CASE("cephes::lgam1p d->d scipy_special_tests", "[cephes::lgam1p][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cephes::lgam1p(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_log1p.cpp b/tests/scipy_special_tests/test_log1p.cpp
new file mode 100644
index 0000000000..71cd6e433a
--- /dev/null
+++ b/tests/scipy_special_tests/test_log1p.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/log.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "log1p"};
+
+TEST_CASE("log1p D->D scipy_special_tests", "[log1p][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::log1p(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_log1pmx.cpp b/tests/scipy_special_tests/test_log1pmx.cpp
new file mode 100644
index 0000000000..6ab870b4fb
--- /dev/null
+++ b/tests/scipy_special_tests/test_log1pmx.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/log.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "log1pmx"};
+
+TEST_CASE("log1pmx d->d scipy_special_tests", "[log1pmx][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::log1pmx(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_log_expit.cpp b/tests/scipy_special_tests/test_log_expit.cpp
new file mode 100644
index 0000000000..008399f3b5
--- /dev/null
+++ b/tests/scipy_special_tests/test_log_expit.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/log_exp.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "log_expit"};
+
+TEST_CASE("log_expit f->f scipy_special_tests", "[log_expit][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::log_expit(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_log_wright_bessel.cpp b/tests/scipy_special_tests/test_log_wright_bessel.cpp
new file mode 100644
index 0000000000..aaa9cda29b
--- /dev/null
+++ b/tests/scipy_special_tests/test_log_wright_bessel.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/wright_bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "log_wright_bessel"};
+
+TEST_CASE("log_wright_bessel ddd->d scipy_special_tests", "[log_wright_bessel][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [a, b, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::log_wright_bessel(a, b, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, b, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_loggamma.cpp b/tests/scipy_special_tests/test_loggamma.cpp
new file mode 100644
index 0000000000..e8ca03fcdd
--- /dev/null
+++ b/tests/scipy_special_tests/test_loggamma.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/gamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "loggamma"};
+
+TEST_CASE("loggamma D->D scipy_special_tests", "[loggamma][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::loggamma(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_logit.cpp b/tests/scipy_special_tests/test_logit.cpp
new file mode 100644
index 0000000000..c12982c546
--- /dev/null
+++ b/tests/scipy_special_tests/test_logit.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/log_exp.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "logit"};
+
+TEST_CASE("logit f->f scipy_special_tests", "[logit][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto p = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::logit(p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_mcm1.cpp b/tests/scipy_special_tests/test_mcm1.cpp
new file mode 100644
index 0000000000..213dc058d8
--- /dev/null
+++ b/tests/scipy_special_tests/test_mcm1.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/mathieu.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "mcm1"};
+
+TEST_CASE("mcm1 ddd->dd scipy_special_tests", "[mcm1][ddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 std::tuple<double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d-d_d.parquet", tables_path / "Out_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [m, q, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::mcm1(m, q, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, q, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, q, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_mcm2.cpp b/tests/scipy_special_tests/test_mcm2.cpp
new file mode 100644
index 0000000000..4ae20d5a38
--- /dev/null
+++ b/tests/scipy_special_tests/test_mcm2.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/mathieu.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "mcm2"};
+
+TEST_CASE("mcm2 ddd->dd scipy_special_tests", "[mcm2][ddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 std::tuple<double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d-d_d.parquet", tables_path / "Out_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [m, q, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::mcm2(m, q, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, q, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, q, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_modified_fresnel_minus.cpp b/tests/scipy_special_tests/test_modified_fresnel_minus.cpp
new file mode 100644
index 0000000000..dbf1a239dd
--- /dev/null
+++ b/tests/scipy_special_tests/test_modified_fresnel_minus.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/fresnel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "modified_fresnel_minus"};
+
+TEST_CASE("modified_fresnel_minus d->DD scipy_special_tests", "[modified_fresnel_minus][d->DD][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 double, std::tuple<std::complex<double>, std::complex<double>, bool>, std::tuple<double, double>>(
+            tables_path / "In_d-cd_cd.parquet", tables_path / "Out_d-cd_cd.parquet",
+            tables_path / ("Err_d-cd_cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, fallback] = output;
+
+    std::complex<double> out0;
+    std::complex<double> out1;
+
+    xsf::modified_fresnel_minus(x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_modified_fresnel_plus.cpp b/tests/scipy_special_tests/test_modified_fresnel_plus.cpp
new file mode 100644
index 0000000000..c1222053a8
--- /dev/null
+++ b/tests/scipy_special_tests/test_modified_fresnel_plus.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/fresnel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "modified_fresnel_plus"};
+
+TEST_CASE("modified_fresnel_plus d->DD scipy_special_tests", "[modified_fresnel_plus][d->DD][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 double, std::tuple<std::complex<double>, std::complex<double>, bool>, std::tuple<double, double>>(
+            tables_path / "In_d-cd_cd.parquet", tables_path / "Out_d-cd_cd.parquet",
+            tables_path / ("Err_d-cd_cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, fallback] = output;
+
+    std::complex<double> out0;
+    std::complex<double> out1;
+
+    xsf::modified_fresnel_plus(x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_msm1.cpp b/tests/scipy_special_tests/test_msm1.cpp
new file mode 100644
index 0000000000..2d47a2f267
--- /dev/null
+++ b/tests/scipy_special_tests/test_msm1.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/mathieu.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "msm1"};
+
+TEST_CASE("msm1 ddd->dd scipy_special_tests", "[msm1][ddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 std::tuple<double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d-d_d.parquet", tables_path / "Out_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [m, q, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::msm1(m, q, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, q, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, q, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_msm2.cpp b/tests/scipy_special_tests/test_msm2.cpp
new file mode 100644
index 0000000000..421be0e8a1
--- /dev/null
+++ b/tests/scipy_special_tests/test_msm2.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/mathieu.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "msm2"};
+
+TEST_CASE("msm2 ddd->dd scipy_special_tests", "[msm2][ddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 std::tuple<double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d-d_d.parquet", tables_path / "Out_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [m, q, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::msm2(m, q, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, q, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, q, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_nbdtr.cpp b/tests/scipy_special_tests/test_nbdtr.cpp
new file mode 100644
index 0000000000..afec118966
--- /dev/null
+++ b/tests/scipy_special_tests/test_nbdtr.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "nbdtr"};
+
+TEST_CASE("nbdtr ddd->d scipy_special_tests", "[nbdtr][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [k, n, p] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::nbdtr(k, n, p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(k, n, p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_nbdtrc.cpp b/tests/scipy_special_tests/test_nbdtrc.cpp
new file mode 100644
index 0000000000..2f99ee574b
--- /dev/null
+++ b/tests/scipy_special_tests/test_nbdtrc.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "nbdtrc"};
+
+TEST_CASE("nbdtrc ddd->d scipy_special_tests", "[nbdtrc][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [k, n, p] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::nbdtrc(k, n, p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(k, n, p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_ndtr.cpp b/tests/scipy_special_tests/test_ndtr.cpp
new file mode 100644
index 0000000000..ea776d8da9
--- /dev/null
+++ b/tests/scipy_special_tests/test_ndtr.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ndtr"};
+
+TEST_CASE("ndtr f->f scipy_special_tests", "[ndtr][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::ndtr(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_ndtri.cpp b/tests/scipy_special_tests/test_ndtri.cpp
new file mode 100644
index 0000000000..8fcde33c52
--- /dev/null
+++ b/tests/scipy_special_tests/test_ndtri.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "ndtri"};
+
+TEST_CASE("ndtri f->f scipy_special_tests", "[ndtri][f->f][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<float, std::tuple<float, bool>, float>(
+        tables_path / "In_f-f.parquet", tables_path / "Out_f-f.parquet",
+        tables_path / ("Err_f-f_" + get_platform_str() + ".parquet")
+    ));
+
+    auto y = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::ndtri(y);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(y, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_oblate_aswfa.cpp b/tests/scipy_special_tests/test_oblate_aswfa.cpp
new file mode 100644
index 0000000000..df8ebda93b
--- /dev/null
+++ b/tests/scipy_special_tests/test_oblate_aswfa.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "oblate_aswfa"};
+
+TEST_CASE("oblate_aswfa ddddd->dd scipy_special_tests", "[oblate_aswfa][ddddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<
+                                         std::tuple<double, double, double, double, double>,
+                                         std::tuple<double, double, bool>, std::tuple<double, double>>(
+        tables_path / "In_d_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d_d-d_d.parquet",
+        tables_path / ("Err_d_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [m, n, c, cv, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::oblate_aswfa(m, n, c, cv, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, cv, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, cv, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_oblate_aswfa_nocv.cpp b/tests/scipy_special_tests/test_oblate_aswfa_nocv.cpp
new file mode 100644
index 0000000000..bbe125b04e
--- /dev/null
+++ b/tests/scipy_special_tests/test_oblate_aswfa_nocv.cpp
@@ -0,0 +1,37 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "oblate_aswfa_nocv"};
+
+TEST_CASE("oblate_aswfa_nocv dddd->dd scipy_special_tests", "[oblate_aswfa_nocv][dddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<double, double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [m, n, c, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::oblate_aswfa_nocv(m, n, c, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_oblate_radial1.cpp b/tests/scipy_special_tests/test_oblate_radial1.cpp
new file mode 100644
index 0000000000..c376bb4a36
--- /dev/null
+++ b/tests/scipy_special_tests/test_oblate_radial1.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "oblate_radial1"};
+
+TEST_CASE("oblate_radial1 ddddd->dd scipy_special_tests", "[oblate_radial1][ddddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<
+                                         std::tuple<double, double, double, double, double>,
+                                         std::tuple<double, double, bool>, std::tuple<double, double>>(
+        tables_path / "In_d_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d_d-d_d.parquet",
+        tables_path / ("Err_d_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [m, n, c, cv, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::oblate_radial1(m, n, c, cv, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, cv, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, cv, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_oblate_radial1_nocv.cpp b/tests/scipy_special_tests/test_oblate_radial1_nocv.cpp
new file mode 100644
index 0000000000..5aaabd5aab
--- /dev/null
+++ b/tests/scipy_special_tests/test_oblate_radial1_nocv.cpp
@@ -0,0 +1,37 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "oblate_radial1_nocv"};
+
+TEST_CASE("oblate_radial1_nocv dddd->dd scipy_special_tests", "[oblate_radial1_nocv][dddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<double, double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [m, n, c, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::oblate_radial1_nocv(m, n, c, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_oblate_radial2.cpp b/tests/scipy_special_tests/test_oblate_radial2.cpp
new file mode 100644
index 0000000000..cc54bf47cd
--- /dev/null
+++ b/tests/scipy_special_tests/test_oblate_radial2.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "oblate_radial2"};
+
+TEST_CASE("oblate_radial2 ddddd->dd scipy_special_tests", "[oblate_radial2][ddddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<
+                                         std::tuple<double, double, double, double, double>,
+                                         std::tuple<double, double, bool>, std::tuple<double, double>>(
+        tables_path / "In_d_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d_d-d_d.parquet",
+        tables_path / ("Err_d_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [m, n, c, cv, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::oblate_radial2(m, n, c, cv, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, cv, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, cv, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_oblate_radial2_nocv.cpp b/tests/scipy_special_tests/test_oblate_radial2_nocv.cpp
new file mode 100644
index 0000000000..e71f8e2740
--- /dev/null
+++ b/tests/scipy_special_tests/test_oblate_radial2_nocv.cpp
@@ -0,0 +1,37 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "oblate_radial2_nocv"};
+
+TEST_CASE("oblate_radial2_nocv dddd->dd scipy_special_tests", "[oblate_radial2_nocv][dddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<double, double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [m, n, c, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::oblate_radial2_nocv(m, n, c, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_owens_t.cpp b/tests/scipy_special_tests/test_owens_t.cpp
new file mode 100644
index 0000000000..1261d2fb51
--- /dev/null
+++ b/tests/scipy_special_tests/test_owens_t.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "owens_t"};
+
+TEST_CASE("owens_t dd->d scipy_special_tests", "[owens_t][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [h, a] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::owens_t(h, a);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(h, a, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_pbdv.cpp b/tests/scipy_special_tests/test_pbdv.cpp
new file mode 100644
index 0000000000..6229ce54db
--- /dev/null
+++ b/tests/scipy_special_tests/test_pbdv.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/par_cyl.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "pbdv"};
+
+TEST_CASE("pbdv dd->dd scipy_special_tests", "[pbdv][dd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<std::tuple<double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d-d_d.parquet", tables_path / "Out_d_d-d_d.parquet",
+            tables_path / ("Err_d_d-d_d_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [v, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::pbdv(v, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(v, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(v, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_pbvv.cpp b/tests/scipy_special_tests/test_pbvv.cpp
new file mode 100644
index 0000000000..f7af9763e3
--- /dev/null
+++ b/tests/scipy_special_tests/test_pbvv.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/par_cyl.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "pbvv"};
+
+TEST_CASE("pbvv dd->dd scipy_special_tests", "[pbvv][dd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<std::tuple<double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d-d_d.parquet", tables_path / "Out_d_d-d_d.parquet",
+            tables_path / ("Err_d_d-d_d_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [v, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::pbvv(v, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(v, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(v, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_pbwa.cpp b/tests/scipy_special_tests/test_pbwa.cpp
new file mode 100644
index 0000000000..f153504fcd
--- /dev/null
+++ b/tests/scipy_special_tests/test_pbwa.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/par_cyl.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "pbwa"};
+
+TEST_CASE("pbwa dd->dd scipy_special_tests", "[pbwa][dd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<std::tuple<double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d-d_d.parquet", tables_path / "Out_d_d-d_d.parquet",
+            tables_path / ("Err_d_d-d_d_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [v, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::pbwa(v, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(v, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(v, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_pdtr.cpp b/tests/scipy_special_tests/test_pdtr.cpp
new file mode 100644
index 0000000000..283d5630bd
--- /dev/null
+++ b/tests/scipy_special_tests/test_pdtr.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "pdtr"};
+
+TEST_CASE("pdtr dd->d scipy_special_tests", "[pdtr][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [k, m] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::pdtr(k, m);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(k, m, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_pdtrc.cpp b/tests/scipy_special_tests/test_pdtrc.cpp
new file mode 100644
index 0000000000..f90f770875
--- /dev/null
+++ b/tests/scipy_special_tests/test_pdtrc.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "pdtrc"};
+
+TEST_CASE("pdtrc dd->d scipy_special_tests", "[pdtrc][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [k, m] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::pdtrc(k, m);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(k, m, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_pdtri.cpp b/tests/scipy_special_tests/test_pdtri.cpp
new file mode 100644
index 0000000000..a2eddc80e0
--- /dev/null
+++ b/tests/scipy_special_tests/test_pdtri.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "pdtri"};
+
+TEST_CASE("pdtri dd->d scipy_special_tests", "[pdtri][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [k, y] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::pdtri(k, y);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(k, y, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_pmv.cpp b/tests/scipy_special_tests/test_pmv.cpp
new file mode 100644
index 0000000000..3ce4c5468b
--- /dev/null
+++ b/tests/scipy_special_tests/test_pmv.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/specfun.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "pmv"};
+
+TEST_CASE("pmv ddd->d scipy_special_tests", "[pmv][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [m, v, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::pmv(m, v, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(m, v, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_poch.cpp b/tests/scipy_special_tests/test_poch.cpp
new file mode 100644
index 0000000000..b1f50cea5d
--- /dev/null
+++ b/tests/scipy_special_tests/test_poch.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/cephes/poch.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "poch"};
+
+TEST_CASE("cephes::poch dd->d scipy_special_tests", "[cephes::poch][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [z, m] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cephes::poch(z, m);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, m, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_prolate_aswfa.cpp b/tests/scipy_special_tests/test_prolate_aswfa.cpp
new file mode 100644
index 0000000000..82d66dbca4
--- /dev/null
+++ b/tests/scipy_special_tests/test_prolate_aswfa.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "prolate_aswfa"};
+
+TEST_CASE("prolate_aswfa ddddd->dd scipy_special_tests", "[prolate_aswfa][ddddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<
+                                         std::tuple<double, double, double, double, double>,
+                                         std::tuple<double, double, bool>, std::tuple<double, double>>(
+        tables_path / "In_d_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d_d-d_d.parquet",
+        tables_path / ("Err_d_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [m, n, c, cv, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::prolate_aswfa(m, n, c, cv, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, cv, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, cv, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_prolate_aswfa_nocv.cpp b/tests/scipy_special_tests/test_prolate_aswfa_nocv.cpp
new file mode 100644
index 0000000000..a459deb5ff
--- /dev/null
+++ b/tests/scipy_special_tests/test_prolate_aswfa_nocv.cpp
@@ -0,0 +1,37 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "prolate_aswfa_nocv"};
+
+TEST_CASE("prolate_aswfa_nocv dddd->dd scipy_special_tests", "[prolate_aswfa_nocv][dddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<double, double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [m, n, c, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::prolate_aswfa_nocv(m, n, c, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_prolate_radial1.cpp b/tests/scipy_special_tests/test_prolate_radial1.cpp
new file mode 100644
index 0000000000..59ffa091bb
--- /dev/null
+++ b/tests/scipy_special_tests/test_prolate_radial1.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "prolate_radial1"};
+
+TEST_CASE("prolate_radial1 ddddd->dd scipy_special_tests", "[prolate_radial1][ddddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<
+                                         std::tuple<double, double, double, double, double>,
+                                         std::tuple<double, double, bool>, std::tuple<double, double>>(
+        tables_path / "In_d_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d_d-d_d.parquet",
+        tables_path / ("Err_d_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [m, n, c, cv, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::prolate_radial1(m, n, c, cv, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, cv, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, cv, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_prolate_radial1_nocv.cpp b/tests/scipy_special_tests/test_prolate_radial1_nocv.cpp
new file mode 100644
index 0000000000..17439503ad
--- /dev/null
+++ b/tests/scipy_special_tests/test_prolate_radial1_nocv.cpp
@@ -0,0 +1,39 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "prolate_radial1_nocv"};
+
+TEST_CASE(
+    "prolate_radial1_nocv dddd->dd scipy_special_tests", "[prolate_radial1_nocv][dddd->dd][scipy_special_tests]"
+) {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<double, double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [m, n, c, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::prolate_radial1_nocv(m, n, c, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_prolate_radial2.cpp b/tests/scipy_special_tests/test_prolate_radial2.cpp
new file mode 100644
index 0000000000..d546ee0c82
--- /dev/null
+++ b/tests/scipy_special_tests/test_prolate_radial2.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "prolate_radial2"};
+
+TEST_CASE("prolate_radial2 ddddd->dd scipy_special_tests", "[prolate_radial2][ddddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<
+                                         std::tuple<double, double, double, double, double>,
+                                         std::tuple<double, double, bool>, std::tuple<double, double>>(
+        tables_path / "In_d_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d_d-d_d.parquet",
+        tables_path / ("Err_d_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [m, n, c, cv, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::prolate_radial2(m, n, c, cv, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, cv, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, cv, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_prolate_radial2_nocv.cpp b/tests/scipy_special_tests/test_prolate_radial2_nocv.cpp
new file mode 100644
index 0000000000..3476c77efc
--- /dev/null
+++ b/tests/scipy_special_tests/test_prolate_radial2_nocv.cpp
@@ -0,0 +1,39 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "prolate_radial2_nocv"};
+
+TEST_CASE(
+    "prolate_radial2_nocv dddd->dd scipy_special_tests", "[prolate_radial2_nocv][dddd->dd][scipy_special_tests]"
+) {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<double, double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d_d-d_d.parquet", tables_path / "Out_d_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        )
+    );
+
+    auto [m, n, c, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::prolate_radial2_nocv(m, n, c, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, n, c, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, n, c, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_prolate_segv.cpp b/tests/scipy_special_tests/test_prolate_segv.cpp
new file mode 100644
index 0000000000..96192a3eb0
--- /dev/null
+++ b/tests/scipy_special_tests/test_prolate_segv.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/sphd_wave.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "prolate_segv"};
+
+TEST_CASE("prolate_segv ddd->d scipy_special_tests", "[prolate_segv][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [m, n, c] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::prolate_segv(m, n, c);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(m, n, c, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_radian.cpp b/tests/scipy_special_tests/test_radian.cpp
new file mode 100644
index 0000000000..38b7522928
--- /dev/null
+++ b/tests/scipy_special_tests/test_radian.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/trig.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "radian"};
+
+TEST_CASE("radian ddd->d scipy_special_tests", "[radian][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [d, m, s] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::radian(d, m, s);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(d, m, s, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_rgamma.cpp b/tests/scipy_special_tests/test_rgamma.cpp
new file mode 100644
index 0000000000..af7a600f32
--- /dev/null
+++ b/tests/scipy_special_tests/test_rgamma.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/loggamma.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "rgamma"};
+
+TEST_CASE("rgamma D->D scipy_special_tests", "[rgamma][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::rgamma(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_riemann_zeta.cpp b/tests/scipy_special_tests/test_riemann_zeta.cpp
new file mode 100644
index 0000000000..805af795fd
--- /dev/null
+++ b/tests/scipy_special_tests/test_riemann_zeta.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/zeta.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "riemann_zeta"};
+
+TEST_CASE("riemann_zeta D->D scipy_special_tests", "[riemann_zeta][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::riemann_zeta(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_round.cpp b/tests/scipy_special_tests/test_round.cpp
new file mode 100644
index 0000000000..338c3f0262
--- /dev/null
+++ b/tests/scipy_special_tests/test_round.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/cephes/round.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "round"};
+
+TEST_CASE("cephes::round d->d scipy_special_tests", "[cephes::round][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cephes::round(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_scaled_exp1.cpp b/tests/scipy_special_tests/test_scaled_exp1.cpp
new file mode 100644
index 0000000000..4b4cdb6ce0
--- /dev/null
+++ b/tests/scipy_special_tests/test_scaled_exp1.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/expint.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "scaled_exp1"};
+
+TEST_CASE("scaled_exp1 d->d scipy_special_tests", "[scaled_exp1][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::scaled_exp1(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_sem.cpp b/tests/scipy_special_tests/test_sem.cpp
new file mode 100644
index 0000000000..b90fc85dda
--- /dev/null
+++ b/tests/scipy_special_tests/test_sem.cpp
@@ -0,0 +1,36 @@
+#include "../testing_utils.h"
+
+#include <xsf/mathieu.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "sem"};
+
+TEST_CASE("sem ddd->dd scipy_special_tests", "[sem][ddd->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<
+                 std::tuple<double, double, double>, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d_d_d-d_d.parquet", tables_path / "Out_d_d_d-d_d.parquet",
+            tables_path / ("Err_d_d_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [m, q, x] = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::sem(m, q, x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(m, q, x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(m, q, x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_sem_cva.cpp b/tests/scipy_special_tests/test_sem_cva.cpp
new file mode 100644
index 0000000000..815fd5e841
--- /dev/null
+++ b/tests/scipy_special_tests/test_sem_cva.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/mathieu.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "sem_cva"};
+
+TEST_CASE("sem_cva dd->d scipy_special_tests", "[sem_cva][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [m, q] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::sem_cva(m, q);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(m, q, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_shichi.cpp b/tests/scipy_special_tests/test_shichi.cpp
new file mode 100644
index 0000000000..e94e2d88c4
--- /dev/null
+++ b/tests/scipy_special_tests/test_shichi.cpp
@@ -0,0 +1,35 @@
+#include "../testing_utils.h"
+
+#include <xsf/sici.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "shichi"};
+
+TEST_CASE("shichi d->dd scipy_special_tests", "[shichi][d->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<double, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d-d_d.parquet", tables_path / "Out_d-d_d.parquet",
+            tables_path / ("Err_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::shichi(x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_sici.cpp b/tests/scipy_special_tests/test_sici.cpp
new file mode 100644
index 0000000000..477f45e61a
--- /dev/null
+++ b/tests/scipy_special_tests/test_sici.cpp
@@ -0,0 +1,35 @@
+#include "../testing_utils.h"
+
+#include <xsf/sici.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "sici"};
+
+TEST_CASE("sici d->dd scipy_special_tests", "[sici][d->dd][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<double, std::tuple<double, double, bool>, std::tuple<double, double>>(
+            tables_path / "In_d-d_d.parquet", tables_path / "Out_d-d_d.parquet",
+            tables_path / ("Err_d-d_d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired0, desired1, fallback] = output;
+
+    double out0;
+    double out1;
+
+    xsf::sici(x, out0, out1);
+    auto [tol0, tol1] = tol;
+
+    auto error0 = xsf::extended_relative_error(out0, desired0);
+    tol0 = adjust_tolerance(tol0);
+    CAPTURE(x, out0, desired0, error0, tol0, fallback);
+    REQUIRE(error0 <= tol0);
+
+    auto error1 = xsf::extended_relative_error(out1, desired1);
+    tol1 = adjust_tolerance(tol1);
+    CAPTURE(x, out1, desired1, error1, tol1, fallback);
+    REQUIRE(error1 <= tol1);
+}
diff --git a/tests/scipy_special_tests/test_sindg.cpp b/tests/scipy_special_tests/test_sindg.cpp
new file mode 100644
index 0000000000..8b45b95e9d
--- /dev/null
+++ b/tests/scipy_special_tests/test_sindg.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/trig.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "sindg"};
+
+TEST_CASE("sindg d->d scipy_special_tests", "[sindg][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::sindg(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_sinpi.cpp b/tests/scipy_special_tests/test_sinpi.cpp
new file mode 100644
index 0000000000..42552d862e
--- /dev/null
+++ b/tests/scipy_special_tests/test_sinpi.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/trig.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "sinpi"};
+
+TEST_CASE("sinpi D->D scipy_special_tests", "[sinpi][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::sinpi(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_smirnov.cpp b/tests/scipy_special_tests/test_smirnov.cpp
new file mode 100644
index 0000000000..7c26aa6b49
--- /dev/null
+++ b/tests/scipy_special_tests/test_smirnov.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "smirnov"};
+
+TEST_CASE("smirnov dd->d scipy_special_tests", "[smirnov][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [n, d] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::smirnov(n, d);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(n, d, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_smirnovc.cpp b/tests/scipy_special_tests/test_smirnovc.cpp
new file mode 100644
index 0000000000..50d922d3d3
--- /dev/null
+++ b/tests/scipy_special_tests/test_smirnovc.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "smirnovc"};
+
+TEST_CASE("smirnovc pd->d scipy_special_tests", "[smirnovc][pd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<std::ptrdiff_t, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_p_d-d.parquet", tables_path / "Out_p_d-d.parquet",
+            tables_path / ("Err_p_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [n, d] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::smirnovc(n, d);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(n, d, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_smirnovci.cpp b/tests/scipy_special_tests/test_smirnovci.cpp
new file mode 100644
index 0000000000..890d890063
--- /dev/null
+++ b/tests/scipy_special_tests/test_smirnovci.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "smirnovci"};
+
+TEST_CASE("smirnovci pd->d scipy_special_tests", "[smirnovci][pd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<std::ptrdiff_t, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_p_d-d.parquet", tables_path / "Out_p_d-d.parquet",
+            tables_path / ("Err_p_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [n, p] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::smirnovci(n, p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(n, p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_smirnovi.cpp b/tests/scipy_special_tests/test_smirnovi.cpp
new file mode 100644
index 0000000000..02edd9f194
--- /dev/null
+++ b/tests/scipy_special_tests/test_smirnovi.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "smirnovi"};
+
+TEST_CASE("smirnovi dd->d scipy_special_tests", "[smirnovi][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [n, p] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::smirnovi(n, p);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(n, p, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_smirnovp.cpp b/tests/scipy_special_tests/test_smirnovp.cpp
new file mode 100644
index 0000000000..fbec93e684
--- /dev/null
+++ b/tests/scipy_special_tests/test_smirnovp.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/stats.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "smirnovp"};
+
+TEST_CASE("smirnovp pd->d scipy_special_tests", "[smirnovp][pd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<std::ptrdiff_t, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_p_d-d.parquet", tables_path / "Out_p_d-d.parquet",
+            tables_path / ("Err_p_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [n, d] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::smirnovp(n, d);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(n, d, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_spence.cpp b/tests/scipy_special_tests/test_spence.cpp
new file mode 100644
index 0000000000..13a982dc1b
--- /dev/null
+++ b/tests/scipy_special_tests/test_spence.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/cephes/spence.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "spence"};
+
+TEST_CASE("cephes::spence d->d scipy_special_tests", "[cephes::spence][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::cephes::spence(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_struve_h.cpp b/tests/scipy_special_tests/test_struve_h.cpp
new file mode 100644
index 0000000000..f56aa21751
--- /dev/null
+++ b/tests/scipy_special_tests/test_struve_h.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/struve.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "struve_h"};
+
+TEST_CASE("struve_h dd->d scipy_special_tests", "[struve_h][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::struve_h(v, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_struve_l.cpp b/tests/scipy_special_tests/test_struve_l.cpp
new file mode 100644
index 0000000000..62c523ac2e
--- /dev/null
+++ b/tests/scipy_special_tests/test_struve_l.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/struve.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "struve_l"};
+
+TEST_CASE("struve_l dd->d scipy_special_tests", "[struve_l][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [v, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::struve_l(v, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(v, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_tandg.cpp b/tests/scipy_special_tests/test_tandg.cpp
new file mode 100644
index 0000000000..2992521662
--- /dev/null
+++ b/tests/scipy_special_tests/test_tandg.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/trig.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "tandg"};
+
+TEST_CASE("tandg d->d scipy_special_tests", "[tandg][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::tandg(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_voigt_profile.cpp b/tests/scipy_special_tests/test_voigt_profile.cpp
new file mode 100644
index 0000000000..4edf3bfe17
--- /dev/null
+++ b/tests/scipy_special_tests/test_voigt_profile.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/erf.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "voigt_profile"};
+
+TEST_CASE("voigt_profile ddd->d scipy_special_tests", "[voigt_profile][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [x, sigma, gamma] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::voigt_profile(x, sigma, gamma);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, sigma, gamma, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_wofz.cpp b/tests/scipy_special_tests/test_wofz.cpp
new file mode 100644
index 0000000000..3d641d4da7
--- /dev/null
+++ b/tests/scipy_special_tests/test_wofz.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/erf.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "wofz"};
+
+TEST_CASE("wofz D->D scipy_special_tests", "[wofz][D->D][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::complex<double>, std::tuple<std::complex<double>, bool>, double>(
+            tables_path / "In_cd-cd.parquet", tables_path / "Out_cd-cd.parquet",
+            tables_path / ("Err_cd-cd_" + get_platform_str() + ".parquet")
+        ));
+
+    auto x = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::wofz(x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_wright_bessel.cpp b/tests/scipy_special_tests/test_wright_bessel.cpp
new file mode 100644
index 0000000000..cf68afdad6
--- /dev/null
+++ b/tests/scipy_special_tests/test_wright_bessel.cpp
@@ -0,0 +1,24 @@
+#include "../testing_utils.h"
+
+#include <xsf/wright_bessel.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "wright_bessel"};
+
+TEST_CASE("wright_bessel ddd->d scipy_special_tests", "[wright_bessel][ddd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] =
+        GENERATE(xsf_test_cases<std::tuple<double, double, double>, std::tuple<double, bool>, double>(
+            tables_path / "In_d_d_d-d.parquet", tables_path / "Out_d_d_d-d.parquet",
+            tables_path / ("Err_d_d_d-d_" + get_platform_str() + ".parquet")
+        ));
+
+    auto [a, b, x] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::wright_bessel(a, b, x);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(a, b, x, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_xlog1py.cpp b/tests/scipy_special_tests/test_xlog1py.cpp
new file mode 100644
index 0000000000..c4cccf1511
--- /dev/null
+++ b/tests/scipy_special_tests/test_xlog1py.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/log.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "xlog1py"};
+
+TEST_CASE("xlog1py dd->d scipy_special_tests", "[xlog1py][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [x, y] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::xlog1py(x, y);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, y, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_xlogy.cpp b/tests/scipy_special_tests/test_xlogy.cpp
new file mode 100644
index 0000000000..07a1cbc8db
--- /dev/null
+++ b/tests/scipy_special_tests/test_xlogy.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/log.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "xlogy"};
+
+TEST_CASE("xlogy dd->d scipy_special_tests", "[xlogy][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [x, y] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::xlogy(x, y);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(x, y, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_zeta.cpp b/tests/scipy_special_tests/test_zeta.cpp
new file mode 100644
index 0000000000..2063f45d12
--- /dev/null
+++ b/tests/scipy_special_tests/test_zeta.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/zeta.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "zeta"};
+
+TEST_CASE("zeta dd->d scipy_special_tests", "[zeta][dd->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<std::tuple<double, double>, std::tuple<double, bool>, double>(
+        tables_path / "In_d_d-d.parquet", tables_path / "Out_d_d-d.parquet",
+        tables_path / ("Err_d_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto [z, q] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::zeta(z, q);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, q, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/scipy_special_tests/test_zetac.cpp b/tests/scipy_special_tests/test_zetac.cpp
new file mode 100644
index 0000000000..1379f570a8
--- /dev/null
+++ b/tests/scipy_special_tests/test_zetac.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+
+#include <xsf/zeta.h>
+
+namespace fs = std::filesystem;
+
+fs::path tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "zetac"};
+
+TEST_CASE("zetac d->d scipy_special_tests", "[zetac][d->d][scipy_special_tests]") {
+    SET_FP_FORMAT()
+    auto [input, output, tol] = GENERATE(xsf_test_cases<double, std::tuple<double, bool>, double>(
+        tables_path / "In_d-d.parquet", tables_path / "Out_d-d.parquet",
+        tables_path / ("Err_d-d_" + get_platform_str() + ".parquet")
+    ));
+
+    auto z = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::zetac(z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    CAPTURE(z, out, desired, error, tol, fallback);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/test_hyp2f1.cpp b/tests/test_hyp2f1.cpp
new file mode 100644
index 0000000000..c3d7e21cd6
--- /dev/null
+++ b/tests/test_hyp2f1.cpp
@@ -0,0 +1,68 @@
+#include <complex>
+#include <filesystem>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <tuple>
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <xsf/fp_error_metrics.h>
+#include <xsf/hyp2f1.h>
+
+#include "testing_utils.h"
+
+namespace fs = std::filesystem;
+
+fs::path hyp2f1_tables_path{fs::path(XSREF_TABLES_PATH) / "scipy_special_tests" / "hyp2f1"};
+
+
+TEST_CASE("hyp2f1 complex scipy.special cases", "[hyp2f1][complex][scipy-special]") {
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<double, double, double, std::complex<double>>, std::tuple<std::complex<double>, bool>, double>(
+            hyp2f1_tables_path / "In_d_d_d_cd-cd.parquet", hyp2f1_tables_path / "Out_d_d_d_cd-cd.parquet",
+            hyp2f1_tables_path / ("Err_d_d_d_cd-cd_" + get_platform_str() + ".parquet")
+        )
+    );
+    auto [a, b, c, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::hyp2f1(a, b, c, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    INFO("a := " << std::setprecision(std::numeric_limits<double>::max_digits10) << a << '\n'
+	 << "b := " << b << '\n'
+	 << "c := " << c << '\n'
+	 << "z := " << z << '\n'
+	 << "out := " << out << '\n'
+	 << "desired := " << desired << '\n'
+	 << "error := " << error << '\n'
+	 << "tolerance := " << tol << '\n'
+	);
+    REQUIRE(error <= tol);
+}
+
+TEST_CASE("hyp2f1 real scipy.special cases", "[hyp2f1][real][scipy-special]") {
+    auto [input, output, tol] = GENERATE(
+        xsf_test_cases<
+            std::tuple<double, double, double, double>, std::tuple<double, bool>, double>(
+            hyp2f1_tables_path / "In_d_d_d_d-d.parquet", hyp2f1_tables_path / "Out_d_d_d_d-d.parquet",
+            hyp2f1_tables_path / ("Err_d_d_d_d-d_" + get_platform_str() + ".parquet")
+        )
+    );
+    auto [a, b, c, z] = input;
+    auto [desired, fallback] = output;
+    auto out = xsf::hyp2f1(a, b, c, z);
+    auto error = xsf::extended_relative_error(out, desired);
+    tol = adjust_tolerance(tol);
+    INFO("a := " << std::setprecision(std::numeric_limits<double>::max_digits10) << a << '\n'
+	 << "b := " << b << '\n'
+	 << "c := " << c << '\n'
+	 << "z := " << z << '\n'
+	 << "out := " << out << '\n'
+	 << "desired := " << desired << '\n'
+	 << "error := " << error << '\n'
+	 << "tolerance := " << tol << '\n'
+	);
+    REQUIRE(error <= tol);
+}
diff --git a/tests/testing_utils.h b/tests/testing_utils.h
new file mode 100644
index 0000000000..49448d1cda
--- /dev/null
+++ b/tests/testing_utils.h
@@ -0,0 +1,151 @@
+#pragma once
+
+#include <algorithm>
+#include <complex>
+#include <filesystem>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <stdexcept>
+#include <tuple>
+#include <type_traits>
+
+#include <arrow/io/file.h>
+#include <parquet/stream_reader.h>
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators.hpp>
+#include <catch2/generators/catch_generators_adapters.hpp>
+
+#include <xsf/fp_error_metrics.h>
+
+namespace {
+
+template <typename T>
+struct remove_complex {
+    using type = T;
+};
+
+template <typename T>
+struct remove_complex<std::complex<T>> {
+    using type = T;
+};
+
+template <typename T>
+using remove_complex_t = typename remove_complex<T>::type;
+
+template <typename T>
+class TableReader {
+  public:
+    TableReader(const std::string &file_path) {
+        PARQUET_ASSIGN_OR_THROW(infile_, arrow::io::ReadableFile::Open(file_path));
+        parquet::StreamReader stream{parquet::ParquetFileReader::Open(infile_)};
+        stream_ = std::make_unique<parquet::StreamReader>(std::move(stream));
+        file_path_ = std::move(file_path);
+    }
+
+    T next() {
+        T row;
+        fill_row(row);
+        return row;
+    }
+
+    bool eof() const { return stream_->eof(); }
+
+  private:
+    void fill_row(T &elements) {
+        if constexpr (std::is_scalar_v<remove_complex_t<T>>) {
+            fill_element(elements);
+        } else {
+            std::apply([this](auto &...x) { (fill_element(x), ...); }, elements);
+        }
+        stream_->EndRow();
+    }
+
+    template <typename U>
+    void fill_element(U &element) {
+        if constexpr (std::is_same_v<U, std::complex<remove_complex_t<U>>>) {
+            using V = remove_complex_t<U>;
+            V real;
+            V imag;
+            *stream_ >> real >> imag;
+            element = U(real, imag);
+        } else {
+            *stream_ >> element;
+        }
+    }
+
+    std::shared_ptr<arrow::io::ReadableFile> infile_;
+    std::unique_ptr<parquet::StreamReader> stream_;
+    std::string file_path_;
+};
+
+template <typename T1, typename T2, typename T3>
+class XsfTestCaseGenerator final : public Catch::Generators::IGenerator<std::tuple<T1, T2, T3>> {
+  public:
+    XsfTestCaseGenerator(
+        std::unique_ptr<TableReader<T1>> input_reader, std::unique_ptr<TableReader<T2>> output_reader,
+        std::unique_ptr<TableReader<T3>> tol_reader
+    )
+        : input_reader_(std::move(input_reader)), output_reader_(std::move(output_reader)),
+          tol_reader_(std::move(tol_reader)) {
+        if (!next()) {
+            throw std::runtime_error("XsfTestCaseGenerator received an empty table\n");
+        }
+    }
+
+    std::tuple<T1, T2, T3> const &get() const override { return current_case_; }
+
+    bool next() override {
+        if (input_reader_->eof() || output_reader_->eof() || tol_reader_->eof()) {
+            return false;
+        }
+        current_case_ = std::make_tuple(input_reader_->next(), output_reader_->next(), tol_reader_->next());
+        return true;
+    }
+
+  private:
+    std::unique_ptr<TableReader<T1>> input_reader_;
+    std::unique_ptr<TableReader<T2>> output_reader_;
+    std::unique_ptr<TableReader<T3>> tol_reader_;
+    std::tuple<T1, T2, T3> current_case_;
+};
+
+template <typename T1, typename T2, typename T3>
+Catch::Generators::GeneratorWrapper<std::tuple<T1, T2, T3>> xsf_test_cases(
+    const std::filesystem::path &input_path, const std::filesystem::path &output_path,
+    const std::filesystem::path &tol_path
+) {
+    auto input_reader = std::make_unique<TableReader<T1>>(input_path.string());
+    auto output_reader = std::make_unique<TableReader<T2>>(output_path.string());
+    auto tol_reader = std::make_unique<TableReader<T3>>(tol_path.string());
+
+    return Catch::Generators::GeneratorWrapper<std::tuple<T1, T2, T3>>(
+        Catch::Detail::make_unique<XsfTestCaseGenerator<T1, T2, T3>>(
+            std::move(input_reader), std::move(output_reader), std::move(tol_reader)
+        )
+    );
+}
+
+
+template <typename T>
+T adjust_tolerance(T tol) {
+    // Add some wiggle room to tolerance from table.
+    return 2 * std::max(std::numeric_limits<T>::epsilon(), tol);
+}
+
+
+std::string get_platform_str() {
+    /* This is going to get a string "<compiler>-<os>-<architecture>" using conditional
+     * compilation, but for now we're just stubbing things out. */
+    return "gcc-linux-x86_64";
+}
+
+
+
+} // namespace
+
+
+#define SET_FP_FORMAT()							\
+    Catch::StringMaker<double>::precision = std::numeric_limits<double>::max_digits10; \
+    Catch::StringMaker<float>::precision = std::numeric_limits<float>::max_digits10;