diff --git a/src/airy.rs b/src/airy.rs
index 42872ea..3dafd91 100644
--- a/src/airy.rs
+++ b/src/airy.rs
@@ -217,17 +217,17 @@ pub fn airyzo(nt: u32, kf: AiryKind) -> (f64, f64, f64, f64) {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_airy_f64() {
-        xsref::test::<(f64, f64, f64, f64), _>("airy", "d-d_d_d_d", |x: &[f64]| airy(x[0]));
+        testing::test::<(f64, f64, f64, f64), _>("airy", "d-d_d_d_d", |x: &[f64]| airy(x[0]));
     }
 
     #[test]
     fn test_airy_c64() {
-        xsref::test::<(Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>), _>(
+        testing::test::<(Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>), _>(
             "airy",
             "cd-cd_cd_cd_cd",
             |x: &[f64]| airy(c64(x[0], x[1])),
@@ -236,12 +236,12 @@ mod tests {
 
     #[test]
     fn test_airye_f64() {
-        xsref::test::<(f64, f64, f64, f64), _>("airye", "d-d_d_d_d", |x: &[f64]| airye(x[0]));
+        testing::test::<(f64, f64, f64, f64), _>("airye", "d-d_d_d_d", |x: &[f64]| airye(x[0]));
     }
 
     #[test]
     fn test_airye_c64() {
-        xsref::test::<(Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>), _>(
+        testing::test::<(Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>), _>(
             "airye",
             "cd-cd_cd_cd_cd",
             |x: &[f64]| airye(c64(x[0], x[1])),
@@ -250,6 +250,6 @@ mod tests {
 
     #[test]
     fn test_itairy() {
-        xsref::test::<(f64, f64, f64, f64), _>("itairy", "d-d_d_d_d", |x: &[f64]| itairy(x[0]));
+        testing::test::<(f64, f64, f64, f64), _>("itairy", "d-d_d_d_d", |x: &[f64]| itairy(x[0]));
     }
 }
diff --git a/src/alg.rs b/src/alg.rs
index 8bd706c..9fa8be4 100644
--- a/src/alg.rs
+++ b/src/alg.rs
@@ -8,10 +8,10 @@ pub fn cbrt(x: f64) -> f64 {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_cbrt() {
-        xsref::test::<f64, _>("cbrt", "d-d", |x: &[f64]| cbrt(x[0]));
+        testing::test::<f64, _>("cbrt", "d-d", |x: &[f64]| cbrt(x[0]));
     }
 }
diff --git a/src/bessel.rs b/src/bessel.rs
index 8b4ade2..cb7d6f9 100644
--- a/src/bessel.rs
+++ b/src/bessel.rs
@@ -234,45 +234,45 @@ pub fn it2i0k0(x: f64) -> (f64, f64) {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     // bessel j
 
     #[test]
     fn test_cyl_bessel_j0() {
-        xsref::test::<f64, _>("cyl_bessel_j0", "d-d", |x: &[f64]| cyl_bessel_j0(x[0]));
+        testing::test::<f64, _>("cyl_bessel_j0", "d-d", |x: &[f64]| cyl_bessel_j0(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_j1() {
-        xsref::test::<f64, _>("cyl_bessel_j1", "d-d", |x: &[f64]| cyl_bessel_j1(x[0]));
+        testing::test::<f64, _>("cyl_bessel_j1", "d-d", |x: &[f64]| cyl_bessel_j1(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_j_f64() {
-        xsref::test::<f64, _>("cyl_bessel_j", "d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("cyl_bessel_j", "d_d-d", |x: &[f64]| {
             cyl_bessel_j(x[0], x[1])
         });
     }
 
     #[test]
     fn test_cyl_bessel_j_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_bessel_j", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_bessel_j", "d_cd-cd", |x: &[f64]| {
             cyl_bessel_j(x[0], c64(x[1], x[2]))
         });
     }
 
     #[test]
     fn test_cyl_bessel_je_f64() {
-        xsref::test::<f64, _>("cyl_bessel_je", "d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("cyl_bessel_je", "d_d-d", |x: &[f64]| {
             cyl_bessel_je(x[0], x[1])
         });
     }
 
     #[test]
     fn test_cyl_bessel_je_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_bessel_je", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_bessel_je", "d_cd-cd", |x: &[f64]| {
             cyl_bessel_je(x[0], c64(x[1], x[2]))
         });
     }
@@ -281,38 +281,38 @@ mod tests {
 
     #[test]
     fn test_cyl_bessel_y0() {
-        xsref::test::<f64, _>("cyl_bessel_y0", "d-d", |x: &[f64]| cyl_bessel_y0(x[0]));
+        testing::test::<f64, _>("cyl_bessel_y0", "d-d", |x: &[f64]| cyl_bessel_y0(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_y1() {
-        xsref::test::<f64, _>("cyl_bessel_y1", "d-d", |x: &[f64]| cyl_bessel_y1(x[0]));
+        testing::test::<f64, _>("cyl_bessel_y1", "d-d", |x: &[f64]| cyl_bessel_y1(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_y_f64() {
-        xsref::test::<f64, _>("cyl_bessel_y", "d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("cyl_bessel_y", "d_d-d", |x: &[f64]| {
             cyl_bessel_y(x[0], x[1])
         });
     }
 
     #[test]
     fn test_cyl_bessel_y_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_bessel_y", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_bessel_y", "d_cd-cd", |x: &[f64]| {
             cyl_bessel_y(x[0], c64(x[1], x[2]))
         });
     }
 
     #[test]
     fn test_cyl_bessel_ye_f64() {
-        xsref::test::<f64, _>("cyl_bessel_ye", "d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("cyl_bessel_ye", "d_d-d", |x: &[f64]| {
             cyl_bessel_ye(x[0], x[1])
         });
     }
 
     #[test]
     fn test_cyl_bessel_ye_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_bessel_ye", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_bessel_ye", "d_cd-cd", |x: &[f64]| {
             cyl_bessel_ye(x[0], c64(x[1], x[2]))
         });
     }
@@ -321,48 +321,48 @@ mod tests {
 
     #[test]
     fn test_cyl_bessel_i0() {
-        xsref::test::<f64, _>("cyl_bessel_i0", "d-d", |x: &[f64]| cyl_bessel_i0(x[0]));
+        testing::test::<f64, _>("cyl_bessel_i0", "d-d", |x: &[f64]| cyl_bessel_i0(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_i0e() {
-        xsref::test::<f64, _>("cyl_bessel_i0e", "d-d", |x: &[f64]| cyl_bessel_i0e(x[0]));
+        testing::test::<f64, _>("cyl_bessel_i0e", "d-d", |x: &[f64]| cyl_bessel_i0e(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_i1() {
-        xsref::test::<f64, _>("cyl_bessel_i1", "d-d", |x: &[f64]| cyl_bessel_i1(x[0]));
+        testing::test::<f64, _>("cyl_bessel_i1", "d-d", |x: &[f64]| cyl_bessel_i1(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_i1e() {
-        xsref::test::<f64, _>("cyl_bessel_i1e", "d-d", |x: &[f64]| cyl_bessel_i1e(x[0]));
+        testing::test::<f64, _>("cyl_bessel_i1e", "d-d", |x: &[f64]| cyl_bessel_i1e(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_i_f64() {
-        xsref::test::<f64, _>("cyl_bessel_i", "d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("cyl_bessel_i", "d_d-d", |x: &[f64]| {
             cyl_bessel_i(x[0], x[1])
         });
     }
 
     #[test]
     fn test_cyl_bessel_i_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_bessel_i", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_bessel_i", "d_cd-cd", |x: &[f64]| {
             cyl_bessel_i(x[0], c64(x[1], x[2]))
         });
     }
 
     #[test]
     fn test_cyl_bessel_ie_f64() {
-        xsref::test::<f64, _>("cyl_bessel_ie", "d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("cyl_bessel_ie", "d_d-d", |x: &[f64]| {
             cyl_bessel_ie(x[0], x[1])
         });
     }
 
     #[test]
     fn test_cyl_bessel_ie_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_bessel_ie", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_bessel_ie", "d_cd-cd", |x: &[f64]| {
             cyl_bessel_ie(x[0], c64(x[1], x[2]))
         });
     }
@@ -371,48 +371,48 @@ mod tests {
 
     #[test]
     fn test_cyl_bessel_k0() {
-        xsref::test::<f64, _>("cyl_bessel_k0", "d-d", |x: &[f64]| cyl_bessel_k0(x[0]));
+        testing::test::<f64, _>("cyl_bessel_k0", "d-d", |x: &[f64]| cyl_bessel_k0(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_k0e() {
-        xsref::test::<f64, _>("cyl_bessel_k0e", "d-d", |x: &[f64]| cyl_bessel_k0e(x[0]));
+        testing::test::<f64, _>("cyl_bessel_k0e", "d-d", |x: &[f64]| cyl_bessel_k0e(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_k1() {
-        xsref::test::<f64, _>("cyl_bessel_k1", "d-d", |x: &[f64]| cyl_bessel_k1(x[0]));
+        testing::test::<f64, _>("cyl_bessel_k1", "d-d", |x: &[f64]| cyl_bessel_k1(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_k1e() {
-        xsref::test::<f64, _>("cyl_bessel_k1e", "d-d", |x: &[f64]| cyl_bessel_k1e(x[0]));
+        testing::test::<f64, _>("cyl_bessel_k1e", "d-d", |x: &[f64]| cyl_bessel_k1e(x[0]));
     }
 
     #[test]
     fn test_cyl_bessel_k_f64() {
-        xsref::test::<f64, _>("cyl_bessel_k", "d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("cyl_bessel_k", "d_d-d", |x: &[f64]| {
             cyl_bessel_k(x[0], x[1])
         });
     }
 
     #[test]
     fn test_cyl_bessel_k_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_bessel_k", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_bessel_k", "d_cd-cd", |x: &[f64]| {
             cyl_bessel_k(x[0], c64(x[1], x[2]))
         });
     }
 
     #[test]
     fn test_cyl_bessel_ke_f64() {
-        xsref::test::<f64, _>("cyl_bessel_ke", "d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("cyl_bessel_ke", "d_d-d", |x: &[f64]| {
             cyl_bessel_ke(x[0], x[1])
         });
     }
 
     #[test]
     fn test_cyl_bessel_ke_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_bessel_ke", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_bessel_ke", "d_cd-cd", |x: &[f64]| {
             cyl_bessel_ke(x[0], c64(x[1], x[2]))
         });
     }
@@ -421,28 +421,28 @@ mod tests {
 
     #[test]
     fn test_cyl_hankel_1_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_hankel_1", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_hankel_1", "d_cd-cd", |x: &[f64]| {
             cyl_hankel_1(x[0], c64(x[1], x[2]))
         });
     }
 
     #[test]
     fn test_cyl_hankel_1e_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_hankel_1e", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_hankel_1e", "d_cd-cd", |x: &[f64]| {
             cyl_hankel_1e(x[0], c64(x[1], x[2]))
         });
     }
 
     #[test]
     fn test_cyl_hankel_2_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_hankel_2", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_hankel_2", "d_cd-cd", |x: &[f64]| {
             cyl_hankel_2(x[0], c64(x[1], x[2]))
         });
     }
 
     #[test]
     fn test_cyl_hankel_2e_c64() {
-        xsref::test::<Complex<f64>, _>("cyl_hankel_2e", "d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("cyl_hankel_2e", "d_cd-cd", |x: &[f64]| {
             cyl_hankel_2e(x[0], c64(x[1], x[2]))
         });
     }
@@ -451,7 +451,7 @@ mod tests {
 
     #[test]
     fn test_besselpoly_f64() {
-        xsref::test::<f64, _>("besselpoly", "d_d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("besselpoly", "d_d_d-d", |x: &[f64]| {
             besselpoly(x[0], x[1], x[2])
         });
     }
@@ -460,21 +460,21 @@ mod tests {
 
     #[test]
     fn test_it1j0y0_f64() {
-        xsref::test::<(f64, f64), _>("it1j0y0", "d-d_d", |x: &[f64]| it1j0y0(x[0]));
+        testing::test::<(f64, f64), _>("it1j0y0", "d-d_d", |x: &[f64]| it1j0y0(x[0]));
     }
 
     #[test]
     fn test_it2j0y0_f64() {
-        xsref::test::<(f64, f64), _>("it2j0y0", "d-d_d", |x: &[f64]| it2j0y0(x[0]));
+        testing::test::<(f64, f64), _>("it2j0y0", "d-d_d", |x: &[f64]| it2j0y0(x[0]));
     }
 
     #[test]
     fn test_it1i0k0_f64() {
-        xsref::test::<(f64, f64), _>("it1i0k0", "d-d_d", |x: &[f64]| it1i0k0(x[0]));
+        testing::test::<(f64, f64), _>("it1i0k0", "d-d_d", |x: &[f64]| it1i0k0(x[0]));
     }
 
     #[test]
     fn test_it2i0k0_f64() {
-        xsref::test::<(f64, f64), _>("it2i0k0", "d-d_d", |x: &[f64]| it2i0k0(x[0]));
+        testing::test::<(f64, f64), _>("it2i0k0", "d-d_d", |x: &[f64]| it2i0k0(x[0]));
     }
 }
diff --git a/src/beta.rs b/src/beta.rs
index 105642c..f2b2fb3 100644
--- a/src/beta.rs
+++ b/src/beta.rs
@@ -6,15 +6,15 @@ xsf_impl!(betaln, (a: f64, b: f64), "Natural log of `|beta|`");
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_beta() {
-        xsref::test::<f64, _>("beta", "d_d-d", |x: &[f64]| beta(x[0], x[1]));
+        testing::test::<f64, _>("beta", "d_d-d", |x: &[f64]| beta(x[0], x[1]));
     }
 
     #[test]
     fn test_betaln() {
-        xsref::test::<f64, _>("betaln", "d_d-d", |x: &[f64]| betaln(x[0], x[1]));
+        testing::test::<f64, _>("betaln", "d_d-d", |x: &[f64]| betaln(x[0], x[1]));
     }
 }
diff --git a/src/binom.rs b/src/binom.rs
index bd8f680..81330c8 100644
--- a/src/binom.rs
+++ b/src/binom.rs
@@ -5,10 +5,10 @@ xsf_impl!(binom, (n: f64, k: f64), "Binomial coefficient");
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_binom() {
-        xsref::test::<f64, _>("binom", "d_d-d", |x: &[f64]| binom(x[0], x[1]));
+        testing::test::<f64, _>("binom", "d_d-d", |x: &[f64]| binom(x[0], x[1]));
     }
 }
diff --git a/src/cdflib.rs b/src/cdflib.rs
index e3821cb..b5b370f 100644
--- a/src/cdflib.rs
+++ b/src/cdflib.rs
@@ -5,10 +5,10 @@ xsf_impl!(gdtrib, (a: f64, p: f64, x: f64), "Inverse of `p = gdtr(a, b, x)` with
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_gdtrib() {
-        xsref::test::<f64, _>("gdtrib", "d_d_d-d", |x: &[f64]| gdtrib(x[0], x[1], x[2]));
+        testing::test::<f64, _>("gdtrib", "d_d_d-d", |x: &[f64]| gdtrib(x[0], x[1], x[2]));
     }
 }
diff --git a/src/digamma.rs b/src/digamma.rs
index c06ab33..0b2872b 100644
--- a/src/digamma.rs
+++ b/src/digamma.rs
@@ -38,16 +38,16 @@ pub fn digamma<T: DigammaArg>(x: T) -> T::Output {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_digamma_f64() {
-        xsref::test::<f64, _>("digamma", "d-d", |x: &[f64]| digamma(x[0]));
+        testing::test::<f64, _>("digamma", "d-d", |x: &[f64]| digamma(x[0]));
     }
 
     #[test]
     fn test_digamma_c64() {
-        xsref::test::<num_complex::Complex<f64>, _>("digamma", "cd-cd", |x: &[f64]| {
+        testing::test::<num_complex::Complex<f64>, _>("digamma", "cd-cd", |x: &[f64]| {
             digamma(num_complex::c64(x[0], x[1]))
         });
     }
diff --git a/src/ellip.rs b/src/ellip.rs
index 72aed5a..806b09a 100644
--- a/src/ellip.rs
+++ b/src/ellip.rs
@@ -37,42 +37,42 @@ pub fn ellipj(u: f64, m: f64) -> (f64, f64, f64, f64) {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     // ellipk*
 
     #[test]
     fn test_ellipk() {
-        xsref::test::<f64, _>("ellipk", "d-d", |x: &[f64]| ellipk(x[0]));
+        testing::test::<f64, _>("ellipk", "d-d", |x: &[f64]| ellipk(x[0]));
     }
 
     #[test]
     fn test_ellipkm1() {
-        xsref::test::<f64, _>("ellipkm1", "d-d", |x: &[f64]| ellipkm1(x[0]));
+        testing::test::<f64, _>("ellipkm1", "d-d", |x: &[f64]| ellipkm1(x[0]));
     }
 
     #[test]
     fn test_ellipkinc() {
-        xsref::test::<f64, _>("ellipkinc", "d_d-d", |x: &[f64]| ellipkinc(x[0], x[1]));
+        testing::test::<f64, _>("ellipkinc", "d_d-d", |x: &[f64]| ellipkinc(x[0], x[1]));
     }
 
     // ellipe*
 
     #[test]
     fn test_ellipe() {
-        xsref::test::<f64, _>("ellipe", "d-d", |x: &[f64]| ellipe(x[0]));
+        testing::test::<f64, _>("ellipe", "d-d", |x: &[f64]| ellipe(x[0]));
     }
 
     #[test]
     fn test_ellipeinc() {
-        xsref::test::<f64, _>("ellipeinc", "d_d-d", |x: &[f64]| ellipeinc(x[0], x[1]));
+        testing::test::<f64, _>("ellipeinc", "d_d-d", |x: &[f64]| ellipeinc(x[0], x[1]));
     }
 
     // ellipj
 
     #[test]
     fn test_ellipj() {
-        xsref::test::<(f64, f64, f64, f64), _>("ellipj", "d_d-d_d_d_d", |x: &[f64]| {
+        testing::test::<(f64, f64, f64, f64), _>("ellipj", "d_d-d_d_d_d", |x: &[f64]| {
             ellipj(x[0], x[1])
         });
     }
diff --git a/src/erf.rs b/src/erf.rs
index 5f60c82..ca9dc4e 100644
--- a/src/erf.rs
+++ b/src/erf.rs
@@ -68,81 +68,81 @@ pub fn voigt_profile(x: f64, sigma: f64, gamma: f64) -> f64 {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     // erf
 
     #[test]
     fn test_erf_f64() {
-        xsref::test::<f64, _>("erf", "d-d", |x: &[f64]| erf(x[0]));
+        testing::test::<f64, _>("erf", "d-d", |x: &[f64]| erf(x[0]));
     }
 
     #[test]
     fn test_erf_c64() {
-        xsref::test::<Complex<f64>, _>("erf", "cd-cd", |x: &[f64]| erf(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("erf", "cd-cd", |x: &[f64]| erf(c64(x[0], x[1])));
     }
 
     // erfc
 
     #[test]
     fn test_erfc_f64() {
-        xsref::test::<f64, _>("erfc", "d-d", |x: &[f64]| erfc(x[0]));
+        testing::test::<f64, _>("erfc", "d-d", |x: &[f64]| erfc(x[0]));
     }
 
     #[test]
     fn test_erfc_c64() {
-        xsref::test::<Complex<f64>, _>("erfc", "cd-cd", |x: &[f64]| erfc(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("erfc", "cd-cd", |x: &[f64]| erfc(c64(x[0], x[1])));
     }
 
     // erfcx
 
     #[test]
     fn test_erfcx_f64() {
-        xsref::test::<f64, _>("erfcx", "d-d", |x: &[f64]| erfcx(x[0]));
+        testing::test::<f64, _>("erfcx", "d-d", |x: &[f64]| erfcx(x[0]));
     }
 
     #[test]
     fn test_erfcx_c64() {
-        xsref::test::<Complex<f64>, _>("erfcx", "cd-cd", |x: &[f64]| erfcx(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("erfcx", "cd-cd", |x: &[f64]| erfcx(c64(x[0], x[1])));
     }
 
     // erfi
 
     #[test]
     fn test_erfi_f64() {
-        xsref::test::<f64, _>("erfi", "d-d", |x: &[f64]| erfi(x[0]));
+        testing::test::<f64, _>("erfi", "d-d", |x: &[f64]| erfi(x[0]));
     }
 
     #[test]
     fn test_erfi_c64() {
-        xsref::test::<Complex<f64>, _>("erfi", "cd-cd", |x: &[f64]| erfi(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("erfi", "cd-cd", |x: &[f64]| erfi(c64(x[0], x[1])));
     }
 
     // dawsn
 
     #[test]
     fn test_dawsn_f64() {
-        xsref::test::<f64, _>("dawsn", "d-d", |x: &[f64]| dawsn(x[0]));
+        testing::test::<f64, _>("dawsn", "d-d", |x: &[f64]| dawsn(x[0]));
     }
 
     #[test]
     fn test_dawsn_c64() {
-        xsref::test::<Complex<f64>, _>("dawsn", "cd-cd", |x: &[f64]| dawsn(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("dawsn", "cd-cd", |x: &[f64]| dawsn(c64(x[0], x[1])));
     }
 
     // wofz
 
     #[test]
     fn test_wofz() {
-        xsref::test::<Complex<f64>, _>("wofz", "cd-cd", |x: &[f64]| wofz(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("wofz", "cd-cd", |x: &[f64]| wofz(c64(x[0], x[1])));
     }
 
     // voigt_profile
 
     #[test]
     fn test_voigt_profile() {
-        xsref::test::<f64, _>("voigt_profile", "d_d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("voigt_profile", "d_d_d-d", |x: &[f64]| {
             voigt_profile(x[0], x[1], x[2])
         });
     }
diff --git a/src/exp.rs b/src/exp.rs
index b2166df..092e1f9 100644
--- a/src/exp.rs
+++ b/src/exp.rs
@@ -48,26 +48,26 @@ pub fn exp10(x: f64) -> f64 {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_expm1_f64() {
-        xsref::test::<f64, _>("expm1", "d-d", |x: &[f64]| expm1(x[0]));
+        testing::test::<f64, _>("expm1", "d-d", |x: &[f64]| expm1(x[0]));
     }
 
     #[test]
     fn test_expm1_c64() {
-        xsref::test::<Complex<f64>, _>("expm1", "cd-cd", |x: &[f64]| expm1(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("expm1", "cd-cd", |x: &[f64]| expm1(c64(x[0], x[1])));
     }
 
     #[test]
     fn test_exp2_f64() {
-        xsref::test::<f64, _>("exp2", "d-d", |x: &[f64]| exp2(x[0]));
+        testing::test::<f64, _>("exp2", "d-d", |x: &[f64]| exp2(x[0]));
     }
 
     #[test]
     fn test_exp10_f64() {
-        xsref::test::<f64, _>("exp10", "d-d", |x: &[f64]| exp10(x[0]));
+        testing::test::<f64, _>("exp10", "d-d", |x: &[f64]| exp10(x[0]));
     }
 }
diff --git a/src/expint.rs b/src/expint.rs
index 53d46f3..23c171f 100644
--- a/src/expint.rs
+++ b/src/expint.rs
@@ -57,37 +57,37 @@ pub fn scaled_exp1(x: f64) -> f64 {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     // expi
 
     #[test]
     fn test_expi_f64() {
-        xsref::test::<f64, _>("expi", "d-d", |x: &[f64]| expi(x[0]));
+        testing::test::<f64, _>("expi", "d-d", |x: &[f64]| expi(x[0]));
     }
 
     #[test]
     fn test_expi_c64() {
-        xsref::test::<Complex<f64>, _>("expi", "cd-cd", |x: &[f64]| expi(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("expi", "cd-cd", |x: &[f64]| expi(c64(x[0], x[1])));
     }
 
     // exp1
 
     #[test]
     fn test_exp1_f64() {
-        xsref::test::<f64, _>("exp1", "d-d", |x: &[f64]| exp1(x[0]));
+        testing::test::<f64, _>("exp1", "d-d", |x: &[f64]| exp1(x[0]));
     }
 
     #[test]
     fn test_exp1_c64() {
-        xsref::test::<Complex<f64>, _>("exp1", "cd-cd", |x: &[f64]| exp1(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("exp1", "cd-cd", |x: &[f64]| exp1(c64(x[0], x[1])));
     }
 
     // scaled_exp1
 
     #[test]
     fn test_scaled_exp1_f64() {
-        xsref::test::<f64, _>("scaled_exp1", "d-d", |x: &[f64]| scaled_exp1(x[0]));
+        testing::test::<f64, _>("scaled_exp1", "d-d", |x: &[f64]| scaled_exp1(x[0]));
     }
 }
diff --git a/src/fresnel.rs b/src/fresnel.rs
index 59e1834..77c4007 100644
--- a/src/fresnel.rs
+++ b/src/fresnel.rs
@@ -130,19 +130,19 @@ pub fn modified_fresnel_minus(x: f64) -> (Complex<f64>, Complex<f64>) {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     // fresnel
 
     #[test]
     fn test_fresnel_f64() {
-        xsref::test::<(f64, f64), _>("fresnel", "d-d_d", |x: &[f64]| fresnel(x[0]));
+        testing::test::<(f64, f64), _>("fresnel", "d-d_d", |x: &[f64]| fresnel(x[0]));
     }
 
     #[test]
     fn test_fresnel_c64() {
-        xsref::test::<(Complex<f64>, Complex<f64>), _>("fresnel", "cd-cd_cd", |x: &[f64]| {
+        testing::test::<(Complex<f64>, Complex<f64>), _>("fresnel", "cd-cd_cd", |x: &[f64]| {
             fresnel(c64(x[0], x[1]))
         });
     }
@@ -151,7 +151,7 @@ mod tests {
 
     #[test]
     fn test_modified_fresnel_plus_c64() {
-        xsref::test::<(Complex<f64>, Complex<f64>), _>(
+        testing::test::<(Complex<f64>, Complex<f64>), _>(
             "modified_fresnel_plus",
             "d-cd_cd",
             |x: &[f64]| modified_fresnel_plus(x[0]),
@@ -162,7 +162,7 @@ mod tests {
 
     #[test]
     fn test_modified_fresnel_minus_c64() {
-        xsref::test::<(Complex<f64>, Complex<f64>), _>(
+        testing::test::<(Complex<f64>, Complex<f64>), _>(
             "modified_fresnel_minus",
             "d-cd_cd",
             |x: &[f64]| modified_fresnel_minus(x[0]),
diff --git a/src/gamma.rs b/src/gamma.rs
index 53b8c41..77de084 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -44,48 +44,48 @@ xsf_impl!(gammasgn, (x: f64), "Sign of the Gamma function");
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_gamma_f64() {
-        xsref::test::<f64, _>("gamma", "d-d", |x: &[f64]| gamma(x[0]));
+        testing::test::<f64, _>("gamma", "d-d", |x: &[f64]| gamma(x[0]));
     }
 
     #[test]
     fn test_gamma_c64() {
-        xsref::test::<Complex<f64>, _>("gamma", "cd-cd", |x: &[f64]| gamma(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("gamma", "cd-cd", |x: &[f64]| gamma(c64(x[0], x[1])));
     }
 
     #[test]
     fn test_gammainc() {
-        xsref::test::<f64, _>("gammainc", "d_d-d", |x: &[f64]| gammainc(x[0], x[1]));
+        testing::test::<f64, _>("gammainc", "d_d-d", |x: &[f64]| gammainc(x[0], x[1]));
     }
 
     #[test]
     fn test_gammaincc() {
-        xsref::test::<f64, _>("gammaincc", "d_d-d", |x: &[f64]| gammaincc(x[0], x[1]));
+        testing::test::<f64, _>("gammaincc", "d_d-d", |x: &[f64]| gammaincc(x[0], x[1]));
     }
 
     #[test]
     fn test_gammaincinv() {
-        xsref::test::<f64, _>("gammaincinv", "d_d-d", |x: &[f64]| gammaincinv(x[0], x[1]));
+        testing::test::<f64, _>("gammaincinv", "d_d-d", |x: &[f64]| gammaincinv(x[0], x[1]));
     }
 
     #[test]
     fn test_gammainccinv() {
-        xsref::test::<f64, _>("gammainccinv", "d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("gammainccinv", "d_d-d", |x: &[f64]| {
             gammainccinv(x[0], x[1])
         });
     }
 
     #[test]
     fn test_gammaln() {
-        xsref::test::<f64, _>("gammaln", "d-d", |x: &[f64]| gammaln(x[0]));
+        testing::test::<f64, _>("gammaln", "d-d", |x: &[f64]| gammaln(x[0]));
     }
 
     #[test]
     fn test_gammasgn() {
-        xsref::test::<f64, _>("gammasgn", "d-d", |x: &[f64]| gammasgn(x[0]));
+        testing::test::<f64, _>("gammasgn", "d-d", |x: &[f64]| gammasgn(x[0]));
     }
 }
diff --git a/src/hyp2f1.rs b/src/hyp2f1.rs
index 319d20c..030c510 100644
--- a/src/hyp2f1.rs
+++ b/src/hyp2f1.rs
@@ -36,19 +36,19 @@ pub fn hyp2f1<T: Hyp2F1Arg>(a: f64, b: f64, c: f64, z: T) -> T::Output {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_hyp2f1_f64() {
-        xsref::test::<f64, _>("hyp2f1", "d_d_d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("hyp2f1", "d_d_d_d-d", |x: &[f64]| {
             hyp2f1(x[0], x[1], x[2], x[3])
         });
     }
 
     #[test]
     fn test_hyp2f1_c64() {
-        xsref::test::<Complex<f64>, _>("hyp2f1", "d_d_d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("hyp2f1", "d_d_d_cd-cd", |x: &[f64]| {
             hyp2f1(x[0], x[1], x[2], c64(x[3], x[4]))
         });
     }
diff --git a/src/iv_ratio.rs b/src/iv_ratio.rs
index 2968fb6..d86c84e 100644
--- a/src/iv_ratio.rs
+++ b/src/iv_ratio.rs
@@ -14,15 +14,15 @@ xsf_impl!(
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_iv_ratio() {
-        xsref::test::<f64, _>("iv_ratio", "d_d-d", |x: &[f64]| iv_ratio(x[0], x[1]));
+        testing::test::<f64, _>("iv_ratio", "d_d-d", |x: &[f64]| iv_ratio(x[0], x[1]));
     }
 
     #[test]
     fn test_iv_ratio_c() {
-        xsref::test::<f64, _>("iv_ratio_c", "d_d-d", |x: &[f64]| iv_ratio_c(x[0], x[1]));
+        testing::test::<f64, _>("iv_ratio_c", "d_d-d", |x: &[f64]| iv_ratio_c(x[0], x[1]));
     }
 }
diff --git a/src/kelvin.rs b/src/kelvin.rs
index a685b37..9f6aeff 100644
--- a/src/kelvin.rs
+++ b/src/kelvin.rs
@@ -36,52 +36,52 @@ pub fn kelvin(x: f64) -> (Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::Complex;
 
     #[test]
     fn test_ber_f64() {
-        xsref::test::<f64, _>("ber", "d-d", |x: &[f64]| ber(x[0]));
+        testing::test::<f64, _>("ber", "d-d", |x: &[f64]| ber(x[0]));
     }
 
     #[test]
     fn test_bei_f64() {
-        xsref::test::<f64, _>("bei", "d-d", |x: &[f64]| bei(x[0]));
+        testing::test::<f64, _>("bei", "d-d", |x: &[f64]| bei(x[0]));
     }
 
     #[test]
     fn test_ker_f64() {
-        xsref::test::<f64, _>("ker", "d-d", |x: &[f64]| ker(x[0]));
+        testing::test::<f64, _>("ker", "d-d", |x: &[f64]| ker(x[0]));
     }
 
     #[test]
     fn test_kei_f64() {
-        xsref::test::<f64, _>("kei", "d-d", |x: &[f64]| kei(x[0]));
+        testing::test::<f64, _>("kei", "d-d", |x: &[f64]| kei(x[0]));
     }
 
     #[test]
     fn test_berp_f64() {
-        xsref::test::<f64, _>("berp", "d-d", |x: &[f64]| berp(x[0]));
+        testing::test::<f64, _>("berp", "d-d", |x: &[f64]| berp(x[0]));
     }
 
     #[test]
     fn test_beip_f64() {
-        xsref::test::<f64, _>("beip", "d-d", |x: &[f64]| beip(x[0]));
+        testing::test::<f64, _>("beip", "d-d", |x: &[f64]| beip(x[0]));
     }
 
     #[test]
     fn test_kerp_f64() {
-        xsref::test::<f64, _>("kerp", "d-d", |x: &[f64]| kerp(x[0]));
+        testing::test::<f64, _>("kerp", "d-d", |x: &[f64]| kerp(x[0]));
     }
 
     #[test]
     fn test_keip_f64() {
-        xsref::test::<f64, _>("keip", "d-d", |x: &[f64]| keip(x[0]));
+        testing::test::<f64, _>("keip", "d-d", |x: &[f64]| keip(x[0]));
     }
 
     #[test]
     fn test_kelvin_c64() {
-        xsref::test::<(Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>), _>(
+        testing::test::<(Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>), _>(
             "kelvin",
             "d-cd_cd_cd_cd",
             |x: &[f64]| kelvin(x[0]),
diff --git a/src/lambertw.rs b/src/lambertw.rs
index f3eb95b..0e682eb 100644
--- a/src/lambertw.rs
+++ b/src/lambertw.rs
@@ -16,12 +16,12 @@ pub fn lambertw(z: Complex<f64>, k: isize, tol: f64) -> Complex<f64> {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_lambertw_c64() {
-        xsref::test::<Complex<f64>, _>("lambertw", "cd_p_d-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("lambertw", "cd_p_d-cd", |x: &[f64]| {
             lambertw(c64(x[0], x[1]), x[2] as isize, x[3])
         });
     }
diff --git a/src/lib.rs b/src/lib.rs
index 88add4d..00dbf06 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,9 +1,7 @@
 mod bindings;
 
 #[cfg(test)]
-mod test_utils;
-#[cfg(test)]
-use test_utils::xsref;
+mod testing;
 
 // airy.h
 mod airy;
diff --git a/src/log.rs b/src/log.rs
index b676b5f..9898a0b 100644
--- a/src/log.rs
+++ b/src/log.rs
@@ -70,44 +70,44 @@ where
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_log1p_f64() {
-        xsref::test::<f64, _>("log1p", "d-d", |x: &[f64]| log1p(x[0]));
+        testing::test::<f64, _>("log1p", "d-d", |x: &[f64]| log1p(x[0]));
     }
 
     #[test]
     fn test_log1p_c64() {
-        xsref::test::<Complex<f64>, _>("log1p", "cd-cd", |x: &[f64]| log1p(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("log1p", "cd-cd", |x: &[f64]| log1p(c64(x[0], x[1])));
     }
 
     #[test]
     fn test_log1pmx_f64() {
-        xsref::test::<f64, _>("log1pmx", "d-d", |x: &[f64]| log1pmx(x[0]));
+        testing::test::<f64, _>("log1pmx", "d-d", |x: &[f64]| log1pmx(x[0]));
     }
 
     #[test]
     fn test_xlogy_f64() {
-        xsref::test::<f64, _>("xlogy", "d_d-d", |x: &[f64]| xlogy(x[0], x[1]));
+        testing::test::<f64, _>("xlogy", "d_d-d", |x: &[f64]| xlogy(x[0], x[1]));
     }
 
     #[test]
     fn test_xlogy_c64() {
-        xsref::test::<Complex<f64>, _>("xlogy", "cd_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("xlogy", "cd_cd-cd", |x: &[f64]| {
             xlogy(c64(x[0], x[1]), c64(x[2], x[3]))
         });
     }
 
     #[test]
     fn test_xlog1py_f64() {
-        xsref::test::<f64, _>("xlog1py", "d_d-d", |x: &[f64]| xlog1py(x[0], x[1]));
+        testing::test::<f64, _>("xlog1py", "d_d-d", |x: &[f64]| xlog1py(x[0], x[1]));
     }
 
     #[test]
     fn test_xlog1py_c64() {
-        xsref::test::<Complex<f64>, _>("xlog1py", "cd_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("xlog1py", "cd_cd-cd", |x: &[f64]| {
             xlog1py(c64(x[0], x[1]), c64(x[2], x[3]))
         });
     }
diff --git a/src/log_exp.rs b/src/log_exp.rs
index bb787e1..6a9acbe 100644
--- a/src/log_exp.rs
+++ b/src/log_exp.rs
@@ -10,25 +10,25 @@ xsf_impl!(log1mexp, (x: f64), "Compute `log(1 - exp(x))`");
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_expit() {
-        xsref::test::<f64, _>("expit", "d-d", |x: &[f64]| expit(x[0]));
+        testing::test::<f64, _>("expit", "d-d", |x: &[f64]| expit(x[0]));
     }
 
     #[test]
     fn test_exprel() {
-        xsref::test::<f64, _>("exprel", "d-d", |x: &[f64]| exprel(x[0]));
+        testing::test::<f64, _>("exprel", "d-d", |x: &[f64]| exprel(x[0]));
     }
 
     #[test]
     fn test_logit() {
-        xsref::test::<f64, _>("logit", "d-d", |x: &[f64]| logit(x[0]));
+        testing::test::<f64, _>("logit", "d-d", |x: &[f64]| logit(x[0]));
     }
 
     #[test]
     fn test_log_expit() {
-        xsref::test::<f64, _>("log_expit", "d-d", |x: &[f64]| log_expit(x[0]));
+        testing::test::<f64, _>("log_expit", "d-d", |x: &[f64]| log_expit(x[0]));
     }
 }
diff --git a/src/loggamma.rs b/src/loggamma.rs
index 9cd4914..61eab80 100644
--- a/src/loggamma.rs
+++ b/src/loggamma.rs
@@ -48,30 +48,32 @@ pub fn rgamma<T: LogGammaArg>(z: T) -> T::Output {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     // loggamma
 
     #[test]
     fn test_loggamma_f64() {
-        xsref::test::<f64, _>("loggamma", "d-d", |x: &[f64]| loggamma(x[0]));
+        testing::test::<f64, _>("loggamma", "d-d", |x: &[f64]| loggamma(x[0]));
     }
 
     #[test]
     fn test_loggamma_c64() {
-        xsref::test::<Complex<f64>, _>("loggamma", "cd-cd", |x: &[f64]| loggamma(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("loggamma", "cd-cd", |x: &[f64]| {
+            loggamma(c64(x[0], x[1]))
+        });
     }
 
     // rgamma
 
     #[test]
     fn test_rgamma_f64() {
-        xsref::test::<f64, _>("rgamma", "d-d", |x: &[f64]| rgamma(x[0]));
+        testing::test::<f64, _>("rgamma", "d-d", |x: &[f64]| rgamma(x[0]));
     }
 
     #[test]
     fn test_rgamma_c64() {
-        xsref::test::<Complex<f64>, _>("rgamma", "cd-cd", |x: &[f64]| rgamma(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("rgamma", "cd-cd", |x: &[f64]| rgamma(c64(x[0], x[1])));
     }
 }
diff --git a/src/mathieu.rs b/src/mathieu.rs
index a8905c1..c5ad9c6 100644
--- a/src/mathieu.rs
+++ b/src/mathieu.rs
@@ -139,45 +139,45 @@ pub fn msm2(m: f64, q: f64, x: f64) -> (f64, f64) {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_cem_cva() {
-        xsref::test::<f64, _>("cem_cva", "d_d-d", |x: &[f64]| cem_cva(x[0], x[1]));
+        testing::test::<f64, _>("cem_cva", "d_d-d", |x: &[f64]| cem_cva(x[0], x[1]));
     }
 
     #[test]
     fn test_sem_cva() {
-        xsref::test::<f64, _>("sem_cva", "d_d-d", |x: &[f64]| sem_cva(x[0], x[1]));
+        testing::test::<f64, _>("sem_cva", "d_d-d", |x: &[f64]| sem_cva(x[0], x[1]));
     }
 
     #[test]
     fn test_cem() {
-        xsref::test::<(f64, f64), _>("cem", "d_d_d-d_d", |x: &[f64]| cem(x[0], x[1], x[2]));
+        testing::test::<(f64, f64), _>("cem", "d_d_d-d_d", |x: &[f64]| cem(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_sem() {
-        xsref::test::<(f64, f64), _>("sem", "d_d_d-d_d", |x: &[f64]| sem(x[0], x[1], x[2]));
+        testing::test::<(f64, f64), _>("sem", "d_d_d-d_d", |x: &[f64]| sem(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_mcm1() {
-        xsref::test::<(f64, f64), _>("mcm1", "d_d_d-d_d", |x: &[f64]| mcm1(x[0], x[1], x[2]));
+        testing::test::<(f64, f64), _>("mcm1", "d_d_d-d_d", |x: &[f64]| mcm1(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_msm1() {
-        xsref::test::<(f64, f64), _>("msm1", "d_d_d-d_d", |x: &[f64]| msm1(x[0], x[1], x[2]));
+        testing::test::<(f64, f64), _>("msm1", "d_d_d-d_d", |x: &[f64]| msm1(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_mcm2() {
-        xsref::test::<(f64, f64), _>("mcm2", "d_d_d-d_d", |x: &[f64]| mcm2(x[0], x[1], x[2]));
+        testing::test::<(f64, f64), _>("mcm2", "d_d_d-d_d", |x: &[f64]| mcm2(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_msm2() {
-        xsref::test::<(f64, f64), _>("msm2", "d_d_d-d_d", |x: &[f64]| msm2(x[0], x[1], x[2]));
+        testing::test::<(f64, f64), _>("msm2", "d_d_d-d_d", |x: &[f64]| msm2(x[0], x[1], x[2]));
     }
 }
diff --git a/src/par_cyl.rs b/src/par_cyl.rs
index 325e6bb..ebc2baa 100644
--- a/src/par_cyl.rs
+++ b/src/par_cyl.rs
@@ -73,20 +73,20 @@ pub fn pbvv(v: f64, x: f64) -> (f64, f64) {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_pbwa() {
-        xsref::test::<(f64, f64), _>("pbwa", "d_d-d_d", |x: &[f64]| pbwa(x[0], x[1]));
+        testing::test::<(f64, f64), _>("pbwa", "d_d-d_d", |x: &[f64]| pbwa(x[0], x[1]));
     }
 
     #[test]
     fn test_pbdv() {
-        xsref::test::<(f64, f64), _>("pbdv", "d_d-d_d", |x: &[f64]| pbdv(x[0], x[1]));
+        testing::test::<(f64, f64), _>("pbdv", "d_d-d_d", |x: &[f64]| pbdv(x[0], x[1]));
     }
 
     #[test]
     fn test_pbvv() {
-        xsref::test::<(f64, f64), _>("pbvv", "d_d-d_d", |x: &[f64]| pbvv(x[0], x[1]));
+        testing::test::<(f64, f64), _>("pbvv", "d_d-d_d", |x: &[f64]| pbvv(x[0], x[1]));
     }
 }
diff --git a/src/sici.rs b/src/sici.rs
index fb5d26b..246fe29 100644
--- a/src/sici.rs
+++ b/src/sici.rs
@@ -119,19 +119,19 @@ pub fn shichi<T: SiciArg>(z: T) -> (T::Output, T::Output) {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     // sici
 
     #[test]
     fn test_sici_f64() {
-        xsref::test::<(f64, f64), _>("sici", "d-d_d", |x: &[f64]| sici(x[0]));
+        testing::test::<(f64, f64), _>("sici", "d-d_d", |x: &[f64]| sici(x[0]));
     }
 
     #[test]
     fn test_sici_c64() {
-        xsref::test::<(Complex<f64>, Complex<f64>), _>("sici", "cd-cd_cd", |x: &[f64]| {
+        testing::test::<(Complex<f64>, Complex<f64>), _>("sici", "cd-cd_cd", |x: &[f64]| {
             sici(c64(x[0], x[1]))
         });
     }
@@ -140,12 +140,12 @@ mod tests {
 
     #[test]
     fn test_shichi_f64() {
-        xsref::test::<(f64, f64), _>("shichi", "d-d_d", |x: &[f64]| shichi(x[0]));
+        testing::test::<(f64, f64), _>("shichi", "d-d_d", |x: &[f64]| shichi(x[0]));
     }
 
     #[test]
     fn test_shichi_c64() {
-        xsref::test::<(Complex<f64>, Complex<f64>), _>("shichi", "cd-cd_cd", |x: &[f64]| {
+        testing::test::<(Complex<f64>, Complex<f64>), _>("shichi", "cd-cd_cd", |x: &[f64]| {
             shichi(c64(x[0], x[1]))
         });
     }
diff --git a/src/specfun.rs b/src/specfun.rs
index c30d669..fad260b 100644
--- a/src/specfun.rs
+++ b/src/specfun.rs
@@ -47,24 +47,24 @@ pub fn pmv(m: i64, v: f64, x: f64) -> f64 {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_hypu() {
         // the table is called "hyperu" instead of "hypu"
-        xsref::test::<f64, _>("hyperu", "d_d_d-d", |x: &[f64]| hypu(x[0], x[1], x[2]));
+        testing::test::<f64, _>("hyperu", "d_d_d-d", |x: &[f64]| hypu(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_hyp1f1() {
-        xsref::test::<Complex<f64>, _>("hyp1f1", "d_d_cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("hyp1f1", "d_d_cd-cd", |x: &[f64]| {
             hyp1f1(x[0], x[1], c64(x[2], x[3]))
         });
     }
 
     #[test]
     fn test_pmv() {
-        xsref::test::<f64, _>("pmv", "d_d_d-d", |x: &[f64]| pmv(x[0] as i64, x[1], x[2]));
+        testing::test::<f64, _>("pmv", "d_d_d-d", |x: &[f64]| pmv(x[0] as i64, x[1], x[2]));
     }
 }
diff --git a/src/sphd_wave.rs b/src/sphd_wave.rs
index e81d46b..81c9abd 100644
--- a/src/sphd_wave.rs
+++ b/src/sphd_wave.rs
@@ -295,95 +295,103 @@ pub fn oblate_radial2(m: u64, n: u64, c: f64, cv: f64, x: f64) -> (f64, f64) {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_prolate_segv() {
-        xsref::test::<f64, _>("prolate_segv", "d_d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("prolate_segv", "d_d_d-d", |x: &[f64]| {
             prolate_segv(x[0] as u64, x[1] as u64, x[2])
         });
     }
 
+    // https://github.com/scipy/xsref/issues/9
+    // #[test]
+    // fn test_oblate_segv() {
+    //     testing::test::<f64, _>("oblate_segv", "d_d_d-d", |x: &[f64]| {
+    //         oblate_segv(x[0] as u64, x[1] as u64, x[2])
+    //     });
+    // }
+
     #[test]
     fn test_prolate_aswfa_nocv() {
-        xsref::test::<(f64, f64), _>("prolate_aswfa_nocv", "d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("prolate_aswfa_nocv", "d_d_d_d-d_d", |x: &[f64]| {
             prolate_aswfa_nocv(x[0] as u64, x[1] as u64, x[2], x[3])
         });
     }
 
     #[test]
     fn test_oblate_aswfa_nocv() {
-        xsref::test::<(f64, f64), _>("oblate_aswfa_nocv", "d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("oblate_aswfa_nocv", "d_d_d_d-d_d", |x: &[f64]| {
             oblate_aswfa_nocv(x[0] as u64, x[1] as u64, x[2], x[3])
         });
     }
 
     #[test]
     fn test_prolate_radial1_nocv() {
-        xsref::test::<(f64, f64), _>("prolate_radial1_nocv", "d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("prolate_radial1_nocv", "d_d_d_d-d_d", |x: &[f64]| {
             prolate_radial1_nocv(x[0] as u64, x[1] as u64, x[2], x[3])
         });
     }
 
     #[test]
     fn test_oblate_radial1_nocv() {
-        xsref::test::<(f64, f64), _>("oblate_radial1_nocv", "d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("oblate_radial1_nocv", "d_d_d_d-d_d", |x: &[f64]| {
             oblate_radial1_nocv(x[0] as u64, x[1] as u64, x[2], x[3])
         });
     }
 
     #[test]
     fn test_prolate_radial2_nocv() {
-        xsref::test::<(f64, f64), _>("prolate_radial2_nocv", "d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("prolate_radial2_nocv", "d_d_d_d-d_d", |x: &[f64]| {
             prolate_radial2_nocv(x[0] as u64, x[1] as u64, x[2], x[3])
         });
     }
 
     #[test]
     fn test_oblate_radial2_nocv() {
-        xsref::test::<(f64, f64), _>("oblate_radial2_nocv", "d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("oblate_radial2_nocv", "d_d_d_d-d_d", |x: &[f64]| {
             oblate_radial2_nocv(x[0] as u64, x[1] as u64, x[2], x[3])
         });
     }
 
     #[test]
     fn test_prolate_aswfa() {
-        xsref::test::<(f64, f64), _>("prolate_aswfa", "d_d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("prolate_aswfa", "d_d_d_d_d-d_d", |x: &[f64]| {
             prolate_aswfa(x[0] as u64, x[1] as u64, x[2], x[3], x[4])
         });
     }
 
     #[test]
     fn test_oblate_aswfa() {
-        xsref::test::<(f64, f64), _>("oblate_aswfa", "d_d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("oblate_aswfa", "d_d_d_d_d-d_d", |x: &[f64]| {
             oblate_aswfa(x[0] as u64, x[1] as u64, x[2], x[3], x[4])
         });
     }
 
     #[test]
     fn test_prolate_radial1() {
-        xsref::test::<(f64, f64), _>("prolate_radial1", "d_d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("prolate_radial1", "d_d_d_d_d-d_d", |x: &[f64]| {
             prolate_radial1(x[0] as u64, x[1] as u64, x[2], x[3], x[4])
         });
     }
 
     #[test]
     fn test_oblate_radial1() {
-        xsref::test::<(f64, f64), _>("oblate_radial1", "d_d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("oblate_radial1", "d_d_d_d_d-d_d", |x: &[f64]| {
             oblate_radial1(x[0] as u64, x[1] as u64, x[2], x[3], x[4])
         });
     }
 
     #[test]
     fn test_prolate_radial2() {
-        xsref::test::<(f64, f64), _>("prolate_radial2", "d_d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("prolate_radial2", "d_d_d_d_d-d_d", |x: &[f64]| {
             prolate_radial2(x[0] as u64, x[1] as u64, x[2], x[3], x[4])
         });
     }
 
     #[test]
     fn test_oblate_radial2() {
-        xsref::test::<(f64, f64), _>("oblate_radial2", "d_d_d_d_d-d_d", |x: &[f64]| {
+        testing::test::<(f64, f64), _>("oblate_radial2", "d_d_d_d_d-d_d", |x: &[f64]| {
             oblate_radial2(x[0] as u64, x[1] as u64, x[2], x[3], x[4])
         });
     }
diff --git a/src/stats.rs b/src/stats.rs
index d4b2df8..967c6d0 100644
--- a/src/stats.rs
+++ b/src/stats.rs
@@ -91,167 +91,167 @@ xsf_impl!(nbdtri, (k: c_int, n: c_int, p: f64), "Negative binomial quantile func
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_ndtr_f64() {
-        xsref::test::<f64, _>("ndtr", "d-d", |x: &[f64]| ndtr(x[0]));
+        testing::test::<f64, _>("ndtr", "d-d", |x: &[f64]| ndtr(x[0]));
     }
 
     // TODO: no log_ndtr xsref tables => needs manual smoketests
 
     #[test]
     fn test_ndtr_c64() {
-        xsref::test::<Complex<f64>, _>("ndtr", "cd-cd", |x: &[f64]| ndtr(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("ndtr", "cd-cd", |x: &[f64]| ndtr(c64(x[0], x[1])));
     }
 
     #[test]
     fn test_ndtri() {
-        xsref::test::<f64, _>("ndtri", "d-d", |x: &[f64]| ndtri(x[0]));
+        testing::test::<f64, _>("ndtri", "d-d", |x: &[f64]| ndtri(x[0]));
     }
 
     #[test]
     fn test_kolmogorov() {
-        xsref::test::<f64, _>("kolmogorov", "d-d", |x: &[f64]| kolmogorov(x[0]));
+        testing::test::<f64, _>("kolmogorov", "d-d", |x: &[f64]| kolmogorov(x[0]));
     }
 
     #[test]
     fn test_kolmogc() {
-        xsref::test::<f64, _>("kolmogc", "d-d", |x: &[f64]| kolmogc(x[0]));
+        testing::test::<f64, _>("kolmogc", "d-d", |x: &[f64]| kolmogc(x[0]));
     }
 
     #[test]
     fn test_kolmogi() {
-        xsref::test::<f64, _>("kolmogi", "d-d", |x: &[f64]| kolmogi(x[0]));
+        testing::test::<f64, _>("kolmogi", "d-d", |x: &[f64]| kolmogi(x[0]));
     }
 
     #[test]
     fn test_kolmogci() {
-        xsref::test::<f64, _>("kolmogci", "d-d", |x: &[f64]| kolmogci(x[0]));
+        testing::test::<f64, _>("kolmogci", "d-d", |x: &[f64]| kolmogci(x[0]));
     }
 
     #[test]
     fn test_kolmogp() {
-        xsref::test::<f64, _>("kolmogp", "d-d", |x: &[f64]| kolmogp(x[0]));
+        testing::test::<f64, _>("kolmogp", "d-d", |x: &[f64]| kolmogp(x[0]));
     }
 
     #[test]
     fn test_smirnov() {
-        xsref::test::<f64, _>("smirnov", "p_d-d", |x: &[f64]| smirnov(x[0] as i32, x[1]));
+        testing::test::<f64, _>("smirnov", "p_d-d", |x: &[f64]| smirnov(x[0] as i32, x[1]));
     }
 
     #[test]
     fn test_smirnovc() {
-        xsref::test::<f64, _>("smirnovc", "p_d-d", |x: &[f64]| smirnovc(x[0] as i32, x[1]));
+        testing::test::<f64, _>("smirnovc", "p_d-d", |x: &[f64]| smirnovc(x[0] as i32, x[1]));
     }
 
     #[test]
     fn test_smirnovi() {
-        xsref::test::<f64, _>("smirnovi", "p_d-d", |x: &[f64]| smirnovi(x[0] as i32, x[1]));
+        testing::test::<f64, _>("smirnovi", "p_d-d", |x: &[f64]| smirnovi(x[0] as i32, x[1]));
     }
 
     #[test]
     fn test_smirnovci() {
-        xsref::test::<f64, _>("smirnovci", "p_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("smirnovci", "p_d-d", |x: &[f64]| {
             smirnovci(x[0] as i32, x[1])
         });
     }
 
     #[test]
     fn test_smirnovp() {
-        xsref::test::<f64, _>("smirnovp", "p_d-d", |x: &[f64]| smirnovp(x[0] as i32, x[1]));
+        testing::test::<f64, _>("smirnovp", "p_d-d", |x: &[f64]| smirnovp(x[0] as i32, x[1]));
     }
 
     #[test]
     fn test_owens_t() {
-        xsref::test::<f64, _>("owens_t", "d_d-d", |x: &[f64]| owens_t(x[0], x[1]));
+        testing::test::<f64, _>("owens_t", "d_d-d", |x: &[f64]| owens_t(x[0], x[1]));
     }
 
     #[test]
     fn test_chdtr() {
-        xsref::test::<f64, _>("chdtr", "d_d-d", |x: &[f64]| chdtr(x[0], x[1]));
+        testing::test::<f64, _>("chdtr", "d_d-d", |x: &[f64]| chdtr(x[0], x[1]));
     }
 
     #[test]
     fn test_chdtrc() {
-        xsref::test::<f64, _>("chdtrc", "d_d-d", |x: &[f64]| chdtrc(x[0], x[1]));
+        testing::test::<f64, _>("chdtrc", "d_d-d", |x: &[f64]| chdtrc(x[0], x[1]));
     }
 
     #[test]
     fn test_chdtri() {
-        xsref::test::<f64, _>("chdtri", "d_d-d", |x: &[f64]| chdtri(x[0], x[1]));
+        testing::test::<f64, _>("chdtri", "d_d-d", |x: &[f64]| chdtri(x[0], x[1]));
     }
 
     #[test]
     fn test_fdtr() {
-        xsref::test::<f64, _>("fdtr", "d_d_d-d", |x: &[f64]| fdtr(x[0], x[1], x[2]));
+        testing::test::<f64, _>("fdtr", "d_d_d-d", |x: &[f64]| fdtr(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_fdtrc() {
-        xsref::test::<f64, _>("fdtrc", "d_d_d-d", |x: &[f64]| fdtrc(x[0], x[1], x[2]));
+        testing::test::<f64, _>("fdtrc", "d_d_d-d", |x: &[f64]| fdtrc(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_fdtri() {
-        xsref::test::<f64, _>("fdtri", "d_d_d-d", |x: &[f64]| fdtri(x[0], x[1], x[2]));
+        testing::test::<f64, _>("fdtri", "d_d_d-d", |x: &[f64]| fdtri(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_gdtr() {
-        xsref::test::<f64, _>("gdtr", "d_d_d-d", |x: &[f64]| gdtr(x[0], x[1], x[2]));
+        testing::test::<f64, _>("gdtr", "d_d_d-d", |x: &[f64]| gdtr(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_gdtrc() {
-        xsref::test::<f64, _>("gdtrc", "d_d_d-d", |x: &[f64]| gdtrc(x[0], x[1], x[2]));
+        testing::test::<f64, _>("gdtrc", "d_d_d-d", |x: &[f64]| gdtrc(x[0], x[1], x[2]));
     }
 
     #[test]
     fn test_pdtr() {
-        xsref::test::<f64, _>("pdtr", "d_d-d", |x: &[f64]| pdtr(x[0], x[1]));
+        testing::test::<f64, _>("pdtr", "d_d-d", |x: &[f64]| pdtr(x[0], x[1]));
     }
 
     #[test]
     fn test_pdtrc() {
-        xsref::test::<f64, _>("pdtrc", "d_d-d", |x: &[f64]| pdtrc(x[0], x[1]));
+        testing::test::<f64, _>("pdtrc", "d_d-d", |x: &[f64]| pdtrc(x[0], x[1]));
     }
 
     #[test]
     fn test_pdtri() {
-        xsref::test::<f64, _>("pdtri", "p_d-d", |x: &[f64]| pdtri(x[0] as i32, x[1]));
+        testing::test::<f64, _>("pdtri", "p_d-d", |x: &[f64]| pdtri(x[0] as i32, x[1]));
     }
 
     #[test]
     fn test_bdtr() {
-        xsref::test::<f64, _>("bdtr", "d_p_d-d", |x: &[f64]| bdtr(x[0], x[1] as i32, x[2]));
+        testing::test::<f64, _>("bdtr", "d_p_d-d", |x: &[f64]| bdtr(x[0], x[1] as i32, x[2]));
     }
 
     #[test]
     fn test_bdtrc() {
-        xsref::test::<f64, _>("bdtrc", "d_p_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("bdtrc", "d_p_d-d", |x: &[f64]| {
             bdtrc(x[0], x[1] as i32, x[2])
         });
     }
 
     #[test]
     fn test_bdtri() {
-        xsref::test::<f64, _>("bdtri", "d_p_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("bdtri", "d_p_d-d", |x: &[f64]| {
             bdtri(x[0], x[1] as i32, x[2])
         });
     }
 
     #[test]
     fn test_nbdtr() {
-        xsref::test::<f64, _>("nbdtr", "p_p_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("nbdtr", "p_p_d-d", |x: &[f64]| {
             nbdtr(x[0] as i32, x[1] as i32, x[2])
         });
     }
 
     #[test]
     fn test_nbdtrc() {
-        xsref::test::<f64, _>("nbdtrc", "p_p_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("nbdtrc", "p_p_d-d", |x: &[f64]| {
             nbdtrc(x[0] as i32, x[1] as i32, x[2])
         });
     }
diff --git a/src/struve.rs b/src/struve.rs
index 751d63c..6e4aca3 100644
--- a/src/struve.rs
+++ b/src/struve.rs
@@ -10,30 +10,30 @@ xsf_impl!(struve_l, (v: f64, x: f64), "Struve `L` function");
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_itstruve0_f64() {
-        xsref::test::<f64, _>("itstruve0", "d-d", |x: &[f64]| itstruve0(x[0]));
+        testing::test::<f64, _>("itstruve0", "d-d", |x: &[f64]| itstruve0(x[0]));
     }
 
     #[test]
     fn test_it2struve0_f64() {
-        xsref::test::<f64, _>("it2struve0", "d-d", |x: &[f64]| it2struve0(x[0]));
+        testing::test::<f64, _>("it2struve0", "d-d", |x: &[f64]| it2struve0(x[0]));
     }
 
     #[test]
     fn test_itmodstruve0_f64() {
-        xsref::test::<f64, _>("itmodstruve0", "d-d", |x: &[f64]| itmodstruve0(x[0]));
+        testing::test::<f64, _>("itmodstruve0", "d-d", |x: &[f64]| itmodstruve0(x[0]));
     }
 
     #[test]
     fn test_struve_h_f64() {
-        xsref::test::<f64, _>("struve_h", "d_d-d", |x: &[f64]| struve_h(x[0], x[1]));
+        testing::test::<f64, _>("struve_h", "d_d-d", |x: &[f64]| struve_h(x[0], x[1]));
     }
 
     #[test]
     fn test_struve_l_f64() {
-        xsref::test::<f64, _>("struve_l", "d_d-d", |x: &[f64]| struve_l(x[0], x[1]));
+        testing::test::<f64, _>("struve_l", "d_d-d", |x: &[f64]| struve_l(x[0], x[1]));
     }
 }
diff --git a/src/test_utils.rs b/src/test_utils.rs
deleted file mode 100644
index 853a99a..0000000
--- a/src/test_utils.rs
+++ /dev/null
@@ -1,436 +0,0 @@
-#[cfg(test)]
-pub mod xsref {
-    use arrow::array::{Array, Float64Array, Int32Array, Int64Array};
-    use num_complex::{Complex, c64};
-    use parquet::arrow::arrow_reader::ParquetRecordBatchReaderBuilder;
-    use std::env;
-    use std::fs::File;
-    use std::path::PathBuf;
-
-    pub trait TestOutputValue: Copy {
-        fn error(self, expected: Self) -> f64;
-        fn magnitude(self) -> f64;
-        fn format(self) -> String;
-    }
-
-    impl TestOutputValue for f64 {
-        fn error(self, expected: Self) -> f64 {
-            relative_error(self, expected)
-        }
-
-        fn magnitude(self) -> f64 {
-            self.abs()
-        }
-
-        fn format(self) -> String {
-            format!("{:.6e}", self)
-        }
-    }
-
-    impl TestOutputValue for Complex<f64> {
-        fn error(self, expected: Self) -> f64 {
-            complex_relative_error(self, expected)
-        }
-
-        fn magnitude(self) -> f64 {
-            self.norm()
-        }
-
-        fn format(self) -> String {
-            format!("{:.6e}+{:.6e}i", self.re, self.im)
-        }
-    }
-
-    pub trait TestOutput: Copy {
-        type Value: TestOutputValue;
-
-        fn values(&self) -> Vec<Self::Value>;
-
-        fn from_parquet_row(row: Vec<f64>) -> Self;
-
-        fn from_parquet_rows(rows: Vec<Vec<f64>>) -> Vec<Self> {
-            rows.into_iter().map(Self::from_parquet_row).collect()
-        }
-
-        fn magnitude(self) -> f64 {
-            let values = self.values();
-            values.iter().map(|x| x.magnitude()).sum::<f64>() / values.len() as f64
-        }
-
-        fn is_huge(self) -> bool {
-            self.magnitude() > f64::EPSILON
-        }
-
-        fn error(self, expected: Self) -> f64 {
-            // max adjusted relative error
-            self.values()
-                .iter()
-                .zip(expected.values().iter())
-                .map(|(&a, &e)| a.error(e))
-                .fold(0.0, |acc, x| acc.max(x))
-        }
-
-        fn format(self) -> String {
-            let values = self.values();
-            if values.len() == 1 {
-                values[0].format()
-            } else {
-                format!(
-                    "({})",
-                    values
-                        .iter()
-                        .map(|x| x.format())
-                        .collect::<Vec<_>>()
-                        .join(", ")
-                )
-            }
-        }
-    }
-
-    impl TestOutput for f64 {
-        type Value = f64;
-
-        fn values(&self) -> Vec<Self::Value> {
-            vec![*self]
-        }
-
-        fn from_parquet_row(row: Vec<f64>) -> Self {
-            row[0]
-        }
-    }
-
-    impl TestOutput for Complex<f64> {
-        type Value = Complex<f64>;
-
-        fn values(&self) -> Vec<Self::Value> {
-            vec![*self]
-        }
-
-        fn from_parquet_row(row: Vec<f64>) -> Self {
-            c64(row[0], row[1])
-        }
-    }
-
-    impl TestOutput for (f64, f64) {
-        type Value = f64;
-
-        fn values(&self) -> Vec<Self::Value> {
-            vec![self.0, self.1]
-        }
-
-        fn from_parquet_row(row: Vec<f64>) -> Self {
-            (row[0], row[1])
-        }
-    }
-
-    impl TestOutput for (Complex<f64>, Complex<f64>) {
-        type Value = Complex<f64>;
-
-        fn values(&self) -> Vec<Self::Value> {
-            vec![self.0, self.1]
-        }
-
-        fn from_parquet_row(row: Vec<f64>) -> Self {
-            (c64(row[0], row[1]), c64(row[2], row[3]))
-        }
-    }
-
-    impl TestOutput for (f64, f64, f64, f64) {
-        type Value = f64;
-
-        fn values(&self) -> Vec<Self::Value> {
-            vec![self.0, self.1, self.2, self.3]
-        }
-
-        fn from_parquet_row(row: Vec<f64>) -> Self {
-            (row[0], row[1], row[2], row[3])
-        }
-    }
-
-    impl TestOutput for (Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>) {
-        type Value = Complex<f64>;
-
-        fn values(&self) -> Vec<Self::Value> {
-            vec![self.0, self.1, self.2, self.3]
-        }
-
-        fn from_parquet_row(row: Vec<f64>) -> Self {
-            (
-                c64(row[0], row[1]),
-                c64(row[2], row[3]),
-                c64(row[4], row[5]),
-                c64(row[6], row[7]),
-            )
-        }
-    }
-
-    #[derive(Debug)]
-    struct TestCase<T: TestOutput> {
-        pub inputs: Vec<f64>,
-        pub expected: T,
-        pub tolerance: f64,
-    }
-
-    #[derive(Debug)]
-    #[allow(dead_code)]
-    pub enum TestError {
-        Io(std::io::Error),
-        DataFormat,
-    }
-
-    impl From<std::io::Error> for TestError {
-        fn from(err: std::io::Error) -> Self {
-            TestError::Io(err)
-        }
-    }
-
-    impl From<arrow::error::ArrowError> for TestError {
-        fn from(_: arrow::error::ArrowError) -> Self {
-            TestError::DataFormat
-        }
-    }
-
-    impl From<parquet::errors::ParquetError> for TestError {
-        fn from(_: parquet::errors::ParquetError) -> Self {
-            TestError::DataFormat
-        }
-    }
-
-    fn read_parquet_rows(file: File) -> Vec<Vec<f64>> {
-        let builder = ParquetRecordBatchReaderBuilder::try_new(file).unwrap();
-        let reader = builder.build().unwrap();
-
-        let mut rows = Vec::new();
-        for batch_result in reader {
-            let batch = batch_result.unwrap();
-            for row_idx in 0..batch.num_rows() {
-                let mut row = Vec::new();
-                for col_idx in 0..batch.num_columns() {
-                    let column = batch.column(col_idx).as_any();
-                    if let Some(f64_array) = column.downcast_ref::<Float64Array>() {
-                        row.push(f64_array.value(row_idx));
-                    } else if let Some(i32_array) = column.downcast_ref::<Int32Array>() {
-                        row.push(i32_array.value(row_idx) as f64);
-                    } else if let Some(i64_array) = column.downcast_ref::<Int64Array>() {
-                        row.push(i64_array.value(row_idx) as f64);
-                    }
-                }
-                rows.push(row);
-            }
-        }
-        rows
-    }
-
-    fn read_parquet_column(file: File) -> Vec<f64> {
-        let builder = ParquetRecordBatchReaderBuilder::try_new(file).unwrap();
-        let reader = builder.build().unwrap();
-
-        let mut values = Vec::new();
-        for batch in reader {
-            if let Some(column) = batch
-                .unwrap()
-                .column(0)
-                .as_any()
-                .downcast_ref::<Float64Array>()
-            {
-                values.extend(column.iter().map(|v| v.unwrap()));
-            }
-        }
-        values
-    }
-
-    fn read_parquet_output<T: TestOutput>(file: File) -> Vec<T> {
-        let builder = ParquetRecordBatchReaderBuilder::try_new(file).unwrap();
-        let reader = builder.build().unwrap();
-
-        let mut rows = Vec::new();
-        for batch_result in reader {
-            let batch = batch_result.unwrap();
-            for row_idx in 0..batch.num_rows() {
-                let mut row = Vec::new();
-                for col_idx in 0..batch.num_columns() {
-                    if let Some(f64_array) = batch
-                        .column(col_idx)
-                        .as_any()
-                        .downcast_ref::<Float64Array>()
-                    {
-                        row.push(f64_array.value(row_idx));
-                    }
-                }
-                rows.push(row);
-            }
-        }
-        T::from_parquet_rows(rows)
-    }
-
-    fn load_testcases<T: TestOutput>(
-        name: &str,
-        signature: &str,
-    ) -> Result<Vec<TestCase<T>>, TestError> {
-        let tables_dir = env::var("XSREF_TABLES_PATH")
-            .map(PathBuf::from)
-            .unwrap_or_else(|_| PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("xsref/tables"))
-            .join("scipy_special_tests")
-            .join(name);
-
-        let get_table = |name: &str| File::open(tables_dir.join(format!("{}.parquet", name)));
-
-        let table_in = read_parquet_rows(get_table(&format!("In_{}", signature))?);
-        let table_out = read_parquet_output::<T>(get_table(&format!("Out_{}", signature))?);
-
-        let platform = if cfg!(all(target_arch = "x86_64", target_os = "linux")) {
-            "gcc-linux-x86_64"
-        } else if cfg!(all(target_arch = "aarch64", target_os = "macos")) {
-            "clang-darwin-aarch64"
-        } else {
-            "other"
-        };
-
-        let table_err_file = get_table(&format!("Err_{}_{}", signature, platform))
-            .or_else(|_| get_table(&format!("Err_{}_other", signature)))?;
-        let table_err = read_parquet_column(table_err_file);
-
-        let test_cases = table_in
-            .into_iter()
-            .zip(table_out)
-            .zip(table_err)
-            .map(|((inputs, expected), tolerance)| TestCase {
-                inputs,
-                expected,
-                tolerance,
-            })
-            .collect();
-
-        Ok(test_cases)
-    }
-
-    /// based on `xsref.float_tools.extended_absolute_error`
-    fn absolute_error(actual: f64, desired: f64) -> f64 {
-        if actual == desired || (actual.is_nan() && desired.is_nan()) {
-            return 0.0;
-        } else if desired.is_nan() || actual.is_nan() {
-            return f64::INFINITY;
-        }
-
-        // f64::MAX + ULP
-        let max1p = f64::MAX * (1.0 + 2.0_f64.powi(-(f64::MANTISSA_DIGITS as i32)));
-        if actual.is_infinite() {
-            (actual.signum() * max1p - desired).abs()
-        } else if desired.is_infinite() {
-            (actual - desired.signum() * max1p).abs()
-        } else {
-            (actual - desired).abs()
-        }
-    }
-
-    /// based on `xsref.float_tools.extended_relative_error`
-    fn relative_error(actual: f64, desired: f64) -> f64 {
-        let abserr0 = if desired == 0.0 {
-            f64::MIN_POSITIVE
-        } else if desired.is_infinite() {
-            f64::MAX
-        } else if desired.is_nan() {
-            1.0
-        } else {
-            desired.abs()
-        };
-
-        let abserr = absolute_error(actual, desired);
-        (abserr / abserr0).min(abserr)
-    }
-
-    /// based on `xsref.float_tools.extended_absolute_error_complex`
-    fn complex_absolute_error(actual: Complex<f64>, expected: Complex<f64>) -> f64 {
-        let real_error = (actual.re - expected.re).abs();
-        let imag_error = (actual.im - expected.im).abs();
-
-        let max_error = real_error.max(imag_error);
-        if max_error == 0.0 {
-            0.0
-        } else if max_error.is_infinite() {
-            f64::INFINITY
-        } else {
-            real_error.hypot(imag_error)
-        }
-    }
-
-    /// based on `xsref.float_tools.extended_relative_error_complex`
-    fn complex_relative_error(actual: Complex<f64>, expected: Complex<f64>) -> f64 {
-        if actual.re.is_nan() || actual.im.is_nan() {
-            if expected.re.is_nan() || expected.im.is_nan() {
-                0.0
-            } else {
-                f64::INFINITY
-            }
-        } else if expected.re.is_infinite() || expected.im.is_infinite() {
-            // If any component of expected is infinite, use component-wise relative error
-            fn rel_error(x: f64, y: f64) -> f64 {
-                if y.is_infinite() {
-                    if x.is_infinite() && x.signum() == y.signum() {
-                        0.0
-                    } else {
-                        f64::INFINITY
-                    }
-                } else if y == 0.0 {
-                    x.abs()
-                } else {
-                    (x / y - 1.0).abs()
-                }
-            }
-
-            let real_rel_err = rel_error(actual.re, expected.re);
-            let imag_rel_err = rel_error(actual.im, expected.im);
-            real_rel_err.max(imag_rel_err)
-        } else {
-            // For finite expected values, use the magnitude-based relative error
-            let abs_error = complex_absolute_error(actual, expected);
-            let expected_magnitude = expected.re.hypot(expected.im);
-            if expected_magnitude == 0.0 {
-                abs_error
-            } else {
-                abs_error / expected_magnitude
-            }
-        }
-    }
-
-    pub fn test<T, F>(name: &str, signature: &str, test_fn: F)
-    where
-        T: TestOutput,
-        F: Fn(&[f64]) -> T,
-    {
-        let cases = load_testcases::<T>(name, signature).unwrap();
-
-        let mut failed = 0;
-        for (i, case) in cases.iter().enumerate() {
-            let actual = test_fn(&case.inputs);
-            let desired = case.expected;
-            let max_error = case.tolerance.max(f64::EPSILON) * 16.0;
-
-            if max_error.is_huge() || desired.is_huge() && actual.is_huge() {
-                continue;
-            }
-
-            let error = T::error(actual, desired);
-            if error > max_error {
-                failed += 1;
-                eprintln!(
-                    "{}: test {}, expected {}, got {}, error {:.2e}, max {:.2e}",
-                    name,
-                    i,
-                    desired.format(),
-                    actual.format(),
-                    error,
-                    max_error
-                );
-            }
-        }
-
-        assert!(
-            failed == 0,
-            "{}: {}/{} tests failed",
-            name,
-            failed,
-            cases.len()
-        );
-    }
-}
diff --git a/src/testing.rs b/src/testing.rs
new file mode 100644
index 0000000..7a50519
--- /dev/null
+++ b/src/testing.rs
@@ -0,0 +1,433 @@
+use arrow::array::{Array, Float64Array, Int32Array, Int64Array};
+use num_complex::{Complex, c64};
+use parquet::arrow::arrow_reader::ParquetRecordBatchReaderBuilder;
+use std::env;
+use std::fs::File;
+use std::path::PathBuf;
+
+pub trait TestOutputValue: Copy {
+    fn error(self, expected: Self) -> f64;
+    fn magnitude(self) -> f64;
+    fn format(self) -> String;
+}
+
+impl TestOutputValue for f64 {
+    fn error(self, expected: Self) -> f64 {
+        relative_error(self, expected)
+    }
+
+    fn magnitude(self) -> f64 {
+        self.abs()
+    }
+
+    fn format(self) -> String {
+        format!("{:.6e}", self)
+    }
+}
+
+impl TestOutputValue for Complex<f64> {
+    fn error(self, expected: Self) -> f64 {
+        complex_relative_error(self, expected)
+    }
+
+    fn magnitude(self) -> f64 {
+        self.norm()
+    }
+
+    fn format(self) -> String {
+        format!("{:.6e}+{:.6e}i", self.re, self.im)
+    }
+}
+
+pub trait TestOutput: Copy {
+    type Value: TestOutputValue;
+
+    fn values(&self) -> Vec<Self::Value>;
+
+    fn from_parquet_row(row: Vec<f64>) -> Self;
+
+    fn from_parquet_rows(rows: Vec<Vec<f64>>) -> Vec<Self> {
+        rows.into_iter().map(Self::from_parquet_row).collect()
+    }
+
+    fn magnitude(self) -> f64 {
+        let values = self.values();
+        values.iter().map(|x| x.magnitude()).sum::<f64>() / values.len() as f64
+    }
+
+    fn is_huge(self) -> bool {
+        self.magnitude() > f64::EPSILON
+    }
+
+    fn error(self, expected: Self) -> f64 {
+        // max adjusted relative error
+        self.values()
+            .iter()
+            .zip(expected.values().iter())
+            .map(|(&a, &e)| a.error(e))
+            .fold(0.0, |acc, x| acc.max(x))
+    }
+
+    fn format(self) -> String {
+        let values = self.values();
+        if values.len() == 1 {
+            values[0].format()
+        } else {
+            format!(
+                "({})",
+                values
+                    .iter()
+                    .map(|x| x.format())
+                    .collect::<Vec<_>>()
+                    .join(", ")
+            )
+        }
+    }
+}
+
+impl TestOutput for f64 {
+    type Value = f64;
+
+    fn values(&self) -> Vec<Self::Value> {
+        vec![*self]
+    }
+
+    fn from_parquet_row(row: Vec<f64>) -> Self {
+        row[0]
+    }
+}
+
+impl TestOutput for Complex<f64> {
+    type Value = Complex<f64>;
+
+    fn values(&self) -> Vec<Self::Value> {
+        vec![*self]
+    }
+
+    fn from_parquet_row(row: Vec<f64>) -> Self {
+        c64(row[0], row[1])
+    }
+}
+
+impl TestOutput for (f64, f64) {
+    type Value = f64;
+
+    fn values(&self) -> Vec<Self::Value> {
+        vec![self.0, self.1]
+    }
+
+    fn from_parquet_row(row: Vec<f64>) -> Self {
+        (row[0], row[1])
+    }
+}
+
+impl TestOutput for (Complex<f64>, Complex<f64>) {
+    type Value = Complex<f64>;
+
+    fn values(&self) -> Vec<Self::Value> {
+        vec![self.0, self.1]
+    }
+
+    fn from_parquet_row(row: Vec<f64>) -> Self {
+        (c64(row[0], row[1]), c64(row[2], row[3]))
+    }
+}
+
+impl TestOutput for (f64, f64, f64, f64) {
+    type Value = f64;
+
+    fn values(&self) -> Vec<Self::Value> {
+        vec![self.0, self.1, self.2, self.3]
+    }
+
+    fn from_parquet_row(row: Vec<f64>) -> Self {
+        (row[0], row[1], row[2], row[3])
+    }
+}
+
+impl TestOutput for (Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>) {
+    type Value = Complex<f64>;
+
+    fn values(&self) -> Vec<Self::Value> {
+        vec![self.0, self.1, self.2, self.3]
+    }
+
+    fn from_parquet_row(row: Vec<f64>) -> Self {
+        (
+            c64(row[0], row[1]),
+            c64(row[2], row[3]),
+            c64(row[4], row[5]),
+            c64(row[6], row[7]),
+        )
+    }
+}
+
+#[derive(Debug)]
+struct TestCase<T: TestOutput> {
+    pub inputs: Vec<f64>,
+    pub expected: T,
+    pub tolerance: f64,
+}
+
+#[derive(Debug)]
+#[allow(dead_code)]
+pub enum TestError {
+    Io(std::io::Error),
+    DataFormat,
+}
+
+impl From<std::io::Error> for TestError {
+    fn from(err: std::io::Error) -> Self {
+        TestError::Io(err)
+    }
+}
+
+impl From<arrow::error::ArrowError> for TestError {
+    fn from(_: arrow::error::ArrowError) -> Self {
+        TestError::DataFormat
+    }
+}
+
+impl From<parquet::errors::ParquetError> for TestError {
+    fn from(_: parquet::errors::ParquetError) -> Self {
+        TestError::DataFormat
+    }
+}
+
+fn read_parquet_rows(file: File) -> Vec<Vec<f64>> {
+    let builder = ParquetRecordBatchReaderBuilder::try_new(file).unwrap();
+    let reader = builder.build().unwrap();
+
+    let mut rows = Vec::new();
+    for batch_result in reader {
+        let batch = batch_result.unwrap();
+        for row_idx in 0..batch.num_rows() {
+            let mut row = Vec::new();
+            for col_idx in 0..batch.num_columns() {
+                let column = batch.column(col_idx).as_any();
+                if let Some(f64_array) = column.downcast_ref::<Float64Array>() {
+                    row.push(f64_array.value(row_idx));
+                } else if let Some(i32_array) = column.downcast_ref::<Int32Array>() {
+                    row.push(i32_array.value(row_idx) as f64);
+                } else if let Some(i64_array) = column.downcast_ref::<Int64Array>() {
+                    row.push(i64_array.value(row_idx) as f64);
+                }
+            }
+            rows.push(row);
+        }
+    }
+    rows
+}
+
+fn read_parquet_column(file: File) -> Vec<f64> {
+    let builder = ParquetRecordBatchReaderBuilder::try_new(file).unwrap();
+    let reader = builder.build().unwrap();
+
+    let mut values = Vec::new();
+    for batch in reader {
+        if let Some(column) = batch
+            .unwrap()
+            .column(0)
+            .as_any()
+            .downcast_ref::<Float64Array>()
+        {
+            values.extend(column.iter().map(|v| v.unwrap()));
+        }
+    }
+    values
+}
+
+fn read_parquet_output<T: TestOutput>(file: File) -> Vec<T> {
+    let builder = ParquetRecordBatchReaderBuilder::try_new(file).unwrap();
+    let reader = builder.build().unwrap();
+
+    let mut rows = Vec::new();
+    for batch_result in reader {
+        let batch = batch_result.unwrap();
+        for row_idx in 0..batch.num_rows() {
+            let mut row = Vec::new();
+            for col_idx in 0..batch.num_columns() {
+                if let Some(f64_array) = batch
+                    .column(col_idx)
+                    .as_any()
+                    .downcast_ref::<Float64Array>()
+                {
+                    row.push(f64_array.value(row_idx));
+                }
+            }
+            rows.push(row);
+        }
+    }
+    T::from_parquet_rows(rows)
+}
+
+fn load_testcases<T: TestOutput>(
+    name: &str,
+    signature: &str,
+) -> Result<Vec<TestCase<T>>, TestError> {
+    let tables_dir = env::var("XSREF_TABLES_PATH")
+        .map(PathBuf::from)
+        .unwrap_or_else(|_| PathBuf::from(env!("CARGO_MANIFEST_DIR")).join("xsref/tables"))
+        .join("scipy_special_tests")
+        .join(name);
+
+    let get_table = |name: &str| File::open(tables_dir.join(format!("{}.parquet", name)));
+
+    let table_in = read_parquet_rows(get_table(&format!("In_{}", signature))?);
+    let table_out = read_parquet_output::<T>(get_table(&format!("Out_{}", signature))?);
+
+    let platform = if cfg!(all(target_arch = "x86_64", target_os = "linux")) {
+        "gcc-linux-x86_64"
+    } else if cfg!(all(target_arch = "aarch64", target_os = "macos")) {
+        "clang-darwin-aarch64"
+    } else {
+        "other"
+    };
+
+    let table_err_file = get_table(&format!("Err_{}_{}", signature, platform))
+        .or_else(|_| get_table(&format!("Err_{}_other", signature)))?;
+    let table_err = read_parquet_column(table_err_file);
+
+    let test_cases = table_in
+        .into_iter()
+        .zip(table_out)
+        .zip(table_err)
+        .map(|((inputs, expected), tolerance)| TestCase {
+            inputs,
+            expected,
+            tolerance,
+        })
+        .collect();
+
+    Ok(test_cases)
+}
+
+/// based on `xsref.float_tools.extended_absolute_error`
+fn absolute_error(actual: f64, desired: f64) -> f64 {
+    if actual == desired || (actual.is_nan() && desired.is_nan()) {
+        return 0.0;
+    } else if desired.is_nan() || actual.is_nan() {
+        return f64::INFINITY;
+    }
+
+    // f64::MAX + ULP
+    let max1p = f64::MAX * (1.0 + 2.0_f64.powi(-(f64::MANTISSA_DIGITS as i32)));
+    if actual.is_infinite() {
+        (actual.signum() * max1p - desired).abs()
+    } else if desired.is_infinite() {
+        (actual - desired.signum() * max1p).abs()
+    } else {
+        (actual - desired).abs()
+    }
+}
+
+/// based on `xsref.float_tools.extended_relative_error`
+fn relative_error(actual: f64, desired: f64) -> f64 {
+    let abserr0 = if desired == 0.0 {
+        f64::MIN_POSITIVE
+    } else if desired.is_infinite() {
+        f64::MAX
+    } else if desired.is_nan() {
+        1.0
+    } else {
+        desired.abs()
+    };
+
+    let abserr = absolute_error(actual, desired);
+    (abserr / abserr0).min(abserr)
+}
+
+/// based on `xsref.float_tools.extended_absolute_error_complex`
+fn complex_absolute_error(actual: Complex<f64>, expected: Complex<f64>) -> f64 {
+    let real_error = (actual.re - expected.re).abs();
+    let imag_error = (actual.im - expected.im).abs();
+
+    let max_error = real_error.max(imag_error);
+    if max_error == 0.0 {
+        0.0
+    } else if max_error.is_infinite() {
+        f64::INFINITY
+    } else {
+        real_error.hypot(imag_error)
+    }
+}
+
+/// based on `xsref.float_tools.extended_relative_error_complex`
+fn complex_relative_error(actual: Complex<f64>, expected: Complex<f64>) -> f64 {
+    if actual.re.is_nan() || actual.im.is_nan() {
+        if expected.re.is_nan() || expected.im.is_nan() {
+            0.0
+        } else {
+            f64::INFINITY
+        }
+    } else if expected.re.is_infinite() || expected.im.is_infinite() {
+        // If any component of expected is infinite, use component-wise relative error
+        fn rel_error(x: f64, y: f64) -> f64 {
+            if y.is_infinite() {
+                if x.is_infinite() && x.signum() == y.signum() {
+                    0.0
+                } else {
+                    f64::INFINITY
+                }
+            } else if y == 0.0 {
+                x.abs()
+            } else {
+                (x / y - 1.0).abs()
+            }
+        }
+
+        let real_rel_err = rel_error(actual.re, expected.re);
+        let imag_rel_err = rel_error(actual.im, expected.im);
+        real_rel_err.max(imag_rel_err)
+    } else {
+        // For finite expected values, use the magnitude-based relative error
+        let abs_error = complex_absolute_error(actual, expected);
+        let expected_magnitude = expected.re.hypot(expected.im);
+        if expected_magnitude == 0.0 {
+            abs_error
+        } else {
+            abs_error / expected_magnitude
+        }
+    }
+}
+
+pub fn test<T, F>(name: &str, signature: &str, test_fn: F)
+where
+    T: TestOutput,
+    F: Fn(&[f64]) -> T,
+{
+    let cases = load_testcases::<T>(name, signature).unwrap();
+
+    let mut failed = 0;
+    for (i, case) in cases.iter().enumerate() {
+        let actual = test_fn(&case.inputs);
+        let desired = case.expected;
+        let max_error = case.tolerance.max(f64::EPSILON) * 16.0;
+
+        if max_error.is_huge() || desired.is_huge() && actual.is_huge() {
+            continue;
+        }
+
+        let error = T::error(actual, desired);
+        if error > max_error {
+            failed += 1;
+            eprintln!(
+                "{}: test {}, expected {}, got {}, error {:.2e}, max {:.2e}",
+                name,
+                i,
+                desired.format(),
+                actual.format(),
+                error,
+                max_error
+            );
+        }
+    }
+
+    assert!(
+        failed == 0,
+        "{}: {}/{} tests failed",
+        name,
+        failed,
+        cases.len()
+    );
+}
diff --git a/src/trig.rs b/src/trig.rs
index f7b6b1b..654186d 100644
--- a/src/trig.rs
+++ b/src/trig.rs
@@ -58,56 +58,56 @@ xsf_impl!(radian, (d: f64, m: f64, s: f64), "Degrees, minutes, seconds to radian
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_sinpi_f64() {
-        xsref::test::<f64, _>("sinpi", "d-d", |x: &[f64]| sinpi(x[0]));
+        testing::test::<f64, _>("sinpi", "d-d", |x: &[f64]| sinpi(x[0]));
     }
 
     #[test]
     fn test_sinpi_c64() {
-        xsref::test::<Complex<f64>, _>("sinpi", "cd-cd", |x: &[f64]| sinpi(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("sinpi", "cd-cd", |x: &[f64]| sinpi(c64(x[0], x[1])));
     }
 
     #[test]
     fn test_cospi_f64() {
-        xsref::test::<f64, _>("cospi", "d-d", |x: &[f64]| cospi(x[0]));
+        testing::test::<f64, _>("cospi", "d-d", |x: &[f64]| cospi(x[0]));
     }
 
     #[test]
     fn test_cospi_c64() {
-        xsref::test::<Complex<f64>, _>("cospi", "cd-cd", |x: &[f64]| cospi(c64(x[0], x[1])));
+        testing::test::<Complex<f64>, _>("cospi", "cd-cd", |x: &[f64]| cospi(c64(x[0], x[1])));
     }
 
     #[test]
     fn test_sindg() {
-        xsref::test::<f64, _>("sindg", "d-d", |x: &[f64]| sindg(x[0]));
+        testing::test::<f64, _>("sindg", "d-d", |x: &[f64]| sindg(x[0]));
     }
 
     #[test]
     fn test_cosdg() {
-        xsref::test::<f64, _>("cosdg", "d-d", |x: &[f64]| cosdg(x[0]));
+        testing::test::<f64, _>("cosdg", "d-d", |x: &[f64]| cosdg(x[0]));
     }
 
     #[test]
     fn test_tandg() {
-        xsref::test::<f64, _>("tandg", "d-d", |x: &[f64]| tandg(x[0]));
+        testing::test::<f64, _>("tandg", "d-d", |x: &[f64]| tandg(x[0]));
     }
 
     #[test]
     fn test_cotdg() {
-        xsref::test::<f64, _>("cotdg", "d-d", |x: &[f64]| cotdg(x[0]));
+        testing::test::<f64, _>("cotdg", "d-d", |x: &[f64]| cotdg(x[0]));
     }
 
     #[test]
     fn test_cosm1() {
-        xsref::test::<f64, _>("cosm1", "d-d", |x: &[f64]| cosm1(x[0]));
+        testing::test::<f64, _>("cosm1", "d-d", |x: &[f64]| cosm1(x[0]));
     }
 
     #[test]
     fn test_radian() {
-        xsref::test::<f64, _>("radian", "d_d_d-d", |x: &[f64]| radian(x[0], x[1], x[2]));
+        testing::test::<f64, _>("radian", "d_d_d-d", |x: &[f64]| radian(x[0], x[1], x[2]));
     }
 }
diff --git a/src/wright_bessel.rs b/src/wright_bessel.rs
index 2a94194..e6349c0 100644
--- a/src/wright_bessel.rs
+++ b/src/wright_bessel.rs
@@ -6,18 +6,18 @@ xsf_impl!(log_wright_bessel, (a: f64, b: f64, x: f64), "Logarithm of `wright_bes
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
 
     #[test]
     fn test_wright_bessel_f64() {
-        xsref::test::<f64, _>("wright_bessel", "d_d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("wright_bessel", "d_d_d-d", |x: &[f64]| {
             wright_bessel(x[0], x[1], x[2])
         });
     }
 
     #[test]
     fn test_log_wright_bessel_f64() {
-        xsref::test::<f64, _>("log_wright_bessel", "d_d_d-d", |x: &[f64]| {
+        testing::test::<f64, _>("log_wright_bessel", "d_d_d-d", |x: &[f64]| {
             log_wright_bessel(x[0], x[1], x[2])
         });
     }
diff --git a/src/zeta.rs b/src/zeta.rs
index 9d1ee1b..879369e 100644
--- a/src/zeta.rs
+++ b/src/zeta.rs
@@ -52,28 +52,28 @@ xsf_impl!(zetac, (x: f64), "Riemann zeta function, minus one");
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::xsref;
+    use crate::testing;
     use num_complex::{Complex, c64};
 
     #[test]
     fn test_riemann_zeta_f64() {
-        xsref::test::<f64, _>("riemann_zeta", "d-d", |x: &[f64]| riemann_zeta(x[0]));
+        testing::test::<f64, _>("riemann_zeta", "d-d", |x: &[f64]| riemann_zeta(x[0]));
     }
 
     #[test]
     fn test_riemann_zeta_c64() {
-        xsref::test::<Complex<f64>, _>("riemann_zeta", "cd-cd", |x: &[f64]| {
+        testing::test::<Complex<f64>, _>("riemann_zeta", "cd-cd", |x: &[f64]| {
             riemann_zeta(c64(x[0], x[1]))
         });
     }
 
     #[test]
     fn test_zeta_f64() {
-        xsref::test::<f64, _>("zeta", "d_d-d", |x: &[f64]| zeta(x[0], x[1]));
+        testing::test::<f64, _>("zeta", "d_d-d", |x: &[f64]| zeta(x[0], x[1]));
     }
 
     #[test]
     fn test_zetac_f64() {
-        xsref::test::<f64, _>("zetac", "d-d", |x: &[f64]| zetac(x[0]));
+        testing::test::<f64, _>("zetac", "d-d", |x: &[f64]| zetac(x[0]));
     }
 }