scipycpp is a header-only C++ library implementing scipy's core APIs with
bit-level precision alignment against Python scipy.
Built on proven C++ libraries:
- numpycpp — numpy primitives
- Eigen3 — linear algebra (solve, inv, det, eig, svd, cholesky)
- pocketfft (bundled) — FFT (same library numpy/scipy use internally)
- ckdtree (bundled) — KDTree from scipy.spatial
#include"scipy/core.h"// Integrationauto [val, err] = scipy::integrate::quad([](double x){ return x*x; }, 0.0, 4.0);
// Stats
std::vector<double> pdf(3);
scipy::stats::norm_pdf(data, pdf.data(), 3);
// Linalg (Eigen3-backed)double A[4] = {2,1,1,3}, b[2] = {5,6}, x[2];
scipy::linalg::solve(A, x, b, 2);
// FFT (pocketfft-backed — same as numpy/scipy)scipy::fft::fft(signal, spectrum, 4);
// Spatial — KDTree with distances
scipy::spatial::KDTree<double> tree(points, n_pts, dim);
double dist; size_t idx;
tree.query(query_point, dist, idx); // returns distance AND index// Spatial — cross-set distance matrix (cdist)scipy::spatial::distance::cdist(XA, mA, XB, mB, dim, dm);
// ndimage — 1D Gaussian filter
std::vector<double> out(n);
scipy::ndimage::gaussian_filter1d(src, out.data(), n, sigma);
// Signal — median filterscipy::signal::medfilt(src, dst.data(), n, kernel_size);
// Transform — Rotation (from_matrix, from_euler, as_euler, as_matrix)auto rot = scipy::spatial::transform::Rotation<double>::from_matrix(R);
auto euler = rot.as_euler_vec("xyz"); // returns [rx, ry, rz]// from_euler + as_matrix — single axis (key use case: ego yaw)double m9[9];
auto rot_z = scipy::spatial::transform::Rotation<double>::from_euler("z", yaw);
rot_z.as_matrix(m9); // 3×3 row-major rotation matrix// from_euler — multi-axis extrinsic / intrinsicdouble angles[3] = {rx, ry, rz};
auto rot_xyz = scipy::spatial::transform::Rotation<double>::from_euler("xyz", angles);
rot_xyz.as_matrix(m9);# Install dependencies
sudo dpkg -i numpycpp-dev-*.deb
sudo apt-get install libeigen3-dev
# Install scipycpp
mkdir build &&cd build
cmake .. && make deb
sudo dpkg -i scipycpp-dev-*.debfind_package(scipycppREQUIRED)
target_link_libraries(myappPRIVATEscipycpp::scipycpp)Under the bitexact build (numpycpp dlsym backend, same math kernels as scipy),
all APIs are 0 ULP bit-identical for both float64/float32, including extreme
values (±inf, NaN, ±0.0, subnormals, saturation inputs).
| Module | Backend | Key APIs | ULP Status |
|---|---|---|---|
stats | Cephes + numpcpp npy_exp (dlsym) | norm.pdf, norm.cdf, norm.ppf | ✅ 0 ULP |
integrate | sequential C++ sum | trapezoid, simpson | |
linalg | Eigen3 partialPivLu | solve | |
spatial | pure C++ / ckdtree | cdist, KDTree | ✅ 0 ULP |
ndimage | numpy kernel + C++ convolution | gaussian_filter1d | ✅ 0 ULP |
signal | sort-based median | medfilt | ✅ 0 ULP |
transform | pure C++ (no scipy delegation) | Rotation.from_matrix, from_euler, as_euler, as_matrix | ✅ 0 ULP single-axis; ≤200 ULP multi-axis |
Full per-test ULP report: doc/ulp_report.csv
(Auto-generated by pytest tests/test_all.py, see doc/ulp_report.md for summary)
| API | Max ULP | Root Cause |
|---|---|---|
linalg.solve | ≤8.9e4 | Eigen3 partialPivLu vs LAPACK gesv: different LU pivots produce different roundoff paths. Well-conditioned small matrices (2×2, identity) are bit-identical. All results within atol=1e-14 (ill-conditioned: atol=1e-10). |
| API | Input | Max ULP | Root Cause |
|---|---|---|---|
trapezoid / simpson | typical scientific data | 0 ULP | Sequential sum matches numpy pairwise for non-uniform data |
trapezoid | float64 uniform array | ≤5 f64-ULP | Sequential vs SIMD pairwise reorder (unavoidable without SIMD intrinsics) |
trapezoid | float32 uniform array | ≤4 f32-ULP | Same, measured in native float32 precision |
simpson | float64 uniform array | ≤6 f64-ULP | Same |
simpson | float32 input (any) | ≤6 f32-ULP | scipy computes internal sum in float32 SIMD; C++ uses sequential float32 |
scipy.integrate.simpson always returns float64 regardless of input dtype
(Python float constants promote the result). scipy::integrate::simpson<T>
matches this: it computes the intermediate sum in T (preserving scipy's
precision path), then multiplies by 1.0/3.0 in double and returns double.
norm.pdf: numpycpp v1.21.2+ resolves exp via dlsym("npy_exp"),
matching scipy's internal numpy math path bit-for-bit.
Bit-level alignment tests against Python scipy. Build flags mirror the numpycpp README bit-exact backend spec exactly.
# Build (cmake with bit-exact flags: -O2 -ffp-contract=off -msse4.1 -mfma, …)
cmake -S tests -B tests/build
cmake --build tests/build -j$(nproc)# Run alignment testscd tests && python3 -m pytest test_all.py -q --tb=short --no-header
# ULP report printed to stderr + exported to doc/ulp_report.csvSummary: doc/ulp_report.md — auto-generated aggregate of doc/ulp_report.csv.
MIT