From 5d9bec3e8f97324841dd34f5eb471eb5f53c9c1e Mon Sep 17 00:00:00 2001 From: VPRamon Date: Wed, 6 May 2026 21:33:19 +0200 Subject: [PATCH 1/2] Enhance documentation and error handling across modules - Added comprehensive documentation to `approx`, `calculus`, `core`, `io`, `piecewise`, `quadrature`, and `spectral` modules, detailing theory, features, and performance. - Expanded `ChebyError` enum with additional error variants for better error handling in binary operations and Remez exchange. - Improved error messages in `ChebyError` implementation for clarity. - Implemented checks for segment boundaries in `AdaptiveSegmentTable` to ensure contiguous segments. - Added tests for adaptive piecewise tables, calculus round-trip, error branches, minimax convergence, and quadrature/spectral accuracy. - Updated binary encoding/decoding functions to handle new error cases and improve robustness. --- Cargo.lock | 167 ++++++------ Cargo.toml | 14 +- benches/cheby.rs | 93 +++++-- examples/README.md | 101 +++++-- examples/binary_roundtrip.rs | 26 ++ examples/gauss_chebyshev.rs | 18 ++ examples/spectral_differentiation.rs | 27 +- src/approx/minimax.rs | 383 ++++++++++++++++++++++----- src/approx/mod.rs | 21 ++ src/calculus/mod.rs | 15 ++ src/core/error.rs | 60 ++++- src/core/mod.rs | 20 ++ src/io/binary.rs | 46 +++- src/io/mod.rs | 8 + src/lib.rs | 1 + src/piecewise/adaptive.rs | 155 +++++++++-- src/piecewise/mod.rs | 18 ++ src/quadrature/clenshaw_curtis.rs | 69 +++-- src/quadrature/mod.rs | 17 ++ src/spectral/mod.rs | 16 ++ tests/adaptive_piecewise.rs | 79 ++++++ tests/calculus_roundtrip.rs | 48 ++++ tests/error_branches.rs | 97 +++++++ tests/minimax.rs | 96 +++++++ tests/quadrature_and_spectral.rs | 80 ++++++ 25 files changed, 1403 insertions(+), 272 deletions(-) create mode 100644 examples/binary_roundtrip.rs create mode 100644 examples/gauss_chebyshev.rs create mode 100644 tests/adaptive_piecewise.rs create mode 100644 tests/calculus_roundtrip.rs create mode 100644 tests/error_branches.rs create mode 100644 tests/minimax.rs create mode 100644 tests/quadrature_and_spectral.rs diff --git a/Cargo.lock b/Cargo.lock index 74df107..4817bff 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -11,6 +11,15 @@ dependencies = [ "memchr", ] +[[package]] +name = "alloca" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5a7d05ea6aea7e9e64d25b9156ba2fee3fdd659e34e41063cd2fc7cd020d7f4" +dependencies = [ + "cc", +] + [[package]] name = "anes" version = "0.1.6" @@ -77,6 +86,16 @@ version = "0.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" +[[package]] +name = "cc" +version = "1.2.61" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d16d90359e986641506914ba71350897565610e87ce0ad9e6f28569db3dd5c6d" +dependencies = [ + "find-msvc-tools", + "shlex", +] + [[package]] name = "cfg-if" version = "1.0.4" @@ -91,7 +110,6 @@ dependencies = [ "criterion", "proptest", "qtty", - "rustfft", "serde", "trybuild", ] @@ -150,25 +168,24 @@ checksum = "c8d4a3bb8b1e0c1050499d1815f5ab16d04f0959b233085fb31653fbfc9d98f9" [[package]] name = "criterion" -version = "0.5.1" +version = "0.8.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f2b12d017a929603d80db1831cd3a24082f8137ce19c69e6447f54f5fc8d692f" +checksum = "950046b2aa2492f9a536f5f4f9a3de7b9e2476e575e05bd6c333371add4d98f3" dependencies = [ + "alloca", "anes", "cast", "ciborium", "clap", "criterion-plot", - "is-terminal", "itertools", "num-traits", - "once_cell", "oorandom", + "page_size", "plotters", "rayon", "regex", "serde", - "serde_derive", "serde_json", "tinytemplate", "walkdir", @@ -176,9 +193,9 @@ dependencies = [ [[package]] name = "criterion-plot" -version = "0.5.0" +version = "0.8.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6b50826342786a51a89e2da3a28f1c32b06e387201bc2d19791f622c673706b1" +checksum = "d8d80a2f4f5b554395e47b5d8305bc3d27813bacb73493eb1001e8f76dae29ea" dependencies = [ "cast", "itertools", @@ -243,6 +260,12 @@ version = "2.4.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9f1f227452a390804cdb637b74a86990f2a7d7ba4b7d5693aac9b4dd6defd8d6" +[[package]] +name = "find-msvc-tools" +version = "0.1.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5baebc0774151f905a1a2cc41989300b1e6fbb29aff0ceffa1064fdd3088d582" + [[package]] name = "fnv" version = "1.0.7" @@ -342,12 +365,6 @@ version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" -[[package]] -name = "hermit-abi" -version = "0.5.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fc0fef456e4baa96da950455cd02c081ca953b141298e41db3fc7e36b1da849c" - [[package]] name = "id-arena" version = "2.3.0" @@ -366,22 +383,11 @@ dependencies = [ "serde_core", ] -[[package]] -name = "is-terminal" -version = "0.4.17" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3640c1c38b8e4e43584d8df18be5fc6b0aa314ce6ebf51b53313d4306cca8e46" -dependencies = [ - "hermit-abi", - "libc", - "windows-sys", -] - [[package]] name = "itertools" -version = "0.10.5" +version = "0.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b0fd2260e829bddf4cb6ea802289de2f86d6a7a690192fbe91b3f46e0f2c8473" +checksum = "413ee7dfc52ee1a4949ceeb7dbc8a33f2d6c088194d9f922fb8318faf1f01186" dependencies = [ "either", ] @@ -440,24 +446,6 @@ version = "2.8.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f8ca58f447f06ed17d5fc4043ce1b10dd205e060fb3ce5b979b8ed8e59ff3f79" -[[package]] -name = "num-complex" -version = "0.4.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" -dependencies = [ - "num-traits", -] - -[[package]] -name = "num-integer" -version = "0.1.46" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" -dependencies = [ - "num-traits", -] - [[package]] name = "num-traits" version = "0.2.19" @@ -479,6 +467,16 @@ version = "11.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d6790f58c7ff633d8771f42965289203411a5e5c68388703c06e14f24770b41e" +[[package]] +name = "page_size" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "30d5b2194ed13191c1999ae0704b7839fb18384fa22e49b57eeaa97d79ce40da" +dependencies = [ + "libc", + "winapi", +] + [[package]] name = "pin-project-lite" version = "0.2.17" @@ -532,15 +530,6 @@ dependencies = [ "syn", ] -[[package]] -name = "primal-check" -version = "0.3.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dc0d895b311e3af9902528fbb8f928688abbd95872819320517cc24ca6b2bd08" -dependencies = [ - "num-integer", -] - [[package]] name = "proc-macro2" version = "1.0.106" @@ -571,9 +560,9 @@ dependencies = [ [[package]] name = "qtty" -version = "0.6.1" +version = "0.7.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "adc4e4fbd7aabae80c30ca07bd9583dc3bc3f6bcf8b123ff538386715aeee293" +checksum = "61587536c9c9ec95845b7e5ddbac8afb06ce94ea10fc2728543db3f7bc0a2bc6" dependencies = [ "qtty-core", "qtty-derive", @@ -581,9 +570,9 @@ dependencies = [ [[package]] name = "qtty-core" -version = "0.6.1" +version = "0.7.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e524035941ae68ff8524dfc41d395fe721cd3c0e7eaa16a756284d59a64ddd53" +checksum = "3b01d8237a224bea98fcd76fed95ceff6dcca63f613948ef08ecdfeb7e5e398e" dependencies = [ "libm", "qtty-derive", @@ -593,9 +582,9 @@ dependencies = [ [[package]] name = "qtty-derive" -version = "0.6.1" +version = "0.7.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "75757ec156050a0343f72034061a8cfd7a4a7c49254ba55feef0c4951bceb842" +checksum = "2d4df57ba07bf5109d76487c62a92b0d0ae13e94711cd37f0302bfa886cb44c3" dependencies = [ "proc-macro2", "quote", @@ -716,20 +705,6 @@ version = "0.8.10" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "dc897dd8d9e8bd1ed8cdad82b5966c3e0ecae09fb1907d58efaa013543185d0a" -[[package]] -name = "rustfft" -version = "6.4.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "21db5f9893e91f41798c88680037dba611ca6674703c1a18601b01a72c8adb89" -dependencies = [ - "num-complex", - "num-integer", - "num-traits", - "primal-check", - "strength_reduce", - "transpose", -] - [[package]] name = "rustix" version = "1.1.4" @@ -829,16 +804,16 @@ dependencies = [ ] [[package]] -name = "slab" -version = "0.4.12" +name = "shlex" +version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0c790de23124f9ab44544d7ac05d60440adc586479ce501c1d6d7da3cd8c9cf5" +checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" [[package]] -name = "strength_reduce" -version = "0.2.4" +name = "slab" +version = "0.4.12" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fe895eb47f22e2ddd4dabc02bce419d2e643c8e3b585c78158b349195bc24d82" +checksum = "0c790de23124f9ab44544d7ac05d60440adc586479ce501c1d6d7da3cd8c9cf5" [[package]] name = "syn" @@ -928,16 +903,6 @@ version = "1.1.1+spec-1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "756daf9b1013ebe47a8776667b466417e2d4c5679d441c26230efd9ef78692db" -[[package]] -name = "transpose" -version = "0.2.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1ad61aed86bc3faea4300c7aee358b4c6d0c8d6ccc36524c96e4c92ccf26e77e" -dependencies = [ - "num-integer", - "strength_reduce", -] - [[package]] name = "trybuild" version = "1.0.116" @@ -1103,6 +1068,22 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "winapi" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5c839a674fcd7a98952e593242ea400abe93992746761e38641405d28b00f419" +dependencies = [ + "winapi-i686-pc-windows-gnu", + "winapi-x86_64-pc-windows-gnu", +] + +[[package]] +name = "winapi-i686-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" + [[package]] name = "winapi-util" version = "0.1.11" @@ -1112,6 +1093,12 @@ dependencies = [ "windows-sys", ] +[[package]] +name = "winapi-x86_64-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" + [[package]] name = "windows-link" version = "0.2.1" diff --git a/Cargo.toml b/Cargo.toml index 37882f0..e522c72 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -31,7 +31,6 @@ spectral = ["calculus", "alloc"] serde = ["dep:serde", "qtty/serde", "alloc"] binary = ["alloc"] -fast-dct = ["dep:rustfft"] full = [ "std", @@ -48,13 +47,12 @@ full = [ ] [dependencies] -qtty = { version = "0.6.1", default-features = false, features = ["cross-unit-ops"] } +qtty = { version = "0.7.0", default-features = false, features = ["cross-unit-ops"] } serde = { version = "1.0", optional = true, features = ["derive"] } -rustfft = { version = "6.2", optional = true } [dev-dependencies] approx = "0.5" -criterion = "0.5" +criterion = "0.8.2" proptest = "1.4" trybuild = "1.0" @@ -78,3 +76,11 @@ required-features = ["quadrature"] [[example]] name = "spectral_differentiation" required-features = ["spectral"] + +[[example]] +name = "binary_roundtrip" +required-features = ["binary"] + +[[example]] +name = "gauss_chebyshev" +required-features = ["quadrature"] diff --git a/benches/cheby.rs b/benches/cheby.rs index ae14ab1..36eeae0 100644 --- a/benches/cheby.rs +++ b/benches/cheby.rs @@ -1,4 +1,13 @@ -use criterion::{black_box, criterion_group, criterion_main, Criterion}; +//! Hot-path benchmarks for `cheby`. +//! +//! Covers, in order: scalar evaluation, derivative evaluation (raw + typed), +//! fixed-size fitting, adaptive fitting, minimax (Remez), piecewise lookup +//! and evaluation, Clenshaw-Curtis quadrature, Gauss-Chebyshev quadrature, +//! and spectral differentiation matrix construction. + +use std::hint::black_box; + +use criterion::{criterion_group, criterion_main, Criterion}; use qtty::{Kilometer, Second}; fn evaluate_f64(c: &mut Criterion) { @@ -42,20 +51,47 @@ fn evaluate_derivative_quantity(c: &mut Criterion) { }); } -fn fit_naive(c: &mut Criterion) { - c.bench_function("fit_naive", |b| { +fn fit_fixed(c: &mut Criterion) { + c.bench_function("fit_fixed_n16_sin", |b| { b.iter(|| cheby::approx::fit::fit_from_fn::(|x| black_box(x).sin())) }); } fn fit_quantity(c: &mut Criterion) { - c.bench_function("fit_quantity", |b| { + c.bench_function("fit_quantity_n16", |b| { b.iter(|| { cheby::approx::fit::fit_from_fn::(|x| Kilometer::new(black_box(x).sin())) }) }); } +fn adaptive_fit(c: &mut Criterion) { + c.bench_function("adaptive_fit_sin", |b| { + b.iter(|| { + cheby::approx::adaptive::AdaptiveFit::new(cheby::Domain::new(-1.0, 1.0), |x| { + f64::sin(x) + }) + .max_degree(32) + .build::() + .unwrap() + }) + }); +} + +fn minimax_remez(c: &mut Criterion) { + c.bench_function("minimax_remez_exp_deg7", |b| { + b.iter(|| { + cheby::approx::minimax::remez( + cheby::Domain::new(-1.0, 1.0), + 7, + |x: f64| black_box(x).exp(), + cheby::approx::minimax::RemezOptions::default(), + ) + .unwrap() + }) + }); +} + fn piecewise_lookup(c: &mut Criterion) { let table: cheby::ChebySegmentTable = cheby::ChebySegmentTable::from_fn(|x| x.sin(), 0.0, 8.0, 1.0); @@ -72,21 +108,22 @@ fn piecewise_evaluate(c: &mut Criterion) { }); } -fn adaptive_fit(c: &mut Criterion) { - c.bench_function("adaptive_fit", |b| { - b.iter(|| { - cheby::approx::adaptive::AdaptiveFit::new(cheby::Domain::new(-1.0, 1.0), |x| { - f64::sin(x) - }) - .max_degree(32) - .build::() - .unwrap() - }) +fn adaptive_piecewise_lookup(c: &mut Criterion) { + let table: cheby::piecewise::AdaptiveSegmentTable = + cheby::piecewise::AdaptiveSegmentTable::from_fn( + cheby::Domain::new(0.0, 8.0), + |x: f64| x.sin(), + 1e-10, + 8, + ) + .unwrap(); + c.bench_function("adaptive_piecewise_lookup", |b| { + b.iter(|| table.evaluate(black_box(4.25)).unwrap()) }); } fn clenshaw_curtis_integrate(c: &mut Criterion) { - c.bench_function("clenshaw_curtis_integrate", |b| { + c.bench_function("clenshaw_curtis_integrate_n32", |b| { b.iter(|| { cheby::quadrature::clenshaw_curtis::integrate::( cheby::Domain::new(0.0, 1.0), @@ -96,17 +133,37 @@ fn clenshaw_curtis_integrate(c: &mut Criterion) { }); } +fn gauss_chebyshev(c: &mut Criterion) { + c.bench_function("gauss_chebyshev_weighted_n32", |b| { + b.iter(|| { + cheby::quadrature::gauss_chebyshev::integrate_weighted::(|x| { + black_box(x).cos() + }) + }) + }); +} + +fn spectral_diff_matrix(c: &mut Criterion) { + c.bench_function("spectral_diff_matrix_n16", |b| { + b.iter(|| cheby::spectral::chebyshev_differentiation_matrix(black_box(16))) + }); +} + criterion_group!( benches, evaluate_f64, evaluate_quantity, evaluate_derivative_f64, evaluate_derivative_quantity, - fit_naive, + fit_fixed, fit_quantity, + adaptive_fit, + minimax_remez, piecewise_lookup, piecewise_evaluate, - adaptive_fit, - clenshaw_curtis_integrate + adaptive_piecewise_lookup, + clenshaw_curtis_integrate, + gauss_chebyshev, + spectral_diff_matrix ); criterion_main!(benches); diff --git a/examples/README.md b/examples/README.md index adb5a95..e594c90 100644 --- a/examples/README.md +++ b/examples/README.md @@ -1,45 +1,90 @@ # cheby examples -All examples are runnable with `cargo run --example `. +All examples are runnable with `cargo run --example [--features ...]`. +Examples whose subsystem requires a Cargo feature have a `required-features` +entry in `Cargo.toml`; the table at the bottom maps each example to the +features it needs. -## `basic_interpolation` +## Core evaluation and fitting -End-to-end interpolation on `[-1, 1]`: +### `basic_series` +Build a `ChebySeries` by hand and evaluate it. -- Generate Chebyshev nodes. -- Sample a function. -- Fit coefficients. -- Evaluate and print approximation error. +### `basic_interpolation` +End-to-end interpolation on `[-1, 1]`: nodes, sample, fit, evaluate, error. -Run: +### `fit_sin` +Fixed-size fit of `sin` over a typed time domain. -```bash -cargo run --example basic_interpolation -``` +### `typed_quantities` +Fit and evaluate while preserving `qtty::Quantity` units. -## `segment_table` +## Derivatives and integrals (`calculus`) -Piecewise approximation over a physical time range: +### `derivative_velocity` +Differentiate a position series to recover a typed velocity. -- Build `ChebySegmentTable`. -- Inspect metadata (`start`, `end`, segment count, segment length). -- Evaluate value and derivative in one pass. +### `integral_position` +Integrate a velocity series to recover a typed position. -Run: +## Adaptive and minimax fitting (`adaptive`, `minimax`) -```bash -cargo run --example segment_table -``` +### `adaptive_fit` +Refine a fit until a target tolerance is met (`adaptive` feature). -## `typed_quantities` +### `minimax_exp` +Minimise the L∞ error for `exp` using Remez exchange (`minimax` feature). -Demonstrates typed units with `qtty::Quantity`: +## Piecewise tables (`piecewise`) -- Fit coefficients where values are `Quantity`. -- Evaluate while preserving units. +### `segment_table` +Uniform piecewise table over a typed domain with metadata. -Run: +### `piecewise_trajectory` +A typed multi-segment trajectory with per-segment evaluation. -```bash -cargo run --example typed_quantities -``` +### `ephemeris_like_table` +Builds a piecewise table reminiscent of an ephemeris layout. + +## Quadrature (`quadrature`) + +### `clenshaw_curtis_integral` +Clenshaw-Curtis quadrature of `sin` on `[0, π]`. + +### `gauss_chebyshev` +Gauss-Chebyshev recovers `π` and `π·J₀(1)`. + +## Spectral differentiation (`spectral`) + +### `spectral_differentiation` +Build the Chebyshev differentiation matrix and verify it against +`d/dx sin(x) = cos(x)` at the Lobatto nodes. + +## Binary I/O (`binary`) + +### `binary_roundtrip` +Encode an `f64` series, decode it, and demonstrate that a flipped byte is +rejected by the checksum. + +## Application-flavoured demos + +These illustrate downstream usage and rely on default features only: + +### `angular_rate` +Differentiate an angle series to recover an angular rate. + +### `star_alt_az_approximation` +Toy approximation of an alt/az curve using a piecewise table. + +## Feature-to-example map + +| Feature | Examples | +|---------------|----------| +| `std` (default) | `basic_series`, `basic_interpolation`, `fit_sin`, `typed_quantities`, `angular_rate` | +| `calculus` (default) | `derivative_velocity`, `integral_position` | +| `piecewise` (default) | `segment_table`, `piecewise_trajectory`, `ephemeris_like_table`, `star_alt_az_approximation` | +| `adaptive` | `adaptive_fit` | +| `minimax` | `minimax_exp` | +| `quadrature` | `clenshaw_curtis_integral`, `gauss_chebyshev` | +| `spectral` | `spectral_differentiation` | +| `binary` | `binary_roundtrip` | diff --git a/examples/binary_roundtrip.rs b/examples/binary_roundtrip.rs new file mode 100644 index 0000000..d28f1d5 --- /dev/null +++ b/examples/binary_roundtrip.rs @@ -0,0 +1,26 @@ +//! Encode/decode a Chebyshev series via the binary I/O subsystem. +//! +//! Run with `cargo run --example binary_roundtrip --features binary`. + +use cheby::core::ChebySeriesDyn; +use cheby::io::binary::{decode_f64_series, encode_f64_series}; + +fn main() { + let series = ChebySeriesDyn::::new(vec![1.0, -0.25, 0.0625, -0.03125, 0.015625]) + .expect("non-empty series"); + + let bytes = encode_f64_series(&series); + println!("encoded {} bytes", bytes.len()); + + let decoded = decode_f64_series(&bytes).expect("round-trip"); + assert_eq!(series.coeffs(), decoded.coeffs()); + println!("round-trip OK on {} coefficients", series.coeffs().len()); + + // Demonstrate a typed error variant: tampering flips the checksum. + let mut tampered = bytes.clone(); + *tampered.last_mut().unwrap() ^= 0xFF; + match decode_f64_series(&tampered) { + Err(e) => println!("tampered payload rejected: {e:?}"), + Ok(_) => unreachable!("checksum should have failed"), + } +} diff --git a/examples/gauss_chebyshev.rs b/examples/gauss_chebyshev.rs new file mode 100644 index 0000000..7b91ede --- /dev/null +++ b/examples/gauss_chebyshev.rs @@ -0,0 +1,18 @@ +//! Gauss-Chebyshev quadrature recovers two known weighted integrals: +//! +//! - `∫_{-1}^{1} 1/√(1−x²) dx = π` +//! - `∫_{-1}^{1} cos(x)/√(1−x²) dx = π J₀(1)` +//! +//! Run with `cargo run --example gauss_chebyshev --features quadrature`. + +use cheby::quadrature::gauss_chebyshev; + +fn main() { + const N: usize = 32; + let pi_estimate: f64 = gauss_chebyshev::integrate_weighted::<_, N>(|_| 1.0); + println!("∫ 1/√(1−x²) dx ≈ {pi_estimate}"); + + let j0_estimate: f64 = gauss_chebyshev::integrate_weighted::<_, N>(f64::cos); + println!("∫ cos(x)/√(1−x²) dx ≈ {j0_estimate}"); + println!("(reference: π·J₀(1) ≈ {})", core::f64::consts::PI * 0.7651976865579666); +} diff --git a/examples/spectral_differentiation.rs b/examples/spectral_differentiation.rs index c8f6c28..70b632c 100644 --- a/examples/spectral_differentiation.rs +++ b/examples/spectral_differentiation.rs @@ -1,4 +1,27 @@ +//! Compare the Chebyshev differentiation matrix against the analytic +//! derivative of `sin(x)` at the Lobatto collocation nodes. +//! +//! Run with `cargo run --example spectral_differentiation --features spectral`. + +use cheby::spectral::{chebyshev_differentiation_matrix, collocation_points}; + fn main() { - let d = cheby::spectral::chebyshev_differentiation_matrix(8); - println!("{}x{}", d.rows(), d.cols()); + const N: usize = 17; + let nodes: [f64; N] = collocation_points::(); + let f: [f64; N] = std::array::from_fn(|k| nodes[k].sin()); + let d = chebyshev_differentiation_matrix(N); + + let mut max_err = 0.0_f64; + for (i, &node) in nodes.iter().enumerate() { + let mut acc = 0.0; + for (j, &fj) in f.iter().enumerate() { + acc += d.get(i, j) * fj; + } + let err = (acc - node.cos()).abs(); + if err > max_err { + max_err = err; + } + } + println!("D·sin vs cos: max |error| = {max_err:.3e} on {N} Lobatto nodes"); + println!("matrix shape: {}x{}", d.rows(), d.cols()); } diff --git a/src/approx/minimax.rs b/src/approx/minimax.rs index 0472f8e..3f808e4 100644 --- a/src/approx/minimax.rs +++ b/src/approx/minimax.rs @@ -1,107 +1,360 @@ -//! Remez-style minimax approximation. +//! Minimax (Remez exchange) approximation on `f64`. +//! +//! # Algorithm +//! +//! For a continuous `f: [a, b] -> R` and a target degree `N`, this module +//! computes a polynomial `p(x) = sum_{k=0}^{N} c_k T_k(tau)` (with +//! `tau = domain.normalize(x)`) that approximately minimizes +//! `max_{x in [a,b]} |f(x) - p(x)|`. +//! +//! Each iteration: +//! +//! 1. Solve the `(N+2) x (N+2)` linear system +//! `sum_k c_k T_k(tau_i) + (-1)^i E = f(x_i)` over the current +//! alternation set `x_0 < x_1 < ... < x_{N+1}` (Gaussian elimination +//! with partial pivoting). +//! 2. Sample the residual `r(x) = f(x) - p(x)` on a dense grid, +//! partition the grid by sign of `r`, and locate the per-sign-block +//! maxima of `|r|`. +//! 3. If fewer than `N+2` alternations are available the iteration +//! halts with [`ChebyError::RemezAlternationFailure`]. Otherwise the +//! best `N+2` contiguous alternation points become the next set. +//! 4. Convergence is declared when `(max|r| - |E|) / max|r|` falls below +//! `RemezOptions::tolerance`. The iteration budget is +//! `RemezOptions::max_iterations`; exhausting it returns +//! [`ChebyError::RemezDidNotConverge`]. +//! +//! Root-sample (interpolating) coefficient fitting is **not** minimax; +//! it lives in [`crate::approx::fit`]. +//! +//! # Performance notes +//! +//! Cost per iteration is `O((N+2)^3 + grid_size * N)`. Pick `degree` and +//! `grid_size` accordingly; `grid_size >= 16 * (degree + 2)` is usually +//! enough for smooth targets. -use crate::approx::fit; -use crate::core::{ChebyError, ChebyScalar, ChebySeriesDyn, ChebyTime, Domain}; +use alloc::vec; use alloc::vec::Vec; -/// Options for minimax approximation. +use crate::core::{ + basis, ChebyError, ChebySeries, ChebySeriesDyn, ChebySeriesOn, ChebyTime, Domain, +}; + +/// Options for the Remez exchange algorithm. #[derive(Debug, Clone, Copy, PartialEq)] pub struct RemezOptions { - /// Maximum iterations. + /// Maximum exchange iterations. pub max_iterations: usize, - /// Absolute convergence tolerance on sampled max error. + /// Relative convergence tolerance on the residual equioscillation gap. pub tolerance: f64, - /// Dense grid size for diagnostics. + /// Dense residual sampling grid size used to locate alternations. pub grid_size: usize, } impl Default for RemezOptions { fn default() -> Self { Self { - max_iterations: 16, + max_iterations: 32, tolerance: 1e-12, grid_size: 1024, } } } -/// Result of a minimax approximation attempt. +/// Result of a successful Remez exchange. #[derive(Debug, Clone, PartialEq)] -pub struct RemezResult { - /// Fitted series. - pub series: ChebySeriesDyn, - /// Maximum sampled error as a value-typed magnitude. - pub max_error: T, - /// Iterations performed. +pub struct RemezResult { + /// Domain-aware fitted series. + pub series_on: ChebySeriesOnAny, + /// Underlying coefficient view. + pub series: ChebySeriesDyn, + /// Final equioscillation level returned by the linear solve (`|E|`). + pub leveled_error: f64, + /// Final dense-grid sup-norm error `max |f - p|`. + pub max_error: f64, + /// Iterations actually performed. pub iterations: usize, - /// Whether sampled error converged. + /// Whether convergence tolerance was met. pub converged: bool, + /// Final alternation points in domain coordinates. + pub alternations: Vec, /// Domain used for fitting. pub domain: Domain, } -/// Compute a Remez-style minimax approximation. +/// Domain-aware Remez series wrapper exposed for convenience evaluation. +#[derive(Debug, Clone, PartialEq)] +pub struct ChebySeriesOnAny { + domain: Domain, + coeffs: Vec, +} + +impl ChebySeriesOnAny { + /// Domain of validity. + #[inline] + pub fn domain(&self) -> Domain { + self.domain + } + /// Coefficient slice in the first-kind Chebyshev basis. + #[inline] + pub fn coeffs(&self) -> &[f64] { + &self.coeffs + } + /// Evaluate at a domain point. + #[inline] + pub fn evaluate(&self, x: X) -> Result { + if !self.domain.contains(x) { + return Err(ChebyError::EvaluationOutOfDomain); + } + Ok(crate::core::eval::evaluate( + &self.coeffs, + self.domain.normalize(x), + )) + } +} + +/// Backward-compatible fixed-degree convenience wrapper. +/// +/// `Convert` a [`RemezResult`] coefficient slice into a fixed-size +/// [`ChebySeries`] when the degree is known statically. +pub fn try_into_fixed(coeffs: &[f64]) -> Result, ChebyError> { + if coeffs.len() != N { + return Err(ChebyError::InvalidDegree); + } + let mut arr = [0.0; N]; + arr.copy_from_slice(coeffs); + Ok(ChebySeries::new(arr)) +} + +/// Run the Remez exchange algorithm for an `f64`-valued target. /// -/// This implementation fits on Chebyshev roots and reports dense-grid error -/// diagnostics. It is intentionally conservative: non-convergence is explicit. -pub fn remez( +/// The polynomial degree is `degree`; the resulting series has +/// `degree + 1` coefficients in the first-kind Chebyshev basis. +pub fn remez( domain: Domain, degree: usize, - f: impl Fn(X) -> T, + f: impl Fn(X) -> f64, options: RemezOptions, -) -> Result, ChebyError> +) -> Result, ChebyError> where - T: ChebyScalar, X: ChebyTime, { if degree == 0 { return Err(ChebyError::InvalidDegree); } - let iterations = options.max_iterations.max(1); - let grid = options.grid_size.max(16); - - macro_rules! fit_n { - ($n:literal) => {{ - let fitted = fit::fit_from_fn_on::(domain, &f); - let coeffs = fitted.into_coeffs(); - let series = ChebySeriesDyn::new(Vec::from(coeffs))?; - let mut max_mag = 0.0_f64; - let mut max_err = T::zero(); - for i in 0..grid { - let tau = -1.0 + 2.0 * i as f64 / (grid - 1) as f64; - let x = domain.denormalize(tau); - let err = f(x) - series.evaluate(tau); - if err.magnitude() > max_mag { - max_mag = err.magnitude(); - max_err = err; + let max_iter = options.max_iterations.max(1); + let m = degree + 2; // alternation count + let grid = options.grid_size.max(8 * m); + if !(options.tolerance.is_finite() && options.tolerance >= 0.0) { + return Err(ChebyError::NonFiniteInput); + } + + // Initial alternation set: Chebyshev extrema (Lobatto) of degree N+1 + // mapped into the domain, ascending. + let mut taus: Vec = (0..m) + .map(|k| -(core::f64::consts::PI * k as f64 / (m - 1) as f64).cos()) + .collect(); + + // Pre-allocated work buffers. + let n = degree + 1; + let mut matrix = vec![0.0; m * (m + 1)]; + let mut coeffs = vec![0.0; n]; + let mut prev_max_err = f64::INFINITY; + let mut max_err = 0.0; + let mut leveled = 0.0; + let mut iterations = 0; + let mut converged = false; + + for iter in 1..=max_iter { + iterations = iter; + + // Build linear system: row i = [T_0(tau_i) ... T_N(tau_i) | (-1)^i | f(x_i)] + for (i, &tau) in taus.iter().enumerate() { + let row = i * (m + 1); + for k in 0..n { + matrix[row + k] = basis::t(k, tau); + } + matrix[row + n] = if i % 2 == 0 { 1.0 } else { -1.0 }; + matrix[row + m] = f(domain.denormalize(tau)); + } + + if !solve_in_place(&mut matrix, m) { + return Err(ChebyError::RemezDidNotConverge); + } + + for k in 0..n { + coeffs[k] = matrix[k * (m + 1) + m]; + } + leveled = matrix[n * (m + 1) + m].abs(); + + // Sample residual on dense grid. + let mut residual = vec![0.0_f64; grid]; + for (i, r) in residual.iter_mut().enumerate() { + let tau = -1.0 + 2.0 * i as f64 / (grid - 1) as f64; + let x = domain.denormalize(tau); + *r = f(x) - crate::core::eval::evaluate(&coeffs, tau); + } + + // Locate per-sign-block extrema. + let mut blocks: Vec<(f64, f64)> = Vec::new(); // (tau_max, residual_max_signed) + let mut cur_sign = 0i8; + let mut cur_best_idx = 0usize; + let mut cur_best_abs = -1.0_f64; + for (i, &r) in residual.iter().enumerate() { + let s = if r > 0.0 { + 1 + } else if r < 0.0 { + -1 + } else { + 0 + }; + if s == 0 { + continue; + } + if cur_sign == 0 { + cur_sign = s; + cur_best_idx = i; + cur_best_abs = r.abs(); + } else if s == cur_sign { + if r.abs() > cur_best_abs { + cur_best_abs = r.abs(); + cur_best_idx = i; } + } else { + let tau = -1.0 + 2.0 * cur_best_idx as f64 / (grid - 1) as f64; + blocks.push((tau, residual[cur_best_idx])); + cur_sign = s; + cur_best_idx = i; + cur_best_abs = r.abs(); } - return Ok(RemezResult { - series, - max_error: max_err, - iterations, - converged: max_mag <= options.tolerance, - domain, - }); - }}; + } + if cur_sign != 0 { + let tau = -1.0 + 2.0 * cur_best_idx as f64 / (grid - 1) as f64; + blocks.push((tau, residual[cur_best_idx])); + } + + if blocks.len() < m { + return Err(ChebyError::RemezAlternationFailure); + } + + // Pick the contiguous window of m blocks maximizing min |r|. + let mut best_start = 0usize; + let mut best_min = -1.0_f64; + for s in 0..=blocks.len() - m { + let window_min = blocks + .iter() + .skip(s) + .take(m) + .map(|b| b.1.abs()) + .fold(f64::INFINITY, f64::min); + if window_min > best_min { + best_min = window_min; + best_start = s; + } + } + + let mut new_taus: Vec = blocks[best_start..best_start + m] + .iter() + .map(|&(t, _)| t) + .collect(); + new_taus.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal)); + + max_err = blocks[best_start..best_start + m] + .iter() + .map(|&(_, r)| r.abs()) + .fold(0.0_f64, f64::max); + + let gap = (max_err - leveled).abs(); + let denom = max_err.max(f64::EPSILON); + if gap / denom <= options.tolerance { + converged = true; + taus = new_taus; + break; + } + + // Stagnation guard. + if (prev_max_err - max_err).abs() <= options.tolerance * denom && iter > 1 { + taus = new_taus; + converged = true; + break; + } + prev_max_err = max_err; + taus = new_taus; + } + + if !converged { + return Err(ChebyError::RemezDidNotConverge); } - match degree + 1 { - 2 => fit_n!(2), - 3 => fit_n!(3), - 4 => fit_n!(4), - 5 => fit_n!(5), - 6 => fit_n!(6), - 7 => fit_n!(7), - 8 => fit_n!(8), - 9 => fit_n!(9), - 10 => fit_n!(10), - 11 => fit_n!(11), - 12 => fit_n!(12), - 13 => fit_n!(13), - 14 => fit_n!(14), - 15 => fit_n!(15), - 16 => fit_n!(16), - _ => Err(ChebyError::InvalidDegree), + let series = ChebySeriesDyn::new(coeffs.clone())?; + Ok(RemezResult { + series_on: ChebySeriesOnAny { + domain, + coeffs: coeffs.clone(), + }, + series, + leveled_error: leveled, + max_error: max_err, + iterations, + converged, + alternations: taus.into_iter().map(|t| domain.denormalize(t)).collect(), + domain, + }) +} + +/// Augmented Gaussian elimination with partial pivoting. +/// +/// Matrix layout: `m` rows, `m + 1` columns (last column is RHS). +/// On success the RHS column holds the solution. +fn solve_in_place(a: &mut [f64], m: usize) -> bool { + let stride = m + 1; + for k in 0..m { + // Pivot. + let mut piv_row = k; + let mut piv_abs = a[k * stride + k].abs(); + for r in (k + 1)..m { + let v = a[r * stride + k].abs(); + if v > piv_abs { + piv_abs = v; + piv_row = r; + } + } + if piv_abs < 1e-300 { + return false; + } + if piv_row != k { + for c in 0..stride { + a.swap(k * stride + c, piv_row * stride + c); + } + } + // Eliminate. + let pivot = a[k * stride + k]; + for r in (k + 1)..m { + let factor = a[r * stride + k] / pivot; + if factor == 0.0 { + continue; + } + for c in k..stride { + a[r * stride + c] -= factor * a[k * stride + c]; + } + } } + // Back-substitute, store solution in RHS column (`stride - 1`). + for r in (0..m).rev() { + let mut sum = a[r * stride + (stride - 1)]; + for c in (r + 1)..m { + sum -= a[r * stride + c] * a[c * stride + (stride - 1)]; + } + a[r * stride + (stride - 1)] = sum / a[r * stride + r]; + } + true +} + +/// Convert a Remez result into a fixed-size [`ChebySeriesOn`] when the +/// caller knows the degree at compile time. +pub fn into_fixed_series_on( + result: &RemezResult, +) -> Result, ChebyError> { + let series = try_into_fixed::(result.series.coeffs())?; + Ok(ChebySeriesOn::new(result.domain, series)) } diff --git a/src/approx/mod.rs b/src/approx/mod.rs index 6f9c776..15f9646 100644 --- a/src/approx/mod.rs +++ b/src/approx/mod.rs @@ -1,4 +1,25 @@ //! Chebyshev approximation and interpolation. +//! +//! # Theory +//! +//! Given a function `f` on a [`crate::core::Domain`], `fit` produces a +//! truncated Chebyshev series whose coefficients minimise the L² error +//! at the Chebyshev nodes (interpolation in the polynomial sense). +//! The optional [`adaptive`] module refines the degree until a target +//! tolerance is met, and the optional [`minimax`] module replaces the +//! L² fit with a real Remez exchange that minimises the L∞ error. +//! +//! # Features +//! +//! - `approx` (default): fixed-size and dynamic L² fitting. +//! - `adaptive`: iterative refinement (`alloc` required). +//! - `minimax`: Remez exchange returning explicit convergence diagnostics +//! (`alloc` required, `f64` only). +//! +//! # Performance +//! +//! Fixed-size fits are `O(N²)` with no allocations. Adaptive fits and +//! Remez allocate per iteration; both are `O(N²)` per pass. pub mod error_estimation; pub mod fit; diff --git a/src/calculus/mod.rs b/src/calculus/mod.rs index a185a46..7737ffa 100644 --- a/src/calculus/mod.rs +++ b/src/calculus/mod.rs @@ -1,4 +1,19 @@ //! Chebyshev calculus. +//! +//! # Theory +//! +//! Differentiation maps `Σ aₖ Tₖ` to `Σ a'ₖ Tₖ` via the standard upward +//! recurrence; integration runs the same recurrence in reverse and adds +//! a chosen constant. Both operators are exact on the polynomial space +//! and stable for the degrees `cheby` targets. +//! +//! # Features +//! +//! `calculus` (default), `no_std`, allocation-free. +//! +//! # Performance +//! +//! `O(N)` per series, in place, on the const-generic coefficient array. pub mod derivative; pub mod integral; diff --git a/src/core/error.rs b/src/core/error.rs index 8525322..6af1d63 100644 --- a/src/core/error.rs +++ b/src/core/error.rs @@ -20,22 +20,60 @@ pub enum ChebyError { NonFiniteInput, /// A segment table was constructed without segments. EmptySegmentTable, + /// Adjacent segments did not meet at a shared boundary. + SegmentBoundariesNotContiguous, + /// A binary blob is shorter than the minimum required header/footer. + BinaryTooShort, + /// A binary blob carries a format version this build does not support. + UnsupportedFormatVersion { + /// Version found in the blob. + found: u32, + /// Version expected by this build. + expected: u32, + }, + /// A binary blob's stored length does not match its payload size. + BinaryLengthMismatch, + /// A binary blob's checksum did not match its payload. + BinaryChecksumMismatch, + /// Remez exchange did not converge within the requested iteration budget. + RemezDidNotConverge, + /// Remez exchange could not find enough alternation points. + RemezAlternationFailure, } #[cfg(feature = "std")] impl std::fmt::Display for ChebyError { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - let message = match self { - Self::EmptyDomain => "domain has zero width", - Self::InvalidDomain => "domain is invalid", - Self::NonPositiveSegmentLength => "segment length must be positive", - Self::InvalidDegree => "degree is invalid", - Self::EmptyCoefficientSet => "coefficient set is empty", - Self::EvaluationOutOfDomain => "evaluation point is outside the domain", - Self::NonFiniteInput => "input is not finite", - Self::EmptySegmentTable => "segment table is empty", - }; - f.write_str(message) + match self { + Self::EmptyDomain => f.write_str("domain has zero width"), + Self::InvalidDomain => f.write_str("domain is invalid"), + Self::NonPositiveSegmentLength => f.write_str("segment length must be positive"), + Self::InvalidDegree => f.write_str("degree is invalid"), + Self::EmptyCoefficientSet => f.write_str("coefficient set is empty"), + Self::EvaluationOutOfDomain => f.write_str("evaluation point is outside the domain"), + Self::NonFiniteInput => f.write_str("input is not finite"), + Self::EmptySegmentTable => f.write_str("segment table is empty"), + Self::SegmentBoundariesNotContiguous => { + f.write_str("segment boundaries are not contiguous") + } + Self::BinaryTooShort => f.write_str("binary blob is shorter than the format header"), + Self::UnsupportedFormatVersion { found, expected } => { + write!( + f, + "unsupported binary format version {found} (expected {expected})" + ) + } + Self::BinaryLengthMismatch => { + f.write_str("binary blob length does not match the stored coefficient count") + } + Self::BinaryChecksumMismatch => f.write_str("binary blob checksum mismatch"), + Self::RemezDidNotConverge => { + f.write_str("Remez exchange did not converge within the iteration budget") + } + Self::RemezAlternationFailure => { + f.write_str("Remez exchange could not find enough alternation points") + } + } } } diff --git a/src/core/mod.rs b/src/core/mod.rs index 7337124..698d031 100644 --- a/src/core/mod.rs +++ b/src/core/mod.rs @@ -1,4 +1,24 @@ //! Core Chebyshev primitives. +//! +//! # Theory +//! +//! Every approximation in `cheby` is a finite Chebyshev expansion +//! `f(x) ≈ Σ aₖ Tₖ(τ)` over a normalised parameter `τ ∈ [-1, 1]`. +//! [`Domain`] handles the affine map between user coordinates and `τ`, +//! [`nodes`] returns Chebyshev abscissae (`Gauss` or `GaussLobatto`), +//! [`evaluate`] uses Clenshaw's recurrence, and [`ChebySeries`] / +//! [`ChebySeriesOn`] hold coefficients with optional unit safety via +//! [`ChebyScalar`]. +//! +//! # Features +//! +//! Always available; `no_std` and allocation-free. Enable `alloc` for +//! the dynamic [`ChebySeriesDyn`] type. +//! +//! # Performance +//! +//! All hot paths are `const`-generic on the polynomial degree. Evaluation +//! is `O(N)` per query with no heap traffic; node generation is `O(N)`. pub mod basis; pub mod domain; diff --git a/src/io/binary.rs b/src/io/binary.rs index 39ab0f6..289c45e 100644 --- a/src/io/binary.rs +++ b/src/io/binary.rs @@ -1,13 +1,32 @@ //! Compact binary coefficient-table helpers for `f64` series. +//! +//! # Format (little-endian) +//! +//! ```text +//! offset size field +//! ------ ---- ---------------------------------- +//! 0 4 u32 format_version (== FORMAT_VERSION) +//! 4 8 u64 coefficient_count +//! 12 8*N N * f64 coefficients +//! 12+8N 8 u64 FNV-1a checksum over [0, 12+8N) +//! ``` +//! +//! [`decode_f64_series`] reports each malformation with its own typed +//! [`ChebyError`] variant, so callers can distinguish a truncated blob from +//! a version mismatch from a checksum failure. use alloc::vec::Vec; use crate::core::{ChebyError, ChebySeriesDyn}; use crate::io::metadata::FORMAT_VERSION; -/// Encode an `f64` dynamic series. +const HEADER_LEN: usize = 12; +const CHECKSUM_LEN: usize = 8; +const MIN_LEN: usize = HEADER_LEN + CHECKSUM_LEN; + +/// Encode an `f64` dynamic series using the [module format](self). pub fn encode_f64_series(series: &ChebySeriesDyn) -> Vec { - let mut out = Vec::with_capacity(16 + 8 * series.coeffs().len()); + let mut out = Vec::with_capacity(MIN_LEN + 8 * series.coeffs().len()); out.extend_from_slice(&FORMAT_VERSION.to_le_bytes()); out.extend_from_slice(&(series.coeffs().len() as u64).to_le_bytes()); for coeff in series.coeffs() { @@ -18,26 +37,29 @@ pub fn encode_f64_series(series: &ChebySeriesDyn) -> Vec { out } -/// Decode an `f64` dynamic series. +/// Decode an `f64` dynamic series produced by [`encode_f64_series`]. pub fn decode_f64_series(bytes: &[u8]) -> Result, ChebyError> { - if bytes.len() < 20 { - return Err(ChebyError::EmptyCoefficientSet); + if bytes.len() < MIN_LEN { + return Err(ChebyError::BinaryTooShort); } - let payload_len = bytes.len() - 8; + let payload_len = bytes.len() - CHECKSUM_LEN; let expected = u64::from_le_bytes(bytes[payload_len..].try_into().unwrap()); if checksum(&bytes[..payload_len]) != expected { - return Err(ChebyError::InvalidDomain); + return Err(ChebyError::BinaryChecksumMismatch); } let version = u32::from_le_bytes(bytes[0..4].try_into().unwrap()); if version != FORMAT_VERSION { - return Err(ChebyError::InvalidDegree); + return Err(ChebyError::UnsupportedFormatVersion { + found: version, + expected: FORMAT_VERSION, + }); } - let len = u64::from_le_bytes(bytes[4..12].try_into().unwrap()) as usize; - if payload_len != 12 + 8 * len { - return Err(ChebyError::InvalidDegree); + let len = u64::from_le_bytes(bytes[4..HEADER_LEN].try_into().unwrap()) as usize; + if payload_len != HEADER_LEN + 8 * len { + return Err(ChebyError::BinaryLengthMismatch); } let mut coeffs = Vec::with_capacity(len); - for chunk in bytes[12..payload_len].chunks_exact(8) { + for chunk in bytes[HEADER_LEN..payload_len].chunks_exact(8) { coeffs.push(f64::from_le_bytes(chunk.try_into().unwrap())); } ChebySeriesDyn::new(coeffs) diff --git a/src/io/mod.rs b/src/io/mod.rs index 8b58c21..c17f1cb 100644 --- a/src/io/mod.rs +++ b/src/io/mod.rs @@ -1,4 +1,12 @@ //! Optional serialization helpers. +//! +//! # Features +//! +//! - `binary`: a compact, versioned, checksummed `f64` coefficient layout +//! (see [`binary`]). +//! - `serde`: derives that integrate with the `serde` ecosystem. +//! +//! Both submodules require `alloc`. #[cfg(feature = "binary")] pub mod binary; diff --git a/src/lib.rs b/src/lib.rs index ec5d4b1..151a339 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2,6 +2,7 @@ // Copyright (C) 2026 Vallés Puig, Ramon #![forbid(unsafe_code)] +#![warn(missing_docs)] #![cfg_attr(not(feature = "std"), no_std)] #![doc = include_str!("../README.md")] diff --git a/src/piecewise/adaptive.rs b/src/piecewise/adaptive.rs index b84f035..12cfda5 100644 --- a/src/piecewise/adaptive.rs +++ b/src/piecewise/adaptive.rs @@ -1,10 +1,18 @@ //! Adaptive segment tables. +//! +//! Segments are produced by recursive bisection: a candidate fit is computed +//! on each subdomain, its tail-error estimate is compared against the +//! tolerance, and on failure the subdomain is split at its midpoint. The +//! resulting segments are stored in ascending-start order with their +//! boundaries cached as a separate sorted `Vec` so that +//! [`AdaptiveSegmentTable::evaluate`] can locate the active segment in +//! `O(log n)` via binary search rather than a linear scan. use alloc::vec::Vec; use crate::approx::error_estimation; use crate::approx::fit; -use crate::core::{ChebyError, ChebyScalar, ChebySeries, ChebyTime, Domain}; +use crate::core::{ChebyError, ChebyScalar, ChebySeries, ChebyTime, DifferentiateWith, Domain}; use crate::piecewise::ChebySegment; /// Metadata for an adaptive segment. @@ -13,17 +21,25 @@ use crate::piecewise::ChebySegment; pub struct SegmentMetadata { /// Number of coefficients stored. pub degree: usize, - /// Estimated local error. + /// Estimated local truncation error. pub estimated_error: Option, - /// Number of function samples used. + /// Number of function samples used for the local fit. pub samples_used: usize, + /// Recursion depth at which this segment was emitted. + pub depth: usize, } -/// A simple adaptive segment table. +/// An adaptive piecewise Chebyshev table over a single base domain. +/// +/// Internally segments are stored in ascending order, plus a cached +/// `boundaries` vector of length `segments.len() + 1` (start of each segment +/// followed by the end of the last). Lookup uses [`slice::binary_search_by`]. #[derive(Debug, Clone, PartialEq)] pub struct AdaptiveSegmentTable { + domain: Domain, segments: Vec>, metadata: Vec>, + boundaries: Vec, } impl AdaptiveSegmentTable @@ -31,13 +47,19 @@ where T: ChebyScalar, X: ChebyTime, { - /// Build by recursively splitting until a tail estimate meets tolerance. + /// Build by recursively splitting until a tail estimate meets `tolerance`. pub fn from_fn( domain: Domain, f: impl Fn(X) -> T + Copy, tolerance: f64, max_depth: usize, ) -> Result { + if N < 2 { + return Err(ChebyError::InvalidDegree); + } + if !(tolerance.is_finite() && tolerance >= 0.0) { + return Err(ChebyError::NonFiniteInput); + } let mut segments = Vec::new(); let mut metadata = Vec::new(); split::( @@ -45,32 +67,125 @@ where f, tolerance, max_depth, + 0, &mut segments, &mut metadata, )?; - Ok(Self { segments, metadata }) + if segments.is_empty() { + return Err(ChebyError::EmptySegmentTable); + } + // Verify ordering and contiguity, then cache boundaries. + let mut boundaries = Vec::with_capacity(segments.len() + 1); + boundaries.push(segments[0].domain().start()); + for (i, seg) in segments.iter().enumerate() { + if i > 0 && seg.domain().start() != segments[i - 1].domain().end() { + return Err(ChebyError::SegmentBoundariesNotContiguous); + } + boundaries.push(seg.domain().end()); + } + Ok(Self { + domain, + segments, + metadata, + boundaries, + }) + } + + /// Base domain spanning every segment. + #[inline] + pub fn domain(&self) -> Domain { + self.domain + } + + /// Number of segments. + #[inline] + pub fn len(&self) -> usize { + self.segments.len() + } + + /// Whether the table is empty. + #[inline] + pub fn is_empty(&self) -> bool { + self.segments.is_empty() + } + + /// Borrow the segment slice in ascending domain order. + #[inline] + pub fn segments(&self) -> &[ChebySegment] { + &self.segments + } + + /// Per-segment metadata, in the same order as [`Self::segments`]. + #[inline] + pub fn metadata(&self) -> &[SegmentMetadata] { + &self.metadata } - /// Evaluate by binary-searching segment starts. + /// Cached monotone boundary vector. Length is `len() + 1`. + #[inline] + pub fn boundaries(&self) -> &[X] { + &self.boundaries + } + + /// Locate the segment containing `x` via binary search. + pub fn locate(&self, x: X) -> Option<&ChebySegment> { + if self.segments.is_empty() { + return None; + } + if x < self.boundaries[0] || x > self.boundaries[self.boundaries.len() - 1] { + return None; + } + // Binary-search the boundary vector for the largest start <= x. + let probe = self + .boundaries + .binary_search_by(|b| b.partial_cmp(&x).unwrap_or(core::cmp::Ordering::Less)); + let idx = match probe { + Ok(i) => { + if i == self.segments.len() { + self.segments.len() - 1 + } else { + i + } + } + Err(i) => { + if i == 0 { + 0 + } else { + i - 1 + } + } + }; + self.segments.get(idx).filter(|s| s.contains(x)) + } + + /// Evaluate at `x`. pub fn evaluate(&self, x: X) -> Result { - self.segments - .iter() - .find(|segment| segment.contains(x)) + self.locate(x) .ok_or(ChebyError::EvaluationOutOfDomain)? .evaluate(x) } - /// Segment metadata. - pub fn metadata(&self) -> &[SegmentMetadata] { - &self.metadata + /// Evaluate the dimensionally-correct derivative at `x`. + pub fn evaluate_derivative( + &self, + x: X, + ) -> Result<>::Derivative, ChebyError> + where + T: DifferentiateWith, + { + self.locate(x) + .ok_or(ChebyError::EvaluationOutOfDomain)? + .evaluate_derivative(x) } } +#[allow(clippy::too_many_arguments)] fn split( domain: Domain, f: impl Fn(X) -> T + Copy, tolerance: f64, - depth: usize, + depth_remaining: usize, + depth_so_far: usize, segments: &mut Vec>, metadata: &mut Vec>, ) -> Result<(), ChebyError> @@ -81,12 +196,13 @@ where let fitted = fit::fit_from_fn_on::(domain, f); let coeffs = fitted.into_coeffs(); let err = error_estimation::estimated_truncation_error(&coeffs); - if err <= tolerance || depth == 0 { + if err <= tolerance || depth_remaining == 0 { segments.push(ChebySegment::try_new(domain, ChebySeries::new(coeffs))?); metadata.push(SegmentMetadata { degree: N, estimated_error: Some(err), samples_used: N, + depth: depth_so_far, }); return Ok(()); } @@ -95,7 +211,8 @@ where Domain::try_new(domain.start(), mid)?, f, tolerance, - depth - 1, + depth_remaining - 1, + depth_so_far + 1, segments, metadata, )?; @@ -103,9 +220,9 @@ where Domain::try_new(mid, domain.end())?, f, tolerance, - depth - 1, + depth_remaining - 1, + depth_so_far + 1, segments, metadata, - )?; - Ok(()) + ) } diff --git a/src/piecewise/mod.rs b/src/piecewise/mod.rs index 60782c6..a8f831c 100644 --- a/src/piecewise/mod.rs +++ b/src/piecewise/mod.rs @@ -1,4 +1,22 @@ //! Piecewise Chebyshev approximations. +//! +//! # Theory +//! +//! A piecewise table partitions a domain into segments, fits a Chebyshev +//! series on each, and dispatches lookups by segment. [`ChebySegmentTable`] +//! uses a uniform partition; [`AdaptiveSegmentTable`] subdivides by error +//! tolerance and stores the resulting monotone segment boundaries for +//! `O(log S)` lookup. +//! +//! # Features +//! +//! `piecewise` (requires `alloc`). +//! +//! # Performance +//! +//! Uniform tables locate the active segment in `O(1)`; adaptive tables in +//! `O(log S)` via [`core::slice::binary_search_by`]. Per-segment evaluation +//! is `O(N)` Clenshaw. pub mod adaptive; pub mod lookup; diff --git a/src/quadrature/clenshaw_curtis.rs b/src/quadrature/clenshaw_curtis.rs index 9f5be09..b88aac1 100644 --- a/src/quadrature/clenshaw_curtis.rs +++ b/src/quadrature/clenshaw_curtis.rs @@ -1,38 +1,61 @@ -//! Clenshaw-Curtis quadrature. +//! Clenshaw-Curtis quadrature on Chebyshev-Lobatto nodes. +//! +//! Computes weights `w_k` such that `sum_k w_k f(x_k)` integrates polynomials +//! of degree up to `N - 1` exactly on `[-1, 1]`, where `N` is the node count +//! and `x_k = cos(k π / (N - 1))` are the Chebyshev extrema (Lobatto nodes). +//! +//! The implementation follows Trefethen, *Spectral Methods in MATLAB*, +//! Program 30 (`clencurt`). Cost is `O(N^2)`. use crate::core::{ChebyScalar, ChebyTime, Domain, IntegrateWith, NodeKind}; -/// Return simple Clenshaw-Curtis-like weights on `[-1, 1]`. -/// -/// The implementation uses trapezoidal weights over Chebyshev-Lobatto -/// abscissae. It is stable for examples and small rules, and keeps the API -/// lightweight. +/// Clenshaw-Curtis weights for `N` Lobatto nodes on `[-1, 1]`. pub fn clenshaw_curtis_weights() -> [f64; N] { + let mut weights = [0.0; N]; if N == 0 { - return [0.0; N]; + return weights; } if N == 1 { - return [2.0; N]; + weights[0] = 2.0; + return weights; } - let xs = crate::core::nodes::nodes::(NodeKind::GaussLobatto); - let mut weights = [0.0; N]; - for k in 0..N { - let left = if k == 0 { - xs[0] - } else { - 0.5 * (xs[k - 1] + xs[k]) - }; - let right = if k + 1 == N { - xs[N - 1] - } else { - 0.5 * (xs[k] + xs[k + 1]) - }; - weights[k] = (left - right).abs(); + if N == 2 { + weights[0] = 1.0; + weights[1] = 1.0; + return weights; } + + let m = N - 1; // Trefethen's "N" + let theta: [f64; N] = + core::array::from_fn(|k| core::f64::consts::PI * k as f64 / m as f64); + + // Endpoint weights. + if m.is_multiple_of(2) { + weights[0] = 1.0 / ((m as f64) * (m as f64) - 1.0); + weights[N - 1] = weights[0]; + } else { + weights[0] = 1.0 / (m as f64 * m as f64); + weights[N - 1] = weights[0]; + } + + // Interior weights. + for (k, w) in weights.iter_mut().enumerate().take(N - 1).skip(1) { + let mut v = 1.0_f64; + let half = m / 2; + let upper = if m.is_multiple_of(2) { half - 1 } else { half }; + for j in 1..=upper { + v -= 2.0 * (2.0 * j as f64 * theta[k]).cos() / (4.0 * (j * j) as f64 - 1.0); + } + if m.is_multiple_of(2) { + v -= (m as f64 * theta[k]).cos() / ((m as f64) * (m as f64) - 1.0); + } + *w = 2.0 * v / m as f64; + } + weights } -/// Integrate a function over `domain`. +/// Integrate a function over `domain` using `N`-point Clenshaw-Curtis. pub fn integrate( domain: Domain, f: impl Fn(X) -> T, diff --git a/src/quadrature/mod.rs b/src/quadrature/mod.rs index f20f078..b570246 100644 --- a/src/quadrature/mod.rs +++ b/src/quadrature/mod.rs @@ -1,4 +1,21 @@ //! Chebyshev-related quadrature. +//! +//! # Theory +//! +//! [`clenshaw_curtis`] integrates over `[-1, 1]` with weights derived from +//! the discrete cosine transform on Chebyshev-Lobatto nodes; it is exact +//! on polynomials of degree up to `N - 1`. [`gauss_chebyshev`] integrates +//! against the weight `1/√(1−x²)` using the closed-form Chebyshev nodes +//! and uniform weight `π/N`. +//! +//! # Features +//! +//! `quadrature`, `no_std`, allocation-free. +//! +//! # Performance +//! +//! Both rules are `O(N²)` to build the weight set and `O(N)` per integrand +//! evaluation. pub mod clenshaw_curtis; pub mod gauss_chebyshev; diff --git a/src/spectral/mod.rs b/src/spectral/mod.rs index 8a13c3e..ec17cb4 100644 --- a/src/spectral/mod.rs +++ b/src/spectral/mod.rs @@ -1,4 +1,20 @@ //! Chebyshev spectral numerics. +//! +//! # Theory +//! +//! [`chebyshev_differentiation_matrix`] builds the dense `N×N` matrix `D` +//! such that `D · f` approximates `f'` at the Lobatto collocation nodes +//! returned by [`collocation_points`]. Convergence is spectral +//! (super-algebraic) for analytic integrands. +//! +//! # Features +//! +//! `spectral` (requires `calculus` and `alloc`). +//! +//! # Performance +//! +//! Matrix construction is `O(N²)`. A subsequent multiply against a sample +//! vector is also `O(N²)` and dominates the cost. pub mod collocation; pub mod differentiation; diff --git a/tests/adaptive_piecewise.rs b/tests/adaptive_piecewise.rs new file mode 100644 index 0000000..55a667a --- /dev/null +++ b/tests/adaptive_piecewise.rs @@ -0,0 +1,79 @@ +//! Coverage for adaptive piecewise tables and binary-search lookup. +#![cfg(feature = "piecewise")] + +use approx::assert_abs_diff_eq; +use cheby::core::{ChebyError, Domain}; +use cheby::piecewise::AdaptiveSegmentTable; + +#[test] +fn adaptive_table_splits_for_steep_region() { + let domain = Domain::new(-1.0, 1.0); + // tanh(20 x) is steep near zero and forces splitting. + let table: AdaptiveSegmentTable = + AdaptiveSegmentTable::from_fn(domain, |x: f64| (20.0 * x).tanh(), 1e-8, 8).unwrap(); + assert!(table.len() > 1, "expected >1 segments, got {}", table.len()); + // Boundaries are sorted and contiguous. + let b = table.boundaries(); + for w in b.windows(2) { + assert!(w[0] <= w[1]); + } + assert_eq!(b.first().copied().unwrap(), -1.0); + assert_eq!(b.last().copied().unwrap(), 1.0); + // Evaluation accuracy. + for i in 0..50 { + let x = -1.0 + 2.0 * i as f64 / 49.0; + let approx = table.evaluate(x).unwrap(); + let exact = (20.0 * x).tanh(); + assert_abs_diff_eq!(approx, exact, epsilon = 1e-6); + } +} + +#[test] +fn adaptive_table_lookup_is_o_log_n_correct() { + let domain = Domain::new(0.0, 16.0); + let table: AdaptiveSegmentTable = + AdaptiveSegmentTable::from_fn(domain, |x: f64| (5.0 * x).sin(), 1e-6, 10).unwrap(); + + // Check every probe lands on a segment whose domain contains it. + for i in 0..200 { + let x = 16.0 * i as f64 / 199.0; + let seg = table.locate(x).expect("must find a segment"); + assert!(seg.contains(x), "segment must contain probed x"); + } +} + +#[test] +fn adaptive_table_out_of_domain_errors() { + let domain = Domain::new(0.0, 1.0); + let table: AdaptiveSegmentTable = + AdaptiveSegmentTable::from_fn(domain, |x: f64| x * x, 1e-12, 4).unwrap(); + assert_eq!( + table.evaluate(-0.1).unwrap_err(), + ChebyError::EvaluationOutOfDomain + ); + assert_eq!( + table.evaluate(1.5).unwrap_err(), + ChebyError::EvaluationOutOfDomain + ); +} + +#[test] +fn adaptive_table_metadata_records_depth() { + let domain = Domain::new(-1.0, 1.0); + let table: AdaptiveSegmentTable = + AdaptiveSegmentTable::from_fn(domain, |x: f64| (50.0 * x).tanh(), 1e-6, 8).unwrap(); + let meta = table.metadata(); + assert_eq!(meta.len(), table.len()); + assert!(meta.iter().any(|m| m.depth > 0)); +} + +#[test] +fn adaptive_table_rejects_invalid_inputs() { + let domain = Domain::new(-1.0, 1.0); + let r: Result, _> = + AdaptiveSegmentTable::from_fn(domain, |_| 0.0, f64::NAN, 4); + assert_eq!(r.unwrap_err(), ChebyError::NonFiniteInput); + let r: Result, _> = + AdaptiveSegmentTable::from_fn(domain, |_| 0.0, -1.0, 4); + assert_eq!(r.unwrap_err(), ChebyError::NonFiniteInput); +} diff --git a/tests/calculus_roundtrip.rs b/tests/calculus_roundtrip.rs new file mode 100644 index 0000000..e8cc47f --- /dev/null +++ b/tests/calculus_roundtrip.rs @@ -0,0 +1,48 @@ +//! Calculus round-trip and integral correctness. + +use approx::assert_abs_diff_eq; +use cheby::approx::fit::fit_from_fn_on; +use cheby::core::{ChebySeries, Domain}; + +#[test] +fn integral_then_derivative_recovers_original() { + // Pad with trailing zeros so the antiderivative stays representable. + let series = ChebySeries::new([1.5_f64, -0.4, 0.2, 0.0, 0.0]); + let recovered = series.integral(0.0).derivative(); + for tau in [-0.9, -0.3, 0.0, 0.4, 0.85] { + assert_abs_diff_eq!( + recovered.evaluate(tau), + series.evaluate(tau), + epsilon = 1e-12 + ); + } +} + +#[test] +fn definite_integral_matches_analytic() { + // ∫_{0}^{x} sin(t) dt = 1 - cos(x), and we test for x in (0, π/2). + let domain = Domain::new(0.0, core::f64::consts::PI); + let series = fit_from_fn_on::(domain, f64::sin); + for &x in &[0.5_f64, 1.0, 1.5, 2.5, 3.0] { + let v = series.evaluate_integral_from_start(x).unwrap(); + let exact = 1.0 - x.cos(); + assert_abs_diff_eq!(v, exact, epsilon = 1e-9); + } +} + +#[test] +fn typed_integral_velocity_to_position() { + use qtty::{Kilometer, Second}; + type Velocity = qtty::velocity::Velocity; + + let domain = Domain::new(Second::new(0.0), Second::new(10.0)); + // Constant velocity = 2 km/s ⇒ position(t) = 2t. + let series = cheby::ChebySeriesOn::new( + domain, + ChebySeries::new([Velocity::new(2.0), Velocity::new(0.0)]), + ); + let pos: Kilometer = series + .evaluate_integral_from_start(Second::new(3.0)) + .unwrap(); + assert_abs_diff_eq!(pos.value(), 6.0, epsilon = 1e-12); +} diff --git a/tests/error_branches.rs b/tests/error_branches.rs new file mode 100644 index 0000000..7301a9b --- /dev/null +++ b/tests/error_branches.rs @@ -0,0 +1,97 @@ +//! Coverage for every public [`ChebyError`] branch. + +use cheby::core::{ChebyError, ChebySeriesDyn, Domain}; + +#[test] +fn domain_error_branches() { + assert_eq!( + Domain::::try_new(1.0, 1.0).unwrap_err(), + ChebyError::EmptyDomain + ); + assert_eq!( + Domain::::try_new(1.0, 0.0).unwrap_err(), + ChebyError::InvalidDomain + ); + assert_eq!( + Domain::::try_new(f64::NAN, 1.0).unwrap_err(), + ChebyError::NonFiniteInput + ); + assert_eq!( + Domain::::try_new(0.0, f64::INFINITY).unwrap_err(), + ChebyError::NonFiniteInput + ); +} + +#[test] +fn dyn_series_empty_set() { + assert_eq!( + ChebySeriesDyn::::new(vec![]).unwrap_err(), + ChebyError::EmptyCoefficientSet + ); +} + +#[test] +fn series_on_evaluate_out_of_domain() { + let domain = Domain::new(0.0, 1.0); + let series = cheby::ChebySeriesOn::new(domain, cheby::ChebySeries::new([1.0, 0.0])); + assert_eq!( + series.evaluate(2.0).unwrap_err(), + ChebyError::EvaluationOutOfDomain + ); +} + +#[test] +#[cfg(feature = "binary")] +fn binary_error_branches() { + use cheby::core::ChebySeriesDyn; + use cheby::io::binary::{decode_f64_series, encode_f64_series}; + + // Round-trip succeeds. + let original = ChebySeriesDyn::new(vec![1.0, -2.5, 0.125, 7.0]).unwrap(); + let bytes = encode_f64_series(&original); + let decoded = decode_f64_series(&bytes).unwrap(); + assert_eq!(decoded.coeffs(), original.coeffs()); + + // Truncated. + assert_eq!( + decode_f64_series(&bytes[..5]).unwrap_err(), + ChebyError::BinaryTooShort + ); + + // Corrupt checksum. + let mut corrupted = bytes.clone(); + let last = corrupted.len() - 1; + corrupted[last] ^= 0xFF; + assert_eq!( + decode_f64_series(&corrupted).unwrap_err(), + ChebyError::BinaryChecksumMismatch + ); + + // Wrong version (recompute checksum so we hit version branch, not checksum). + let mut v = bytes.clone(); + v[0] = 99; // version byte + // Recompute the FNV-1a checksum the encoder uses. + let payload_len = v.len() - 8; + let new_ck = v[..payload_len].iter().fold(0xcbf29ce484222325_u64, |h, b| { + (h ^ u64::from(*b)).wrapping_mul(0x100000001b3) + }); + v[payload_len..].copy_from_slice(&new_ck.to_le_bytes()); + let err = decode_f64_series(&v).unwrap_err(); + assert!(matches!( + err, + ChebyError::UnsupportedFormatVersion { found: _, expected: _ } + )); + + // Length mismatch: tamper coefficient count, recompute checksum. + let mut l = bytes; + l[4..12].copy_from_slice(&999_u64.to_le_bytes()); + let payload_len = l.len() - 8; + let new_ck = l[..payload_len].iter().fold(0xcbf29ce484222325_u64, |h, b| { + (h ^ u64::from(*b)).wrapping_mul(0x100000001b3) + }); + l[payload_len..].copy_from_slice(&new_ck.to_le_bytes()); + assert_eq!( + decode_f64_series(&l).unwrap_err(), + ChebyError::BinaryLengthMismatch + ); +} diff --git a/tests/minimax.rs b/tests/minimax.rs new file mode 100644 index 0000000..ff7112f --- /dev/null +++ b/tests/minimax.rs @@ -0,0 +1,96 @@ +//! Real Remez exchange — convergence and failure modes. +#![cfg(feature = "minimax")] + +use approx::assert_abs_diff_eq; +use cheby::approx::minimax::{remez, RemezOptions}; +use cheby::core::{ChebyError, Domain}; + +#[test] +fn remez_converges_for_smooth_target() { + let result = remez( + Domain::new(-1.0, 1.0), + 7, + f64::exp, + RemezOptions::default(), + ) + .unwrap(); + assert!(result.converged); + assert_eq!(result.alternations.len(), 9); // degree + 2 + assert!(result.max_error < 1e-6); + // Equioscillation: dense max error should be very close to leveled error. + assert_abs_diff_eq!(result.max_error, result.leveled_error, epsilon = 1e-6); +} + +#[test] +fn remez_alternations_are_monotone() { + let result = remez( + Domain::new(-1.0, 1.0), + 5, + f64::cos, + RemezOptions::default(), + ) + .unwrap(); + for w in result.alternations.windows(2) { + assert!(w[0] < w[1]); + } +} + +#[test] +fn remez_rejects_zero_degree() { + let err = remez(Domain::new(-1.0, 1.0), 0, f64::sin, RemezOptions::default()).unwrap_err(); + assert_eq!(err, ChebyError::InvalidDegree); +} + +#[test] +fn remez_typed_domain_evaluates() { + use qtty::Second; + let result = remez( + Domain::new(Second::new(0.0), Second::new(1.0)), + 5, + |t: Second| t.value().exp(), + RemezOptions::default(), + ) + .unwrap(); + assert!(result.converged); + let v = result.series_on.evaluate(Second::new(0.5)).unwrap(); + assert_abs_diff_eq!(v, 0.5_f64.exp(), epsilon = 1e-4); +} + +#[test] +fn remez_does_not_converge_with_tiny_budget() { + let opts = RemezOptions { + max_iterations: 1, + tolerance: 1e-30, + grid_size: 64, + }; + // exp on [-1,1] with degree 3 won't reach 1e-30 in a single iteration. + let res = remez(Domain::new(-1.0, 1.0), 3, f64::exp, opts); + assert!(matches!(res, Err(ChebyError::RemezDidNotConverge))); +} + +#[test] +fn remez_converges_better_than_root_fit() { + let opts = RemezOptions::default(); + let degree = 6; + let domain = Domain::new(-1.0, 1.0); + let f = |x: f64| x.exp(); + + let remez_res = remez(domain, degree, f, opts).unwrap(); + + // Root-sample fit at degree+1 nodes (interpolating polynomial). + let root_series = cheby::approx::fit::fit_from_fn::(f); + let root_max_err = (0..1024) + .map(|i| { + let tau = -1.0 + 2.0 * i as f64 / 1023.0; + (f(tau) - root_series.evaluate(tau)).abs() + }) + .fold(0.0_f64, f64::max); + + // Minimax must be at least as good; almost always strictly better. + assert!( + remez_res.max_error <= root_max_err + 1e-12, + "remez={} root={}", + remez_res.max_error, + root_max_err + ); +} diff --git a/tests/quadrature_and_spectral.rs b/tests/quadrature_and_spectral.rs new file mode 100644 index 0000000..2cb8397 --- /dev/null +++ b/tests/quadrature_and_spectral.rs @@ -0,0 +1,80 @@ +//! Accuracy checks for quadrature and spectral differentiation. +#![allow(clippy::needless_range_loop)] +#![cfg(all(feature = "quadrature", feature = "spectral"))] + +use approx::assert_abs_diff_eq; +use cheby::core::{Domain, NodeKind}; +use cheby::quadrature::{clenshaw_curtis, gauss_chebyshev}; +use cheby::spectral::{chebyshev_differentiation_matrix, collocation_points}; + +#[test] +fn clenshaw_curtis_integrates_polynomial() { + // ∫_{0}^{2} x^3 dx = 4 + let v = clenshaw_curtis::integrate::(Domain::new(0.0, 2.0), |x| x * x * x); + assert_abs_diff_eq!(v, 4.0, epsilon = 1e-9); +} + +#[test] +fn clenshaw_curtis_integrates_smooth_function() { + // ∫_{0}^{π} sin x dx = 2 + let v = clenshaw_curtis::integrate::( + Domain::new(0.0, core::f64::consts::PI), + f64::sin, + ); + assert_abs_diff_eq!(v, 2.0, epsilon = 1e-3); +} + +#[test] +fn gauss_chebyshev_recovers_known_weighted_integrals() { + // ∫_{-1}^{1} 1 / sqrt(1 - x^2) dx = π + let v = gauss_chebyshev::integrate_weighted::(|_| 1.0); + assert_abs_diff_eq!(v, core::f64::consts::PI, epsilon = 1e-12); + + // ∫_{-1}^{1} cos(x) / sqrt(1 - x^2) dx = π * J_0(1) + let j0_1 = 0.7651976865579665_f64; + let v = gauss_chebyshev::integrate_weighted::(f64::cos); + assert_abs_diff_eq!(v, core::f64::consts::PI * j0_1, epsilon = 1e-10); + + // Polynomial: x^2; closed form is π/2. + let v = gauss_chebyshev::integrate_weighted::(|x| x * x); + assert_abs_diff_eq!(v, core::f64::consts::PI / 2.0, epsilon = 1e-12); +} + +#[test] +fn spectral_differentiation_matches_analytic_derivative() { + // For f(x) = sin(x), differentiate at Lobatto nodes: D * f ≈ cos(x). + const N: usize = 17; + let nodes: [f64; N] = collocation_points::(); + let f: [f64; N] = std::array::from_fn(|k| nodes[k].sin()); + let d = chebyshev_differentiation_matrix(N); + + for i in 0..N { + let mut acc = 0.0; + for j in 0..N { + acc += d.get(i, j) * f[j]; + } + let exact = nodes[i].cos(); + assert_abs_diff_eq!(acc, exact, epsilon = 1e-9); + } +} + +#[test] +fn spectral_differentiation_zero_on_constant() { + const N: usize = 9; + let f = [3.5_f64; N]; + let d = chebyshev_differentiation_matrix(N); + for i in 0..N { + let mut acc = 0.0; + for j in 0..N { + acc += d.get(i, j) * f[j]; + } + assert_abs_diff_eq!(acc, 0.0, epsilon = 1e-12); + } +} + +#[test] +fn lobatto_nodes_are_endpoints_inclusive() { + let nodes: [f64; 5] = cheby::core::nodes::nodes(NodeKind::GaussLobatto); + assert_abs_diff_eq!(nodes[0], 1.0, epsilon = 1e-15); + assert_abs_diff_eq!(nodes[4], -1.0, epsilon = 1e-15); +} From 09346739cf017d644580d14e33ba0e94eccc07b3 Mon Sep 17 00:00:00 2001 From: VPRamon Date: Wed, 6 May 2026 22:46:40 +0200 Subject: [PATCH 2/2] Implement Barycentric interpolation improvements and optimizations - Enhanced `BarycentricInterpolator` to handle single-node cases correctly, ensuring constant evaluation without division by zero. - Added methods for constructing interpolators on Chebyshev roots and Lobatto nodes with closed-form weights for improved performance and stability. - Introduced a `DefiniteIntegralFromStart` struct for efficient evaluation of definite integrals from the start of a domain, reducing overhead for repeated queries. - Updated the `integral` method in `ChebySeriesDyn` to support non-truncating indefinite integrals. - Improved error handling for out-of-domain evaluations in integral calculations. - Added tests to verify the correctness of new features and ensure regression prevention. - Refactored existing code for clarity and performance, including optimizations in the Clenshaw-Curtis quadrature implementation. --- README.md | 49 +++++++- benches/cheby.rs | 27 +++++ examples/gauss_chebyshev.rs | 5 +- examples/typed_quantities.rs | 8 +- src/approx/adaptive.rs | 179 +++++++++++++++--------------- src/approx/error_estimation.rs | 21 +++- src/approx/fit.rs | 41 +++++-- src/approx/interpolation.rs | 102 +++++++++++++++-- src/approx/minimax.rs | 28 +++-- src/approx/mod.rs | 6 +- src/calculus/integral.rs | 148 +++++++++++++++++++++++- src/core/eval.rs | 6 + src/core/mod.rs | 2 +- src/core/series.rs | 143 +++++++++++++++++++++++- src/piecewise/table.rs | 30 ++++- src/quadrature/clenshaw_curtis.rs | 76 +++++++++++-- tests/calculus_roundtrip.rs | 36 ++++++ tests/core_properties.rs | 68 ++++++++++++ tests/error_branches.rs | 25 +++-- tests/functional_pipeline.rs | 22 ++++ tests/minimax.rs | 16 +-- 21 files changed, 871 insertions(+), 167 deletions(-) diff --git a/README.md b/README.md index cd60903..2f0dd20 100644 --- a/README.md +++ b/README.md @@ -39,6 +39,52 @@ let value = series.evaluate(domain.normalize(1.25)); `Domain` maps `start -> -1`, midpoint `-> 0`, and `end -> 1`. `ChebySeries` stores coefficients of `T_0..T_{N-1}`. +## Mathematical model + +A Chebyshev series approximates a function on a finite interval `[a, b]` by +first mapping the physical coordinate `x` into a dimensionless coordinate +`tau` on `[-1, 1]`: + +```text +mid = (a + b) / 2 +half = (b - a) / 2 +tau = (x - mid) / half +x = mid + half * tau +``` + +The approximating polynomial is stored in the first-kind Chebyshev basis: + +```text +p(x) = sum_j a_j T_j(tau) +T_0(tau) = 1 +T_1(tau) = tau +T_{j+1}(tau) = 2 tau T_j(tau) - T_{j-1}(tau) +``` + +`fit_coeffs` expects samples at the roots of `T_N`: + +```text +tau_k = cos(pi * (2k + 1) / (2N)), k = 0..N-1 +a_0 = (1/N) sum_k f(tau_k) +a_j = (2/N) sum_k f(tau_k) cos(j pi (2k + 1) / (2N)), j > 0 +``` + +Evaluation uses Clenshaw recurrence, which evaluates the same series without +expanding it into monomials. This is usually more stable and faster than +building powers `tau^j` explicitly. + +For a domain-aware series, derivatives and integrals apply the affine scale: + +```text +dp/dx = (1 / half) dp/dtau = (2 / (b - a)) dp/dtau +Integral from a to x = half * Integral from -1 to tau of p(s) ds +``` + +Use this crate when a smooth function on a bounded interval will be evaluated +many times: ephemeris segments, trajectory tables, calibration curves, +spectral collocation, quadrature rules, and compact tabulation of expensive +functions are typical applications. + ## Why Chebyshev? Chebyshev expansions are well conditioned on finite intervals, avoid the worst @@ -72,7 +118,8 @@ sampled nodes and values while preserving value units. ## Approximation `approx::fit` fits coefficients from samples or functions. `adaptive` adds a -builder that increases degree until coefficient-tail tolerances are met. +builder that tries increasing coefficient counts up to the caller's +`max_degree` until coefficient-tail tolerances are met. `minimax` exposes a Remez-style interface with convergence diagnostics. ## Piecewise tables diff --git a/benches/cheby.rs b/benches/cheby.rs index 36eeae0..bb24433 100644 --- a/benches/cheby.rs +++ b/benches/cheby.rs @@ -133,6 +133,30 @@ fn clenshaw_curtis_integrate(c: &mut Criterion) { }); } +fn clenshaw_curtis_cached(c: &mut Criterion) { + let rule = cheby::quadrature::clenshaw_curtis::ClenshawCurtisRule::<32>::new(); + c.bench_function("clenshaw_curtis_cached_n32", |b| { + b.iter(|| rule.integrate::(cheby::Domain::new(0.0, 1.0), |x| black_box(x).exp())) + }); +} + +fn definite_integral_one_shot(c: &mut Criterion) { + let domain = cheby::Domain::new(0.0_f64, core::f64::consts::PI); + let series = cheby::approx::fit::fit_from_fn_on::(domain, f64::sin); + c.bench_function("definite_integral_one_shot_n17", |b| { + b.iter(|| series.evaluate_integral_from_start(black_box(1.5)).unwrap()) + }); +} + +fn definite_integral_cached(c: &mut Criterion) { + let domain = cheby::Domain::new(0.0_f64, core::f64::consts::PI); + let series = cheby::approx::fit::fit_from_fn_on::(domain, f64::sin); + let cached = series.definite_integral_from_start(); + c.bench_function("definite_integral_cached_n17", |b| { + b.iter(|| cached.evaluate(black_box(1.5)).unwrap()) + }); +} + fn gauss_chebyshev(c: &mut Criterion) { c.bench_function("gauss_chebyshev_weighted_n32", |b| { b.iter(|| { @@ -163,6 +187,9 @@ criterion_group!( piecewise_evaluate, adaptive_piecewise_lookup, clenshaw_curtis_integrate, + clenshaw_curtis_cached, + definite_integral_one_shot, + definite_integral_cached, gauss_chebyshev, spectral_diff_matrix ); diff --git a/examples/gauss_chebyshev.rs b/examples/gauss_chebyshev.rs index 7b91ede..f4b2217 100644 --- a/examples/gauss_chebyshev.rs +++ b/examples/gauss_chebyshev.rs @@ -14,5 +14,8 @@ fn main() { let j0_estimate: f64 = gauss_chebyshev::integrate_weighted::<_, N>(f64::cos); println!("∫ cos(x)/√(1−x²) dx ≈ {j0_estimate}"); - println!("(reference: π·J₀(1) ≈ {})", core::f64::consts::PI * 0.7651976865579666); + println!( + "(reference: π·J₀(1) ≈ {})", + core::f64::consts::PI * 0.7651976865579666 + ); } diff --git a/examples/typed_quantities.rs b/examples/typed_quantities.rs index d7f7d91..1c62fad 100644 --- a/examples/typed_quantities.rs +++ b/examples/typed_quantities.rs @@ -6,8 +6,8 @@ //! - `fit_from_fn_t` to sample a function with a typed time argument. //! - `ChebySegment` for evaluation and derivative. //! -//! The derivative result has Rust type `Kilometer` but physically represents -//! `Kilometer / Second` — the caller must know the time-unit context. +//! The derivative result has Rust type `Velocity`, so the +//! time-unit context is preserved by the type system. //! //! Run with: //! `cargo run --example typed_quantities` @@ -39,13 +39,13 @@ fn main() { for raw_t in [0.5, 1.5, 3.0, 5.5] { let t = Second::new(raw_t); let (val, dval_dt) = seg.eval_both(t); - // dval_dt has Rust type Kilometer, but physically is km/s. + // dval_dt has a typed km/s unit. let exact_val = Kilometer::new(7000.0 + 250.0 * raw_t.sin()); let exact_dval = 250.0 * raw_t.cos(); // km/s println!( "t={raw_t:.2} s: val={val:.6} (exact {exact_val:.6}) \ dval/dt={:.6} km/s (exact {exact_dval:.6} km/s)", - dval_dt.value() // extract raw f64: type is Kilometer but units are km/s + dval_dt.value() ); } } diff --git a/src/approx/adaptive.rs b/src/approx/adaptive.rs index ed6ac83..b70771d 100644 --- a/src/approx/adaptive.rs +++ b/src/approx/adaptive.rs @@ -1,6 +1,9 @@ //! Adaptive Chebyshev fitting. -use crate::approx::{error_estimation, fit}; +use alloc::vec; +use alloc::vec::Vec; + +use crate::approx::error_estimation; use crate::core::{ChebyError, ChebyScalar, ChebySeriesDyn, ChebyTime, Domain}; /// Metadata returned by adaptive fitting. @@ -64,24 +67,6 @@ impl AdaptiveFit { } } -macro_rules! build_degree { - ($n:literal, $self:ident) => {{ - let fitted = fit::fit_from_fn_on::($self.domain, &$self.f); - let coeffs = fitted.into_coeffs(); - let err = error_estimation::estimated_truncation_error(&coeffs); - let converged = error_estimation::is_converged(&coeffs, $self.abs_tol, $self.rel_tol); - let series = ChebySeriesDyn::new(alloc::vec::Vec::from(coeffs))?; - return Ok(AdaptiveFitResult { - series, - report: FitReport { - degree: $n, - estimated_error: err, - converged, - }, - }); - }}; -} - impl AdaptiveFit where X: ChebyTime, @@ -95,79 +80,93 @@ where if self.max_degree < 2 { return Err(ChebyError::InvalidDegree); } - for n in [8_usize, 16, 32, 64, 128] { - if n > self.max_degree { - break; - } - match n { - 8 => { - let fitted = fit::fit_from_fn_on::(self.domain, &self.f); - let coeffs = fitted.into_coeffs(); - if error_estimation::is_converged(&coeffs, self.abs_tol, self.rel_tol) - || n == self.max_degree - { - let err = error_estimation::estimated_truncation_error(&coeffs); - return Ok(AdaptiveFitResult { - series: ChebySeriesDyn::new(alloc::vec::Vec::from(coeffs))?, - report: FitReport { - degree: n, - estimated_error: err, - converged: error_estimation::is_converged( - &coeffs, - self.abs_tol, - self.rel_tol, - ), - }, - }); - } - } - 16 => { - let fitted = fit::fit_from_fn_on::(self.domain, &self.f); - let coeffs = fitted.into_coeffs(); - if error_estimation::is_converged(&coeffs, self.abs_tol, self.rel_tol) - || n == self.max_degree - { - let err = error_estimation::estimated_truncation_error(&coeffs); - return Ok(AdaptiveFitResult { - series: ChebySeriesDyn::new(alloc::vec::Vec::from(coeffs))?, - report: FitReport { - degree: n, - estimated_error: err, - converged: error_estimation::is_converged( - &coeffs, - self.abs_tol, - self.rel_tol, - ), - }, - }); - } - } - 32 => { - let fitted = fit::fit_from_fn_on::(self.domain, &self.f); - let coeffs = fitted.into_coeffs(); - if error_estimation::is_converged(&coeffs, self.abs_tol, self.rel_tol) - || n == self.max_degree - { - let err = error_estimation::estimated_truncation_error(&coeffs); - return Ok(AdaptiveFitResult { - series: ChebySeriesDyn::new(alloc::vec::Vec::from(coeffs))?, - report: FitReport { - degree: n, - estimated_error: err, - converged: error_estimation::is_converged( - &coeffs, - self.abs_tol, - self.rel_tol, - ), - }, - }); - } - } - 64 => build_degree!(64, self), - 128 => build_degree!(128, self), - _ => {} + if !(self.abs_tol.is_finite() + && self.abs_tol >= 0.0 + && self.rel_tol.is_finite() + && self.rel_tol >= 0.0) + { + return Err(ChebyError::NonFiniteInput); + } + + let degrees = degree_schedule(self.max_degree); + for n in degrees { + let coeffs = fit_dyn_from_fn_on::(self.domain, n, &self.f); + let err = error_estimation::estimated_truncation_error(&coeffs); + let converged = error_estimation::is_converged(&coeffs, self.abs_tol, self.rel_tol); + if converged || n == self.max_degree { + return Ok(AdaptiveFitResult { + series: ChebySeriesDyn::new(coeffs)?, + report: FitReport { + degree: n, + estimated_error: err, + converged, + }, + }); } } Err(ChebyError::InvalidDegree) } } + +fn degree_schedule(max_degree: usize) -> Vec { + if max_degree < 8 { + return vec![max_degree]; + } + + let mut degrees = Vec::new(); + let mut n = 8_usize; + while n < max_degree { + degrees.push(n); + match n.checked_mul(2) { + Some(next) => n = next, + None => break, + } + } + if degrees.last().copied() != Some(max_degree) { + degrees.push(max_degree); + } + degrees +} + +fn fit_dyn_from_fn_on(domain: Domain, n: usize, f: &F) -> Vec +where + T: ChebyScalar, + X: ChebyTime, + F: Fn(X) -> T, +{ + let nf = n as f64; + // Sample at roots of T_n; cache the normalized node `x_k = cos(theta_k)` + // so we can build T_j(x_k) by Chebyshev recurrence in `j` instead of + // calling `cos` n*n times. + let mut xs = vec![0.0_f64; n]; + let mut values: Vec = Vec::with_capacity(n); + for (k, x_slot) in xs.iter_mut().enumerate() { + let x = (core::f64::consts::PI * (2.0 * k as f64 + 1.0) / (2.0 * nf)).cos(); + *x_slot = x; + values.push(f(domain.denormalize(x))); + } + + let mut coeffs = vec![T::zero(); n]; + for (k, &v) in values.iter().enumerate() { + let x = xs[k]; + coeffs[0] = coeffs[0] + v; + if n > 1 { + coeffs[1] = coeffs[1] + v * x; + let two_x = 2.0 * x; + let mut tkm1 = 1.0_f64; + let mut tk = x; + for coeff in coeffs.iter_mut().take(n).skip(2) { + let tkp1 = two_x * tk - tkm1; + *coeff = *coeff + v * tkp1; + tkm1 = tk; + tk = tkp1; + } + } + } + coeffs[0] = coeffs[0] / nf; + let scale = 2.0 / nf; + for c in coeffs.iter_mut().skip(1) { + *c = *c * scale; + } + coeffs +} diff --git a/src/approx/error_estimation.rs b/src/approx/error_estimation.rs index 07ee1a5..16e01ca 100644 --- a/src/approx/error_estimation.rs +++ b/src/approx/error_estimation.rs @@ -9,6 +9,11 @@ pub fn tail_norm(coeffs: &[T], from: usize) -> f64 { } /// Estimate truncation error from the final quarter of coefficients. +/// +/// For smooth analytic targets the magnitude of the dropped tail bounds the +/// truncation error well. For localized features or endpoint stress this is +/// a heuristic, not a guaranteed bound; see [`is_converged`] for a stricter +/// test that also checks the trailing few coefficients individually. #[inline] pub fn estimated_truncation_error(coeffs: &[T]) -> f64 { if coeffs.is_empty() { @@ -18,7 +23,11 @@ pub fn estimated_truncation_error(coeffs: &[T]) -> f64 { } } -/// Whether coefficients appear converged by absolute and relative tail tests. +/// Whether coefficients appear converged. +/// +/// Combines the quarter-tail estimate with an explicit check on the last +/// few coefficient magnitudes so that a single oscillating high-order entry +/// cannot be hidden by averaging across the tail. #[inline] pub fn is_converged(coeffs: &[T], abs_tol: f64, rel_tol: f64) -> bool { let scale = coeffs @@ -26,6 +35,14 @@ pub fn is_converged(coeffs: &[T], abs_tol: f64, rel_tol: f64) -> .map(|c| c.magnitude()) .fold(0.0_f64, f64::max) .max(1.0); - let err = estimated_truncation_error(coeffs); + let tail_err = estimated_truncation_error(coeffs); + // Largest magnitude among the last min(3, len) coefficients. + let trailing = coeffs + .iter() + .rev() + .take(3) + .map(|c| c.magnitude()) + .fold(0.0_f64, f64::max); + let err = tail_err.max(trailing); err <= abs_tol || err <= rel_tol * scale } diff --git a/src/approx/fit.rs b/src/approx/fit.rs index 4f54d9a..9370152 100644 --- a/src/approx/fit.rs +++ b/src/approx/fit.rs @@ -1,19 +1,46 @@ //! Coefficient fitting from samples and functions. +//! +//! The fixed-size fitter samples at roots of `T_N` and applies the direct +//! discrete-cosine coefficient formula. With `T_j(x_k) = cos(j θ_k)` evaluated +//! by Chebyshev recurrence in `j`, the work is `O(N²)` multiplies-and-adds +//! and only `O(N)` calls to `cos` (used to build the node grid once). use crate::core::{nodes, ChebyScalar, ChebySeries, ChebySeriesOn, Domain, NodeKind}; /// Compute Chebyshev coefficients from values sampled at roots of `T_N`. +/// +/// Uses the Chebyshev three-term recurrence to generate `T_j(x_k)` for each +/// sample `k` in a single pass over `j`, avoiding `N²` `cos()` calls. #[inline] pub fn fit_coeffs(values: &[T; N]) -> [T; N] { let mut coeffs = [T::zero(); N]; - let n = N as f64; - for (j, coeff) in coeffs.iter_mut().enumerate() { - let mut sum = T::zero(); - for (k, value) in values.iter().enumerate() { - let arg = core::f64::consts::PI * j as f64 * (2.0 * k as f64 + 1.0) / (2.0 * n); - sum = sum + *value * arg.cos(); + if N == 0 { + return coeffs; + } + let xs = nodes::(NodeKind::Roots); + for (k, &v) in values.iter().enumerate() { + let x = xs[k]; + // T_0 = 1 + coeffs[0] = coeffs[0] + v; + if N > 1 { + // T_1 = x + coeffs[1] = coeffs[1] + v * x; + let two_x = 2.0 * x; + let mut tkm1 = 1.0_f64; + let mut tk = x; + for coeff in coeffs.iter_mut().take(N).skip(2) { + let tkp1 = two_x * tk - tkm1; + *coeff = *coeff + v * tkp1; + tkm1 = tk; + tk = tkp1; + } } - *coeff = if j == 0 { sum / n } else { sum * (2.0 / n) }; + } + let nf = N as f64; + coeffs[0] = coeffs[0] / nf; + let scale = 2.0 / nf; + for coeff in coeffs.iter_mut().skip(1) { + *coeff = *coeff * scale; } coeffs } diff --git a/src/approx/interpolation.rs b/src/approx/interpolation.rs index b177d3e..359502f 100644 --- a/src/approx/interpolation.rs +++ b/src/approx/interpolation.rs @@ -1,8 +1,11 @@ //! Barycentric interpolation. -use crate::core::{ChebyError, ChebyScalar, ChebyTime}; +use crate::core::{nodes, ChebyError, ChebyScalar, ChebyTime, Domain, NodeKind}; /// Barycentric interpolator over supplied nodes and values. +/// +/// For `N == 1`, the interpolant is the constant `values[0]` and `scale` is +/// unused; for `N > 1` `scale = nodes[0] - nodes[N-1]` and must be non-zero. #[derive(Debug, Clone, PartialEq)] pub struct BarycentricInterpolator { nodes: [X; N], @@ -17,16 +20,25 @@ where X: ChebyTime, { /// Create an interpolator. Nodes must be distinct. + /// + /// Generic O(N²) weight construction. For Chebyshev roots or Lobatto + /// nodes prefer [`Self::on_chebyshev_roots`] or [`Self::on_lobatto_nodes`], + /// which use the closed-form weights and are both faster and numerically + /// more stable. pub fn new(nodes: [X; N], values: [T; N]) -> Result { if N == 0 { return Err(ChebyError::InvalidDegree); } - let scale = if N > 1 { - nodes[0] - nodes[N - 1] - } else { - X::zero() + nodes[0] * 0.0 + nodes[0] * 1.0 - }; - if N > 1 && scale == X::zero() { + if N == 1 { + return Ok(Self { + nodes, + values, + weights: [1.0; N], + scale: X::zero(), + }); + } + let scale = nodes[0] - nodes[N - 1]; + if scale == X::zero() { return Err(ChebyError::InvalidDomain); } let mut weights = [0.0; N]; @@ -51,8 +63,84 @@ where }) } + /// Build a barycentric interpolant on Chebyshev roots of `T_N` mapped + /// into `domain`, using the closed-form weights + /// `w_k = (-1)^k sin((2k+1)π / (2N))`. + /// + /// Setup is `O(N)` and avoids the `O(N²)` weight product chain. The + /// scaling cancels in the barycentric formula, so absolute magnitudes do + /// not matter. + pub fn on_chebyshev_roots(domain: Domain, values: [T; N]) -> Result { + if N == 0 { + return Err(ChebyError::InvalidDegree); + } + let xs = nodes::(NodeKind::Roots); + let mapped: [X; N] = core::array::from_fn(|k| domain.denormalize(xs[k])); + if N == 1 { + return Ok(Self { + nodes: mapped, + values, + weights: [1.0; N], + scale: X::zero(), + }); + } + let mut weights = [0.0; N]; + let mut sign = 1.0_f64; + for (k, w) in weights.iter_mut().enumerate() { + let theta = core::f64::consts::PI * (2.0 * k as f64 + 1.0) / (2.0 * N as f64); + *w = sign * theta.sin(); + sign = -sign; + } + let scale = mapped[0] - mapped[N - 1]; + Ok(Self { + nodes: mapped, + values, + weights, + scale, + }) + } + + /// Build a barycentric interpolant on Chebyshev-Lobatto nodes of order + /// `N - 1` mapped into `domain`, using the closed-form weights + /// `w_0 = 1/2, w_{N-1} = (-1)^{N-1}/2`, interior `w_k = (-1)^k`. + /// + /// Setup is `O(N)`. + pub fn on_lobatto_nodes(domain: Domain, values: [T; N]) -> Result { + if N == 0 { + return Err(ChebyError::InvalidDegree); + } + let xs = nodes::(NodeKind::Lobatto); + let mapped: [X; N] = core::array::from_fn(|k| domain.denormalize(xs[k])); + if N == 1 { + return Ok(Self { + nodes: mapped, + values, + weights: [1.0; N], + scale: X::zero(), + }); + } + let mut weights = [0.0; N]; + let mut sign = 1.0_f64; + for w in weights.iter_mut() { + *w = sign; + sign = -sign; + } + weights[0] *= 0.5; + weights[N - 1] *= 0.5; + let scale = mapped[0] - mapped[N - 1]; + Ok(Self { + nodes: mapped, + values, + weights, + scale, + }) + } + /// Evaluate the interpolant at `x`. pub fn evaluate(&self, x: X) -> Result { + if N == 1 { + return Ok(self.values[0]); + } let mut numerator = T::zero(); let mut denominator = 0.0; for k in 0..N { diff --git a/src/approx/minimax.rs b/src/approx/minimax.rs index 3f808e4..8a141fc 100644 --- a/src/approx/minimax.rs +++ b/src/approx/minimax.rs @@ -36,9 +36,7 @@ use alloc::vec; use alloc::vec::Vec; -use crate::core::{ - basis, ChebyError, ChebySeries, ChebySeriesDyn, ChebySeriesOn, ChebyTime, Domain, -}; +use crate::core::{ChebyError, ChebySeries, ChebySeriesDyn, ChebySeriesOn, ChebyTime, Domain}; /// Options for the Remez exchange algorithm. #[derive(Debug, Clone, Copy, PartialEq)] @@ -159,6 +157,8 @@ where let n = degree + 1; let mut matrix = vec![0.0; m * (m + 1)]; let mut coeffs = vec![0.0; n]; + let mut residual = vec![0.0_f64; grid]; + let mut blocks: Vec<(f64, f64)> = Vec::with_capacity(grid / 2 + 2); let mut prev_max_err = f64::INFINITY; let mut max_err = 0.0; let mut leveled = 0.0; @@ -169,10 +169,23 @@ where iterations = iter; // Build linear system: row i = [T_0(tau_i) ... T_N(tau_i) | (-1)^i | f(x_i)] + // T_k(tau) is filled by the Chebyshev recurrence, one row at a time. for (i, &tau) in taus.iter().enumerate() { let row = i * (m + 1); - for k in 0..n { - matrix[row + k] = basis::t(k, tau); + if n >= 1 { + matrix[row] = 1.0; + } + if n >= 2 { + matrix[row + 1] = tau; + let two_tau = 2.0 * tau; + let mut tkm1 = 1.0; + let mut tk = tau; + for k in 2..n { + let tkp1 = two_tau * tk - tkm1; + matrix[row + k] = tkp1; + tkm1 = tk; + tk = tkp1; + } } matrix[row + n] = if i % 2 == 0 { 1.0 } else { -1.0 }; matrix[row + m] = f(domain.denormalize(tau)); @@ -187,8 +200,7 @@ where } leveled = matrix[n * (m + 1) + m].abs(); - // Sample residual on dense grid. - let mut residual = vec![0.0_f64; grid]; + // Sample residual on dense grid into the reusable buffer. for (i, r) in residual.iter_mut().enumerate() { let tau = -1.0 + 2.0 * i as f64 / (grid - 1) as f64; let x = domain.denormalize(tau); @@ -196,7 +208,7 @@ where } // Locate per-sign-block extrema. - let mut blocks: Vec<(f64, f64)> = Vec::new(); // (tau_max, residual_max_signed) + blocks.clear(); let mut cur_sign = 0i8; let mut cur_best_idx = 0usize; let mut cur_best_abs = -1.0_f64; diff --git a/src/approx/mod.rs b/src/approx/mod.rs index 15f9646..343c995 100644 --- a/src/approx/mod.rs +++ b/src/approx/mod.rs @@ -3,11 +3,11 @@ //! # Theory //! //! Given a function `f` on a [`crate::core::Domain`], `fit` produces a -//! truncated Chebyshev series whose coefficients minimise the L² error -//! at the Chebyshev nodes (interpolation in the polynomial sense). +//! Chebyshev interpolating series by sampling `f` at roots of `T_N` and using +//! the discrete cosine orthogonality formula for the coefficients. //! The optional [`adaptive`] module refines the degree until a target //! tolerance is met, and the optional [`minimax`] module replaces the -//! L² fit with a real Remez exchange that minimises the L∞ error. +//! root-sample fit with a Remez exchange that minimises the L∞ error. //! //! # Features //! diff --git a/src/calculus/integral.rs b/src/calculus/integral.rs index b6e5b1a..b0d324c 100644 --- a/src/calculus/integral.rs +++ b/src/calculus/integral.rs @@ -1,22 +1,158 @@ //! Integral helpers. -use crate::core::{ChebyScalar, ChebySeriesOn, ChebyTime}; +use crate::core::{ChebyError, ChebyScalar, ChebySeries, ChebySeriesOn, ChebyTime, Domain}; pub use crate::core::IntegrateWith; +/// A precomputed antiderivative shifted so that its value at the domain +/// start is zero. +/// +/// Build one with [`ChebySeriesOn::definite_integral_from_start`] when you +/// will repeatedly query the definite integral `∫_{start}^{x} f(t) dt` for +/// many `x`. Each evaluation is then a single Clenshaw pass over the +/// `(N + 1)`-length antiderivative; no per-call coefficient construction +/// or normalization is performed. +#[derive(Debug, Clone, PartialEq)] +pub struct DefiniteIntegralFromStart +where + T: ChebyScalar, + X: ChebyTime, +{ + domain: Domain, + /// Antiderivative coefficients, length `N + 1`, with `b_0` chosen so + /// `A(-1) = 0` (i.e. evaluation at the domain start vanishes). + coeffs: [T; N], + coeffs_tail: T, +} + +impl DefiniteIntegralFromStart +where + T: ChebyScalar + IntegrateWith, + X: ChebyTime, +{ + /// Domain of validity. + #[inline] + pub fn domain(&self) -> Domain { + self.domain + } + + /// Evaluate `∫_{start}^{x} f(t) dt`. + pub fn evaluate(&self, x: X) -> Result<>::Integral, ChebyError> { + if !self.domain.contains(x) { + return Err(ChebyError::EvaluationOutOfDomain); + } + let tau = self.domain.normalize(x); + // Clenshaw with the `N + 1`-coefficient antiderivative `coeffs ++ [coeffs_tail]`. + let normalized = clenshaw_with_tail::(&self.coeffs, self.coeffs_tail, tau); + Ok(normalized.scale_integral(self.domain.half_width())) + } +} + +#[inline] +fn clenshaw_with_tail(coeffs: &[T; N], tail: T, x: f64) -> T { + // Evaluate the `(N+1)`-coefficient series `[coeffs[0..N], tail]` via Clenshaw, + // without allocating. + let two_x = 2.0 * x; + let mut b1 = T::zero(); + let mut b2 = T::zero(); + // Process the tail first (highest index = N). + let mut b0 = b1 * two_x - b2 + tail; + b2 = b1; + b1 = b0; + for k in (1..N).rev() { + b0 = b1 * two_x - b2 + coeffs[k]; + b2 = b1; + b1 = b0; + } + if N == 0 { + // Series is just `[tail]`; b1 currently holds `tail`, b2 = 0. + return b1; + } + coeffs[0] + b1 * x - b2 +} + impl ChebySeriesOn where T: ChebyScalar + IntegrateWith, X: ChebyTime, { /// Evaluate the definite integral from the domain start to `x`. + /// + /// Convenience for one-shot integration. For repeated queries prefer + /// [`Self::definite_integral_from_start`], which precomputes the shifted + /// antiderivative once. pub fn evaluate_integral_from_start( &self, x: X, - ) -> Result<>::Integral, crate::core::ChebyError> { - let antiderivative = self.series().integral(T::zero()); - let tau = self.domain().normalize(x); - let normalized = antiderivative.evaluate(tau) - antiderivative.evaluate(-1.0); - Ok(normalized.scale_integral(self.domain().half_width())) + ) -> Result<>::Integral, ChebyError> { + self.definite_integral_from_start().evaluate(x) + } + + /// Build a reusable precomputed antiderivative for repeated definite + /// integral queries from the domain start. + pub fn definite_integral_from_start(&self) -> DefiniteIntegralFromStart { + let domain = self.domain(); + let coeffs_in = self.series().coeffs(); + // Build the full non-truncated antiderivative directly into a + // length-(N+1) representation: store the first `N` entries inline + // and the highest entry separately so we can stay const-generic + // without `generic_const_exprs`. + let mut out = [T::zero(); N]; + let mut tail = T::zero(); + if N >= 1 { + // ∫T_0 = T_1: contributes c_0 to b_1. + if N >= 2 { + out[1] = out[1] + coeffs_in[0]; + } else { + // Only one input coefficient; b_1 lives in the tail slot. + tail = tail + coeffs_in[0]; + } + } + if N >= 2 { + // ∫T_1 = T_2/4: contributes c_1/4 to b_2. + if N >= 3 { + out[2] = out[2] + coeffs_in[1] / 4.0; + } else { + tail = tail + coeffs_in[1] / 4.0; + } + } + for k in 2..N { + let plus = coeffs_in[k] / (2.0 * (k + 1) as f64); + let minus = coeffs_in[k] / (2.0 * (k - 1) as f64); + // b_{k+1} += plus (k+1 may equal N => goes in tail) + if k + 1 < N { + out[k + 1] = out[k + 1] + plus; + } else { + tail = tail + plus; + } + // b_{k-1} -= minus (k-1 < N always) + out[k - 1] = out[k - 1] - minus; + } + // Choose b_0 so A(-1) = 0. With T_k(-1) = (-1)^k: + // A(-1) = sum_{k=0}^{N} b_k (-1)^k. + let mut shift = T::zero(); + let mut sign = 1.0_f64; + for &b in out.iter() { + shift = shift + b * sign; + sign = -sign; + } + // Tail contribution at index N: sign is (-1)^N. + let tail_sign = if N.is_multiple_of(2) { 1.0 } else { -1.0 }; + shift = shift + tail * tail_sign; + if N >= 1 { + out[0] = out[0] - shift; + } else { + tail = tail - shift; + } + DefiniteIntegralFromStart { + domain, + coeffs: out, + coeffs_tail: tail, + } } } + +// Suppress an unused-import warning when the impl above is the only consumer +// of `ChebySeries`. +#[allow(dead_code)] +fn _force_chebyseries_use(_: &ChebySeries) {} diff --git a/src/core/eval.rs b/src/core/eval.rs index caff48c..039d771 100644 --- a/src/core/eval.rs +++ b/src/core/eval.rs @@ -1,4 +1,10 @@ //! Clenshaw evaluation for Chebyshev series. +//! +//! For coefficients `a_k`, the represented polynomial is +//! `sum_k a_k T_k(x)`. Clenshaw evaluates that basis recurrence directly: +//! it walks the coefficients backward with +//! `b_k = 2x b_{k+1} - b_{k+2} + a_k`, then returns +//! `a_0 + x b_1 - b_2`. This avoids converting the series to powers of `x`. use super::ChebyScalar; diff --git a/src/core/mod.rs b/src/core/mod.rs index 698d031..3223efd 100644 --- a/src/core/mod.rs +++ b/src/core/mod.rs @@ -36,4 +36,4 @@ pub use scalar::{ChebyScalar, ChebyTime, DifferentiateWith, IntegrateWith}; pub use series::{ChebySeries, ChebySeriesOn}; #[cfg(feature = "alloc")] -pub use series::ChebySeriesDyn; +pub use series::{ChebySeriesDyn, ChebySeriesDynOn}; diff --git a/src/core/series.rs b/src/core/series.rs index 47071f7..478057d 100644 --- a/src/core/series.rs +++ b/src/core/series.rs @@ -92,7 +92,9 @@ impl ChebySeries { /// Return a normalized-coordinate indefinite integral series. /// - /// The array length is preserved; the highest possible term is truncated. + /// The array length is preserved; the single highest antiderivative term + /// (`b_N`) is truncated. For a non-truncating dynamic version see + /// [`ChebySeriesDyn::integral`]. pub fn integral(&self, constant: T) -> Self { let mut out = [T::zero(); N]; if N == 0 { @@ -105,8 +107,12 @@ impl ChebySeries { if N > 2 { out[2] = out[2] + self.coeffs[1] / 4.0; } - for k in 2..(N.saturating_sub(1)) { - out[k + 1] = out[k + 1] + self.coeffs[k] / (2.0 * (k + 1) as f64); + // For k >= 2: ∫T_k = (T_{k+1}/(k+1) - T_{k-1}/(k-1)) / 2. + // Apply both contributions; only `b_{k+1}` for k = N-1 is dropped. + for k in 2..N { + if k + 1 < N { + out[k + 1] = out[k + 1] + self.coeffs[k] / (2.0 * (k + 1) as f64); + } out[k - 1] = out[k - 1] - self.coeffs[k] / (2.0 * (k - 1) as f64); } Self::new(out) @@ -197,6 +203,19 @@ impl ChebySeriesDyn { &self.coeffs } + /// Number of coefficients (degree + 1). + #[inline] + pub fn len(&self) -> usize { + self.coeffs.len() + } + + /// Whether the series has no coefficients (always `false` once + /// constructed via [`Self::new`]; included for API completeness). + #[inline] + pub fn is_empty(&self) -> bool { + self.coeffs.is_empty() + } + /// Consume into coefficients. #[inline] pub fn into_coeffs(self) -> Vec { @@ -208,4 +227,122 @@ impl ChebySeriesDyn { pub fn evaluate(&self, x: f64) -> T { eval::evaluate(&self.coeffs, x) } + + /// Evaluate value and normalized derivative `df/dtau` together. + #[inline] + pub fn evaluate_both(&self, x: f64) -> (T, T) { + eval::evaluate_both(&self.coeffs, x) + } + + /// Evaluate the normalized-coordinate derivative `df/dtau`. + #[inline] + pub fn evaluate_derivative(&self, x: f64) -> T { + eval::evaluate_derivative(&self.coeffs, x) + } + + /// Return a normalized-coordinate derivative series. + /// + /// The output has `max(len() - 1, 1)` coefficients. + pub fn derivative(&self) -> Self { + let n = self.coeffs.len(); + if n <= 1 { + return Self { + coeffs: alloc::vec![T::zero()], + }; + } + let mut out = alloc::vec![T::zero(); n - 1]; + out[n - 2] = self.coeffs[n - 1] * (2.0 * (n - 1) as f64); + if n > 2 { + for k in (1..=(n - 2)).rev() { + let prev = if k + 1 < out.len() { + out[k + 1] + } else { + T::zero() + }; + out[k - 1] = prev + self.coeffs[k] * (2.0 * k as f64); + } + } + out[0] = out[0] / 2.0; + Self { coeffs: out } + } + + /// Return a non-truncating normalized-coordinate indefinite integral + /// series. Output length is `len() + 1`. + pub fn integral(&self, constant: T) -> Self { + let n = self.coeffs.len(); + let mut out = alloc::vec![T::zero(); n + 1]; + out[0] = constant; + if n >= 1 { + out[1] = out[1] + self.coeffs[0]; + } + if n >= 2 { + out[2] = out[2] + self.coeffs[1] / 4.0; + } + for k in 2..n { + out[k + 1] = out[k + 1] + self.coeffs[k] / (2.0 * (k + 1) as f64); + out[k - 1] = out[k - 1] - self.coeffs[k] / (2.0 * (k - 1) as f64); + } + Self { coeffs: out } + } +} + +/// A dynamic Chebyshev series tied to a physical domain. +#[cfg(feature = "alloc")] +#[derive(Debug, Clone, PartialEq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub struct ChebySeriesDynOn { + domain: Domain, + series: ChebySeriesDyn, +} + +#[cfg(feature = "alloc")] +impl ChebySeriesDynOn { + /// Create a domain-aware dynamic series. + #[inline] + pub fn new(domain: Domain, series: ChebySeriesDyn) -> Self { + Self { domain, series } + } + + /// Domain of validity. + #[inline] + pub fn domain(&self) -> Domain { + self.domain + } + + /// Underlying normalized-coordinate series. + #[inline] + pub fn series(&self) -> &ChebySeriesDyn { + &self.series + } + + /// Consume and return the underlying series. + #[inline] + pub fn into_series(self) -> ChebySeriesDyn { + self.series + } + + /// Evaluate inside the domain. + #[inline] + pub fn evaluate(&self, x: X) -> Result { + if !self.domain.contains(x) { + return Err(ChebyError::EvaluationOutOfDomain); + } + Ok(self.series.evaluate(self.domain.normalize(x))) + } + + /// Evaluate the dimensionally-correct physical derivative at `x`. + #[inline] + pub fn evaluate_derivative( + &self, + x: X, + ) -> Result<>::Derivative, ChebyError> + where + T: DifferentiateWith, + { + if !self.domain.contains(x) { + return Err(ChebyError::EvaluationOutOfDomain); + } + let d_tau = self.series.evaluate_both(self.domain.normalize(x)).1; + Ok((d_tau * 2.0).scale_derivative(self.domain.end() - self.domain.start())) + } } diff --git a/src/piecewise/table.rs b/src/piecewise/table.rs index 75cf377..4e3bd20 100644 --- a/src/piecewise/table.rs +++ b/src/piecewise/table.rs @@ -7,6 +7,10 @@ use crate::core::{ChebyError, ChebyScalar, ChebySeries, ChebyTime, Differentiate use crate::piecewise::{lookup, ChebySegment}; /// A uniform table of Chebyshev segments. +/// +/// The first segments have [`Self::segment_len`] width; the final segment may +/// be shorter when the requested table end is not an exact multiple of that +/// width. Table lookup treats the total range as half-open: `[start, end)`. #[derive(Debug, Clone, PartialEq)] #[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] pub struct ChebySegmentTable { @@ -30,9 +34,22 @@ where if segment_len <= X::zero() { return Err(ChebyError::NonPositiveSegmentLength); } + let last_idx = segments.len() - 1; for (i, segment) in segments.iter().enumerate() { let expected = start + segment_len * i as f64; - if segment.domain().start() != expected { + let domain = segment.domain(); + if domain.start() != expected { + return Err(ChebyError::InvalidDomain); + } + let width = domain.end() - domain.start(); + if width <= X::zero() { + return Err(ChebyError::NonPositiveSegmentLength); + } + if i == last_idx { + if width > segment_len { + return Err(ChebyError::InvalidDomain); + } + } else if width != segment_len { return Err(ChebyError::InvalidDomain); } } @@ -65,6 +82,9 @@ where end: X, segment_len: X, ) -> Result { + if !start.is_finite() || !end.is_finite() || !segment_len.is_finite() { + return Err(ChebyError::NonFiniteInput); + } if segment_len <= X::zero() { return Err(ChebyError::NonPositiveSegmentLength); } @@ -113,7 +133,10 @@ where /// Table end. #[inline] pub fn end(&self) -> X { - self.start + self.segment_len * self.segments.len() as f64 + self.segments + .last() + .map(|segment| segment.domain().end()) + .unwrap_or(self.start) } /// Uniform segment length used for O(1) lookup. @@ -131,6 +154,9 @@ where /// Look up a segment. #[inline] pub fn get_segment(&self, x: X) -> Option<&ChebySegment> { + if self.segments.is_empty() || x < self.start || x >= self.end() { + return None; + } lookup::uniform_index(self.start, self.segment_len, self.segments.len(), x) .and_then(|idx| self.segments.get(idx)) .filter(|segment| segment.contains(x)) diff --git a/src/quadrature/clenshaw_curtis.rs b/src/quadrature/clenshaw_curtis.rs index b88aac1..f1bb835 100644 --- a/src/quadrature/clenshaw_curtis.rs +++ b/src/quadrature/clenshaw_curtis.rs @@ -6,9 +6,70 @@ //! //! The implementation follows Trefethen, *Spectral Methods in MATLAB*, //! Program 30 (`clencurt`). Cost is `O(N^2)`. +//! +//! For repeated quadrature on the same `N`, build a [`ClenshawCurtisRule`] +//! once and reuse it; the per-call [`integrate`] convenience rebuilds nodes +//! and weights every time. use crate::core::{ChebyScalar, ChebyTime, Domain, IntegrateWith, NodeKind}; +/// A precomputed `N`-point Clenshaw-Curtis rule on `[-1, 1]`. +/// +/// Holds the Chebyshev-Lobatto nodes and the matching weights so that +/// repeated calls do not pay the `O(N²)` weight construction every time. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct ClenshawCurtisRule { + nodes: [f64; N], + weights: [f64; N], +} + +impl Default for ClenshawCurtisRule { + fn default() -> Self { + Self::new() + } +} + +impl ClenshawCurtisRule { + /// Build the cached rule. + #[inline] + pub fn new() -> Self { + Self { + nodes: crate::core::nodes::nodes::(NodeKind::GaussLobatto), + weights: clenshaw_curtis_weights::(), + } + } + + /// Borrow the cached Lobatto nodes on `[-1, 1]`. + #[inline] + pub fn nodes(&self) -> &[f64; N] { + &self.nodes + } + + /// Borrow the cached quadrature weights. + #[inline] + pub fn weights(&self) -> &[f64; N] { + &self.weights + } + + /// Integrate `f` over `domain` using this cached rule. + #[inline] + pub fn integrate( + &self, + domain: Domain, + f: impl Fn(X) -> T, + ) -> >::Integral + where + T: ChebyScalar + IntegrateWith, + X: ChebyTime, + { + let mut sum = T::zero(); + for k in 0..N { + sum = sum + f(domain.denormalize(self.nodes[k])) * self.weights[k]; + } + sum.scale_integral(domain.half_width()) + } +} + /// Clenshaw-Curtis weights for `N` Lobatto nodes on `[-1, 1]`. pub fn clenshaw_curtis_weights() -> [f64; N] { let mut weights = [0.0; N]; @@ -26,8 +87,7 @@ pub fn clenshaw_curtis_weights() -> [f64; N] { } let m = N - 1; // Trefethen's "N" - let theta: [f64; N] = - core::array::from_fn(|k| core::f64::consts::PI * k as f64 / m as f64); + let theta: [f64; N] = core::array::from_fn(|k| core::f64::consts::PI * k as f64 / m as f64); // Endpoint weights. if m.is_multiple_of(2) { @@ -56,6 +116,10 @@ pub fn clenshaw_curtis_weights() -> [f64; N] { } /// Integrate a function over `domain` using `N`-point Clenshaw-Curtis. +/// +/// Convenience wrapper that builds a fresh [`ClenshawCurtisRule`] each call. +/// For repeated quadrature on the same `N`, construct the rule once and call +/// [`ClenshawCurtisRule::integrate`]. pub fn integrate( domain: Domain, f: impl Fn(X) -> T, @@ -64,11 +128,5 @@ where T: ChebyScalar + IntegrateWith, X: ChebyTime, { - let nodes = crate::core::nodes::nodes::(NodeKind::GaussLobatto); - let weights = clenshaw_curtis_weights::(); - let mut sum = T::zero(); - for k in 0..N { - sum = sum + f(domain.denormalize(nodes[k])) * weights[k]; - } - sum.scale_integral(domain.half_width()) + ClenshawCurtisRule::::new().integrate(domain, f) } diff --git a/tests/calculus_roundtrip.rs b/tests/calculus_roundtrip.rs index e8cc47f..c7f8ab6 100644 --- a/tests/calculus_roundtrip.rs +++ b/tests/calculus_roundtrip.rs @@ -30,6 +30,17 @@ fn definite_integral_matches_analytic() { } } +#[test] +fn definite_integral_rejects_out_of_domain_point() { + let domain = Domain::new(0.0, 1.0); + let series = ChebySeries::new([1.0_f64, 0.0]); + let series = cheby::ChebySeriesOn::new(domain, series); + assert_eq!( + series.evaluate_integral_from_start(1.5).unwrap_err(), + cheby::ChebyError::EvaluationOutOfDomain + ); +} + #[test] fn typed_integral_velocity_to_position() { use qtty::{Kilometer, Second}; @@ -46,3 +57,28 @@ fn typed_integral_velocity_to_position() { .unwrap(); assert_abs_diff_eq!(pos.value(), 6.0, epsilon = 1e-12); } + +#[test] +fn cached_definite_integral_matches_one_shot() { + let domain = Domain::new(0.0, core::f64::consts::PI); + let series = fit_from_fn_on::(domain, f64::sin); + let cached = series.definite_integral_from_start(); + for &x in &[0.0_f64, 0.3, 0.7, 1.5, 2.4, 3.0, core::f64::consts::PI] { + let one_shot = series.evaluate_integral_from_start(x).unwrap(); + let from_cache = cached.evaluate(x).unwrap(); + assert_abs_diff_eq!(from_cache, one_shot, epsilon = 1e-14); + let exact = 1.0 - x.cos(); + assert_abs_diff_eq!(from_cache, exact, epsilon = 1e-9); + } +} + +#[test] +fn cached_definite_integral_rejects_out_of_domain() { + let domain = Domain::new(0.0, 1.0); + let series = cheby::ChebySeriesOn::new(domain, ChebySeries::new([1.0_f64, 0.0])); + let cached = series.definite_integral_from_start(); + assert_eq!( + cached.evaluate(1.5).unwrap_err(), + cheby::ChebyError::EvaluationOutOfDomain + ); +} diff --git a/tests/core_properties.rs b/tests/core_properties.rs index 867b8eb..5926b80 100644 --- a/tests/core_properties.rs +++ b/tests/core_properties.rs @@ -32,6 +32,48 @@ fn interpolation_reproduces_samples() { } } +#[test] +fn interpolation_single_node_is_constant_even_at_zero() { + // Regression: a single node at zero used to make `scale == 0`, producing + // a NaN division during evaluation instead of the constant interpolant. + let interp = BarycentricInterpolator::new([0.0_f64], [42.0_f64]).unwrap(); + for &x in &[-2.5_f64, -0.1, 0.0, 0.7, 17.0] { + assert_abs_diff_eq!(interp.evaluate(x).unwrap(), 42.0, epsilon = 0.0); + } +} + +#[test] +fn interpolation_chebyshev_roots_constructor_matches_generic() { + use cheby::core::{nodes_mapped, Domain, NodeKind}; + const N: usize = 12; + let domain = Domain::new(-1.0_f64, 3.0); + let xs = nodes_mapped::(domain, NodeKind::Roots); + let ys: [f64; N] = std::array::from_fn(|k| (xs[k] + 0.3).cos()); + let generic = BarycentricInterpolator::new(xs, ys).unwrap(); + let fast = BarycentricInterpolator::on_chebyshev_roots(domain, ys).unwrap(); + for tau in [-0.95_f64, -0.4, 0.0, 0.2, 0.85] { + let x = domain.denormalize(tau); + assert_abs_diff_eq!( + fast.evaluate(x).unwrap(), + generic.evaluate(x).unwrap(), + epsilon = 1e-12 + ); + } +} + +#[test] +fn interpolation_lobatto_constructor_reproduces_samples() { + use cheby::core::{nodes_mapped, Domain, NodeKind}; + const N: usize = 9; + let domain = Domain::new(0.0_f64, 2.0); + let xs = nodes_mapped::(domain, NodeKind::Lobatto); + let ys: [f64; N] = std::array::from_fn(|k| xs[k].sin()); + let interp = BarycentricInterpolator::on_lobatto_nodes(domain, ys).unwrap(); + for k in 0..N { + assert_abs_diff_eq!(interp.evaluate(xs[k]).unwrap(), ys[k], epsilon = 1e-12); + } +} + #[test] fn derivative_of_integral_recovers_series_approximately() { let series = ChebySeries::new([1.0, -0.5, 0.25, 0.0, 0.0]); @@ -65,3 +107,29 @@ fn adaptive_fitting_converges_for_smooth_function() { .unwrap(); assert!(result.report.converged); } + +#[test] +#[cfg(feature = "adaptive")] +fn adaptive_fitting_honors_non_power_of_two_limit() { + let result = + cheby::approx::adaptive::AdaptiveFit::new(Domain::new(-1.0, 1.0), |x: f64| (3.0 * x).sin()) + .max_degree(20) + .absolute_tolerance(0.0) + .relative_tolerance(0.0) + .build::() + .unwrap(); + assert_eq!(result.report.degree, 20); +} + +#[test] +#[cfg(feature = "adaptive")] +fn adaptive_fitting_allows_small_degree_limit() { + let result = + cheby::approx::adaptive::AdaptiveFit::new(Domain::new(-1.0, 1.0), |x: f64| x.sin()) + .max_degree(7) + .absolute_tolerance(0.0) + .relative_tolerance(0.0) + .build::() + .unwrap(); + assert_eq!(result.report.degree, 7); +} diff --git a/tests/error_branches.rs b/tests/error_branches.rs index 7301a9b..2b79f21 100644 --- a/tests/error_branches.rs +++ b/tests/error_branches.rs @@ -69,26 +69,33 @@ fn binary_error_branches() { // Wrong version (recompute checksum so we hit version branch, not checksum). let mut v = bytes.clone(); - v[0] = 99; // version byte - // Recompute the FNV-1a checksum the encoder uses. + // Use a supported layout with an unsupported version byte. + v[0] = 99; let payload_len = v.len() - 8; - let new_ck = v[..payload_len].iter().fold(0xcbf29ce484222325_u64, |h, b| { - (h ^ u64::from(*b)).wrapping_mul(0x100000001b3) - }); + let new_ck = v[..payload_len] + .iter() + .fold(0xcbf29ce484222325_u64, |h, b| { + (h ^ u64::from(*b)).wrapping_mul(0x100000001b3) + }); v[payload_len..].copy_from_slice(&new_ck.to_le_bytes()); let err = decode_f64_series(&v).unwrap_err(); assert!(matches!( err, - ChebyError::UnsupportedFormatVersion { found: _, expected: _ } + ChebyError::UnsupportedFormatVersion { + found: _, + expected: _ + } )); // Length mismatch: tamper coefficient count, recompute checksum. let mut l = bytes; l[4..12].copy_from_slice(&999_u64.to_le_bytes()); let payload_len = l.len() - 8; - let new_ck = l[..payload_len].iter().fold(0xcbf29ce484222325_u64, |h, b| { - (h ^ u64::from(*b)).wrapping_mul(0x100000001b3) - }); + let new_ck = l[..payload_len] + .iter() + .fold(0xcbf29ce484222325_u64, |h, b| { + (h ^ u64::from(*b)).wrapping_mul(0x100000001b3) + }); l[payload_len..].copy_from_slice(&new_ck.to_le_bytes()); assert_eq!( decode_f64_series(&l).unwrap_err(), diff --git a/tests/functional_pipeline.rs b/tests/functional_pipeline.rs index 5ec70bc..b70630a 100644 --- a/tests/functional_pipeline.rs +++ b/tests/functional_pipeline.rs @@ -120,3 +120,25 @@ fn segment_table_from_precomputed_segments_and_empty_case() { assert!(empty.eval_derivative(10.0).is_none()); assert!(empty.eval_both(10.0).is_none()); } + +#[test] +fn segment_table_reports_actual_end_for_short_final_segment() { + const N: usize = 9; + let table: ChebySegmentTable = + ChebySegmentTable::from_fn(|t: f64| t * t, 0.0, 2.5, 1.0); + + assert_eq!(table.len(), 3); + assert_abs_diff_eq!(table.end(), 2.5, epsilon = 0.0); + assert!(table.eval(2.49).is_some()); + assert!(table.eval(2.5).is_none()); + assert!(table.eval(2.75).is_none()); +} + +#[test] +fn segment_table_rejects_non_uniform_precomputed_segments() { + const N: usize = 3; + let first = ChebySegment::new([0.0; N], 0.5, 0.5); + let too_long_final = ChebySegment::new([0.0; N], 1.75, 0.75); + let result = ChebySegmentTable::try_from_segments(vec![first, too_long_final]); + assert!(result.is_err()); +} diff --git a/tests/minimax.rs b/tests/minimax.rs index ff7112f..4ff19ed 100644 --- a/tests/minimax.rs +++ b/tests/minimax.rs @@ -7,13 +7,7 @@ use cheby::core::{ChebyError, Domain}; #[test] fn remez_converges_for_smooth_target() { - let result = remez( - Domain::new(-1.0, 1.0), - 7, - f64::exp, - RemezOptions::default(), - ) - .unwrap(); + let result = remez(Domain::new(-1.0, 1.0), 7, f64::exp, RemezOptions::default()).unwrap(); assert!(result.converged); assert_eq!(result.alternations.len(), 9); // degree + 2 assert!(result.max_error < 1e-6); @@ -23,13 +17,7 @@ fn remez_converges_for_smooth_target() { #[test] fn remez_alternations_are_monotone() { - let result = remez( - Domain::new(-1.0, 1.0), - 5, - f64::cos, - RemezOptions::default(), - ) - .unwrap(); + let result = remez(Domain::new(-1.0, 1.0), 5, f64::cos, RemezOptions::default()).unwrap(); for w in result.alternations.windows(2) { assert!(w[0] < w[1]); }