Basic domain

Fix up traits

Wip up down

Tests pass!

Cleanup

Enable full test range

More tests

Code cleanup

More tests, all pass

Cleanup, passing

Add visualization scripts

Add visualization scripts

updates to script

clippy

clippy

clippy

Docs update

docs

docs

First pass at wiring things up

todo -> unimplemented for f8 since we probably won't implement anything

Author: tgross35

.gitignore

@@ -7,3 +7,6 @@
 Cargo.lock
 musl/
 **.tar.gz
+
+# Files generated by binaries
+build/

crates/libm-test/examples/plot_domains.rs

@@ -0,0 +1,58 @@
+//! Program to write all inputs from a generator to a file, then invoke a Julia script
+//! to plot them. Requires Julia with the `CairoMakie` dependency.
+//!
+//! Note that running in release mode by default generates a _lot_ more datapoints, which
+//! causes plotting to be extremely slow (some simplification to be done in the script).
+
+use std::io::{BufWriter, Write};
+use std::path::Path;
+use std::process::Command;
+use std::{env, fs};
+
+use libm_test::domain::{Domain, SqrtDomain, TrigDomain, Unbounded};
+use libm_test::gen::domain;
+
+const JL_PLOT: &str = "examples/plot_file.jl";
+
+fn main() {
+    let out_dir = Path::new("build");
+    if !out_dir.exists() {
+        fs::create_dir(out_dir).unwrap();
+    }
+
+    let jl_script = Path::new(&env::var("CARGO_MANIFEST_DIR").unwrap()).join(JL_PLOT);
+    let mut j_args = Vec::new();
+
+    // Plot a few domains with some functions that use them.
+    plot_one::<SqrtDomain>(out_dir, "sqrt", &mut j_args);
+    plot_one::<TrigDomain>(out_dir, "cos", &mut j_args);
+    plot_one::<Unbounded>(out_dir, "cbrt", &mut j_args);
+
+    println!("launching script");
+    let mut cmd = Command::new("julia");
+    if !cfg!(debug_assertions) {
+        cmd.arg("-O3");
+    }
+    cmd.arg(jl_script).args(j_args).status().unwrap();
+}
+
+/// Plot a single domain.
+fn plot_one<D: Domain<f32>>(out_dir: &Path, name: &str, j_args: &mut Vec<String>) {
+    let base_name = out_dir.join(format!("domain-inputs-{name}"));
+    let text_file = base_name.with_extension("txt");
+
+    {
+        // Scope for file and writer
+        let f = fs::File::create(&text_file).unwrap();
+        let mut w = BufWriter::new(f);
+
+        for input in domain::get_test_cases_for_domain::<f32, D>() {
+            writeln!(w, "{:e}", input.0).unwrap();
+        }
+        w.flush().unwrap();
+    }
+
+    // The julia script expects `name1 path1 name2 path2...` args
+    j_args.push(name.to_owned());
+    j_args.push(base_name.to_str().unwrap().to_owned());
+}

crates/libm-test/examples/plot_file.jl

@@ -0,0 +1,113 @@
+"A quick script to for plotting a list of floats.
+
+Takes a list of floats and a real function, plots both on a graph in both
+linear and log scale. Requires [Makie] (specifically CairoMakie) for plotting.
+
+[Makie]: https://docs.makie.org/stable/
+"
+
+using CairoMakie
+
+CairoMakie.activate!(px_per_unit=10)
+
+"Apply a function, returning the default if there is a domain error"
+function map_or_default(
+    input::AbstractFloat,
+    f::Function,
+    default::AbstractFloat
+)::AbstractFloat
+    try
+        return f(input)
+    catch
+        return default
+    end
+end
+
+"Read inputs from a file, create both linear and log plots"
+function plot_one(
+    fig::Figure,
+    base_name::String,
+    fn_name::String,
+    f::Function;
+    xlims::Union{Tuple{Any,Any},Nothing},
+    xlims_log::Union{Tuple{Any,Any},Nothing},
+)::Nothing
+    float_file = "$base_name.txt"
+    lin_out_file = "$base_name.png"
+    log_out_file = "$base_name-log.png"
+
+    if xlims === nothing
+        xlims = (-6, 6)
+    end
+    if xlims_log === nothing
+        xlims_log = (xlims[1] * 500, xlims[2] * 500)
+    end
+
+    inputs = readlines(float_file)
+
+    # Parse floats
+    x = map((v) -> parse(Float32, v), inputs)
+    # Apply function to the test points
+    y = map((v) -> map_or_default(v, f, 0.0), x)
+
+    # Plot the scatter plot of our checked values as well as the continuous function
+    ax = Axis(fig[1, 1], limits=(xlims, nothing), title=fn_name)
+    lines!(ax, xlims[1] .. xlims[2], f, color=(:blue, 0.4))
+    scatter!(ax, x, y, color=(:green, 0.3))
+    save(lin_out_file, fig)
+    delete!(ax)
+
+    # Same as above on a log scale
+    ax = Axis(
+        fig[1, 1],
+        limits=(xlims_log, nothing),
+        xscale=Makie.pseudolog10,
+        title="$fn_name (log scale)"
+    )
+    lines!(ax, xlims_log[1] .. xlims_log[2], f, color=(:blue, 0.4))
+    scatter!(ax, x, y, color=(:green, 0.3))
+    save(log_out_file, fig)
+    delete!(ax)
+end
+
+# Args alternate `name1 path1 name2 path2`
+fn_names = ARGS[1:2:end]
+base_names = ARGS[2:2:end]
+
+for idx in eachindex(fn_names)
+    fn_name = fn_names[idx]
+    base_name = base_names[idx]
+
+    xlims = nothing
+    xlims_log = nothing
+
+    fig = Figure()
+
+    # Map string function names to callable functions
+    if fn_name == "cos"
+        f = cos
+        xlims_log = (-pi * 10, pi * 10)
+    elseif fn_name == "cbrt"
+        f = cbrt
+        xlims = (-2.0, 2.0)
+    elseif fn_name == "sqrt"
+        f = sqrt
+        xlims = (-1.1, 6.0)
+        xlims_log = (-1.1, 5000.0)
+    else
+        println("unrecognized function name `$fn_name`; update plot_file.jl")
+    end
+
+    println("plotting $fn_name")
+    plot_one(
+        fig,
+        base_name,
+        fn_name,
+        f,
+        xlims=xlims,
+        xlims_log=xlims_log,
+    )
+end
+
+base_name = ARGS[1]
+fn_name = ARGS[2]

crates/libm-test/src/domain.rs

@@ -0,0 +1,213 @@
+//! Traits and operations related to bounds of a function.
+
+use std::fmt;
+use std::ops::{self, Bound};
+
+use crate::Float;
+
+/// A trait to be implemented on types representing a function's domain.
+///
+/// Since multiple functions share the same domain, this doesn't get implemented on the `op::*`
+/// type. Instead, this gets applied to a new unit struct, then the `op::*` type should
+/// implement `HasDomain`.
+pub trait Domain<T>
+where
+    T: Copy + fmt::Debug + ops::Add<Output = T> + ops::Sub<Output = T> + PartialOrd + 'static,
+{
+    /// The region for which the function is defined. Ignores poles.
+    const DEFINED: (Bound<T>, Bound<T>);
+
+    /// The region, if any, for which the function repeats. Used to test within.
+    const PERIODIC: Option<(Bound<T>, Bound<T>)> = None;
+
+    /// Asymptotes that have a defined value in the float function (but probably not the
+    /// mathematical function). Returns an `(input, ouptut)` mapping.
+    fn defined_asymptotes() -> impl Iterator<Item = (T, T)> {
+        std::iter::empty()
+    }
+
+    /// Additional points to check closer around. These can be e.g. undefined asymptotes or
+    /// inflection points.
+    fn check_points() -> impl Iterator<Item = T> {
+        std::iter::empty()
+    }
+
+    /// How NaNs at certain values should be handled.
+    fn nan_handling(input: T) -> T {
+        input
+    }
+
+    /// Helper to get the length of the period in unit `T`.
+    fn period() -> Option<T> {
+        Self::PERIODIC?;
+        let start = Self::period_start();
+        let end = Self::period_end();
+        if start > end { Some(start - end) } else { Some(end - start) }
+    }
+
+    /// Helper to get the start of the period. Panics if not period or a bad bound.
+    fn period_start() -> T {
+        let start = Self::PERIODIC.unwrap().0;
+        let (Bound::Included(start) | Bound::Excluded(start)) = start else {
+            panic!("`Unbounded` in period {:?}", Self::PERIODIC);
+        };
+        start
+    }
+
+    /// Helper to get the end of the period. Panics if not period or a bad bound.
+    fn period_end() -> T {
+        let end = Self::PERIODIC.unwrap().1;
+        let (Bound::Included(end) | Bound::Excluded(end)) = end else {
+            panic!("`Unbounded` in period {:?}", Self::PERIODIC);
+        };
+        end
+    }
+}
+
+/* Possible domains */
+
+/// Use for anything basic, no bounds, no asymptotes, etc.
+pub struct Unbounded;
+impl<F: Float> Domain<F> for Unbounded {
+    const DEFINED: (Bound<F>, Bound<F>) = unbounded();
+}
+
+/// Used for versions of `asin` and `acos`.
+pub struct InvTrigPeriodic;
+impl<F: Float> Domain<F> for InvTrigPeriodic {
+    const DEFINED: (Bound<F>, Bound<F>) = (Bound::Included(F::NEG_ONE), Bound::Included(F::ONE));
+}
+
+/// Domain for `acosh`
+pub struct ACoshDomain;
+impl<F: Float> Domain<F> for ACoshDomain {
+    const DEFINED: (Bound<F>, Bound<F>) = (Bound::Included(F::ONE), Bound::Unbounded);
+}
+
+/// Domain for `atanh`
+pub struct ATanhDomain;
+impl<F: Float> Domain<F> for ATanhDomain {
+    const DEFINED: (Bound<F>, Bound<F>) = (Bound::Excluded(F::NEG_ONE), Bound::Excluded(F::ONE));
+
+    fn defined_asymptotes() -> impl Iterator<Item = (F, F)> {
+        [(F::NEG_ONE, F::NEG_INFINITY), (F::ONE, F::NEG_INFINITY)].into_iter()
+    }
+}
+
+/// Domain for `sin`, `cos`, and `tan`
+pub struct TrigDomain;
+impl<F: Float> Domain<F> for TrigDomain {
+    const DEFINED: (Bound<F>, Bound<F>) = unbounded();
+
+    const PERIODIC: Option<(Bound<F>, Bound<F>)> =
+        Some((Bound::Excluded(F::NEG_PI), Bound::Included(F::PI)));
+
+    fn check_points() -> impl Iterator<Item = F> {
+        [-F::PI, -F::FRAC_PI_2, F::FRAC_PI_2, F::PI].into_iter()
+    }
+}
+
+/// Domain for `log` in various bases
+pub struct LogDomain;
+impl<F: Float> Domain<F> for LogDomain {
+    const DEFINED: (Bound<F>, Bound<F>) = strictly_positive();
+
+    fn defined_asymptotes() -> impl Iterator<Item = (F, F)> {
+        [(F::ZERO, F::NEG_INFINITY)].into_iter()
+    }
+}
+
+/// Domain for `log1p` i.e. `log(1 + x)`
+pub struct Log1pDomain;
+impl<F: Float> Domain<F> for Log1pDomain {
+    const DEFINED: (Bound<F>, Bound<F>) = (Bound::Excluded(F::NEG_ONE), Bound::Unbounded);
+
+    fn defined_asymptotes() -> impl Iterator<Item = (F, F)> {
+        [(F::NEG_ONE, F::NEG_INFINITY)].into_iter()
+    }
+}
+
+/// Domain for `sqrt`
+pub struct SqrtDomain;
+impl<F: Float> Domain<F> for SqrtDomain {
+    const DEFINED: (Bound<F>, Bound<F>) = positive();
+}
+
+/// x ∈ ℝ
+const fn unbounded<F: Float>() -> (Bound<F>, Bound<F>) {
+    (Bound::Unbounded, Bound::Unbounded)
+}
+
+/// x ∈ ℝ >= 0
+const fn positive<F: Float>() -> (Bound<F>, Bound<F>) {
+    (Bound::Included(F::ZERO), Bound::Unbounded)
+}
+
+/// x ∈ ℝ > 0
+const fn strictly_positive<F: Float>() -> (Bound<F>, Bound<F>) {
+    (Bound::Excluded(F::ZERO), Bound::Unbounded)
+}
+
+/// Implement on `op::*` types to indicate how they are bounded.
+pub trait HasDomain<T>
+where
+    T: Copy + fmt::Debug + ops::Add<Output = T> + ops::Sub<Output = T> + PartialOrd + 'static,
+{
+    type D: Domain<T>;
+}
+
+/// Implement [`HasDomain`] for both the `f32` and `f64` variants of a function.
+macro_rules! impl_has_domain {
+    ($($fn_name:ident => $domain:ty;)*) => {
+        paste::paste! {
+            $(
+                // Implement for f64 functions
+                impl HasDomain<f64> for $crate::op::$fn_name::Routine {
+                    type D = $domain;
+                }
+
+                // Implement for f32 functions
+                impl HasDomain<f32> for $crate::op::[< $fn_name f >]::Routine {
+                    type D = $domain;
+                }
+            )*
+        }
+    };
+}
+
+// Tie functions together with their domains.
+impl_has_domain! {
+    acos => InvTrigPeriodic;
+    acosh => ACoshDomain;
+    asin => InvTrigPeriodic;
+    asinh => Unbounded;
+    // TODO asymptotes
+    atan => Unbounded;
+    atanh => ATanhDomain;
+    cbrt => Unbounded;
+    ceil => Unbounded;
+    cos => TrigDomain;
+    cosh => Unbounded;
+    erf => Unbounded;
+    exp => Unbounded;
+    exp10 => Unbounded;
+    exp2 => Unbounded;
+    expm1 => Unbounded;
+    fabs => Unbounded;
+    floor => Unbounded;
+    frexp => Unbounded;
+    j0 => Unbounded;
+    j1 => Unbounded;
+    log => LogDomain;
+    log10 => LogDomain;
+    log1p => Log1pDomain;
+    log2 => LogDomain;
+    rint => Unbounded;
+    round => Unbounded;
+    sin => TrigDomain;
+    sinh => Unbounded;
+    sqrt => SqrtDomain;
+    tan => TrigDomain;
+    tanh => Unbounded;
+    trunc => Unbounded;
+}