From 91fd61aabe4f126e8e615a0cdd047c19c9b9bb41 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Tue, 13 Oct 2026 09:10:31 -0700 Subject: [PATCH] Draw the Python layer's ulps plots with the C++ ulps_plot and ulps_heatmap Reference: cpp_bin_float_50 or a Python callable; ulps_plot gains threads. Co-Authored-By: Claude Opus 5.5 --- include/boost/math/tools/ulps_heatmap.hpp | 84 +++---- include/boost/math/tools/ulps_plot.hpp | 244 ++++++++++++++------ nanobind/CMakeLists.txt | 22 ++ nanobind/boost/math/ulps.py | 257 ++++++++++------------ nanobind/src/bindings.hpp | 2 + nanobind/src/module.cpp | 3 + nanobind/src/ulps.cpp | 233 ++++++++++++++++++++ nanobind/src/ulps.hpp | 59 +++++ nanobind/src/ulps_elliptic.cpp | 56 +++++ nanobind/src/ulps_gamma.cpp | 38 ++++ nanobind/src/ulps_special.cpp | 86 ++++++++ nanobind/tests/test_ulps.py | 74 ++++++- test/test_ulps_plot.cpp | 25 ++- 13 files changed, 932 insertions(+), 251 deletions(-) create mode 100644 nanobind/src/ulps.cpp create mode 100644 nanobind/src/ulps.hpp create mode 100644 nanobind/src/ulps_elliptic.cpp create mode 100644 nanobind/src/ulps_gamma.cpp create mode 100644 nanobind/src/ulps_special.cpp diff --git a/include/boost/math/tools/ulps_heatmap.hpp b/include/boost/math/tools/ulps_heatmap.hpp index 540ef460aa..b42b75fa1c 100644 --- a/include/boost/math/tools/ulps_heatmap.hpp +++ b/include/boost/math/tools/ulps_heatmap.hpp @@ -32,7 +32,6 @@ #ifndef BOOST_MATH_BUILD_MODULE #include #include -#include #include #include #include @@ -43,7 +42,6 @@ #include #include #include -#include #include #endif #include @@ -138,36 +136,6 @@ std::vector heatmap_abscissas(CoarseReal a, CoarseReal b, std::size_ return xs; } -// Runs body(j) for j in [0, n) on several threads. -template -void parallel_for(std::size_t n, unsigned threads, Body body) -{ - if (n == 0) - { - return; - } - if (threads == 0) - { - threads = (std::max)(1u, std::thread::hardware_concurrency()); - } - std::atomic next{0}; - std::vector pool; - for (unsigned t = 0; t < threads; ++t) - { - pool.emplace_back([&]() { - for (std::size_t j = next++; j < n; j = next++) - { - body(j); - } - }); - } - for (auto & th : pool) - { - th.join(); - } -} - - // f_x and f_y at (qx, qy), between the nodes of the grid x (columns) by y (rows), from the bicubic through the 4 x 4 nearest // grid derivatives dx, dy (row-major, like the values z), which must all be trusted (tx, ty). As in interpolate_derivative, // the bicubic through z must predict f(qx, qy) = qz well enough that the condition number is good to about fd_tolerance. @@ -339,6 +307,11 @@ class ulps_heatmap { void write_ulps_envelope(std::string const & filename) const; + // The same SVGs as write_ulps and write_ulps_envelope, as strings. + std::string ulps_svg() const; + + std::string ulps_envelope_svg() const; + // For each panel, counts of pixels without a value and outside 1/2 ulp and the envelope, // and the worst errors and where they occur. std::string summary() const; @@ -366,9 +339,12 @@ class ulps_heatmap { return {probe_x_[source - 1], probe_y_[source - 1]}; } - void write_svg(std::string const & filename, std::vector panel::* values, double floor, + void write_svg(std::ostream & fs, std::vector panel::* values, double floor, std::string const & quantity) const; + void write_file(std::string const & filename, std::vector panel::* values, double floor, + std::string const & quantity) const; + std::vector xs_; std::vector ys_; // Row-major, with row 0 at y_max (the top of the image): @@ -669,13 +645,44 @@ ulps_heatmap& ulps_heatmap void ulps_heatmap::write_ulps(std::string const & filename) const { - write_svg(filename, &panel::ulps, 0.5, "|ulps|"); + write_file(filename, &panel::ulps, 0.5, "|ulps|"); } template void ulps_heatmap::write_ulps_envelope(std::string const & filename) const { - write_svg(filename, &panel::ulps_envelope, 1, "|ulps| / max(1/2, cond/2)"); + write_file(filename, &panel::ulps_envelope, 1, "|ulps| / max(1/2, cond/2)"); +} + +template +std::string ulps_heatmap::ulps_svg() const +{ + std::ostringstream os; + write_svg(os, &panel::ulps, 0.5, "|ulps|"); + return os.str(); +} + +template +std::string ulps_heatmap::ulps_envelope_svg() const +{ + std::ostringstream os; + write_svg(os, &panel::ulps_envelope, 1, "|ulps| / max(1/2, cond/2)"); + return os.str(); +} + +template +void ulps_heatmap::write_file(std::string const & filename, std::vector panel::* values, + double floor, std::string const & quantity) const +{ + if (!detail::ends_with(filename, ".svg")) + { + throw std::logic_error("Only svg files are supported at this time."); + } + // Render first, so that a failure leaves no partial file: + std::ostringstream os; + write_svg(os, values, floor, quantity); + std::ofstream fs(filename); + fs << os.str(); } template @@ -736,17 +743,13 @@ std::string ulps_heatmap::summary() const } template -void ulps_heatmap::write_svg(std::string const & filename, std::vector panel::* values, +void ulps_heatmap::write_svg(std::ostream & fs, std::vector panel::* values, double floor, std::string const & quantity) const { if (panels_.empty()) { throw std::domain_error("No function added; call add_fn first."); } - if (!detail::ends_with(filename, ".svg")) - { - throw std::logic_error("Only svg files are supported at this time."); - } if (!png_encoder_) { throw std::logic_error("No PNG encoder: set one with png_encoder, or define BOOST_MATH_USE_LODEPNG and link lodepng."); @@ -930,7 +933,6 @@ void ulps_heatmap::write_svg(std::string const & fil fn_png_ptr = &fn_png; } - std::ofstream fs(filename); fs << "\n" << "\n" << "\n" diff --git a/include/boost/math/tools/ulps_plot.hpp b/include/boost/math/tools/ulps_plot.hpp index 8c0c8670ae..55798e0ec3 100644 --- a/include/boost/math/tools/ulps_plot.hpp +++ b/include/boost/math/tools/ulps_plot.hpp @@ -21,6 +21,11 @@ #include #include #include +#include +#include +#include +#include +#include #endif #include #include @@ -92,6 +97,103 @@ inline std::string base64(std::vector const & bytes) return out; } +// Runs body(j) for j in [0, n) on several threads; all of the hardware's if threads is 0. If body throws, no further +// work is handed out, and the first exception is rethrown once the threads have finished. +template +void parallel_for(std::size_t n, unsigned threads, Body body) +{ + if (n == 0) + { + return; + } + if (threads == 0) + { + threads = (std::max)(1u, std::thread::hardware_concurrency()); + } + if (threads == 1 || n == 1) + { + for (std::size_t j = 0; j < n; ++j) + { + body(j); + } + return; + } + threads = static_cast((std::min)(static_cast(threads), n)); + std::atomic next{0}; + std::atomic failed{false}; + std::exception_ptr error; + std::mutex error_mutex; + auto work = [&]() { + try + { + for (std::size_t j = next++; j < n && !failed; j = next++) + { + body(j); + } + } + catch (...) + { + std::lock_guard lock(error_mutex); + if (!error) + { + error = std::current_exception(); + } + failed = true; + } + }; + std::vector pool; + try + { + for (unsigned t = 0; t < threads; ++t) + { + pool.emplace_back(work); + } + } + catch (...) + { + // Couldn't start a thread: stop the ones already running and report it. + failed = true; + for (auto & th : pool) + { + th.join(); + } + throw; + } + for (auto & th : pool) + { + th.join(); + } + if (error) + { + std::rethrow_exception(error); + } +} + +#ifdef BOOST_MATH_STANDALONE +// Uniform random numbers in [lo, hi) of a type that std::uniform_real_distribution may not support. +template::value> +class standalone_uniform +{ +public: + standalone_uniform(Real const & lo, Real const & hi) : lo_(lo), hi_(hi) {} + + template + Real operator()(Gen & gen) { return lo_ + (hi_ - lo_)*static_cast(unit_(gen)); } + +private: + Real lo_; + Real hi_; + std::uniform_real_distribution unit_{0.0L, 1.0L}; +}; + +template +class standalone_uniform : public std::uniform_real_distribution +{ +public: + using std::uniform_real_distribution::uniform_real_distribution; +}; +#endif + // ---- Condition numbers by finite differences on a grid ---- // How well a finite-difference derivative must agree with a lower-order one to be trusted, as a fraction of @@ -328,53 +430,59 @@ bool interpolate_derivative(Real const * x, Real const * f, Real const * d, unsi // px (anywhere in the grid's range) with values fpx, returned in pcond. Inf where f = 0, NaN where f is not finite. // The grid's come from finite differences of fx (grid_derivatives), and the probes' from interpolate_derivative; where these // can't be trusted, local_condition_number evaluates f at nearby points. If fallbacks isn't null, it receives the number of those. +// With threads other than 1, f must be thread safe. template std::vector grid_condition_numbers(F const & f, std::vector const & x, std::vector const & fx, std::vector const & px, std::vector const & fpx, std::vector & pcond, - std::size_t * fallbacks = nullptr) + std::size_t * fallbacks = nullptr, unsigned threads = 1) { using std::abs; using std::isfinite; std::size_t const n = x.size(); - std::size_t local = 0; + std::atomic local{0}; std::vector d(n); std::vector trusted(n); grid_derivatives(x.data(), fx.data(), 1, n, d.data(), trusted.data()); std::vector cond(n, std::numeric_limits::quiet_NaN()); - for (std::size_t j = 0; j < n; ++j) - { - if (!isfinite(fx[j])) - { - continue; - } - if (trusted[j]) - { - cond[j] = abs(x[j]*d[j]/fx[j]); - } - else + std::size_t const chunk = 64; + parallel_for((n + chunk - 1)/chunk, threads, [&](std::size_t c) { + for (std::size_t j = c*chunk; j < (std::min)(n, (c + 1)*chunk); ++j) { - cond[j] = local_condition_number(f, x[j], fx[j]); - ++local; + if (!isfinite(fx[j])) + { + continue; + } + if (trusted[j]) + { + cond[j] = abs(x[j]*d[j]/fx[j]); + } + else + { + cond[j] = local_condition_number(f, x[j], fx[j]); + ++local; + } } - } + }); pcond.assign(px.size(), std::numeric_limits::quiet_NaN()); - for (std::size_t j = 0; j < px.size(); ++j) - { - if (!isfinite(fpx[j])) - { - continue; - } - Real dj; - if (interpolate_derivative(x.data(), fx.data(), d.data(), trusted.data(), n, px[j], fpx[j], dj)) - { - pcond[j] = abs(px[j]*dj/fpx[j]); - } - else + parallel_for((px.size() + chunk - 1)/chunk, threads, [&](std::size_t c) { + for (std::size_t j = c*chunk; j < (std::min)(px.size(), (c + 1)*chunk); ++j) { - pcond[j] = local_condition_number(f, px[j], fpx[j]); - ++local; + if (!isfinite(fpx[j])) + { + continue; + } + Real dj; + if (interpolate_derivative(x.data(), fx.data(), d.data(), trusted.data(), n, px[j], fpx[j], dj)) + { + pcond[j] = abs(px[j]*dj/fpx[j]); + } + else + { + pcond[j] = local_condition_number(f, px[j], fpx[j]); + ++local; + } } - } + }); if (fallbacks) { *fallbacks = local; @@ -937,8 +1045,11 @@ class ulps_plot { // (seeded with random_seed, or randomly if it is -1). With perturb_abscissas, f is evaluated at the sample and the // implementation at the nearest CoarseReal; otherwise both at the CoarseReal. // With log_abscissas, the abscissas are sampled uniformly in log(x), and the x-axis is logarithmic; requires 0 < a. + // f, and the implementations passed to add_fn, are evaluated on this many threads (all of the hardware's if 0); with + // more than one, they must be thread safe. ulps_plot(F hi_acc_impl, CoarseReal a, CoarseReal b, - size_t samples = 1000, bool perturb_abscissas = false, int random_seed = -1, bool log_abscissas = false); + size_t samples = 1000, bool perturb_abscissas = false, int random_seed = -1, bool log_abscissas = false, + unsigned threads = 1); ulps_plot& clip(PreciseReal clip); @@ -1418,6 +1529,7 @@ class ulps_plot { CoarseReal a_; CoarseReal b_; bool log_abscissas_; + unsigned threads_; PreciseReal clip_; int width_; std::string envelope_color_; @@ -1622,7 +1734,8 @@ void ulps_plot::write(std::string const & filename) template ulps_plot::ulps_plot(F hi_acc_impl, CoarseReal a, CoarseReal b, - size_t samples, bool perturb_abscissas, int random_seed, bool log_abscissas) : crop_color_("red") + size_t samples, bool perturb_abscissas, int random_seed, bool log_abscissas, unsigned threads) + : threads_(threads), crop_color_("red") { // Use digits10 for this comparison in case the two types have differeing radixes: static_assert(std::numeric_limits::digits10 >= std::numeric_limits::digits10, "PreciseReal must have higher precision that CoarseReal"); @@ -1681,9 +1794,9 @@ ulps_plot::ulps_plot(F hi_acc_impl, CoarseReal a, Co #ifndef BOOST_MATH_STANDALONE boost::random::uniform_real_distribution dis(t_lo, t_hi); #else - // Use std::random in standalone mode if it is a type that the standard library can support (float, double, or long double) - static_assert(std::numeric_limits::digits10 <= std::numeric_limits::digits10, "Standalone mode does not support types with precision that exceeds long double"); - std::uniform_real_distribution dis(t_lo, t_hi); + // std's can only generate float, double and long double; other types get a long double in [0, 1), which leaves the + // last bits of a multiprecision probe zero (or those of t_lo), but the probes needn't be uniform to the last bit: + detail::standalone_uniform dis(t_lo, t_hi); #endif std::vector probe_x(probes); std::vector probe_coarse(probes); @@ -1707,17 +1820,16 @@ ulps_plot::ulps_plot(F hi_acc_impl, CoarseReal a, Co #endif }; std::vector grid_y(grid_size); - for (size_t i = 0; i < grid_size; ++i) - { - grid_y[i] = reference(grid_x[i]); - } std::vector probe_y(probes); - for (size_t i = 0; i < probes; ++i) - { - probe_y[i] = reference(probe_x[i]); - } + size_t const chunk = 64; + detail::parallel_for((samples + chunk - 1)/chunk, threads_, [&](size_t c) { + for (size_t i = c*chunk; i < (std::min)(samples, (c + 1)*chunk); ++i) + { + if (i < grid_size) { grid_y[i] = reference(grid_x[i]); } else { probe_y[i - grid_size] = reference(probe_x[i - grid_size]); } + } + }); std::vector probe_cond; - std::vector grid_cond = detail::grid_condition_numbers(reference, grid_x, grid_y, probe_x, probe_y, probe_cond); + std::vector grid_cond = detail::grid_condition_numbers(reference, grid_x, grid_y, probe_x, probe_y, probe_cond, nullptr, threads_); // Merge the grid and the probes, in order: std::vector order(samples); @@ -1765,28 +1877,36 @@ ulps_plot& ulps_plot::ad size_t samples = precise_abscissas_.size(); std::vector ulps(samples, std::numeric_limits::quiet_NaN()); std::vector classes(samples, detail::ulps_class::ok); + std::vector messages(samples); + size_t const chunk = 64; + detail::parallel_for((samples + chunk - 1)/chunk, threads_, [&](size_t ch) { + for (size_t i = ch*chunk; i < (std::min)(samples, (ch + 1)*chunk); ++i) + { + PreciseReal y_hi_acc = precise_ordinates_[i]; + CoarseReal w = 0; + // The implementation is evaluated everywhere, even where the true value is out of range, to check that + // overflow and underflow are handled correctly: + detail::ulps_class c = detail::classify_result(y_hi_acc, [&]() { return g(coarse_abscissas_[i]); }, w, &messages[i]); + classes[i] = c; + if (c != detail::ulps_class::ok) + { + continue; + } + PreciseReal y_lo_acc = static_cast(w); + PreciseReal absy = abs(y_hi_acc); + PreciseReal dist = static_cast(nextafter(static_cast(absy), (std::numeric_limits::max)()) - static_cast(absy)); + ulps[i] = static_cast((y_lo_acc - y_hi_acc)/dist); + } + }); std::string first_message; for (size_t i = 0; i < samples; ++i) { - PreciseReal y_hi_acc = precise_ordinates_[i]; - CoarseReal w = 0; - std::string message; - // The implementation is evaluated everywhere, even where the true value is out of range, to check that - // overflow and underflow are handled correctly: - detail::ulps_class c = detail::classify_result(y_hi_acc, [&]() { return g(coarse_abscissas_[i]); }, w, &message); - classes[i] = c; - if (c != detail::ulps_class::ok) + detail::ulps_class const c = classes[i]; + if (!messages[i].empty() && (c == detail::ulps_class::failed || c == detail::ulps_class::spurious_overflow || c == detail::ulps_class::mishandled)) { - if (first_message.empty() && !message.empty() && (c == detail::ulps_class::failed || c == detail::ulps_class::spurious_overflow || c == detail::ulps_class::mishandled)) - { - first_message = message; - } - continue; + first_message = messages[i]; + break; } - PreciseReal y_lo_acc = static_cast(w); - PreciseReal absy = abs(y_hi_acc); - PreciseReal dist = static_cast(nextafter(static_cast(absy), (std::numeric_limits::max)()) - static_cast(absy)); - ulps[i] = static_cast((y_lo_acc - y_hi_acc)/dist); } class_list_.emplace_back(classes); first_message_.emplace_back(first_message); diff --git a/nanobind/CMakeLists.txt b/nanobind/CMakeLists.txt index 2f492ecfb4..69c096565f 100644 --- a/nanobind/CMakeLists.txt +++ b/nanobind/CMakeLists.txt @@ -56,6 +56,28 @@ find_package(Eigen3 CONFIG QUIET) if(Eigen3_FOUND) target_link_libraries(_ext PRIVATE Eigen3::Eigen) endif() +# Optional: boost.math.ulps needs Boost.Multiprecision (standalone, which needs Boost.Config) and zlib. Only the +# sibling checkouts (or the cache variables) are searched, never a system Boost of another release. +find_path(BOOST_MATH_NB_MULTIPRECISION_INCLUDE boost/multiprecision/cpp_bin_float.hpp + HINTS "${CMAKE_CURRENT_SOURCE_DIR}/../../multiprecision/include" "$ENV{HOME}/boost/libs/multiprecision/include" + DOC "Boost.Multiprecision's include directory, for boost.math.ulps" NO_DEFAULT_PATH) +find_path(BOOST_MATH_NB_CONFIG_INCLUDE boost/config.hpp + HINTS "${CMAKE_CURRENT_SOURCE_DIR}/../../config/include" "$ENV{HOME}/boost/libs/config/include" + DOC "Boost.Config's include directory, for boost.math.ulps" NO_DEFAULT_PATH) +find_package(ZLIB) +if(BOOST_MATH_NB_MULTIPRECISION_INCLUDE AND BOOST_MATH_NB_CONFIG_INCLUDE AND ZLIB_FOUND) + set(ulps_sources src/ulps.cpp src/ulps_gamma.cpp src/ulps_special.cpp src/ulps_elliptic.cpp) + target_sources(_ext PRIVATE ${ulps_sources}) + # Only these see Multiprecision, so that its presence changes nothing else: + set_source_files_properties(${ulps_sources} PROPERTIES + INCLUDE_DIRECTORIES "${BOOST_MATH_NB_MULTIPRECISION_INCLUDE};${BOOST_MATH_NB_CONFIG_INCLUDE}" + COMPILE_DEFINITIONS BOOST_MP_STANDALONE) + set_source_files_properties(src/module.cpp PROPERTIES COMPILE_DEFINITIONS BMPY_HAVE_ULPS) + target_link_libraries(_ext PRIVATE ZLIB::ZLIB) +else() + message(STATUS "boost.math.ulps disabled: needs Boost.Multiprecision, Boost.Config and zlib " + "(set BOOST_MATH_NB_MULTIPRECISION_INCLUDE and BOOST_MATH_NB_CONFIG_INCLUDE)") +endif() if(NOT MSVC) # Match CI and test/Makefile: contraction allowed everywhere. target_compile_options(_ext PRIVATE -ffp-contract=fast) diff --git a/nanobind/boost/math/ulps.py b/nanobind/boost/math/ulps.py index 1d4736d1ed..efd5ca0edc 100644 --- a/nanobind/boost/math/ulps.py +++ b/nanobind/boost/math/ulps.py @@ -1,151 +1,130 @@ # Copyright 2026 Nick Thompson # Distributed under the Boost Software License, Version 1.0. # https://www.boost.org/LICENSE_1_0.txt -"""ulps plots against an mpmath reference. - - import boost.math as bm, mpmath as mp - from boost.math import ulps - ulps.plot(bm.digamma, mp.digamma, np.linspace(-5, 5, 2000)) - ulps.heatmap(bm.lgamma_p, lambda a, x: mp.log(mp.gammainc(a, 0, x, regularized=True)), - np.logspace(0, 7, 300), np.linspace(0.01, 0.99, 300), y_of=lambda a, r: a * r) - -Each point is classified as one of: - ulps |f - ref| / ulp(ref), both finite; - overflow f overflowed (raised or returned +-inf) where ref is finite: shown cyan; - nan f raised a domain error or returned NaN where ref is finite: shown green; - handled ref isn't a finite double and f agrees (inf, or underflow to 0): shown gray. +"""ulps plots and heatmaps from boost::math::tools::ulps_plot and ulps_heatmap; the reference is the same +function at cpp_bin_float_50, or reference= any Python callable (an mpmath function, say). For example: + + import boost.math as bm; from boost.math import ulps + ulps.plot(bm.digamma, -5, 5).save("digamma.svg") + ulps.render_to_png(ulps.heatmap(bm.lgamma_p, (1, 1e7), (0.01, 1e7), log_x=True, log_y=True)[0], "lgamma_p.png") """ import math +import os +import re +import shutil +import subprocess +import tempfile + +from . import _ext + + +class Svg: + """An SVG: str() is its text, summary the C++ summary(); shown inline in Jupyter.""" + + def __init__(self, text, summary): + self.text = text + self.summary = summary -import mpmath as mp -import numpy as np + def __str__(self): + return self.text -OVERFLOW, NAN, HANDLED = -1.0, -2.0, -3.0 -SPECIAL = {OVERFLOW: ("#2ec4d6", "spurious overflow"), NAN: ("#4caf50", "NaN / domain error"), - HANDLED: ("#8a8a85", "out of range, handled")} -SURFACE, INK, MUTED = "#1a1a19", "#f0efec", "#c3c2b7" + def _repr_svg_(self): + return self.text + def save(self, path): + with open(path, "w", encoding="utf-8") as f: + f.write(self.text) + return path -def error(f, ref, *args, dps=50): - """The error of f(*args) against the mpmath reference ref(*args): ulps, or a SPECIAL code.""" - with mp.workdps(dps): + +def _name(f): + if not hasattr(_ext, "_ulps_plot"): + raise RuntimeError("boost.math.ulps was built without Boost.Multiprecision or zlib; see nanobind/CMakeLists.txt") + if isinstance(f, str): + return f + if getattr(_ext, getattr(f, "__name__", ""), None) is not f: + raise TypeError("ulps plots need a boost.math function or its name; arbitrary Python callables are not supported") + return f.__name__ + + +def _python_reference(reference): + """Wraps reference for the C++ side: arguments come as floats, or as 55-digit strings (made mpmath numbers here) + where they aren't doubles, and the result goes back as a float or a decimal string; NaN if it raises or isn't real.""" + try: + import mpmath as mp + except ImportError: + raise ImportError("reference= needs mpmath: pip install mpmath") from None + + def call(*args): try: - r = ref(*[mp.mpf(a) for a in args]) + r = reference(*[mp.mpf(a) if isinstance(a, str) else a for a in args]) + if isinstance(r, mp.mpf): + return mp.nstr(r, 55) if mp.isfinite(r) else float(r) + return float(r) except Exception: - return math.nan # No reference here; left blank in the plots. - try: - got = f(*args) - except OverflowError: - got = math.inf - except (ValueError, ArithmeticError): - got = math.nan - if isinstance(r, mp.mpc) or not mp.isfinite(r): - return HANDLED if (math.isinf(got) or math.isnan(got)) else float("inf") - rf = float(r) - if math.isinf(rf) or (rf == 0 and r != 0): - # The answer isn't representable: overflow to inf or underflow to 0 is right. - return HANDLED if (got == rf or (rf == 0 and abs(got) <= np.finfo(float).tiny)) else float( - abs(mp.mpf(got) - r) / np.spacing(np.finfo(float).tiny)) - if math.isinf(got): - return OVERFLOW - if math.isnan(got): - return NAN - return float(abs(mp.mpf(got) - r) / np.spacing(abs(rf)) if rf != 0 else abs(got) / np.spacing(0.0)) - - -def _summary(err): - e = err[err >= 0] - missing = int(np.sum(np.isnan(err))) - counts = {name: int(np.sum(err == code)) for code, (_, name) in SPECIAL.items() if np.any(err == code)} - parts = [f"max {e.max():.3g} ulps, median {np.median(e):.2g}" if e.size else "no finite errors"] - parts += [f"{n} {name}" for name, n in counts.items()] - if missing: - parts.append(f"{missing} without a reference") - return "; ".join(parts) - - -def _style(ax, fig): - fig.set_facecolor(SURFACE) - ax.set_facecolor(SURFACE) - ax.tick_params(colors=MUTED) - for s in ax.spines.values(): - s.set_color(MUTED) - - -def _decades(cb): - from matplotlib.ticker import FuncFormatter - cb.ax.yaxis.set_major_formatter(FuncFormatter(lambda v, _: "$10^{%d}$" % round(np.log10(v)))) - cb.ax.tick_params(colors=MUTED) - cb.outline.set_edgecolor(MUTED) - cb.set_label("error in ulps", color=INK) - - -def plot(f, ref, x, title=None, xlabel="x", xscale="linear", out=None, dps=50, reference=None): - """1D: error in ulps of f(x) against ref(x). Returns (fig, err).""" - import matplotlib.pyplot as plt - x = np.asarray(x, dtype=float) - err = np.array([error(f, ref, xi, dps=dps) for xi in x]) - fig, ax = plt.subplots(figsize=(8, 4.5)) - _style(ax, fig) - ok = err >= 0 - ax.scatter(x[ok], np.maximum(err[ok], 0.01), s=3, color="#f5a524", rasterized=True) - top = np.max(err[ok]) if ok.any() else 1 - for code, (color, name) in SPECIAL.items(): - m = err == code - if m.any(): - ax.scatter(x[m], np.full(m.sum(), max(top, 1) * 2), s=6, color=color, label=name, rasterized=True) - ax.set_xscale(xscale) - ax.set_yscale("log") - ax.set_xlabel(xlabel, color=INK) - ax.set_ylabel("error in ulps", color=INK) - ax.set_title(title or getattr(f, "__name__", ""), color=INK) - if any(np.any(err == c) for c in SPECIAL): - ax.legend(frameon=False, labelcolor=INK) - fig.text(0.5, 0.01, f"Reference: {reference or f'mpmath at {dps} digits'}; {_summary(err)}.", ha="center", - color=MUTED, fontsize=8, wrap=True) - fig.tight_layout(rect=(0, 0.03, 1, 1)) - if out: - fig.savefig(out, dpi=300, facecolor=SURFACE) - return fig, err - - -def heatmap(f, ref, x, y, title=None, xlabel="x", ylabel="y", xscale="log", yscale="log", - y_of=None, out=None, dps=50, reference=None): - """2D: error in ulps of f(x, y) over the grid x by y. With y_of, the second argument is - y_of(x, y) instead of y (e.g. y_of=lambda a, r: a * r to plot against x/a). Returns (fig, err).""" - import matplotlib.pyplot as plt - from matplotlib.colors import ListedColormap, LogNorm - x = np.asarray(x, dtype=float) - y = np.asarray(y, dtype=float) - err = np.empty((y.size, x.size)) - for i, yi in enumerate(y): - for j, xj in enumerate(x): - err[i, j] = error(f, ref, xj, y_of(xj, yi) if y_of else yi, dps=dps) - fig, ax = plt.subplots(figsize=(7.2, 6)) - _style(ax, fig) - finite = np.where(err >= 0, np.maximum(err, 0.1), np.nan) - vmax = max(np.nanmax(finite) if np.isfinite(finite).any() else 1, 1) - im = ax.pcolormesh(x, y, finite, norm=LogNorm(vmin=0.1, vmax=vmax), cmap="inferno", shading="nearest", - rasterized=True) - handles = [] - for code, (color, name) in SPECIAL.items(): - m = err == code - if m.any(): - ax.pcolormesh(x, y, np.where(m, 1.0, np.nan), cmap=ListedColormap([color]), shading="nearest", - rasterized=True) - handles.append(plt.Rectangle((0, 0), 1, 1, color=color, label=name)) - ax.set_xscale(xscale) - ax.set_yscale(yscale) - ax.set_xlabel(xlabel, color=INK, fontsize=12) - ax.set_ylabel(ylabel, color=INK, fontsize=12) - ax.set_title(title or getattr(f, "__name__", ""), color=INK, fontsize=12) - if handles: - ax.legend(handles=handles, frameon=False, labelcolor=INK, loc="upper left", fontsize=8) - _decades(fig.colorbar(im, ax=ax)) - fig.text(0.5, 0.01, f"Reference: {reference or f'mpmath at {dps} digits'}; {y.size} x {x.size} grid.\n" - f"{_summary(err)}; errors below 0.1 ulp shown as 0.1.", ha="center", color=MUTED, fontsize=8) - fig.tight_layout(rect=(0, 0.05, 1, 1)) - if out: - fig.savefig(out, dpi=300, facecolor=SURFACE) - return fig, err + return math.nan + + return call + + +def _run(fn, reference, threads): + if reference is None: + return fn(None, threads) + import mpmath as mp + call = _python_reference(reference) + # The Python calls are serialized anyway; more threads only add GIL churn: + with mp.workdps(60): + return fn(call, 1) + + +def plot(f, a, b, *, samples=10000, log_x=False, title=None, reference=None, clip=None, show_function=False, + threads=0): + """The ulps plot of f on [a, b]; threads=0 uses them all. reference= an mpmath-based callable replaces + the cpp_bin_float_50 reference, and runs on one thread.""" + name = _name(f) + svg, summary = _run(lambda ref, th: _ext._ulps_plot(name, a, b, samples, log_x, name if title is None else title, ref, + clip or -1, show_function, th), reference, threads) + return Svg(svg, summary) + + +def heatmap(f, x_range, y_range, *, width=800, height=500, log_x=False, log_y=False, title=None, x_label="x", + y_label="y", log_scale=False, clip=None, reference=None, show_function=False, threads=0): + """The ulps heatmap of f(x, y) over x_range by y_range: (ulps, ulps/envelope), two Svgs. reference= as for plot.""" + name = _name(f) + ulps, envelope, summary = _run( + lambda ref, th: _ext._ulps_heatmap(name, *x_range, *y_range, width, height, log_x, log_y, + name if title is None else title, x_label, y_label, log_scale, clip or -1, + ref, show_function, th), reference, threads) + return Svg(ulps, summary), Svg(envelope, summary) + + +def _chrome(): + candidates = [os.environ.get("CHROME"), "/Applications/Google Chrome.app/Contents/MacOS/Google Chrome"] + candidates += [shutil.which(c) for c in ("google-chrome", "google-chrome-stable", "chromium", "chromium-browser")] + return next((c for c in candidates if c and os.path.isfile(c)), None) + + +def render_to_png(svg, path, *, width=None, height=None): + """Renders an Svg (or SVG text) to a PNG with headless Chrome, or cairosvg. Returns path.""" + text = str(svg) + size = lambda attr: int(float(re.search(r"]*?\b" + attr + r"=['\"]([0-9.]+)", text).group(1))) + width = width or size("width") + height = height or size("height") + path = os.path.abspath(path) + chrome = _chrome() + if chrome: + with tempfile.TemporaryDirectory() as tmp: + src = os.path.join(tmp, "plot.svg") + with open(src, "w", encoding="utf-8") as f: + f.write(text) + subprocess.run([chrome, "--headless", "--disable-gpu", "--hide-scrollbars", f"--screenshot={path}", + f"--window-size={width},{height}", "file://" + src], + check=True, capture_output=True, timeout=120) + return path + try: + import cairosvg + except ImportError: + raise RuntimeError("render_to_png needs Google Chrome (or set $CHROME) or `pip install cairosvg`") from None + cairosvg.svg2png(bytestring=text.encode(), write_to=path, output_width=width, output_height=height) + return path diff --git a/nanobind/src/bindings.hpp b/nanobind/src/bindings.hpp index 44f65ace00..b57873f331 100644 --- a/nanobind/src/bindings.hpp +++ b/nanobind/src/bindings.hpp @@ -26,6 +26,8 @@ void bind_fft(nb::module_& m); void bind_quadrature(nb::module_& m); void bind_optimization(nb::module_& m); void bind_polynomial_roots(nb::module_& m); +// Only with Boost.Multiprecision and zlib; see CMakeLists.txt. +void bind_ulps(nb::module_& m); } // namespace bmpy diff --git a/nanobind/src/module.cpp b/nanobind/src/module.cpp index 22b9447e1c..75b8d10ab2 100644 --- a/nanobind/src/module.cpp +++ b/nanobind/src/module.cpp @@ -23,4 +23,7 @@ NB_MODULE(_ext, m) { bmpy::bind_quadrature(m); bmpy::bind_optimization(m); bmpy::bind_polynomial_roots(m); +#ifdef BMPY_HAVE_ULPS + bmpy::bind_ulps(m); +#endif } diff --git a/nanobind/src/ulps.cpp b/nanobind/src/ulps.cpp new file mode 100644 index 0000000000..1d1e266099 --- /dev/null +++ b/nanobind/src/ulps.cpp @@ -0,0 +1,233 @@ +// Copyright 2026 Nick Thompson +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt +// +// ulps plots and heatmaps for boost.math.ulps: the double implementation against the same +// function at cpp_bin_float_50, or against a Python reference. + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "bindings.hpp" +#include "ulps.hpp" + +namespace nb = bmpy::nb; + +namespace { + +using bmpy::ulps::Precise; +using plot1 = boost::math::tools::ulps_plot, Precise, double>; +using heatmap2 = boost::math::tools::ulps_heatmap, Precise, double>; + +bmpy::ulps::registry const& functions() { + static bmpy::ulps::registry const r = [] { + bmpy::ulps::registry r; + bmpy::ulps::register_gamma(r); + bmpy::ulps::register_special(r); + bmpy::ulps::register_elliptic(r); + return r; + }(); + return r; +} + +template +auto const& lookup(Map const& m, std::string const& name) { + auto it = m.find(name); + if (it == m.end()) { + throw nb::value_error(("no ulps reference for '" + name + "'").c_str()); + } + return it->second; +} + +// An RGB8 PNG: one IDAT, deflated by zlib. +std::vector encode_png(std::vector const& rgb, std::size_t width, std::size_t height) { + std::vector raw; + raw.reserve((3 * width + 1) * height); + for (std::size_t j = 0; j < height; ++j) { + raw.push_back(0); // filter: none + raw.insert(raw.end(), rgb.begin() + 3 * width * j, rgb.begin() + 3 * width * (j + 1)); + } + uLongf size = compressBound(static_cast(raw.size())); + std::vector deflated(size); + if (compress2(deflated.data(), &size, raw.data(), static_cast(raw.size()), Z_DEFAULT_COMPRESSION) != Z_OK) { + throw std::runtime_error("zlib: compress2 failed"); + } + deflated.resize(size); + + std::vector png = {0x89, 'P', 'N', 'G', '\r', '\n', 0x1a, '\n'}; + auto be32 = [](std::vector& v, std::uint32_t n) { + for (int s = 24; s >= 0; s -= 8) { + v.push_back(static_cast(n >> s)); + } + }; + auto chunk = [&](char const* type, std::vector const& data) { + be32(png, static_cast(data.size())); + std::size_t const start = png.size(); + png.insert(png.end(), type, type + 4); + png.insert(png.end(), data.begin(), data.end()); + be32(png, static_cast(crc32(0, png.data() + start, static_cast(png.size() - start)))); + }; + std::vector header; + be32(header, static_cast(width)); + be32(header, static_cast(height)); + header.insert(header.end(), {8, 2, 0, 0, 0}); // 8 bits, RGB, deflate, adaptive filtering, no interlace + chunk("IHDR", header); + chunk("IDAT", deflated); + chunk("IEND", {}); + return png; +} + +// A Python reference, called one at a time: not merely one holder of the GIL at a time, since the interpreter switches +// threads in the middle of a call, and mpmath keeps its working precision in a global, which its functions raise and +// restore. Arguments go as floats where they are doubles, and otherwise (the condition number's finite differences) as +// 55-digit strings, which boost.math.ulps turns into mpmath numbers. It returns a float or a decimal string; anything +// else, or an exception, is NaN. KeyboardInterrupt, SystemExit and a failed signal handler stop the run instead: every +// later call returns NaN at once, and rethrow() raises the error once the computation is over. +class python_reference { +public: + explicit python_reference(nb::handle f) : f_(f) {} + + template + Precise operator()(Args const&... args) { + std::lock_guard lock(mutex_); + if (abort_) { + return std::numeric_limits::quiet_NaN(); + } + nb::gil_scoped_acquire acquire; + bool signal = false; + try { + signal = PyErr_CheckSignals() != 0; + if (signal) { + throw nb::python_error(); + } + return from_python(f_(to_python(args)...)); + } catch (nb::python_error& e) { + if (!signal && !e.matches(PyExc_KeyboardInterrupt) && !e.matches(PyExc_SystemExit)) { + return std::numeric_limits::quiet_NaN(); + } + abort_ = true; + error_ = std::make_unique(std::move(e)); + return std::numeric_limits::quiet_NaN(); + } catch (std::exception const&) { + return std::numeric_limits::quiet_NaN(); + } + } + + // With the GIL held. + void rethrow() { + if (error_) { + nb::python_error e = std::move(*error_); + error_.reset(); + throw e; + } + } + +private: + static nb::object to_python(Precise const& x) { + double const d = static_cast(x); + if (Precise(d) == x) { + return nb::float_(d); + } + return nb::str(x.str(55, std::ios_base::scientific).c_str()); + } + + static Precise from_python(nb::handle r) { + if (nb::isinstance(r) || nb::isinstance(r)) { + return Precise(nb::cast(r)); + } + if (nb::isinstance(r)) { + return Precise(nb::cast(r).c_str()); + } + return std::numeric_limits::quiet_NaN(); + } + + nb::handle f_; + std::mutex mutex_; + bool abort_ = false; + std::unique_ptr error_; +}; + +std::tuple plot(std::string const& name, double a, double b, std::size_t samples, bool log_x, + std::string const& title, nb::handle reference, double clip, bool show_function, + unsigned threads) { + auto const& f = lookup(functions().unary, name); + std::function precise = f.precise; + python_reference py(reference); + if (!reference.is_none()) { + precise = [&py](Precise const& x) { return py(x); }; + } + std::string svg; + std::string summary; + { + nb::gil_scoped_release release; + plot1 p(precise, a, b, samples, false, 20261013, log_x, threads); + p.add_fn(f.coarse).title(title).show_function(show_function).png_encoder(encode_png); + if (clip > 0) { + p.clip(clip); + } + std::ostringstream os; + os << p; + svg = os.str(); + summary = p.summary(); + } + py.rethrow(); + return {svg, summary}; +} + +std::tuple heatmap(std::string const& name, double x_min, double x_max, double y_min, + double y_max, std::size_t width, std::size_t height, bool log_x, + bool log_y, std::string const& title, std::string const& x_label, + std::string const& y_label, bool log_scale, double clip, + nb::handle reference, bool show_function, unsigned threads) { + auto const& f = lookup(functions().binary, name); + std::function precise = f.precise; + python_reference py(reference); + if (!reference.is_none()) { + precise = [&py](Precise const& x, Precise const& y) { return py(x, y); }; + } + std::tuple result; + { + nb::gil_scoped_release release; + heatmap2 h(precise, x_min, x_max, y_min, y_max, width, height, log_x, log_y, threads); + h.add_fn(f.coarse).title(title).x_label(x_label).y_label(y_label).log_scale(log_scale).show_function(show_function) + .png_encoder(encode_png); + if (clip > 0) { + h.clip(clip); + } + result = {h.ulps_svg(), h.ulps_envelope_svg(), h.summary()}; + } + py.rethrow(); + return result; +} + +} // namespace + +void bmpy::bind_ulps(nb::module_& m) { + m.def("_ulps_plot", &plot, nb::arg("name"), nb::arg("a"), nb::arg("b"), nb::arg("samples"), nb::arg("log_x"), + nb::arg("title"), nb::arg("reference").none(), nb::arg("clip"), nb::arg("show_function"), nb::arg("threads"), + "(svg, summary) of an ulps_plot; see boost.math.ulps.plot"); + m.def("_ulps_heatmap", &heatmap, nb::arg("name"), nb::arg("x_min"), nb::arg("x_max"), nb::arg("y_min"), + nb::arg("y_max"), nb::arg("width"), nb::arg("height"), nb::arg("log_x"), nb::arg("log_y"), nb::arg("title"), + nb::arg("x_label"), nb::arg("y_label"), nb::arg("log_scale"), nb::arg("clip"), nb::arg("reference").none(), + nb::arg("show_function"), nb::arg("threads"), + "(ulps svg, ulps/envelope svg, summary) of an ulps_heatmap; see boost.math.ulps.heatmap"); + m.def("_ulps_functions", [] { + std::vector unary, binary; + for (auto const& kv : functions().unary) unary.push_back(kv.first); + for (auto const& kv : functions().binary) binary.push_back(kv.first); + return std::make_tuple(unary, binary); + }, "The names boost.math.ulps.plot and heatmap accept."); +} diff --git a/nanobind/src/ulps.hpp b/nanobind/src/ulps.hpp new file mode 100644 index 0000000000..dd73efd99a --- /dev/null +++ b/nanobind/src/ulps.hpp @@ -0,0 +1,59 @@ +// Copyright 2026 Nick Thompson +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt +// +// The functions boost.math.ulps can plot: each at double (the implementation) and at +// cpp_bin_float_50 (the reference), type-erased so that ulps_plot and ulps_heatmap are +// instantiated once per arity. + +#ifndef BOOST_MATH_NANOBIND_ULPS_HPP +#define BOOST_MATH_NANOBIND_ULPS_HPP + +#include + +#include +#include +#include + +namespace bmpy::ulps { + +using Precise = boost::multiprecision::cpp_bin_float_50; + +struct fn1 { + std::function coarse; + std::function precise; +}; + +struct fn2 { + std::function coarse; + std::function precise; +}; + +struct registry { + std::map unary; + std::map binary; +}; + +template +fn1 make1(F f) { + return {[f](double x) { return static_cast(f(x)); }, + [f](Precise const& x) { return static_cast(f(x)); }}; +} + +template +fn2 make2(F f) { + return {[f](double x, double y) { return static_cast(f(x, y)); }, + [f](Precise const& x, Precise const& y) { return static_cast(f(x, y)); }}; +} + +// One entry per function; fn is called through bm::. +#define BMPY_ULPS1(fn) r.unary.emplace(#fn, bmpy::ulps::make1([](auto x) { return bm::fn(x); })) +#define BMPY_ULPS2(fn) r.binary.emplace(#fn, bmpy::ulps::make2([](auto x, auto y) { return bm::fn(x, y); })) + +void register_gamma(registry& r); +void register_special(registry& r); +void register_elliptic(registry& r); + +} // namespace bmpy::ulps + +#endif diff --git a/nanobind/src/ulps_elliptic.cpp b/nanobind/src/ulps_elliptic.cpp new file mode 100644 index 0000000000..cd9b976a61 --- /dev/null +++ b/nanobind/src/ulps_elliptic.cpp @@ -0,0 +1,56 @@ +// Copyright 2026 Nick Thompson +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "ulps.hpp" + +namespace bm = boost::math; + +void bmpy::ulps::register_elliptic(registry& r) { + BMPY_ULPS1(ellint_1); + BMPY_ULPS2(ellint_1); + BMPY_ULPS1(ellint_2); + BMPY_ULPS2(ellint_2); + BMPY_ULPS2(ellint_3); + BMPY_ULPS1(ellint_d); + BMPY_ULPS2(ellint_d); + BMPY_ULPS2(ellint_rc); + BMPY_ULPS2(jacobi_zeta); + BMPY_ULPS2(heuman_lambda); + + BMPY_ULPS2(jacobi_sn); + BMPY_ULPS2(jacobi_cn); + BMPY_ULPS2(jacobi_dn); + BMPY_ULPS2(jacobi_cd); + BMPY_ULPS2(jacobi_dc); + BMPY_ULPS2(jacobi_ns); + BMPY_ULPS2(jacobi_sd); + BMPY_ULPS2(jacobi_ds); + BMPY_ULPS2(jacobi_nc); + BMPY_ULPS2(jacobi_nd); + BMPY_ULPS2(jacobi_sc); + BMPY_ULPS2(jacobi_cs); + + BMPY_ULPS2(jacobi_theta1); + BMPY_ULPS2(jacobi_theta2); + BMPY_ULPS2(jacobi_theta3); + BMPY_ULPS2(jacobi_theta4); + BMPY_ULPS2(jacobi_theta3m1); + BMPY_ULPS2(jacobi_theta4m1); + BMPY_ULPS2(jacobi_theta1tau); + BMPY_ULPS2(jacobi_theta2tau); + BMPY_ULPS2(jacobi_theta3tau); + BMPY_ULPS2(jacobi_theta4tau); + BMPY_ULPS2(jacobi_theta3m1tau); + BMPY_ULPS2(jacobi_theta4m1tau); +} diff --git a/nanobind/src/ulps_gamma.cpp b/nanobind/src/ulps_gamma.cpp new file mode 100644 index 0000000000..6e886f949d --- /dev/null +++ b/nanobind/src/ulps_gamma.cpp @@ -0,0 +1,38 @@ +// Copyright 2026 Nick Thompson +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include +#include +#include + +#include "ulps.hpp" + +namespace bm = boost::math; + +void bmpy::ulps::register_gamma(registry& r) { + BMPY_ULPS1(tgamma); + BMPY_ULPS1(tgamma1pm1); + BMPY_ULPS1(lgamma); + BMPY_ULPS1(digamma); + BMPY_ULPS1(trigamma); + + BMPY_ULPS2(tgamma); // upper incomplete + BMPY_ULPS2(tgamma_ratio); + BMPY_ULPS2(tgamma_delta_ratio); + BMPY_ULPS2(tgamma_lower); + BMPY_ULPS2(gamma_p); + BMPY_ULPS2(gamma_q); + BMPY_ULPS2(gamma_p_inv); + BMPY_ULPS2(gamma_q_inv); + BMPY_ULPS2(gamma_p_inva); + BMPY_ULPS2(gamma_q_inva); + BMPY_ULPS2(gamma_p_derivative); + BMPY_ULPS2(lgamma_p); + BMPY_ULPS2(lgamma_q); + BMPY_ULPS2(ligamma); + BMPY_ULPS2(ligamma_lower); + BMPY_ULPS2(beta); + BMPY_ULPS2(lbeta); +} diff --git a/nanobind/src/ulps_special.cpp b/nanobind/src/ulps_special.cpp new file mode 100644 index 0000000000..335ed59482 --- /dev/null +++ b/nanobind/src/ulps_special.cpp @@ -0,0 +1,86 @@ +// Copyright 2026 Nick Thompson +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "ulps.hpp" + +namespace bm = boost::math; + +void bmpy::ulps::register_special(registry& r) { + BMPY_ULPS1(erf); + BMPY_ULPS1(erfc); + BMPY_ULPS1(erf_inv); + BMPY_ULPS1(erfc_inv); + BMPY_ULPS1(expint); + BMPY_ULPS1(zeta); + BMPY_ULPS1(lambert_w0); + BMPY_ULPS1(lambert_wm1); + BMPY_ULPS1(lambert_w0_prime); + BMPY_ULPS1(lambert_wm1_prime); + BMPY_ULPS2(owens_t); + + BMPY_ULPS1(log1p); + BMPY_ULPS1(log1pmx); + BMPY_ULPS1(expm1); + BMPY_ULPS1(cbrt); + BMPY_ULPS1(sqrt1pm1); + BMPY_ULPS1(sin_pi); + BMPY_ULPS1(cos_pi); + BMPY_ULPS1(sinc_pi); + BMPY_ULPS1(sinhc_pi); + BMPY_ULPS1(asinh); + BMPY_ULPS1(acosh); + BMPY_ULPS1(atanh); + BMPY_ULPS1(rsqrt); + BMPY_ULPS1(logit); + BMPY_ULPS1(logistic_sigmoid); + BMPY_ULPS2(powm1); + BMPY_ULPS2(pow1p); + BMPY_ULPS2(hypot); + BMPY_ULPS2(logaddexp); + + BMPY_ULPS2(cyl_bessel_j); + BMPY_ULPS2(cyl_neumann); + BMPY_ULPS2(cyl_bessel_i); + BMPY_ULPS2(cyl_bessel_k); + BMPY_ULPS2(cyl_bessel_j_prime); + BMPY_ULPS2(cyl_neumann_prime); + BMPY_ULPS2(cyl_bessel_i_prime); + BMPY_ULPS2(cyl_bessel_k_prime); + BMPY_ULPS1(airy_ai); + BMPY_ULPS1(airy_bi); + BMPY_ULPS1(airy_ai_prime); + BMPY_ULPS1(airy_bi_prime); + + BMPY_ULPS2(hypergeometric_1F0); + BMPY_ULPS2(hypergeometric_0F1); +} diff --git a/nanobind/tests/test_ulps.py b/nanobind/tests/test_ulps.py index c0761cbec0..b1302c94b9 100644 --- a/nanobind/tests/test_ulps.py +++ b/nanobind/tests/test_ulps.py @@ -1,18 +1,76 @@ # Copyright 2026 Nick Thompson # Distributed under the Boost Software License, Version 1.0. # https://www.boost.org/LICENSE_1_0.txt -import math +import importlib.util import mpmath as mp +import pytest import boost.math as bm from boost.math import ulps +pytestmark = pytest.mark.skipif(not hasattr(bm._ext, "_ulps_plot"), reason="built without boost.math.ulps") -def test_error_classification(): - assert ulps.error(bm.tgamma, mp.gamma, 4.5) < 4 - assert ulps.error(lambda x: 7.0, lambda x: mp.mpf(7) + mp.mpf(2) ** -50, 0.0) == 1 # ulp(7) = 2^-50 - assert ulps.error(lambda x: math.inf, lambda x: mp.mpf(1), 0.0) == ulps.OVERFLOW - assert ulps.error(lambda x: math.nan, lambda x: mp.mpf(1), 0.0) == ulps.NAN - assert ulps.error(bm.tgamma, mp.gamma, 200.0) == ulps.HANDLED - assert math.isnan(ulps.error(bm.tgamma, lambda x: 1 / 0, 2.0)) + +def check(svg): + assert str(svg).startswith((" 10000 + assert svg.summary + + +def test_plot(tmp_path): + svg = ulps.plot(bm.digamma, -5, 5, samples=2000) + check(svg) + assert "samples=2000" in svg.summary + assert svg._repr_svg_() == str(svg) + path = tmp_path / "digamma.svg" + svg.save(path) + assert path.read_text() == str(svg) + + +def test_heatmap(): + u, e = ulps.heatmap("lgamma_p", (1, 1e3), (0.01, 1e3), width=40, height=30, log_x=True, log_y=True) + check(u) + check(e) + assert str(u) != str(e) + assert "1200 pixels" in u.summary + + +def test_python_reference(): + svg = ulps.plot(bm.erf, -3, 3, samples=200, reference=mp.erf) + check(svg) + assert "worst_ulps=" in svg.summary + # A reference that raises leaves the points without a value: + svg = ulps.plot(bm.erf, -3, 3, samples=50, reference=lambda x: 1 / 0) + assert "undefined=50" in svg.summary + + +def test_reference_interrupt(): + def ref(x): + raise KeyboardInterrupt + + with pytest.raises(KeyboardInterrupt): + ulps.plot(bm.erf, -3, 3, samples=50, reference=ref) + + +def test_only_bound_functions(): + def digamma(x): + return bm.digamma(x) + + with pytest.raises(TypeError, match="arbitrary Python callables"): + ulps.plot(digamma, 1, 2) + + +def test_unknown_name(): + with pytest.raises(ValueError, match="no ulps reference for 'nonesuch'"): + ulps.plot("nonesuch", 0, 1) + with pytest.raises(ValueError, match="no ulps reference for 'digamma'"): + ulps.heatmap(bm.digamma, (0, 1), (0, 1)) + + +@pytest.mark.skipif(ulps._chrome() is None and importlib.util.find_spec("cairosvg") is None, + reason="needs Chrome or cairosvg") +def test_render_to_png(tmp_path): + path = ulps.render_to_png(ulps.plot(bm.digamma, 1, 2, samples=100), tmp_path / "digamma.png") + with open(path, "rb") as f: + assert f.read(8) == b"\x89PNG\r\n\x1a\n" diff --git a/test/test_ulps_plot.cpp b/test/test_ulps_plot.cpp index 3ff812076e..129dc0d411 100644 --- a/test/test_ulps_plot.cpp +++ b/test/test_ulps_plot.cpp @@ -392,6 +392,28 @@ void test_files() // The summary counts the grid and the probes: CHECK_LE(std::size_t(1), count(plot.summary(), "samples=500 ")); + // Threads change nothing: + { + ulps_plot serial(f, 0.5f, 6.0f, 500, false, 3, false, 1); + ulps_plot threaded(f, 0.5f, 6.0f, 500, false, 3, false, 4); + serial.add_fn(g).add_fn(h); + threaded.add_fn(g).add_fn(h); + CHECK_TRUE(serial.summary() == threaded.summary()); + } + // An exception on a worker thread reaches the caller: + { + bool thrown = false; + try + { + detail::parallel_for(100, 4, [](std::size_t j) { if (j == 37) { throw std::runtime_error("37"); } }); + } + catch (std::runtime_error const & e) + { + thrown = std::string(e.what()) == "37"; + } + CHECK_TRUE(thrown); + } + // A logarithmic axis: ulps_plot log_plot(f, 0.01f, 100.0f, 200, true, 3, true); log_plot.add_fn(g).png_encoder(nullptr); @@ -409,7 +431,8 @@ void test_files() std::string const heat_name = "test_ulps_heatmap.svg"; heat.write_ulps_envelope(heat_name); check_svg(read_file(heat_name)); - CHECK_EQUAL(widths.size(), std::size_t(1)); + CHECK_TRUE(read_file(heat_name) == heat.ulps_envelope_svg()); + CHECK_EQUAL(widths.size(), std::size_t(2)); CHECK_EQUAL(widths[0], std::size_t(40)); CHECK_EQUAL(heights[0], std::size_t(30)); std::remove(heat_name.c_str());