From 913c735165485a036da56270eaa7dc1ed4abc290 Mon Sep 17 00:00:00 2001 From: Jeremy Howard Date: Tue, 10 Nov 2026 19:16:15 +0000 Subject: [PATCH] Share native math execution, selector geometry and conversions --- Cargo.toml | 4 +- DEV.md | 68 ++++----- README.md | 2 +- buffer/src/selection.rs | 47 +++--- core/src/agreement.rs | 34 ++--- core/src/arc.rs | 5 +- core/src/buffer.rs | 19 +-- core/src/buffer/output.rs | 3 +- core/src/convert.rs | 62 +++----- core/src/distribution.rs | 111 +++++++------- core/src/facts.rs | 92 +++++------ core/src/float.rs | 9 -- core/src/integer.rs | 37 ++--- core/src/lib.rs | 4 +- core/src/ln.rs | 33 +++- core/src/parallel.rs | 9 +- core/src/parts.rs | 10 +- core/src/pow.rs | 10 +- core/src/reduce.rs | 8 - core/src/reduce/run.rs | 10 +- core/src/slice.rs | 220 +-------------------------- core/src/special.rs | 27 ++-- core/src/structure.rs | 57 ++++++- core/src/tolerant.rs | 105 ++++--------- kernels/src/operation.rs | 9 +- kernels/src/operation/map.rs | 22 +-- kernels/src/operation/map/complex.rs | 8 +- nbs/index.ipynb | 2 +- src/array.rs | 22 +-- src/lib.md | 45 ++++-- src/lib.rs | 2 +- tests/accuracy.rs | 6 +- tests/linalg.rs | 5 +- tests/numeric.rs | 101 +++++++----- tests/slices.rs | 139 +++++++---------- tests/speed_python.txt | 2 + tests/support/mod.rs | 9 +- tests/test_speed.py | 25 ++- tests/test_structure.py | 7 + tests/test_ufuncs.py | 2 + timing/src/lib.rs | 27 ++-- 41 files changed, 599 insertions(+), 820 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 77dd7a1..928f6b9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,11 +7,11 @@ version = "0.1.60" edition = "2024" [workspace.dependencies] -rustymath = { path = ".", version = ">=0.1.59, <1", default-features = false } +rustymath = { path = ".", version = ">=0.1.60, <1", default-features = false } rustymath-buffer = { path = "buffer", version = ">=0.1.39, <1" } pyo3 = ">=0.29" pyo3-build-config = ">=0.29" -rustymath-core = { path = "core", version = ">=0.1.59, <1", default-features = false } +rustymath-core = { path = "core", version = ">=0.1.60, <1", default-features = false } rustymath-kernels = { path = "kernels", version = ">=0.1.59, <1" } rustymath-timing = { path = "timing", version = ">=0.1.56, <1" } half = { version = ">=2.7.1, <3", features = ["num-traits"] } diff --git a/DEV.md b/DEV.md index fa14e21..a374564 100644 --- a/DEV.md +++ b/DEV.md @@ -92,10 +92,10 @@ The dependency layers are `rustymath-core` → `rustymath-kernels` → the publi | `lib.rs`, `float.rs` | Public functions and objects; `Float`, `Native`, fused/unfused arithmetic and sealed implementation traits. | | `exp.rs`, `ln.rs`, `trig.rs`, `arc.rs`, `hyperbolic.rs`, `pow.rs`, `parts.rs`, `sqrts.rs` | Elementary scalar functions and specialized forms. | | `table.rs`, `dd.rs` | Table-based f64 logarithms/exponentials and double-word arithmetic. | -| `special.rs`, `lambert.rs`, `complex.rs` | Special functions; complex functions, which the facts table's complex formulas call, and the forms of `Abs` and `Angle`, the only complex forms left. | -| `slice.rs`, `parallel.rs` | Scalar/slice object interfaces, output storage and elementwise scheduling. | +| `special.rs`, `lambert.rs`, `complex.rs` | Special and complex functions, called by the facts table's formulas and native bulk routines. | +| `slice.rs`, `parallel.rs` | Shared item/slice loops, native bulk arguments, output storage and elementwise scheduling. | | `reduce.rs`, `stats.rs` | Folds, scans, axis traversal and statistical summaries. | -| `facts.rs` | What a planner knows about each function object before it reads an item: the `Range` of its results, from its arguments' ranges by interval rules, whether any result fails, and which argument a constant keeps (`Facts`). Also the one table of which reduction folds each function and which objects take characters, which rustymath, its Python crate and basedpl each expand. | +| `facts.rs` | The `Formula` table: supported kinds, executable numerical bodies, validity, result ranges, reduction and character rules, native bulk callbacks and partial derivatives. `Op` selects a formula at runtime. | | `structure.rs` | Concrete shape rules and positional plans for indexing, take/drop, rotation, windows, transpose and catenation. Shape functions take `usize` dimensions; `take_steps`, `rotate_steps` and related functions build the `Steps` that read each result. `Walk` and `Slide` recognize shared structural geometry. Symbolic lengths, trace decisions and Program construction belong to BPL. | | `arith.rs`, `integer.rs`, `convert.rs` | Wrapped/checked arithmetic, integer functions and vectorizable conversions. | | `distribution.rs`, `random.rs` | Checked distributions, function objects, sampling, rolls and deals. | @@ -113,55 +113,45 @@ The dependency layers are `rustymath-core` → `rustymath-kernels` → the publi | `buffer/` (`rustymath-buffer`) | The Python buffer protocol for compact buffers, selector decoding, dict keys and data frames, and the Python exception for each error kind. Basedpl's extension shares it. | | `timing/` (`rustymath-timing`) | What the speed tests share: the timing lock, the timing loop, instruction and byte counts, copy times and budget files. It depends on core, and core, kernels and the facade don't depend on it, so editing it rebuilds only the tests, the Python extensions and basedpl's `testkit`. | -## Functions and slice forms +## Functions, formulas and kernels -Free functions are generic over `Float`: f16 computes in f32 and rounds its result, while f32 and f64 use their own working type. `Native` names f32 and f64. `functions!` declares the free wrappers at the crate root and in `fused`/`unfused`; an `as` entry also declares the function object. `FUSED` selects hardware FMA where available. Don't introduce software FMA into targets without it. +Free functions are generic over `Float`: f16 computes in f32 and rounds its result, while f32 and f64 use their own working type. `Native` names f32 and f64. `functions!` declares the free wrappers at the crate root and in `fused`/`unfused`; an `as` entry also declares the formula object. `FUSED` selects hardware FMA where available. Don't introduce software FMA into targets without it. -Objects implement `Unary`, `Binary` or `Ternary` at concrete item types. They can hold settings. `one` computes a scalar, `checked` also reports failure, and `many` appends results to a vector and returns whether every item succeeded. Arguments are `Arg::One` or `Arg::Many`. Generic code calling an object states its bound, such as `Exp: Unary`; free functions need no object bound. +`facts::Formula` owns each object's supported dtypes, numerical bodies, validity, ranges, reductions and partial derivatives. `Kernel::resolve` selects a computing dtype and result dtype for known formula objects; `Kernel::of` selects from a runtime `Op`. `Kernel::{real,integers,complex,exact}` evaluates scalars. `map`, folds, scans, writes and products use the same resolved definition over `Read` plans. `Variant` selects checked or wrapping arithmetic. `Settings` supplies the execution tolerance outside item loops. -`Forms` defines the scalar and slice forms. Its unsafe contract requires `slice` to initialize every output slot, even when an item fails. `forms!` writes its implementations. `at real` writes f32/f64 forms and an f16 adapter through f32. Explicit type lists cover complex, integer and Boolean forms. `=> T` names an output type different from the input. Attributes before the object, such as `#[inline]`, go on its scalar forms, `scalar` and `judged_scalar`, which gives the result with whether it failed, as the slice form judges it. `checked_scalar` comes from `judged_scalar`. The integer arithmetic forms use `#[inline]`, so that basedpl's loops over items inline them, because LTO is off. The item body returns one of three outcomes: +`Formula::native_slices` supplies optional typed bulk callbacks in `NativeSlices`. Kernels choose them at the map boundary and otherwise run the formula's item body. Keep specialized logarithm, power/root, error-function, complex-comparison and distribution routines behind this table. A native callback writes every result slot even when its verdict reports failure. Real and integer callbacks return a Boolean verdict; complex callbacks leave the final validity check to the map boundary. f16 arguments widen to f32 at that boundary. -- A plain float/complex value fails on NaN from arguments that aren't NaN. Plain integer and Boolean values have no NaN failure. -- `Option` fails on `None`. -- `Total` never fails, including where NaN is the intended result. +`Formula::emit_valid` and `FAILS_NAN` define checked formula failures; IEEE NaNs are valid results for arithmetic that accepts them. `Kernel::may_raise_on` bounds those failures from argument ranges. BPL uses that query to admit deferred calls. Native bulk loops retain `Outcome` and `Total` for their verdicts: plain floating results fail on a new NaN, `Option` fails on `None`, a wrapped-result pair carries a fit flag, and `Total` never fails. Test the chosen formula or native verdict rather than imposing one NaN rule on every function. -Arithmetic objects deliberately accept IEEE NaNs. An integer result that doesn't fit its type fails: the integer forms of `Add`, `Subtract`, `Multiply` and `Abs` return `Option`, as the other integer functions do. A `: total` entry in `functions!` wraps the body in `Total`. Test an object's actual verdict rather than applying the ordinary NaN rule to every object. - -Each arity trait's `TOTAL` says whether an object never fails. It is true for an object whose body gives a `Total` output, or a plain integer or Boolean output. `forms!` derives it from the body's type. Basedpl's fusion runs only total objects, because a deferred call mustn't fail. - -A tolerant object's scalar forms compute at its `Tolerance` itself, which gives the exact result at 0 without a test, so that loop has no branch for each item. Only the slice forms choose `Exactly` through `with_tolerance!`. - -`arities!` gives the public traits one implementation through `Forms`. Keep objects concrete: their forms compile in rustymath instead of being instantiated in every caller. Forms that differ from the free functions, including `Sin`, `Cos`, `Tan`, `Pow` and `Root`, live beside their implementations. - -`kernels!` selects an object's function pointers at each dtype. `Fold`, `Scan`, `RightFold` and `LeftScan` share `AccumulationKernel` and `.run`, with the probe alone choosing the operation. Its typed kernel reads an optional right tail or scalar start in the total or computing dtype. +Tolerant formula objects retain their declared tolerance. A kernel's `Settings` overrides it for the execution. Native bulk routines select exact or tolerant implementations once per call through `with_tolerance!`. ### Loops and codegen -`Forms::slice` writes into provided slots. `many` and `pairs` reserve spare capacity through `appended`, without zeroing it. `written` views initialized slots for later passes. Keep that unsafe storage handling in these helpers. `many_into` and `Slots::pairs` run the same loops into a `parallel::Slots`, which `extended` and `Slots::append` hand out. A caller such as basedpl's fused chains then writes its last step straight into its result. Each writer fills every slot it gets, and a writer that doesn't panics before the vector's length changes. +Native slice callbacks write every provided output slot, including when they report failure. `appended` reserves spare capacity without zeroing it; `written` views initialized slots for later passes. Kernel outputs use `MapOutput`, `Output` and `parallel::Slots` to write allocated or caller-provided storage. Those helpers centralize storage initialization and the associated safety checks. -`each` runs the ordinary item loop. It checks outputs first and inspects arguments only for NaN outputs. `two_pass` computes a branch-free fast path, then recomputes its rare misses. `split` separates complex parts into real buffers, calls real slice forms and redoes exceptional items with the scalar function. `CHUNK` governs the cache-sized pieces used by these forms and their adapters. +`each` runs the ordinary item loop. It checks outputs first and inspects arguments only for NaN outputs. `two_pass` computes a branch-free fast path, then recomputes its rare misses. `CHUNK` governs cache-sized pieces in native bulk routines and conversion adapters. In two loops that ANDed a compound flag for each item, x86-64's loop vectorizer left the loop scalar, where arm64's vectorized it. On CI's x86-64 runner in October 2026, the log-normal density took about 2.5 times as long per item, and `3|v` about 3.7 times. `distribution.rs`'s `composed` flagged `!r.is_nan() | v.is_nan()` in its `finish` loop. `each`'s loops flag only `!r.is_nan()`, check `kept` afterwards where a result is NaN, and vectorize on both targets, and `composed` now does the same. The exact `Residue` flagged `!far(x, y)`, whose `x != 0` test shares its comparison with the residue's own `x == 0` select. x86-64 kept that select as a branch. Its loop now flags `near(x, y)`, which leaves out the test of `x`, and vectorizes on both targets, with 3% fewer instructions for `3|v` on arm64. `two_pass!` also vectorized it, but its NaN select and NaN check cost arm64 21% more instructions for `3|v`. basedpl's `bpl-tests disasm --x86` shows whether a loop vectorized. -Keep these performance rules when changing loops: +These measured loop patterns are useful starting points; their reasons matter more than preserving a particular implementation: -- Pass item bodies through `each!` or `two_pass!`. They mark closures `#[inline(always)]`; leaving a large body out of line can stop vectorization. -- Write separate closures for functions selected at runtime. A shared closure can leave the selector inside every iteration. `trig.rs` shows the pattern. -- Keep `pairs` at `#[inline(never)]`. Inlining a function whose loop runs behind `parallel::split`'s `&dyn` only repeats its setup at every call site. -- Check code in basedpl as well as in this crate when changing a generic function that it instantiates. The small test crate is not the same compilation context. -- Fill a large per-call buffer through `appended`'s spare capacity, or with `Vec::with_capacity` and `resize` as `linalg`'s `filled` does. On macOS, `vec![0; n]` allocates through calloc. Calloc's `madvise` costs more than writing the items. Its cost grows with the number of threads that touch the buffer. -- Write each hot loop as a plain loop over slices, with each value's representation chosen before the loop, or as a call to one of our own `#[inline]` helpers, nested helpers included. Local builds of the `test` profile are incremental, in up to 256 codegen units, with thin LTO across each crate's units. In an incremental build, a generic function without `#[inline]` gets one copy per crate, in one unit. std's iterator plumbing, such as `Map::fold` and `Vec`'s `spec_extend`, has no `#[inline]`. Thin LTO inlines such a copy into a caller in another unit only when the copy is below LLVM's import limit of 100 instructions. A longer one, such as a fused iterator `fold` over several adapters, stays out of line, and the loop runs as calls, sometimes one for each item. Which unit holds a copy changes with unrelated edits. A default trait method called on a concrete type also stays out of line unless it has `#[inline]`. CI builds the `test` profile with `CARGO_INCREMENTAL=0`. Its builds split code into units differently from local builds. A loop through `Steps::offsets` also runs `get`'s variant match and reloads its fields for every item, because the compiler can't rule out that the output buffer aliases the `Steps`. In basedpl, `BASEDPL_CALLS=pattern cargo test --test units -- --nocapture` lists each call from a matching function into another unit, in basedpl's build and in the rustymath rlib it links. - - Keep `Steps::offsets` a `Range` mapped through `get`. That iterator is `TrustedLen`, and `Vec::extend` and `zip` use their fast loops only for `TrustedLen` iterators, which stable Rust doesn't let a new type implement. A custom iterator type in its place took a write such as `x.[i]←w` from about 9 copies to 21. A hot loop over a table loops over its slice directly. - - Keep 64-bit integer multiplies out of loops that should vectorize. NEON has no 64-bit vector multiply. Multiplying each table entry by a stride as it was read stopped `Steps::combined`'s loop vectorizing, and a write such as `x.[i]←w` cost about a copy more. A table holds finished offsets, which `Table::scaled` multiplies once. - - Loop over an axis's positions through `Steps::for_each`, which matches the kind of step once, before its loop. LLVM can keep `get` out of line despite its `#[inline]`, even within the caller's unit. A loop over rows through `get` or `offsets` then makes a call for each row. Before `for_each`, each short row cost about 3.5 ns in `scatter_into` and `walk_into`: `[∞;0 1]⌷m` on 65,536 rows of 4 floats cost 15 copies, against 2.2 with it. +- Prefer `each!` or `two_pass!` for item bodies where they fit. They mark closures `#[inline(always)]`; leaving a large body out of line can stop vectorization. +- Separate closures can move runtime function selection outside the item loop. A shared closure can leave the selector inside every iteration. `trig.rs` shows the pattern. +- `pairs` currently uses `#[inline(never)]` because inlining a loop behind `parallel::split`'s `&dyn` repeated its setup at every call site. Revisit the choice if the calling structure or measurements change. +- Measure a generic change in its basedpl caller as well as the small Rust test crate when compilation context can affect the result. +- Filling a large per-call buffer through `appended`'s spare capacity, or `Vec::with_capacity` and `resize` as `linalg`'s `filled` does, avoids a measured macOS calloc cost. `vec![0; n]` uses calloc there; its `madvise` cost exceeded writing the items and grew with the number of threads touching the buffer. +- Plain slice loops with representations chosen before the loop, or calls to small `#[inline]` helpers, have avoided per-item dispatch and cross-unit calls in measured hot paths. Local builds of the `test` profile are incremental, in up to 256 codegen units, with thin LTO across each crate's units. In an incremental build, a generic function without `#[inline]` gets one copy per crate, in one unit. std's iterator plumbing, such as `Map::fold` and `Vec`'s `spec_extend`, has no `#[inline]`. Thin LTO inlines such a copy into a caller in another unit only when the copy is below LLVM's import limit of 100 instructions. A longer one, such as a fused iterator `fold` over several adapters, stays out of line, and the loop runs as calls, sometimes one for each item. Which unit holds a copy changes with unrelated edits. A default trait method called on a concrete type also stays out of line unless it has `#[inline]`. CI builds the `test` profile with `CARGO_INCREMENTAL=0`. Its builds split code into units differently from local builds. A loop through `Steps::offsets` also runs `get`'s variant match and reloads its fields for every item, because the compiler can't rule out that the output buffer aliases the `Steps`. In basedpl, `BASEDPL_CALLS=pattern cargo test --test units -- --nocapture` lists each call from a matching function into another unit, in basedpl's build and in the rustymath rlib it links. + - `Steps::offsets` currently uses a `Range` mapped through `get` for its `TrustedLen` behavior. That iterator is `TrustedLen`, and `Vec::extend` and `zip` use their fast loops only for `TrustedLen` iterators, which stable Rust doesn't let a new type implement. A custom iterator type in its place took a write such as `x.[i]←w` from about 9 copies to 21. A hot loop over a table loops over its slice directly. + - Moving 64-bit integer multiplies outside vector loops helps NEON, which has no 64-bit vector multiply. Multiplying each table entry by a stride as it was read stopped `Steps::combined`'s loop vectorizing, and a write such as `x.[i]←w` cost about a copy more. A table holds finished offsets, which `Table::scaled` multiplies once. + - `Steps::for_each` matches the kind of step once before its loop, avoiding a measured call per row. LLVM can keep `get` out of line despite its `#[inline]`, even within the caller's unit. A loop over rows through `get` or `offsets` then makes a call for each row. Before `for_each`, each short row cost about 3.5 ns in `scatter_into` and `walk_into`: `[∞;0 1]⌷m` on 65,536 rows of 4 floats cost 15 copies, against 2.2 with it. Elementary kernels aim for branch-free vectorizable loops. Incomplete gamma/beta and their inverses also use data-dependent series, continued fractions and iterations, with precision-appropriate stopping. Table gathers help AVX2 but are unavailable in NEON/wasm. Short dependency chains and vectorizable conversions matter as much as instruction count. -`Pow`/`Root` use table-based f64 arithmetic and narrow-width vector paths. Their scalar object forms have separate fast paths; f32 scalar calls use tables that the vector paths avoid. A whole exponent from 2 to 8 squares and multiplies, in the scalar form and in a slice form with a constant exponent. A degree of 2 or 3 takes `sqrt` or `cbrt`. A table-free f64 `Pow` measured about 30% slower: with `ln` and `exp` in one loop, their polynomial constants don't fit in registers. The loop then rebuilds them for every item. `cbrt` computes at each width's own precision. Trig slice forms redo large arguments after their fast reduction. f64 inverse-error slice forms group inputs by approximation region. `Erfc` chooses a single formula for blocks wholly on one side of its split. Normal/log-normal CDF and quantile objects reuse those optimized forms through `composed`. +`Pow`/`Root` use table-based f64 arithmetic and narrow-width vector paths. Their scalar formulas have separate fast paths; f32 scalar calls use tables that the vector paths avoid. A whole exponent from 2 to 8 squares and multiplies, in the scalar form and in a slice form with a constant exponent. A degree of 2 or 3 takes `sqrt` or `cbrt`. A table-free f64 `Pow` measured about 30% slower: with `ln` and `exp` in one loop, their polynomial constants don't fit in registers. The loop then rebuilds them for every item. `cbrt` computes at each width's own precision. Trig slice forms redo large arguments after their fast reduction. f64 inverse-error slice forms group inputs by approximation region. `Erfc` chooses a single formula for blocks wholly on one side of its split. Normal/log-normal bulk CDF and quantile routines reuse those optimized loops through `composed`. ### Parallel execution -`parallel::ranges` schedules independent ranges. `split` applies it to disjoint output slices; `many`, `pairs` and `extended` use that interface. A long call times its first unit, then cuts the rest into pieces of at least `PIECE_TIME` of work, at most four for each thread. `pieces` runs fewer pieces than threads on a pool of the largest power of two up to their number: memory-bound calls of a million items ran faster on small pools than on all 18 threads of an M5 Max. Chunked forms keep their temporary buffers in cache. Work and result initialization finish for every piece even if another piece reports failure. `parallel::extended` lets a caller such as basedpl's fusion fill a vector's spare capacity in pieces, through `Slots`. +`parallel::ranges` schedules independent ranges. `split` applies it to disjoint output slices; `pairs` and `extended` use that interface. A long call times its first unit, then cuts the rest into pieces of at least `PIECE_TIME` of work, at most four for each thread. `pieces` runs fewer pieces than threads on a pool of the largest power of two up to their number: memory-bound calls of a million items ran faster on small pools than on all 18 threads of an M5 Max. Chunked forms keep their temporary buffers in cache. Work and result initialization finish for every piece even if another piece reports failure. `parallel::extended` lets a caller such as basedpl's fusion fill a vector's spare capacity in pieces, through `Slots`. `parallel::stack` runs a function on each cell of a stack of matrices, for basedpl's `f⍤2` and the binding's matrix functions. Each cell writes its part of each of the stack's outputs, a tuple of one to three `Vec`s (`parallel::Outputs`), as `Stacked::qr` writes Q and R and `Stacked::svd` U, S and V. `parallel::each` gives a result of any type for each cell. Both schedule cells through `each_cell`, which times the first cells on the calling thread, then splits the rest across threads in rounds of about `ROUND_TIME` for each thread, with an interrupt check between rounds and one scratch value, such as an `LuWork` or `QrWork`, for each piece. A cell on another thread doesn't split again. @@ -187,11 +177,11 @@ This scheduler covers elementwise forms, folds and matrix products. Scans don't ### Complex functions and distributions -Complex functions live in `complex.rs`. The facts table's complex formulas call them, so kernels run them on each item. Only `Abs` and `Angle` keep scalar forms at complex arguments, which basedpl's polar `∧` calls. Exceptional magnitudes and branch-cut cases take their own paths inside each function. Inverse circle/hyperbolic functions use Kahan's formulas; complex ln Gamma uses Stirling, recurrence and reflection with SciPy's branch convention. +Complex functions live in `complex.rs`. The facts table's complex formulas call them, so kernels run them on each item. `Abs` and `Angle` are formula objects for complex magnitude and angle. Exceptional magnitudes and branch-cut cases take their own paths inside each function. Inverse circle/hyperbolic functions use Kahan's formulas; complex ln Gamma uses Stirling, recurrence and reflection with SciPy's branch convention. Complex errors are measured against a magnitude-based scale, with function-specific scales for ln Gamma, powers, Gamma and binomial. Quarter-turn directions are deliberately exact. See `complex.rs` for that contract rather than applying a separate ulp bound to each component. -`Distribution` validates parameters and holds a `Law`. `Pdf`, `Cdf` and `Quantile` select the law outside their item loops. Each law gives its density, cdf and quantile as slice forms, its `_many` methods, which by default run a map that holds its constants in a loop. `Pdf`, `Cdf` and `Quantile` run the slice forms for one item as for many. `Quantile` gives NaN for a p outside 0 to 1 before any law's formulas run (`quantiles`). The closed-form laws, normal, log-normal, uniform, exponential, Cauchy, Laplace, logistic and Weibull, use f32 working arithmetic for f32 inputs, with parameters and precomputed coefficients narrowed before the item loop. Other laws use f64 arithmetic and narrow their outputs. A closed-form law takes `erfc`, `erfcinv`, `ln`, `log1p`, `atan`, `tan`, `pow` and `root` from their slice forms through `composed`, because those functions vectorize in their slice forms but not inline in a law's loop. The Weibull law takes its powers from `Pow` and its quantile's root from `Root`, which keep its cdf and quantile within 4 ulps where e^(shape·ln z) didn't. Branches in a formula select between values computed beforehand. A division or an `exp` inside one arm of a branch keeps the loop scalar. The accuracy page lists the cases past 4 ulps, and `closed_forms_at_f32` in `tests/numeric.rs` checks the f32 forms against f64 to a relative error of 1e-4. Sampling, rolls and deals take a caller-owned RNG. Keep RNG ownership separate from the numerical algorithms. +`Distribution` validates parameters and holds a `Law`. `Pdf`, `Cdf` and `Quantile` are formula objects whose native bulk callbacks select the law outside their item loops. Each law gives its density, cdf and quantile as slice forms, its `_many` methods, which by default run a map that holds its constants in a loop. `Pdf`, `Cdf` and `Quantile` run the slice forms for one item as for many. `Quantile` gives NaN for a p outside 0 to 1 before any law's formulas run (`quantiles`). The closed-form laws, normal, log-normal, uniform, exponential, Cauchy, Laplace, logistic and Weibull, use f32 working arithmetic for f32 inputs, with parameters and precomputed coefficients narrowed before the item loop. Other laws use f64 arithmetic and narrow their outputs. A closed-form law takes `erfc`, `erfcinv`, `ln`, `log1p`, `atan`, `tan`, `pow` and `root` from their slice forms through `composed`, because those functions vectorize in their slice forms but not inline in a law's loop. The Weibull law takes its powers from `Pow` and its quantile's root from `Root`, which keep its cdf and quantile within 4 ulps where e^(shape·ln z) didn't. Branches in a formula select between values computed beforehand. A division or an `exp` inside one arm of a branch keeps the loop scalar. The accuracy page lists the cases past 4 ulps, and `closed_forms_at_f32` in `tests/numeric.rs` checks the f32 forms against f64 to a relative error of 1e-4. Sampling, rolls and deals take a caller-owned RNG. Keep RNG ownership separate from the numerical algorithms. ### Integers, exact numbers and polynomials @@ -217,7 +207,7 @@ Sorting, partitioning, rank selection and quantiles share exact ordering keys. F `gather` copies through a caller's conversion for cycle, repeat, indexed take and scatter. `positions` lists the positions a mask sets: it skips blocks of 64 counts that are all 0 and writes blocks that are all 1 in one loop. In a mixed block, it and the split into two classes in `search` write only the selected items through the 8-bit `SET` table. `walk` copies the items that one `Steps` per axis reaches: strides, clipped ranges, rotations and offset tables. Build a clipped range with `Steps::clipped` and a rotation with `Steps::rotated`. Each gives a `Stride` where nothing clips or wraps. Every `Clipped` entry reads at least one fill. Every `Rotated` entry wraps. Positions outside the source take the caller's fill. It merges neighbouring runs, copies long runs whole, transposes through `linalg::transpose_into`'s tiles and writes short windows through `windows`. `scatter_walk` writes through the same fast paths: transposed tiles, reversed and strided runs, index values converted inside the write loop (`Table::put_into`, as `Table::take_into` reads them), and cells at a table of starts (`placed`, as `picked` reads them). `walk` and `scatter_walk` take a `Convert`. With `Same`, they move items straight between source and result. A walk with `Same` of at least `SPLIT` items splits across threads in whole `CHUNK`s, through `parallel::extended`, unless the calling thread's calls run serially. Each thread walks its range a `CHUNK` at a time into a buffer of its own, then copies the buffer into its part of the result. On 4 threads it gains about 1.5 to 2 times. Memory bandwidth limits the gain, as does the `PIECE_TIME` of work that `parallel::ranges` gives each piece at least. Any other conversion runs as a separate pass over blocks of about `CHUNK` items: a walk gathers a block of source items, then converts it, and a scatter converts a block of items, then writes it. Their walks therefore compile once for each item type (`Item`), and their fixed-width loops, `windows` and `windows_into`, once for each `Unit`: an unsigned integer, or an array of two, with an item type's size and alignment. Two slice casts, `units` and `written_as`, move items as units. Both casts are sound by the promises of two unsafe traits. An `Item` type has no padding and has its `Unit`'s size and alignment, which `units!` checks at compile time. Every bit pattern of a `Unit` is a value of it, so items read as units are valid units. `written_as` writes only units that `units` read from items of the same type, so each slot it fills holds a valid item. That matters for `bool` and `char`, which have bit patterns that aren't values. A conversion's own pass costs about 0.25 ns per item on this Mac, measured on 20 October 2026. A conversion inside the walk cost nothing extra, because the walk waits on memory. So `scatter_walk` writes short cells along a stride from Booleans and integers into `f64` through fixed-width loops that convert as they write (`Item::cells_from`). `units!` gives those five source types alone such loops, which add about 17 KB of basedpl's machine code. `Steps::through` composes two walks into one only where the result stays regular, or where either walk is already an offset table. Elsewhere callers walk the inner copy into a buffer first. That costs one pass, where building and reading a table of every offset costs two. Basedpl reuses these paths for structural operations. Measure search/gather changes there as well. -`agreement` pairs two layouts by BPL's broadcasting rule. Each argument's two pairings compose into one mapping through `Mapping::composed`, which builds a table of offsets where they don't compose structurally. A kernel reads each argument through a `Read`, a `ReadPlan` over the argument's bound sources. Readers take items `in_place` where they can: one item, a contiguous run, items within one repeat, or a stretch of one cycle. Other items copy through `copied`, whose one walk, `copy_walked`, fills a repeated item, copies a cycle through `gather::cycled` and walks `Steps` through `walk_plan`. `Read::copy_items` walks a `Steps` map straight into its output where the source starts at 0, the map reads the type's default fill and every position lies inside the data. Its other reads go item by item. `tiled_cycle` writes a cycle shorter than `LONG_CYCLE` out once, so that a run reads it in place at any phase. `Mapping::line_block` sizes a block as whole rows where an argument's innermost axis skips items, as a transpose's does. `Agreement` is generic over its layout type, which basedpl uses to carry its renderer. `Mapping::walk` turns a structural walk's tables into a `Mapping`, and `Mapping::tables` gives tables back for code that walks or classifies them, such as `Walk::new`. basedpl's pending copies, programs and `∂` record hold that one map type, and compose maps through `Mapping::through`. A mapping may carry a checked one-item native `Fill`; numerical readers retain it, and `to_walk` returns it beside the base and steps. Geometry-only code uses `tables`, while source conversions retype the fill through `cast_fill`. Maps needing distinct fills do not compose. +`agreement` pairs two layouts by BPL's broadcasting rule. `Agreement` is generic over its layout type, which BPL uses to retain its renderer. `Mapping` describes positions as `Single`, `Linear`, `Tiled` or `Steps`; it carries no fill value. `ReadPlan` names source inputs and their logical lengths, maps, joins, index values and default or invalid positions. Maps compose through `ReadPlan::remapped`, retaining their missing-position plan. `Missing::Scalar` reads a separate one-item source, and `Missing::Input` supplies another plan. `Binding::full` binds a whole source; `Binding::window` binds a range within its producer's logical extent. A kernel reads the bound sources through `Read`. `Read::part` borrows contiguous or repeated items where possible and copies other requested ranges into typed scratch; `part_view` lends a contiguous run at its storage dtype or fills a scratch buffer. `ReadPlan::runs` exposes source, fill, index and invalid runs, and `RunSink::walk` accepts a structural walk whole. The retained `Walked` helpers also serve mapped linear-algebra and reduction paths: `in_place`, `copied`, `tiled_cycle` and `copy_walked`. `Mapping::{entries,tables,to_walk}` exposes positional geometry without data or fill conversion. ### Optimization @@ -253,13 +243,13 @@ The native SVD, Hermitian and general solvers compile everywhere. On macOS the S Expose useful mathematical and array capabilities through a thin Python interface to the shared Rust algorithms. A public Rust helper does not require a matching Python binding. Keep Python argument decoding, defaults, scalar-return policy and result wrappers in `py/src/` or `python/rustymath/`; keep generic mathematics and array algorithms in RM even when a Python convenience is their only direct caller. Related bindings share the core's names, argument order and dtypes. A function with NumPy's name needn't behave as NumPy's does. `remainder` and `floor_divide` round through the comparison tolerance where one is set, set functions keep first-occurrence order, orderings put NaN last, and `diagonal` and `trace` take the last two axes, as the matrix functions do. -Python operations are eager. `python/rustymath/create.py` gives fill constructors their names/defaults over one native fill primitive, and names the dtype facts that Rust returns. Sequence generation stays in `create`. Structural functions and bracket indexing in `py/src/structure.rs` delegate selector decoding to `buffer/src/selection.rs`, shared with BasedPL, for the selection and writeback paths. Brackets plan through `structure::plan_parts`, the planner that BasedPL's `⌷` and Python indexing share, with broadcasting, as BasedPL's `.[]` does: index arrays pair by broadcasting, by arithmetic's rule for names. A selector with one axis name selects along the axis with that name (`structure::addressed`). A selection by strides, padding or rotations is a view that shares storage, and a selection by index arrays copies its items (`Array::selected`). `reshape` keeps a view's mapping, which reads the same items in order under any layout. A view whose mapping isn't affine, such as padded windows or a reshape of windows, keeps the items that its first whole read gathers, and later reads use them (`ArrayView::keeps`). An affine view, such as a slice or an affine reshape of one, keeps nothing, and a reader that can't use its strides gathers it again. A small selection that doesn't compose with a view's mapping lists its offsets, and a larger one gathers the view first (`agreement::small_read`, which BasedPL's reads of pending arrays share). `take`, `pad` and `windows` read positions past an end by `structure::Edge`'s modes, which BasedPL's `•edge` names: `take` through `Edge::positions`, and `pad` and padded windows through `Steps::padded`. `structure::Windows` plans windows for `windows` as it does for BasedPL's `↕`. A `fill=` is a `Fill` checked at the array's dtype (`Fill::new`), which the selection's mapping reads (`Mapping::filled`). Assignment detaches shared storage, preserving independent values. `rm.at` combines through the same selection with `Array::combine_at` and its function's `Scatter` kernels, sharing `Array::assign`'s checks. A view retains its whole base until `compact()` copies its logical items, beside any items it has gathered, and even when it is a contiguous subrange. Methods call the same binding functions as package-level operations. `python/rustymath/structure.py` composes conveniences such as trace and repeated differences. Differences run subtraction kernels on two mappings of the input, not gathered copies. `python/rustymath/random.py` subclasses the native generator for distribution conveniences and multivariate-normal sampling. `py/src/random.rs` owns its RNG. Sampling, rolls and deals stay in the shared core. +Python operations are eager. `python/rustymath/create.py` gives fill constructors their names/defaults over one native fill primitive, and names the dtype facts that Rust returns. Sequence generation stays in `create`. Structural functions and bracket indexing in `py/src/structure.rs` delegate selector decoding to `buffer/src/selection.rs`, shared with BasedPL, for the selection and writeback paths. Brackets plan through `structure::plan_parts`, the planner that BasedPL's `⌷` and Python indexing share, with broadcasting, as BasedPL's `.[]` does: index arrays pair by broadcasting, by arithmetic's rule for names. A selector with one axis name selects along the axis with that name (`structure::addressed`). Every selection remains a view of shared sources, including index arrays and padding (`Array::selected`). `reshape` keeps a view's mapping, which reads the same items in order under any layout. Readers borrow runs they can use in place and gather other requested ranges into scratch storage. A view keeps its `ReadPlan` and source buffers, without a materialization cache. `take`, `pad` and `windows` read positions past an end by `structure::Edge`'s modes, which BasedPL's `•edge` names: `take` through `Edge::positions`, and `pad` and padded windows through `Steps::padded`. `structure::Windows` plans windows for `windows` as it does for BasedPL's `↕`. A `fill=` is checked at the array's dtype through `Buffer::kept` and retained as a separate one-item source, selected by `Missing::Scalar`. Assignment detaches shared storage, preserving independent values. `rm.at` combines through the same selection with `Array::combine_at` and its function's `Scatter` kernels, sharing `Array::assign`'s checks. A view retains its whole base until `compact()` copies its logical items, even when it is a contiguous subrange. Methods call the same binding functions as package-level operations. `python/rustymath/structure.py` composes conveniences such as trace and repeated differences. Differences run subtraction kernels on two mappings of the input, not gathered copies. `python/rustymath/random.py` subclasses the native generator for distribution conveniences and multivariate-normal sampling. `py/src/random.rs` owns its RNG. Sampling, rolls and deals stay in the shared core. `py/src/array.rs` holds `rustymath.Array`, which wraps a rustymath `array::Array`. `numpy.asarray` receives a read-only snapshot through `rustymath_buffer::export`. Each export pins its exact backing allocation and signed-stride geometry until `__releasebuffer__`, independently of later assignment to the Python owner. `rm.array(x)` shares existing Rustymath storage with independent values and copies foreign inputs. `py/src/convert.rs` retains the call's ownership, scalar classification and NEP-50 promotion policy. Its arguments expose core `ArrayView`s. Dense callers prepare once through core `Items`. Elementwise callers keep their mappings. Arrays of different dtypes promote through `DType::common`. Nested lists, tuples, `range`, `dict` and `str` follow the same converter. `None` reads as NaN as an item. A bare `None` doesn't convert, so an operator given one returns `NotImplemented`. A pandas `DataFrame` reads through `rustymath_buffer::frame`: its columns convert together, as a list of them would, and transpose into a matrix of rows by columns keyed by its labels. `rustymath_buffer::labels` is the one rule for a dict's keys and a frame's column labels, shared with BasedPL. `give` returns a rank-0 result as a Python number when every data argument was a Python scalar. `rustymath_buffer::read_array` reads sibling snapshots with keys and names. The shared buffer reader normalizes signed bytes, wider unsigned integers and foreign byte order, borrows contiguous aligned items, and uses CPython's `PyBuffer_ToContiguous` to copy other layouts into C order. It doesn't use PyO3's `PyUntypedBuffer`, which rejects rank-0 exports. The crate has no Rust tests, because PyO3's `extension-module` feature unifies across the workspace, and a test that embeds Python can't link. `tests/test_array.py` tests it through the class. `py/src/functions.rs` holds the elementwise functions, reductions and scans, one entry each in the `functions!` and `reductions!` macros. A function converts its data arguments through `convert`, finds its object's kernels with `kernels!`, and runs them at the dtype that `Kernels::computes` gives. `convert` promotes the arrays' dtypes with `DType::common`. It then applies `DType::common_weak` for the Python values, as NEP 50 does for weak scalars. BPL shares both rules. An argument whose dtype has no kernel computes at the first dtype that has one: `int64`, then `float64`, then `complex128` for Booleans and integers, and a wider float, then a complex dtype at least as wide, for floats. `Kernels` casts the argument to that dtype first. The object's slice form then reads the items in place and splits independent ranges across threads, as it does for Rust callers. `checked` turns a kernel's result into the call's: the `missing` error for a dtype that the function doesn't take, the wrapped integer result where `Kernels::wraps` says that the kernel wrote one, as NumPy does, and `OverflowError` for any other integer result that doesn't fit its dtype. Comparisons, and the functions that round or divide through them, take one scalar `ct`. Where `ct` is None, `functions::tolerance` reads the process's tolerance, which `set_tolerance` sets and which starts at 0, exact. The Rust core takes every tolerance as an argument, so the setting lives in the binding, and BPL's own tolerance setting never changes Python's results. A tolerance of 0 runs the exact algorithms, through `Exactly` (`with_tolerance!`). -With no `axis`, a reduction reduces the whole array and gives a Python number, and a scan runs over the flattened array. Reductions take axis positions, names or tuples of them. `accumulate(f, x, axis=None, start=None)` scans any native binary function from the left, with an optional scalar seed. The same `functions!` entries supply its callbacks and `inner`'s. Known reductions retain their optimized scans. Comparisons keep the computing dtype, reading later results as 0 or 1. A sort, `argsort` or `lexsort`, takes one axis and works along the last axis by default. Grading shares native-width key construction and stable sorting with BPL, including heterogeneous lexsort keys. Folds and scans run through `reduce::along`, `scan_along` and `left_scan`. A view that reads a block of its storage, such as a transpose or a reversal, lends the block, which folds in place (`ArrayView::fold_items`). Any other view folds through its mapping with `reduce::along_mapped` at its own dtype, which folds runs of at least `reduce::RUN` items where they lie and copies shorter runs a piece at a time. A view that keeps what it gathers (`ArrayView::keeps`) and reads shorter runs gathers once instead, and later folds read the gathered items. A view at a dtype that the reduction doesn't fold gathers its items, converted. `minmax` over the whole array takes `reduce::extremes`' single pass for Booleans and integers, and two folds otherwise. `logical_scan` takes a 2×2 truth table. +With no `axis`, a reduction reduces the whole array and gives a Python number, and a scan runs over the flattened array. Reductions take axis positions, names or tuples of them. `accumulate(f, x, axis=None, start=None)` scans any native binary function from the left, with an optional scalar seed. The same `functions!` entries supply its callbacks and `inner`'s. Known reductions retain their optimized scans. Comparisons keep the computing dtype, reading later results as 0 or 1. A sort, `argsort` or `lexsort`, takes one axis and works along the last axis by default. Grading shares native-width key construction and stable sorting with BPL, including heterogeneous lexsort keys. Folds and scans call `Kernel::fold_read` and `Kernel::scan` over shared `Read` plans. The kernel selects accumulation dtypes and reads or converts each requested block through those plans. Views retain no gathered-item cache. `minmax` over the whole array takes `reduce::extremes`' single pass for Booleans and integers, and two folds otherwise. `logical_scan` takes a 2×2 truth table. Statistical bindings in `py/src/functions.rs` share conversion and axis preparation through `stats::Axes`. Numerical loops live in `stats`, including two-pass variance, weighted averages, covariance and histogram counting. Float16 statistics accumulate and return float32. The 2-norm folds through `reduce::Squares`, which computes and sums the squares in float64, for float32 data too, with no buffer of squares. `linalg::norm`'s Frobenius norm is that 2-norm over every item. Complex variance and norms return real values. @@ -287,7 +277,7 @@ Python tests are grouped by topic in `tests/test_ufuncs.py`, `test_linalg.py`, ` To see a new or changed Python test fail on the build before your change, run it on the baseline build that “Commands and builds” notes: `run_pytest(test, rms, build=base)`, where `test` is a node ID such as `tests/test_structure.py::test_edge_modes_pad_take_and_windows`. `run_pytest(test, rms)` runs it on your build. In the shell, `bpl-tests run_pytest TEST --path PATH --build BASE` runs it on the baseline build. -Elementary real accuracy checks use a four-ulp bound with an additional allowance for f64 libm references. Gamma and incomplete functions use scaled errors and function-specific bounds; f64 iterative checks have a wider allowance. `tests/accuracy.rs` defines those bounds. Real checks cover fused/unfused arithmetic, f16 inputs, sampled wider inputs and reference files. Complex checks exercise public scalar/slice forms at the build's FMA mode, including branch cuts. `testdata/references.py` generates the mpmath references; preserve their notices and `THIRD-PARTY.md`. +Elementary real accuracy checks use a four-ulp bound with an additional allowance for f64 libm references. Gamma and incomplete functions use scaled errors and function-specific bounds; f64 iterative checks have a wider allowance. `tests/accuracy.rs` defines those bounds. Real checks cover fused/unfused arithmetic, f16 inputs, sampled wider inputs and reference files. Complex checks exercise scalar formulas and bulk kernels at the build's FMA mode, including branch cuts. `testdata/references.py` generates the mpmath references; preserve their notices and `THIRD-PARTY.md`. `nbs/index.ipynb` generates the README and home page, with short examples linking to feature notebooks. `nbs/sidebar.yml` lists those pages. `src/lib.md` supplies the Rust documentation; `nbs/rust.qmd` symlinks it as a standalone Quarto page. All library manifests set `doctest = false`. Keep documentation examples, and put meaningful example assertions in ordinary Rust tests. Doc comments in `.rs` files show examples as `text` blocks. Edit those sources, not the generated README. `nbdev-test` checks notebook code; `nbdev-preview` serves saved outputs because `_quarto.yml` disables execution. `.github/workflows/docs.yml` deploys the site. diff --git a/README.md b/README.md index a4a1f78..77d6622 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ Float functions give results within four units in the last place (ulps) of the c pip install rustymath ``` -Rust users can add the crate to `Cargo.toml`. See [Using rustymath from Rust](https://answerdotai.github.io/rustymath/rust.html) for scalar and slice forms, build settings, accuracy and speed. +Rust users can add the crate to `Cargo.toml`. See [Using rustymath from Rust](https://answerdotai.github.io/rustymath/rust.html) for function objects and native bulk routines, build settings, accuracy and speed. ``` python import rustymath as rm diff --git a/buffer/src/selection.rs b/buffer/src/selection.rs index 9d82405..fda3246 100644 --- a/buffer/src/selection.rs +++ b/buffer/src/selection.rs @@ -1,15 +1,30 @@ use crate::py_err; use pyo3::{exceptions::{PyIndexError, PyKeyError, PyTypeError}, prelude::*, types::{PyInt, PyList, PySlice, PyString, PyTuple}}; -use rustymath::{buffer::View, gather::{IndexSource, Table}, layout::{Layout, position}, structure::{self, IndexPart, Plan, Selector}}; +use rustymath::{buffer::View, gather::{IndexSource, Table}, layout::{Layout, position}, structure::{self, IndexPart, Plan, Selector, SliceBound, SliceSpec}}; use std::sync::Arc; -/// The position of key `name` on `axis` of `layout`, as a slice's bound reads it: a missing bound is a `KeyError`, since a slice needs a position on the axis. -fn key(layout: &Layout, axis: usize, name: &str) -> PyResult { - layout.keys(axis).and_then(|k| k.position(name)).ok_or_else(|| PyKeyError::new_err(name.to_owned())) +/// Decode a Python slice without binding it to an array's current length. +fn slice_spec<'py>(slice: &Bound<'py, PySlice>) -> PyResult { + let py = slice.py(); + let (start, stop, step) = (slice.getattr("start")?, slice.getattr("stop")?, slice.getattr("step")?); + let numeric = |v: &Bound<'py, PyAny>| if v.is_instance_of::() { py.None().into_bound(py) } else { v.clone() }; + let unpack = py.get_type::().call1((numeric(&start), numeric(&stop), step))?; + let (mut first, mut last, mut stride) = (0, 0, 0); + // SAFETY: unpack is a slice and the output pointers refer to local isize values. + let result = unsafe { pyo3::ffi::PySlice_Unpack(unpack.as_ptr(), &mut first, &mut last, &mut stride) }; + if result < 0 { return Err(PyErr::fetch(py)); } + let bound = |v: &Bound<'_, PyAny>, at| -> PyResult> { + Ok(if v.is_none() { None } + else if let Ok(s) = v.cast::() { Some(SliceBound::Key(s.to_str()?.into())) } + else { Some(SliceBound::Position(at)) }) + }; + Ok(SliceSpec { start: bound(&start, first)?, stop: bound(&stop, last)?, step: stride }) } -/// Whether `selector` decodes `obj` itself, with no conversion: a slice, a key, an integer or a list of keys. Other selectors go to its `array` callback. -pub fn native(obj: &Bound<'_, PyAny>) -> bool { obj.is_instance_of::() || obj.is_instance_of::() || key_names(obj).is_some() } +/// Selection errors, including a missing key used as a slice bound. +pub fn error(error: rustymath::error::Error) -> PyErr { + if let Some(key) = error.message.strip_prefix("missing key: ") { PyKeyError::new_err(key.to_owned()) } else { py_err(error) } +} /// The keys that `obj` names: one for a string, and one for each item of a list or tuple of strings. `None` for any other selector. pub fn key_names(obj: &Bound<'_, PyAny>) -> Option> { @@ -20,21 +35,8 @@ pub fn key_names(obj: &Bound<'_, PyAny>) -> Option> { /// One Python selector, independent of the selectors on the other axes. Keys read by `structure::key_selector`, so a missing key reads the fill, as in BPL. pub fn selector(layout: &Layout, axis: usize, obj: &Bound<'_, PyAny>, array: impl FnOnce() -> PyResult<(Selector, bool)>) -> PyResult<(Selector, bool)> { let size = *layout.shape().get(axis).ok_or_else(|| PyIndexError::new_err("too many axis selectors"))?; - if !native(obj) { return array(); } if let Ok(slice) = obj.cast::() { - let py = obj.py(); - let bound = |name| -> PyResult> { - let value = slice.getattr(name)?; - if let Ok(s) = value.cast::() { Ok(key(layout, axis, s.to_str()?)?.into_pyobject(py)?.into_any()) } else { Ok(value) } - }; - let (start, mut stop, step) = (bound("start")?, bound("stop")?, slice.getattr("step")?); - if slice.getattr("stop")?.is_instance_of::() { - let reverse = !step.is_none() && step.extract::()? < 0; - let at = stop.extract::()?; - stop = if reverse && at == 0 { py.None().into_bound(py) } else { (at + if reverse { -1 } else { 1 }).into_pyobject(py)?.into_any() }; - } - let s = py.get_type::().call1((start, stop, step))?.cast_into::()?.indices(size as isize)?; - return Ok((Selector::Range { start: s.start.max(0) as usize, len: s.slicelength, step: s.step }, false)); + return slice_spec(slice)?.selector(layout, axis).map(|s| (s, false)).map_err(error); } if let Some(names) = key_names(obj) { let one = obj.is_instance_of::(); @@ -45,7 +47,7 @@ pub fn selector(layout: &Layout, axis: usize, obj: &Bound<'_, PyAny>, array: imp let i = obj.extract::()?; return Ok((Selector::At(position(i, size).ok_or_else(|| PyIndexError::new_err(format!("position {i} is outside an axis of length {size}")))?), true)); } - unreachable!("`native` admits only the kinds decoded above") + array() } /// The selector that an array's items give: a character string picks a key, or the fill where the axis lacks it, and a table reads Booleans and other numbers in place as positions, as `⌷` does. @@ -74,6 +76,7 @@ pub fn parts<'py>(index: &Bound<'py, PyAny>) -> PyResult() && ["start", "stop", "step"].iter().all(|n| p.getattr(*n).is_ok_and(|v| v.is_none())) { IndexPart::All } + else if let Ok(slice) = p.cast::() { IndexPart::Slice(Arc::new(slice_spec(slice)?)) } else { IndexPart::Axis(p) }) }).collect() } @@ -92,7 +95,7 @@ fn select_parts( structure::plan_parts(layout, &claims(&written, &names), None, true, |i, axis| { let IndexPart::Axis(p) = &written[i] else { unreachable!("an axis refers to its written payload") }; selector(layout, axis, p, || array(layout, axis, p)) - }, py_err) + }, error) } /// For each axis of `layout`, the keys that a write by `index` names with Python's own strings and lists of strings and the axis lacks (`structure::missing_keys`). The write adds them before it plans (`structure::extended`), as BPL's assignment does. `name` gives a part's one axis name, as for `selection`. diff --git a/core/src/agreement.rs b/core/src/agreement.rs index d16c0d0..6aaf7d7 100644 --- a/core/src/agreement.rs +++ b/core/src/agreement.rs @@ -389,6 +389,11 @@ impl<'s> Read<'s> { if range.start > range.end || range.end > self.plan.len() { return Err(Error::new(Kind::Length, "read range exceeds its result")); } self.bindings() } + /// Append the converted range without changing existing output items. + pub(crate) fn append(&self, range: Range, checked: bool, output: &mut Vec) -> Result<(), Error> { + self.range(&range)?; + self.read_onto(range, checked, output) + } /// The items that `range` reads as one run of a source's storage, at whatever dtype it holds: the borrowed run and its stride, 1 for consecutive items and 0 for one item that every position reads (`ReadPlan::run`). fn run_view(&self, range: &Range) -> Option<(crate::buffer::View<'s>, usize)> { let (source, start, stride) = self.plan.run(range)?; @@ -526,7 +531,7 @@ impl<'s> Read<'s> { let binding = self.sources.get(source)?; let len = binding.view.len(); (binding.start == 0 && len == input.len() && match *map { - Mapping::Single => true, Mapping::Linear(r) => r > 0 && end.saturating_sub(1) / r < len, Mapping::Tiled(m) => m > 0 && m <= len, + Mapping::Single => len > 0, Mapping::Linear(r) => r > 0 && end.saturating_sub(1) / r < len, Mapping::Tiled(m) => m > 0 && m <= len, Mapping::Steps { base, ref steps } => Steps::furthest(base, &steps[..]).is_none_or(|furthest| furthest < len) }).then_some(binding) } @@ -979,7 +984,7 @@ impl Mapping { let base = base + outer.iter().map(|s| if let Steps::Stride { start, .. } = *s { start } else { 0 }).sum::(); Some(if base == 0 && Steps::in_order(inner) { Period::Cycles(p) } else { Period::Walks { len: p, base, steps: inner } }) } - /// Whether `expanded` copies most blocks of `block` items read through this mapping, rather than reading them where they lie: steps whose runs of adjacent offsets (`Steps::contiguous`) are shorter than a block, as a transpose's, a reversal's, windows' and a column slice's are. Repeats and cycles don't count: a run that stays inside one repeat or cycle reads in place, and the caller decides where runs end. + /// Whether `expanded_items` copies most blocks of `block` items read through this mapping, rather than reading them where they lie: steps whose runs of adjacent offsets (`Steps::contiguous`) are shorter than a block, as a transpose's, a reversal's, windows' and a column slice's are. Repeats and cycles don't count: a run that stays inside one repeat or cycle reads in place, and the caller decides where runs end. pub fn copies(&self, block: usize) -> bool { match self { Self::Steps { steps, .. } => Steps::contiguous(steps) < block, _ => false } } /// Whether a walk of blocks of `block` items through this mapping reads each item from its own cache line, where a walk of the whole copy reads them through the tiled transpose (`Steps::transposing`): after `simplify`, the last two axes form a transpose whose inner axis holds a block or more, so that each block walks one strided run. A copy made whole then costs much less than reads through it a block at a time. pub fn tiles_whole(&self, block: usize) -> bool { @@ -1218,23 +1223,17 @@ pub(crate) const LINE: usize = 128; /// The most items in a block of whole rows from [`Mapping::line_block`]. const LINE_BLOCK: usize = 1 << 16; -/// The one item that every position in `start..start + n` reads, where `mapping` repeats a single item over all of them. `None` where the positions read more than one item. -fn single(mapping: &Mapping, data: &[A], start: usize, n: usize) -> Option { - mapping.run(start, n).and_then(|(o, step)| (step == 0).then(|| data[o])) -} - -/// Items `start..start + n` of an argument that `mapping` lays out, read as `expanded` reads them but always as a slice of all `n` items. A repeated item fills `buf`. +/// Items `start..start + n` of an argument that `mapping` lays out. Borrowed runs stay in place; other reads fill `buf` with all `n` items. pub fn expanded_items<'a, A: Element>(mapping: &Walked, data: &'a [A], start: usize, n: usize, buf: &'a mut Vec) -> &'a [A] { - if let Some(x) = single(mapping, data, start, n) { - buf.clear(); - buf.resize(n, x); - return buf; + match in_place(mapping, data, start, n) { + Some(Arg::Many(items)) => items, + Some(Arg::One(x)) => { buf.clear(); buf.resize(n, x); buf }, + None => copied(mapping, data, start, n, buf) } - match expanded(mapping, data, start, n, buf) { Arg::Many(items) => items, Arg::One(_) => unreachable!("`single` gives every repeated item") } } /// Items `start..start + n` of an argument that `mapping` lays out, where they need no copy: a single item, a contiguous run, items within one repeat, or items of a cycle: within one period of a mapping that cycles (`Period::Cycles`), and as far as `data` holds a `Tiled` cycle, once or written out past its end (`tiled_cycle`). A `Steps` mapping reads a forward run in place through `run`. -pub fn in_place<'a, A: Copy>(mapping: &Walked, data: &'a [A], start: usize, n: usize) -> Option> { +fn in_place<'a, A: Copy>(mapping: &Walked, data: &'a [A], start: usize, n: usize) -> Option> { match mapping.run(start, n) { Some((o, 0)) => return data.get(o).copied().map(Arg::One), Some((o, 1)) => return data.get(o..o.checked_add(n)?).map(Arg::Many), @@ -1250,11 +1249,6 @@ pub fn in_place<'a, A: Copy>(mapping: &Walked, data: &'a [A], start: usize, n: u data.get(from..from.checked_add(n)?).map(Arg::Many) } -/// Items `start..start + n` of an argument that `mapping` lays out over the frame, read `in_place` where they can be. Others copy into `buf`. A repeated item fills each of its periods. A cycled argument copies through `gather::cycled`, and a `Steps` mapping through `walk_plan`, over the walk that its reader keeps (`Walked`). -pub fn expanded<'a, A: Element>(mapping: &Walked, data: &'a [A], start: usize, n: usize, buf: &'a mut Vec) -> Arg<'a, A> { - if let Some(arg) = in_place(mapping, data, start, n) { return arg; } - Arg::Many(copied(mapping, data, start, n, buf)) -} /// The cycle that `mapping`'s positions read from `data`, written into `buf` as `Period::tiling` says for runs of up to `run` positions, with the `Tiled` mapping that then reads it in place: a cycle read in order copies its first period, and one that walks entries (`Period::Walks`) gathers its period through them. A cycle shorter than its cover is written out past its end, so that a run reads it in place at any phase. `None` where `tiling` gives none, leaving `buf` alone. pub fn tiled_cycle(data: &[A], mapping: &Mapping, run: usize, buf: &mut Vec) -> Option { @@ -1272,7 +1266,7 @@ pub fn tiled_cycle(data: &[A], mapping: &Mapping, run: usize, buf: & Some(Mapping::Tiled(p)) } -/// Items `start..start + n` that `mapping` lays out, copied into `buf`, cleared first: what `expanded` gives where `in_place` refuses them. A caller that has just had `in_place` refuse them copies here directly, and the copy records that refusal. +/// Items `start..start + n` that `mapping` lays out, copied into `buf`, cleared first: what `expanded_items` gives where `in_place` refuses them. A caller that has just had `in_place` refuse them copies here directly, and the copy records that refusal. pub fn copied<'a, A: Element>(mapping: &Walked, data: &[A], start: usize, n: usize, buf: &'a mut Vec) -> &'a [A] { buf.clear(); let reason = copy_walked(mapping, data, start, n, A::FILL, buf); diff --git a/core/src/arc.rs b/core/src/arc.rs index f94c35d..7e80cae 100644 --- a/core/src/arc.rs +++ b/core/src/arc.rs @@ -70,8 +70,9 @@ pub(crate) fn atan(x: R) -> R { (hi + ((u * z).mla::(p, u) + lo)).copysign(x) } -// The Cauchy distribution's cdf runs this slice form (`distribution::through`). -forms!(Atan at real: |_, [x]| crate::slice::Total(atan::<_, { crate::FUSED }>(x))); +pub(crate) fn atan_slice(_: &Atan, _: &crate::recipe::NativeOps<{ crate::FUSED }>, x: [crate::Arg<'_, R>; 1], n: usize, out: &mut [std::mem::MaybeUninit]) -> bool { + each!(x, n, out, |[x]| crate::slice::Total(atan::<_, { crate::FUSED }>(x))) +} /// [`atan2`](crate::atan2), as atan(y/x) in the point's quadrant. It takes the arctangent of the smaller magnitude over the larger, which never overflows. #[inline] diff --git a/core/src/buffer.rs b/core/src/buffer.rs index 8033e9b..4cfeb62 100644 --- a/core/src/buffer.rs +++ b/core/src/buffer.rs @@ -303,8 +303,8 @@ impl DType { /// The float width that a Boolean or integer computes in under rustymath's own policy, NumPy's: the float that [`DType::inexact`] gives it under `f64`. `f64` for any other dtype. pub fn numpy_float(self) -> FloatWidth { match self.inexact(FloatWidth::F64) { Some((DType::Float(w), _)) => w, _ => FloatWidth::F64 } } - /// The dtype that an argument of this dtype computes in, for a function with a kernel at each dtype where `has` holds. That's this dtype where it has a kernel. Otherwise it's the first with a kernel of these: for a Boolean or an integer, integers of `int` or of this dtype's width, whichever is wider, then floats of `float`; for a float or a complex number, the dtype that [`DType::inexact`] gives it; then each wider float, and each complex dtype whose parts are at least as wide; and last, for a float, integers of `int`, which compute it as the whole number it holds (`operation::whole_numbers`). `int` and `float` are the caller's policy for exact numbers. NumPy's is `i64` and the float that `inexact` gives under `f64`. `None` where no dtype on the way has a kernel. - pub fn computes_in(self, int: Width, float: FloatWidth, has: impl Fn(DType) -> bool) -> Option { + /// The dtype that an argument of this dtype computes in, for a function with a kernel at each dtype where `has` holds. That's this dtype where it has a kernel. Otherwise it's the first with a kernel of these: for a Boolean or an integer, integers of `int` or of this dtype's width, whichever is wider, then floats of `float`; for a float or a complex number, the dtype that [`DType::inexact`] gives it; then each wider float, and each complex dtype whose parts are at least as wide; and last, for a float, integers of `int`, which compute it as the whole number it holds (`operation::whole_numbers`). `int` and `float` are the caller's policy for exact numbers. NumPy's is `i64` and the float that `inexact` gives under `f64`. `None` where no dtype on the way has a kernel. Calls to `has` stop at the first success. + pub fn computes_in(self, int: Width, float: FloatWidth, mut has: impl FnMut(DType) -> bool) -> Option { if has(self) { return Some(self); } let (integer, float, whole) = match self { Self::Boolean => (Some(Self::Integer(int)), Self::Float(float), None), @@ -837,22 +837,13 @@ impl Buffer { } } /// The items of `view` at `dtype`. Matching dtypes copy. Otherwise each item converts as [`Source`] reads it, and conversions to a narrower dtype follow these rules. A float converts to an integer by truncating toward zero. A value that doesn't fit its target gives an Overflow error, as do NaN, infinities and anything but exactly 0 or 1 for Booleans. A complex number converts to a real dtype only when its imaginary part is zero, and otherwise gives a Domain error. Characters convert only to characters. - pub fn cast(view: View<'_>, dtype: DType) -> Result { Self::cast_mapped(view, &crate::agreement::Mapping::Linear(1), view.len(), dtype) } + pub fn cast(view: View<'_>, dtype: DType) -> Result { + crate::convert::Conversion { input: view.dtype(), output: dtype, kind: crate::convert::ConversionKind::Cast }.buffer(view, true, crate::tolerant::Tolerance::EXACT) + } /// The items of `view` at `dtype`, where every item keeps its value. Items convert as [`Buffer::cast`] converts them, and a float that becomes a Boolean or an integer must be whole. A float rounds to a float of another width. A Domain error where an item doesn't keep its value. pub fn kept(view: View<'_>, dtype: DType) -> Result { crate::convert::Conversion { input: view.dtype(), output: dtype, kind: crate::convert::ConversionKind::Keep }.buffer(view, true, crate::tolerant::Tolerance::EXACT) } - /// Shared conversion dispatch for compact and mapped logical items. Unselected storage never participates in checked conversion. - pub fn cast_mapped(view: View<'_>, mapping: &crate::agreement::Mapping, n: usize, dtype: DType) -> Result { - use crate::agreement::{Binding, Missing, Read, ReadPlan}; - let mut sources = vec![Binding::full(view)]; - let missing = Missing::prototype(view.dtype(), &mut sources); - let plan = ReadPlan::copy(0, view.len(), mapping.clone(), n, missing); - let mut output = Self::with_capacity(dtype, n); - crate::convert::Conversion { input: view.dtype(), output: dtype, kind: crate::convert::ConversionKind::Cast } - .map(Read { sources: &sources, plan: &plan }, 0..n, true, &mut output)?; - Ok(output) - } } impl Buffer { diff --git a/core/src/buffer/output.rs b/core/src/buffer/output.rs index 9afb20e..de201d3 100644 --- a/core/src/buffer/output.rs +++ b/core/src/buffer/output.rs @@ -101,7 +101,6 @@ fn converted( apply: &(dyn Fn(std::ops::Range, &mut Slots<'_, O>, Option<&mut Scratch>) -> bool + Sync) ) -> Result<(), Error> { let mut values = work.as_deref_mut().map_or_else(Vec::new, |work| std::mem::take(work.result::())); - let mut local = Scratch::default(); let conversion = crate::convert::Conversion { input: O::DTYPE, output: target.dtype(), kind: crate::convert::ConversionKind::Truncate }; let result = (|| { for at in (range.start..range.end).step_by(crate::CHUNK) { @@ -109,7 +108,7 @@ fn converted( values.clear(); Slots::append(&mut values, n, |slots| apply(at..at + n, slots, work.as_deref_mut())); if let Err(error) = conversion.append(O::view(&values), checked, crate::tolerant::Tolerance::EXACT, - work.as_deref_mut().unwrap_or(&mut local), target.reborrow()) { + work.as_deref_mut(), target.reborrow()) { target.fill(range.end - at - n); return Err(error); } diff --git a/core/src/convert.rs b/core/src/convert.rs index 91abcb8..a63ccd9 100644 --- a/core/src/convert.rs +++ b/core/src/convert.rs @@ -30,11 +30,11 @@ macro_rules! conversion { with_width!(to, I => with_float_width!(width, F => $write!(F, I, |x| Ok(saturated(x))))) } (ConversionKind::Cast | ConversionKind::Truncate, DType::Integer(from), DType::Integer(to)) => with_width!(to, I => with_width!(from, S => { - $write!(S, I, |x| { + $write!(S, I, move |x| { if $checked { I::try_from(x).map_err(|_| Error::new(Kind::Overflow, "integer doesn't fit its converted dtype")) } else { Ok(x as I) } }) })), - (ConversionKind::Truncate, input_type, DType::Character) if input_type.is_integral() => $write!(i64, char, |x| { + (ConversionKind::Truncate, input_type, DType::Character) if input_type.is_integral() => $write!(i64, char, move |x| { match u32::try_from(x).ok().and_then(char::from_u32) { Some(x) => Ok(x), None if $checked => Err(Error::new(Kind::Domain, "result is not a Unicode character")), None => Ok(char::FILL) } @@ -90,44 +90,36 @@ impl Conversion { .validity::(&crate::recipe::RangeOps(self.input), crate::domain::Range::of_dtype(self.input)) .is_some_and(|(_, valid)| valid != Some(true))) } - /// Converts contiguous storage through the shared typed conversion rules, without a read plan. + /// Converts contiguous storage through the mapped conversion executor. pub(crate) fn buffer(self, view: crate::buffer::View<'_>, checked: bool, tolerance: crate::tolerant::Tolerance) -> Result { - use crate::buffer::Cast; - macro_rules! write { - ($s:ty, $t:ty, $convert:expr) => { Ok(converted_buffer::<$s, $t>(view, checked, $convert)) }; - (@keep $t:ty) => { Ok(<$t>::buffer((0..view.len()).map(|i| converted_item(crate::buffer::kept_as::<$t>(view, i)).0).collect())) }; - (@cast $t:ty) => {{ - let mut values = Vec::with_capacity(view.len()); - if let Err(error) = <$t>::extend_cast(view, view.len(), checked, &mut values) { crate::error::record(error); values.resize(view.len(), <$t>::FILL); } - Ok(<$t>::buffer(values)) - }}; - (@unsupported) => { Err(crate::error::Error::new(crate::error::Kind::Domain, "unsupported resolved conversion")) }; - } - crate::error::checked(|| conversion!(self, checked, tolerance, write)) + let mut result = crate::buffer::Buffer::with_capacity(self.output, view.len()); + self.append(view, checked, tolerance, None, crate::buffer::MapOutput::Buffer(&mut result))?; + Ok(result) } pub(crate) fn append<'a, 'b>( self, view: crate::buffer::View<'_>, checked: bool, tolerance: crate::tolerant::Tolerance, - scratch: &'a mut crate::buffer::Scratch, output: crate::buffer::MapOutput<'a, 'b> + scratch: Option<&'a mut crate::buffer::Scratch>, output: crate::buffer::MapOutput<'a, 'b> ) -> Result<(), crate::error::Error> { let sources = [crate::agreement::Binding::full(view)]; let plan = crate::agreement::ReadPlan::input(0, view.len()); self.map_in(crate::agreement::Read { sources: &sources, plan: &plan }, 0..view.len(), checked, tolerance, scratch, output) } - pub(crate) fn map( - self, input: crate::agreement::Read<'_>, range: std::ops::Range, checked: bool, output: &mut crate::buffer::Buffer - ) -> Result<(), crate::error::Error> { - if output.dtype() != self.output { *output = crate::buffer::Buffer::with_capacity(self.output, range.len()); } else { output.clear(); } - self.map_in(input, range, checked, crate::tolerant::Tolerance::EXACT, &mut crate::buffer::Scratch::default(), crate::buffer::MapOutput::Buffer(output)) - } - pub fn map_in<'a, 'b>( self, input: crate::agreement::Read<'_>, range: std::ops::Range, checked: bool, tolerance: crate::tolerant::Tolerance, - scratch: &'a mut crate::buffer::Scratch, target: crate::buffer::MapOutput<'a, 'b> + scratch: Option<&'a mut crate::buffer::Scratch>, target: crate::buffer::MapOutput<'a, 'b> ) -> Result<(), crate::error::Error> { - let mut output = crate::buffer::Output { target, result_dtype: self.output, checked, scratch: Some(scratch) }; + let mut output = crate::buffer::Output { target, result_dtype: self.output, checked, scratch }; macro_rules! write { ($s:ty, $t:ty, $convert:expr) => { written::<$s, $t>(input, range, checked, &mut output, $convert) }; - (@cast $t:ty) => { write!($t, $t, Ok) }; + (@cast $t:ty) => {{ + if let crate::buffer::MapOutput::Buffer(buffer) = &mut output.target && let Some(values) = <$t>::vec_mut(buffer) { + let (at, n) = (values.len(), range.len()); + if let Err(error) = input.append(range, checked, values) { + crate::error::record(error); values.truncate(at); values.resize(at + n, <$t>::FILL); + } + Ok(()) + } else { write!($t, $t, Ok) } + }}; (@keep $t:ty) => { kept_mapped::<$t>(self.input, input, range, checked, &mut output) }; (@unsupported) => {{ output.target.fill(range.len()); Err(Error::new(Kind::Domain, "unsupported resolved conversion")) }}; } @@ -140,24 +132,6 @@ fn converted_item(value: Result (value, true), Err(error) => { crate::error::record(error); (T::FILL, false) } } } -fn converted_buffer( - view: crate::buffer::View<'_>, checked: bool, convert: impl Fn(S) -> Result -) -> crate::buffer::Buffer { - let (mut values, mut scratch) = (Vec::with_capacity(view.len()), Vec::new()); - for from in (0..view.len()).step_by(crate::CHUNK) { - let part = view.slice(from..(from + crate::CHUNK).min(view.len())); - let items = if let Some(items) = S::of(part) { items } else { - scratch.clear(); - if let Err(error) = S::extend_cast(part, part.len(), checked, &mut scratch) { - crate::error::record(error); values.resize(values.len() + part.len(), T::FILL); continue; - } - &scratch - }; - values.extend(items.iter().map(|&x| converted_item(convert(x)).0)); - } - T::buffer(values) -} - #[inline] fn written( input: crate::agreement::Read<'_>, range: std::ops::Range, checked: bool, output: &mut crate::buffer::Output<'_, '_>, diff --git a/core/src/distribution.rs b/core/src/distribution.rs index 6287a19..4c5c5b9 100644 --- a/core/src/distribution.rs +++ b/core/src/distribution.rs @@ -1,6 +1,6 @@ -//! Probability distributions. Each has a density, a cumulative distribution, a quantile function and a sampler. The first three are also function objects, `Pdf`, `Cdf` and `Quantile`, which hold the distribution. Their slice forms choose the distribution's formulas once, outside the loop, so the loops of the closed-form distributions vectorize. Their one-item forms run the slice forms on one item. An item of a slice form fails where its result is NaN and its argument isn't, as for a quantile of p outside 0 to 1. A discrete distribution's density is its mass at a whole x, and its quantile function and sampler give whole numbers. +//! Probability distributions. Each has a density, a cumulative distribution, a quantile function and a sampler. The first three are also function objects, `Pdf`, `Cdf` and `Quantile`, which hold the distribution. Their native bulk routines choose the distribution's formulas outside the item loop. Closed-form distributions then vectorize. Scalar evaluation runs the same routine on one item. An item fails where its result is NaN and its argument isn't, as for a quantile of p outside 0 to 1. A discrete distribution's density is its mass at a whole x, and its quantile function and sampler give whole numbers. use crate::{ - Arg, Atan, Erfc, Erfcinv, Float, Ln, Log1p, Unary, betainc, betaincinv, cos, erfcinv, exp, expm1, gammainc, gammaincc, gammaincinv, lgamma, ln, log1p, pow, + Arg, Atan, Erfc, Erfcinv, Float, Ln, Log1p, Native, betainc, betaincinv, cos, erfcinv, exp, expm1, gammainc, gammaincc, gammaincinv, lgamma, ln, log1p, pow, random::Items, sin, slice::{CHUNK, kept, written} }; use rand::{Rng, distr::{Distribution as _, Open01, StandardUniform}}; @@ -50,7 +50,7 @@ macro_rules! each_law { }; } -/// One distribution's formulas. `quantile`, `quantile_map` and `quantile_many` take p from 0 to 1. The `_many` methods are the slice forms at float type `F`, which `Pdf`, `Cdf` and `Quantile` run for one item as for many. By default a slice form runs its `_map` in a loop. A map gives one item's result at `F` and holds its constants, so that the loop computes them once. By default it computes in `f64` and rounds to `F`. A law whose formulas need nothing wider computes in `F` itself, so that its slice forms run at the arguments' width. A formula that takes a function from that function's slice form overrides its `_many` method, through [`composed`]. +/// One distribution's formulas. `quantile`, `quantile_map` and `quantile_many` take p from 0 to 1. The `_many` methods are the native bulk routines at float type `F`, which `Pdf`, `Cdf` and `Quantile` run for one item as for many. By default a native bulk routine runs its `_map` in a loop. A map gives one item's result at `F` and holds its constants, so that the loop computes them once. By default it computes in `f64` and rounds to `F`. A law whose formulas need nothing wider computes in `F` itself, so that its native bulk routines run at the arguments' width. A formula that takes a function from that function's native bulk routine overrides its `_many` method, through [`composed`]. trait Formulas { fn pdf(&self, x: f64) -> f64; fn cdf(&self, x: f64) -> f64; @@ -66,21 +66,22 @@ trait Formulas { #[inline(always)] fn quantile_map(&self) -> impl Fn(F) -> F + Copy { move |p| narrowed(p, |p| self.quantile(p)) } #[inline(always)] - fn pdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Ln: Unary { + fn pdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let pdf = self.pdf_map(); each!([x], n, out, |[x]| pdf(x)) } #[inline(always)] - fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Erfc: Unary, Atan: Unary { + fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let cdf = self.cdf_map(); each!([x], n, out, |[x]| cdf(x)) } #[inline(always)] - fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Erfcinv: Unary, Log1p: Unary, - Ln: Unary { let q = self.quantile_map(); each!([p], n, out, |[p]| q(p)) } + fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { + let q = self.quantile_map(); each!([p], n, out, |[p]| q(p)) + } } -/// `Formulas`' `f64` methods for a law whose generic methods hold its formulas: each runs its slice form on one item. +/// `Formulas`' `f64` methods for a law whose generic methods hold its formulas. Each runs its native bulk routine on one item. macro_rules! through_slices { () => { fn pdf(&self, x: f64) -> f64 { @@ -95,7 +96,7 @@ macro_rules! through_slices { }; } -/// The result that `many`, a slice form, writes into its one slot. +/// The result that `many`, a native bulk routine, writes into its one slot. #[inline(always)] fn one_item(many: impl FnOnce(&mut [MaybeUninit]) -> bool) -> F { let mut out = [MaybeUninit::uninit()]; @@ -107,35 +108,41 @@ fn one_item(many: impl FnOnce(&mut [MaybeUninit]) -> bool) -> F { #[inline(always)] fn narrowed(x: F, f: impl Fn(f64) -> f64) -> F { F::narrow(f(x.wide())) } -/// [`Distribution::pdf`] as a function object. The log-normal distribution's slice form takes ln from [`Ln`]'s slice form, which vectorizes where a loop of whole densities doesn't. The Weibull distribution uses the shared [`Pow`] formula. +/// [`Distribution::pdf`] as a function object. The log-normal distribution's native bulk routine takes ln from [`Ln`]'s native bulk routine, which vectorizes where a loop of whole densities doesn't. The Weibull distribution uses the shared [`Pow`] formula. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Pdf(pub Distribution); -forms!(Pdf at real: |k, [x]| one_item(|out| k.slice([Arg::One(x)], 1, out)), many |k, [x], n, out| each_law!(&k.0 .0, d => d.pdf_many(x, n, out))); +fn pdf_slice(k: &Pdf, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x]: [Arg<'_, F>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + each_law!(&k.0 .0, d => d.pdf_many(x, n, out)) +} -/// [`Distribution::cdf`] as a function object. The normal and log-normal distributions' slice forms take erfc from [`Erfc`]'s slice form, which computes only the formulas that each block of arguments needs. The Cauchy distribution's takes atan from [`Atan`]'s, and the Weibull distribution uses the shared [`Pow`] formula. +/// [`Distribution::cdf`] as a function object. The normal and log-normal distributions' native bulk routines take erfc from [`Erfc`]'s native bulk routine, which computes only the formulas that each block of arguments needs. The Cauchy distribution's takes atan from [`Atan`]'s, and the Weibull distribution uses the shared [`Pow`] formula. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Cdf(pub Distribution); -forms!(Cdf at real: |k, [x]| one_item(|out| k.slice([Arg::One(x)], 1, out)), many |k, [x], n, out| each_law!(&k.0 .0, d => d.cdf_many(x, n, out))); +fn cdf_slice(k: &Cdf, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x]: [Arg<'_, F>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + each_law!(&k.0 .0, d => d.cdf_many(x, n, out)) +} -/// [`Distribution::quantile`] as a function object, through [`quantiles`]. The normal and log-normal distributions' slice forms take erfcinv(2p) from [`Erfcinv`]'s slice form, whose loops don't branch on each item's region. Logarithms use [`Log1p`] and [`Ln`] slice forms; tangent and root use the shared [`Tan`] and [`Root`] formulas. +/// [`Distribution::quantile`] as a function object, through [`quantiles`]. The normal and log-normal distributions' native bulk routines take erfcinv(2p) from [`Erfcinv`]'s native bulk routine, whose loops don't branch on each item's region. Logarithms use [`Log1p`] and [`Ln`] native bulk routines; tangent and root use the shared [`Tan`] and [`Root`] formulas. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Quantile(pub Distribution); -forms!(Quantile at real: |k, [p]| one_item(|out| k.slice([Arg::One(p)], 1, out)), many |k, [p], n, out| quantiles(p, n, out, |p, n, out| each_law!(&k.0 .0, d => d.quantile_many(p, n, out)))); +fn quantile_slice(k: &Quantile, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [p]: [Arg<'_, F>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + quantiles(p, n, out, |p, n, out| each_law!(&k.0 .0, d => d.quantile_many(p, n, out))) +} macro_rules! native_law { - ($( $object:ident $(valid $valid:ident)? ),+) => {$( + ($( $object:ident => $slice:ident $(valid $valid:ident)? ),+) => {$( impl crate::facts::Formula<1> for $object { fn supports(&self, dtype: crate::buffer::DType) -> bool { matches!(dtype, crate::buffer::DType::Float(_)) } const FAILS_NAN: bool = true; $(fn emit_valid + ?Sized>(&self, _: crate::domain::Variant, o: &O, [p]: [O::V; 1], _: O::V) -> Option<(crate::error::Kind, O::B)> { Some((crate::error::Kind::Domain, $valid(o, p))) })? - fn real(&self, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x]: [F; 1]) -> Option { - Some(F::narrow(if F::WIDE { self.one(x.wide()) } else { f64::from(self.one(x.wide() as f32)) })) + fn real(&self, ops: &crate::recipe::NativeOps<{ crate::FUSED }>, [x]: [F; 1]) -> Option { + Some(one_item(|out| $slice(self, ops, [Arg::One(x)], 1, out))) } - fn native_slices(&self) -> Option> { Some(crate::facts::NativeSlices::unary()) } + fn native_slices(&self) -> Option> { Some(crate::facts::NativeSlices::unary($slice, $slice)) } } )+}; } -native_law!(Pdf, Cdf, Quantile valid probability); +native_law!(Pdf => pdf_slice, Cdf => cdf_slice, Quantile => quantile_slice valid probability); fn probability + ?Sized>(o: &O, p: O::V) -> O::B { o.and(o.less_equal(o.number(F::zero()), p), o.less_equal(p, o.number(F::one()))) @@ -151,7 +158,7 @@ fn quantiles(p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit], law: (0..n).filter(|&i| !inside(p.nth(i))).fold(ok, |_, i| { out[i].write(F::nan()); false }) } -/// `finish(x, v)` for each of the `n` items x of `x`, written into `out`. `values` writes a block's values v into the block's slots of `out`, from the block's arguments `inner(x)`, usually through one slice form, as [`through`] gives it. It may overwrite the arguments. A `One` argument computes its result once. Each step runs over up to `CHUNK` items at a time, so the arguments and values between the steps stay in cache, and the arguments' buffer holds no more items than `x`. An item fails where its result is NaN and its argument isn't, as in every slice form. Gives whether every item succeeded. +/// `finish(x, v)` for each of the `n` items x of `x`, written into `out`. `values` writes a block's values v into the block's slots of `out`, from the block's arguments `inner(x)`, usually through one native bulk routine, as [`through`] gives it. It may overwrite the arguments. A `One` argument computes its result once. Each step runs over up to `CHUNK` items at a time, so the arguments and values between the steps stay in cache, and the arguments' buffer holds no more items than `x`. An item fails where its result is NaN and its argument isn't, as in every native bulk routine. Gives whether every item succeeded. fn composed( x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit], inner: impl Fn(F) -> F, values: impl Fn(&mut [F], &mut [MaybeUninit]), finish: impl Fn(F, F) -> F ) -> bool { @@ -179,8 +186,10 @@ fn composed( ok } -/// The values of `g`'s slice form at a block's arguments, for [`composed`]. -fn through>(g: G) -> impl Fn(&mut [F], &mut [MaybeUninit]) { move |a, v| { g.slice([Arg::Many(a)], a.len(), v); } } +/// The values of `g`'s native bulk routine at a block's arguments, for [`composed`]. +fn through(g: K, slice: crate::facts::Slice) -> impl Fn(&mut [F], &mut [MaybeUninit]) { + move |a, v| { slice(&g, &crate::recipe::NativeOps::DEFAULT, [Arg::Many(a)], a.len(), v); } +} /// Each argument of a block raised to `k`, with `k`'s route chosen once (`pow::powers`), for [`composed`]. fn raised(k: F) -> impl Fn(&mut [F], &mut [MaybeUninit]) { @@ -190,11 +199,11 @@ fn raised(k: F) -> impl Fn(&mut [F], &mut [MaybeUninit]) { impl Distribution { /// The density at x, or for a discrete distribution the probability of x. - pub fn pdf(&self, x: f64) -> f64 { Pdf(*self).one(x) } + pub fn pdf(&self, x: f64) -> f64 { one_item(|out| pdf_slice(&Pdf(*self), &crate::recipe::NativeOps::DEFAULT, [Arg::One(x)], 1, out)) } /// The probability of a value at most x. - pub fn cdf(&self, x: f64) -> f64 { Cdf(*self).one(x) } + pub fn cdf(&self, x: f64) -> f64 { one_item(|out| cdf_slice(&Cdf(*self), &crate::recipe::NativeOps::DEFAULT, [Arg::One(x)], 1, out)) } /// The least x whose `cdf` is at least p. p outside 0 to 1 gives NaN. - pub fn quantile(&self, p: f64) -> f64 { Quantile(*self).one(p) } + pub fn quantile(&self, p: f64) -> f64 { one_item(|out| quantile_slice(&Quantile(*self), &crate::recipe::NativeOps::DEFAULT, [Arg::One(p)], 1, out)) } /// The words that one attempt of a draw reads, which each item of `samples` packs (`random::Items`). pub fn words(&self) -> usize { each_law!(&self.0, d => d.words()) } /// Fills `out` with values drawn from the distribution, one for each item of `items`, which packs this distribution's `words`. A normal distribution's item gives two values, from one pair of uniform words. @@ -367,14 +376,14 @@ impl Formulas for Normal { #[inline(always)] fn pdf_map(&self) -> impl Fn(F) -> F + Copy { self.density() } #[inline(always)] - fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Erfc: Unary { + fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let half = F::narrow(0.5); - composed(x, n, out, self.erfc_line(), through(Erfc), #[inline(always)] move |_, e| half * e) + composed(x, n, out, self.erfc_line(), through(Erfc, crate::special::erfc_slice), #[inline(always)] move |_, e| half * e) } #[inline(always)] - fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Erfcinv: Unary { + fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let (two, from) = (F::narrow(2.0), self.from_erfcinv()); - composed(p, n, out, #[inline(always)] move |p| two * p, through(Erfcinv), #[inline(always)] move |_, e| from(e)) + composed(p, n, out, #[inline(always)] move |p| two * p, through(Erfcinv, crate::special::erfcinv_slice), #[inline(always)] move |_, e| from(e)) } fn sample(&self, rng: &mut R) -> f64 { self.mean + self.sd * standard_normal(rng) } /// An item, a pair of uniform words, gives two values, from the cosine and the sine of one angle. @@ -447,22 +456,22 @@ impl Formulas for Cauchy { let (location, rscale, height) = (F::narrow(self.location), F::narrow(self.scale.recip()), F::narrow((PI * self.scale).recip())); move |x| { let z = (x - location) * rscale; height / (F::one() + z * z) } } - /// The arctangent, from `Atan`'s slice form, is of z near the centre and of -1/z beyond |z| = 1. From that arctangent a, the cdf is 1/2 + a/π near the centre, a/π below -1 and 1 + a/π above 1. There 1/2 + atan(z)/π would lose the lower tail's digits. + /// The arctangent, from `Atan`'s native bulk routine, is of z near the centre and of -1/z beyond |z| = 1. From that arctangent a, the cdf is 1/2 + a/π near the centre, a/π below -1 and 1 + a/π above 1. There 1/2 + atan(z)/π would lose the lower tail's digits. #[inline(always)] - fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Atan: Unary { + fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let (location, rscale, rpi, half) = (F::narrow(self.location), F::narrow(self.scale.recip()), F::narrow(PI.recip()), F::narrow(0.5)); composed(x, n, out, move |x| { let z = (x - location) * rscale; let (num, den) = if z.abs() <= F::one() { (z, F::one()) } else { (-F::one(), z) }; num / den - }, through(Atan), move |x, a| { + }, through(Atan, crate::arc::atan_slice), move |x, a| { let (z, a) = ((x - location) * rscale, a * rpi); if z.abs() <= F::one() { half + a } else if z < F::zero() { a } else { F::one() + a } }) } /// The tangent, from the shared `Tan` formula, is of π·a, where a is p below p = 1/4, 1 - p above 3/4, and p - 1/2 between, so its argument stays within ±π/4. Near 0 and 1 the quantile is location ∓ scale/tan(π·a), which keeps p's digits there. #[inline(always)] - fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { + fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let (location, scale, pi, low, high, half) = (F::narrow(self.location), F::narrow(self.scale), F::narrow(PI), F::narrow(0.25), F::narrow(0.75), F::narrow(0.5)); composed(p, n, out, move |p| { @@ -494,9 +503,9 @@ impl Formulas for Exponential { move |x| { let c = -expm1(-rate * x); if x <= F::zero() { F::zero() } else { c } } } #[inline(always)] - fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Log1p: Unary { + fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let mean = F::narrow(self.rate.recip()); - composed(p, n, out, |p| -p, through(Log1p), move |_, l| -l * mean) + composed(p, n, out, |p| -p, through(Log1p, crate::ln::log1p_slice), move |_, l| -l * mean) } fn sample(&self, rng: &mut R) -> f64 { -ln(open(rng)) / self.rate } } @@ -575,12 +584,12 @@ impl Formulas for Laplace { } } #[inline(always)] - fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Ln: Unary { + fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let (location, scale, half) = (F::narrow(self.location), F::narrow(self.scale), F::narrow(0.5)); composed(p, n, out, move |p| { let r = if p < half { p } else { F::one() - p }; r + r - }, through(Ln), move |p, l| if p < half { location + scale * l } else { location - scale * l }) + }, through(Ln, crate::ln::ln_slice), move |p, l| if p < half { location + scale * l } else { location - scale * l }) } } @@ -610,9 +619,9 @@ impl Formulas for Logistic { } } #[inline(always)] - fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Ln: Unary { + fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let (location, scale) = (F::narrow(self.location), F::narrow(self.scale)); - composed(p, n, out, |p| p / (F::one() - p), through(Ln), move |_, l| location + scale * l) + composed(p, n, out, |p| p / (F::one() - p), through(Ln, crate::ln::ln_slice), move |_, l| location + scale * l) } } @@ -628,23 +637,23 @@ impl Formulas for LogNormal { through_slices!(); /// A nonpositive or NaN x goes to `Ln` as 1. Its density is then 0, or NaN for NaN. #[inline(always)] - fn pdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Ln: Unary { + fn pdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let pdf = self.normal().density(); - composed(x, n, out, #[inline(always)] |x| if x > F::zero() { x } else { F::one() }, through(Ln), #[inline(always)] move |x, l| { + composed(x, n, out, #[inline(always)] |x| if x > F::zero() { x } else { F::one() }, through(Ln, crate::ln::ln_slice), #[inline(always)] move |x, l| { let p = pdf(l) / x; if x <= F::zero() { F::zero() } else { p } }) } #[inline(always)] - fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Erfc: Unary { + fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let (line, half) = (self.normal().erfc_line(), F::narrow(0.5)); composed( - x, n, out, #[inline(always)] move |x| line(ln(x)), through(Erfc), #[inline(always)] move |x, e| if x <= F::zero() { F::zero() } else { half * e }) + x, n, out, #[inline(always)] move |x| line(ln(x)), through(Erfc, crate::special::erfc_slice), #[inline(always)] move |x, e| if x <= F::zero() { F::zero() } else { half * e }) } #[inline(always)] - fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Erfcinv: Unary { + fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let (two, from) = (F::narrow(2.0), self.normal().from_erfcinv()); - composed(p, n, out, #[inline(always)] move |p| two * p, through(Erfcinv), #[inline(always)] move |_, e| exp(from(e))) + composed(p, n, out, #[inline(always)] move |p| two * p, through(Erfcinv, crate::special::erfcinv_slice), #[inline(always)] move |_, e| exp(from(e))) } fn words(&self) -> usize { 2 } fn sample(&self, rng: &mut R) -> f64 { exp(self.normal().sample(rng)) } @@ -682,7 +691,7 @@ impl Formulas for Weibull { through_slices!(); /// shape/scale·z^(shape-1)·e^(-z^shape) for z = x/scale. The shared `Pow` formula gives z^(shape-1), from a base of 1 for a negative x. The density is 0 where e^(-z^shape) underflows, even where z^(shape-1) overflows, and for a negative or infinite x. #[inline(always)] - fn pdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { + fn pdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let (rscale, height) = (F::narrow(self.scale.recip()), F::narrow(self.shape / self.scale)); composed(x, n, out, move |x| if x < F::zero() { F::one() } else { x * rscale }, raised(F::narrow(self.shape - 1.0)), move |x, w| { let e = exp(-times(x * rscale, w)); @@ -692,19 +701,19 @@ impl Formulas for Weibull { } /// 1 - e^(-z^shape) for z = x/scale. The shared `Pow` formula gives z^shape, from a base of 1 for a nonpositive x, whose cdf is 0. #[inline(always)] - fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { + fn cdf_many(&self, x: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let rscale = F::narrow(self.scale.recip()); composed(x, n, out, move |x| if x <= F::zero() { F::one() } else { x * rscale }, raised(F::narrow(self.shape)), |x, w| { let c = -expm1(-w); if x <= F::zero() { F::zero() } else { c } }) } - /// scale·(-log1p(-p))^(1/shape), from `Log1p` slice form and the shared `Root` formula. `Root` takes shape itself as its degree, where `Pow` would need 1/shape rounded. + /// scale·(-log1p(-p))^(1/shape), from `Log1p`'s native bulk routine and the shared `Root` formula. `Root` takes shape itself as its degree, where `Pow` would need 1/shape rounded. #[inline(always)] - fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool where Log1p: Unary { + fn quantile_many(&self, p: Arg<'_, F>, n: usize, out: &mut [MaybeUninit]) -> bool { let (degree, scale) = (F::narrow(self.shape), F::narrow(self.scale)); composed(p, n, out, |p| -p, move |a, v| { - Log1p.slice([Arg::Many(a)], a.len(), v); + crate::ln::log1p_slice(&Log1p, &crate::recipe::NativeOps::DEFAULT, [Arg::Many(a)], a.len(), v); for (x, &l) in a.iter_mut().zip(unsafe { written(v) }.iter()) { *x = -l; } crate::pow::roots(degree, a, v); }, move |_, w| scale * w) diff --git a/core/src/facts.rs b/core/src/facts.rs index ad8f561..1595bd2 100644 --- a/core/src/facts.rs +++ b/core/src/facts.rs @@ -6,7 +6,7 @@ use crate::pow::parity; use crate::recipe::{NativeOps, Ops}; use crate::reduce::{Maximum, Minimum, Product, Reduce, Reduction}; use crate::tolerant::{Residue, Tolerance, Tolerant, TolerantComplex, lcm_from_gcd}; -use crate::{Arg, Betainc, BetaincAt, Binary, FUSED, Flipped, Float, Native, Unary}; +use crate::{Arg, Betainc, BetaincAt, FUSED, Flipped, Float, Native}; use num_complex::Complex; use num_traits::{One, Signed, Zero}; use std::f64::consts::PI; @@ -186,40 +186,30 @@ impl Op { } /// A formula's native slice bodies. Real and integer variants run the body at their dtype (f32 for f16). Complex variants keep the real-input shortcut and final validity check at the map boundary; an absent body uses the scalar formula. pub struct NativeSlices { - pub(crate) single: Option>, pub(crate) double: Option>, pub(crate) byte: Option>, - pub(crate) short: Option>, pub(crate) int: Option>, pub(crate) long: Option>, - pub(crate) complex_single: Option>, pub(crate) complex_double: Option>, - pub(crate) real_complex_single: Option>, pub(crate) real_complex_double: Option> + pub single: Option>, pub double: Option>, pub byte: Option>, + pub short: Option>, pub int: Option>, pub long: Option>, + pub complex_single: Option>, pub complex_double: Option>, + pub real_complex_single: Option>, pub real_complex_double: Option>, + pub complex_test_single: Option>, pub complex_test_double: Option> } /// A slice form writes all `n` results at `U` from arguments at `T`. The caller reports errors from its verdict `R`: `bool` for real/integer forms, `()` for complex forms. pub type Slice = fn(&K, &NativeOps<{ FUSED }>, [crate::Arg<'_, T>; N], usize, &mut [std::mem::MaybeUninit]) -> R; impl NativeSlices { - pub fn single(&self) -> Option> { self.single } - pub fn double(&self) -> Option> { self.double } - pub fn byte(&self) -> Option> { self.byte } - pub fn short(&self) -> Option> { self.short } - pub fn int(&self) -> Option> { self.int } - pub fn long(&self) -> Option> { self.long } - pub fn complex_single(&self) -> Option> { self.complex_single } - pub fn complex_double(&self) -> Option> { self.complex_double } - pub fn real_complex_single(&self) -> Option> { self.real_complex_single } - pub fn real_complex_double(&self) -> Option> { self.real_complex_double } - const NONE: Self = Self { single: None, double: None, byte: None, short: None, int: None, long: None, complex_single: None, complex_double: None, - real_complex_single: None, real_complex_double: None }; + real_complex_single: None, real_complex_double: None, complex_test_single: None, complex_test_double: None }; fn with_complex(self, single: Slice, double: Slice) -> Self { Self { complex_single: Some(single), complex_double: Some(double), ..self } } fn with_real_complex(self, single: Slice, double: Slice) -> Self { Self { real_complex_single: Some(single), real_complex_double: Some(double), ..self } } -} -impl NativeSlices where K: Unary + Unary { - pub(crate) fn unary() -> Self { - Self { single: Some(|k, _, x, n, out| Unary::::slice(k, x, n, out)), double: Some(|k, _, x, n, out| Unary::::slice(k, x, n, out)), - ..Self::NONE } + fn with_complex_tests(self, single: Slice, double: Slice) -> Self { + Self { complex_test_single: Some(single), complex_test_double: Some(double), ..self } } } +impl NativeSlices { + pub(crate) fn unary(single: Slice, double: Slice) -> Self { Self { single: Some(single), double: Some(double), ..Self::NONE } } +} impl NativeSlices { fn power() -> Self { Self { single: Some(|k, _, x, n, out| crate::pow::power_slice::(k, x, n, out)), @@ -234,10 +224,9 @@ impl NativeSlices { } impl NativeSlices { fn residue() -> Self { - Self { single: Some(|k, o, x, n, out| Binary::::slice(&Residue(native_tolerance(o.tolerance, k.0)), x, n, out)), - double: Some(|k, o, x, n, out| Binary::::slice(&Residue(native_tolerance(o.tolerance, k.0)), x, n, out)), - byte: Some(|k, _, x, n, out| Binary::::slice(k, x, n, out)), short: Some(|k, _, x, n, out| Binary::::slice(k, x, n, out)), - int: Some(|k, _, x, n, out| Binary::::slice(k, x, n, out)), long: Some(|k, _, x, n, out| Binary::::slice(k, x, n, out)), ..Self::NONE } + Self { single: Some(crate::tolerant::residue_slice), double: Some(crate::tolerant::residue_slice), + byte: Some(crate::integer::residue_byte), short: Some(crate::integer::residue_short), + int: Some(crate::integer::residue_int), long: Some(crate::integer::residue_long), ..Self::NONE } } } @@ -373,10 +362,11 @@ macro_rules! objects { $(native_integer |$nik:pat_param, $nio:tt, $nix:tt| $native_integer:expr;)? $(native_complex |$nfk:pat_param, $nfo:tt, $nfx:tt| $native_complex:expr;)? $(complex_slices $complex_slices:ident;)? + $(complex_tests $complex_tests:ident;)? $(real_complex_slices $real_complex_slices:ident;)? $(native_real |$nrk:pat_param, $nro:tt, $nrx:tt| $native_real:expr;)? $(real_complex |$rck:pat_param, $rco:tt, $rcx:tt| $real_complex:expr;)? - $(slices $slices:ident;)? + $(slices $slices:expr;)? $(exact |$ek:pat_param, $ex:tt| $exact:expr;)? $(checks $checks:ident;)? $(valid |$vk:pat_param, $vv:tt, $vo:tt, $vx:tt, $vr:tt| $valid:expr;)? @@ -418,9 +408,10 @@ macro_rules! objects { $(#[inline(always)] fn real(&self, $nro: &crate::recipe::NativeOps<{ crate::FUSED }>, $nrx: [F; $n]) -> Option { let $nrk = self; $native_real })? $(const REAL_COMPLEX: bool = true; #[inline(always)] fn real_complex(&self, $rco: &crate::recipe::NativeOps<{ crate::FUSED }>, $rcx: [F; $n]) -> Option> { let $rck = self; $real_complex })? fn native_slices(&self) -> Option> { - None $(.or_else(|| Some(NativeSlices::$slices())))? + None $(.or_else(|| Some($slices)))? $(.or_else(|| Some(NativeSlices::NONE)).map(|s| s.with_complex(crate::complex::bulk::$complex_slices, crate::complex::bulk::$complex_slices)))? $(.or_else(|| Some(NativeSlices::NONE)).map(|s| s.with_real_complex(crate::complex::bulk::$real_complex_slices, crate::complex::bulk::$real_complex_slices)))? + $(.or_else(|| Some(NativeSlices::NONE)).map(|s| s.with_complex_tests(crate::tolerant::$complex_tests, crate::tolerant::$complex_tests)))? } } )+}; @@ -445,8 +436,8 @@ objects!(1: crate::Exp { range [x] => increasing(x, f64::exp).meet(Range::NONNEGATIVE); formula |_, o, [x]| Some(crate::exp::formula(o, x)); slope |_, _, r| r; native_complex |_, _, [x]| Some(crate::complex::exp::(x)); complex_slices exp; } crate::Expm1 { range [x] => increasing(x, f64::exp_m1); formula |_, o, [x]| Some(crate::exp::minus_one_formula(o, x, o.number(F::zero()))); slope |o, _, r| o.add(r, o.number(F::one())); native_complex |_, _, [x]| Some(crate::complex::exp::(x) - Complex::new(F::one(), F::zero())); } crate::Ln { range [x] => increasing(x, |v| if v < 0.0 { f64::NEG_INFINITY } else { v.ln() }); domain _, [x] => inside([x], [Range::NONNEGATIVE]); formula |_, o, [x]| Some(crate::ln::formula(o, x)); slope |o, x, _| o.reciprocal(x); native_complex |_, _, [x]| Some(crate::complex::ln::(x)); complex_slices ln; real_complex |_, o, [x]| Some(num_complex::Complex::new(crate::ln::formula(o, o.magnitude(x)), if x < F::zero() { F::narrow(PI) } else { F::zero() })); checks nan; } - crate::Log2 { range [x] => increasing(x, f64::log2); domain _, [x] => inside([x], [Range::NONNEGATIVE]); formula |_, o, [x]| Some(crate::ln::log_with(o, x, crate::ln::base2(o))); slope |o, x, _| o.divide(o.number(F::narrow(std::f64::consts::LOG2_E)), x); slices unary; checks nan; } - crate::Log10 { range [x] => increasing(x, f64::log10); domain _, [x] => inside([x], [Range::NONNEGATIVE]); formula |_, o, [x]| Some(crate::ln::log_with(o, x, crate::ln::base10(o))); slope |o, x, _| o.divide(o.number(F::narrow(std::f64::consts::LOG10_E)), x); slices unary; checks nan; } + crate::Log2 { range [x] => increasing(x, f64::log2); domain _, [x] => inside([x], [Range::NONNEGATIVE]); formula |_, o, [x]| Some(crate::ln::log_with(o, x, crate::ln::base2(o))); slope |o, x, _| o.divide(o.number(F::narrow(std::f64::consts::LOG2_E)), x); slices NativeSlices::unary(crate::ln::log2_slice, crate::ln::log2_slice); checks nan; } + crate::Log10 { range [x] => increasing(x, f64::log10); domain _, [x] => inside([x], [Range::NONNEGATIVE]); formula |_, o, [x]| Some(crate::ln::log_with(o, x, crate::ln::base10(o))); slope |o, x, _| o.divide(o.number(F::narrow(std::f64::consts::LOG10_E)), x); slices NativeSlices::unary(crate::ln::log10_slice, crate::ln::log10_slice); checks nan; } crate::RealLn { range [x] => increasing(x.magnitude(), f64::ln); output real; formula |_, o, [x]| Some(crate::ln::formula(o, o.magnitude(x))); slope |o, x, _| o.reciprocal(x); native_complex |_, _, [x]| Some(crate::complex::ln::(x)); } crate::RealSqrt { range [x] => exact(x.lo.max(0.0).sqrt(), x.hi.max(0.0).sqrt()); output real; @@ -558,8 +549,8 @@ objects!(1: native_real |_, _, [x]| Some(crate::special::gamma::(x)); checks nan; } crate::Sign { formula |_, o, [x]| Some(crate::recipe::sign(o, x)); slope |o, _, _| o.number(F::zero()); integers arithmetic; native_complex |_, _, [x]| Some(crate::complex::sign::(x)); exact |_, [x]| Ok(Some(x.signum())); } - crate::tolerant::Floor { formula |k, o, [x]| Some(crate::with_tolerance!(k.0, |t| t.floor_in(o, x))); slope |o, _, _| o.number(F::zero()); integers formula; native_complex |k, o, [x]| Some(crate::tolerant::TolerantComplex::complex_floor(native_tolerance(o.tolerance, k.0), x)); exact |_, [x]| Ok(Some(x.floor())); } - crate::tolerant::Ceil { formula |k, o, [x]| Some(crate::with_tolerance!(k.0, |t| t.ceil_in(o, x))); slope |o, _, _| o.number(F::zero()); integers formula; native_complex |k, o, [x]| Some(crate::tolerant::TolerantComplex::complex_ceil(native_tolerance(o.tolerance, k.0), x)); exact |_, [x]| Ok(Some(x.ceil())); } + crate::tolerant::Floor { formula |k, o, [x]| Some(crate::with_tolerance!(k.0, |t| t.floor_in(o, x))); slope |o, _, _| o.number(F::zero()); integers formula; native_complex |k, o, [x]| Some(crate::tolerant::TolerantComplex::complex_floor(k.0.overridden(o.tolerance), x)); exact |_, [x]| Ok(Some(x.floor())); } + crate::tolerant::Ceil { formula |k, o, [x]| Some(crate::with_tolerance!(k.0, |t| t.ceil_in(o, x))); slope |o, _, _| o.number(F::zero()); integers formula; native_complex |k, o, [x]| Some(crate::tolerant::TolerantComplex::complex_ceil(k.0.overridden(o.tolerance), x)); exact |_, [x]| Ok(Some(x.ceil())); } crate::Not { table 1; test |_, o, [x]| Some(o.equal(x, o.number(F::zero()))); integers formula; valid |_, v, o, x, _| truth_valid(v, o, x); } crate::Cis { output complex; native_complex |_, _, [x]| Some(crate::complex::exp::(num_complex::Complex::new(-x.im, x.re))); complex_slices cis; real_complex_slices real_cis; real_complex |_, o, [x]| (!crate::trig::Kind::Cos.predicate(o, x)).then(|| { let (c, s) = crate::complex::quarter(x).unwrap_or_else(|| (crate::trig::formula(o, x, crate::trig::Kind::Cos), crate::trig::formula(o, x, crate::trig::Kind::Sin))); num_complex::Complex::new(c, s) }); } crate::LambertW { domain dtype, [x] => crate::lambert::real_domain(dtype, x); native_complex |_, _, [x]| Some(crate::lambert::complex(x)); native_real |_, _, [x]| Some(crate::lambert::lambertw::(x)); checks nan; } @@ -568,10 +559,10 @@ objects!(1: crate::Digamma { slope |o, x, _| o.polygamma(o.number(F::one()), x); native_real |_, _, [x]| Some(crate::special::digamma::(x)); } // 2/√π e^-x², and its negative for erfc. crate::Erf { slope |o, x, _| o.multiply(o.number(F::narrow(std::f64::consts::FRAC_2_SQRT_PI)), o.exp(o.negate(o.multiply(x, x)))); native_real |_, _, [x]| Some(crate::special::erf::(x)); } - crate::Erfc { slope |o, x, _| o.multiply(o.number(F::narrow(-std::f64::consts::FRAC_2_SQRT_PI)), o.exp(o.negate(o.multiply(x, x)))); native_real |_, _, [x]| Some(crate::special::erfc::(x)); slices unary; } + crate::Erfc { slope |o, x, _| o.multiply(o.number(F::narrow(-std::f64::consts::FRAC_2_SQRT_PI)), o.exp(o.negate(o.multiply(x, x)))); native_real |_, _, [x]| Some(crate::special::erfc::(x)); slices NativeSlices::unary(crate::special::erfc_slice, crate::special::erfc_slice); } // √π/2 e^(erfinv x)², from the result. - crate::Erfinv { slope |o, _, r| o.multiply(o.number(F::narrow(std::f64::consts::FRAC_2_SQRT_PI.recip())), o.exp(o.multiply(r, r))); native_real |_, _, [x]| Some(crate::special::erfinv::(x)); slices unary; checks nan; } - crate::Erfcinv { slope |o, _, r| o.multiply(o.number(F::narrow(-std::f64::consts::FRAC_2_SQRT_PI.recip())), o.exp(o.multiply(r, r))); native_real |_, _, [x]| Some(crate::special::erfcinv::(x)); slices unary; checks nan; } + crate::Erfinv { slope |o, _, r| o.multiply(o.number(F::narrow(std::f64::consts::FRAC_2_SQRT_PI.recip())), o.exp(o.multiply(r, r))); native_real |_, _, [x]| Some(crate::special::erfinv::(x)); slices NativeSlices::unary(crate::special::erfinv_slice, crate::special::erfinv_slice); checks nan; } + crate::Erfcinv { slope |o, _, r| o.multiply(o.number(F::narrow(-std::f64::consts::FRAC_2_SQRT_PI.recip())), o.exp(o.multiply(r, r))); native_real |_, _, [x]| Some(crate::special::erfcinv::(x)); slices NativeSlices::unary(crate::special::erfcinv_slice, crate::special::erfcinv_slice); checks nan; } ); objects!(2: crate::Add { range [x, y] => exact(x.lo + y.lo, x.hi + y.hi); identity [Some(0.0), Some(0.0)]; fold Sum; truth Sum; characters Shift; formula |_, o, [x, y]| Some(o.add(x, y)); slopes |o, _, _, _| [o.number(F::one()), o.number(F::one())]; integers arithmetic; native_complex |_, _, [x, y]| Some(x+y); exact |_, [x, y]| Ok(Some(x + y)); valid |_, v, o, [x, y], _| arithmetic_valid(v, o, crate::recipe::Binary::Add, x, y); } @@ -606,14 +597,14 @@ objects!(2: range [x, y] => powered(x, y); domain _, [x, y] => if whole(y) { Some([x, y]) } else { inside([x, y], [Range::NONNEGATIVE, Range::ANY]) }; identity [None, Some(1.0)]; formula |_, o, [x, y]| Some(crate::pow::formula(o, x, y)); shared power; native_integer |_, o, [x, y]| Some(I::low(crate::integer::native_power(o.variant, x.into(), y.into()))); - native_complex |_, _, [x, y]| Some(crate::complex::pow::(x, y)); complex_slices pow; slices power; exact |_, [x, y]| crate::exact::power(x, y); checks nan; + native_complex |_, _, [x, y]| Some(crate::complex::pow::(x, y)); complex_slices pow; slices NativeSlices::power(); exact |_, [x, y]| crate::exact::power(x, y); checks nan; valid |_, v, o, x, r| power_valid(v, o, x, r); } crate::RealPow { range [x, _] => signs([x]); identity [None, Some(1.0)]; output real; formula |_, o, [x, y]| Some(crate::pow::real_formula(o, x, y)); shared real_power; native_complex |_, _, [x, y]| Some(crate::complex::pow::(x, y)); } // The yth root of x, x^(1/y). - crate::Root { range [_, y] => signs([y]); domain _, [n, y] => root_domain(n, y); formula |_, o, [n, y]| crate::pow::root_formula(o, n, y); slopes |o, x, y, r| [o.negate(o.divide(o.multiply(r, o.ln(y)), o.multiply(x, x))), o.divide(r, o.multiply(x, y))]; native_complex |_, _, [n, y]| Some(crate::complex::nth_root::(n, y)); native_real |_, _, [n, y]| Some(crate::pow::root::(n, y)); slices root; exact |_, [n, y]| { if num_traits::Zero::is_zero(n) { return Err(crate::error::Error::new(crate::error::Kind::Domain, "root degree is zero")); } crate::exact::power(y, &n.recip()) }; checks nan; } - crate::Logb { domain _, [x, y] => inside([x, y], [Range::NONNEGATIVE; 2]); shared logb; native_complex |_, _, [b, x]| Some(crate::complex::divide(crate::complex::ln::(x), crate::complex::ln::(b))); native_real |_, _, [b, x]| Some(crate::ln::logb::(b, x)); slices logarithm; valid |_, _, o, x, r| logarithm_valid(o, x, r); } - crate::MathishLogb { domain _, [x, y] => inside([x, y], [Range::NONNEGATIVE; 2]); shared logb; native_complex |_, _, [b, x]| { let one = num_complex::Complex::new(F::one(), F::zero()); Some(if b == one && x == one { num_complex::Complex::new(F::zero(), F::zero()) } else { crate::complex::divide(crate::complex::ln::(x), crate::complex::ln::(b)) }) }; native_real |_, _, [b, x]| Some(crate::ln::logb_mode::(b, x)); slices mathish_logarithm; valid |_, _, o, x, r| logarithm_valid(o, x, r); } + crate::Root { range [_, y] => signs([y]); domain _, [n, y] => root_domain(n, y); formula |_, o, [n, y]| crate::pow::root_formula(o, n, y); slopes |o, x, y, r| [o.negate(o.divide(o.multiply(r, o.ln(y)), o.multiply(x, x))), o.divide(r, o.multiply(x, y))]; native_complex |_, _, [n, y]| Some(crate::complex::nth_root::(n, y)); native_real |_, _, [n, y]| Some(crate::pow::root::(n, y)); slices NativeSlices::root(); exact |_, [n, y]| { if num_traits::Zero::is_zero(n) { return Err(crate::error::Error::new(crate::error::Kind::Domain, "root degree is zero")); } crate::exact::power(y, &n.recip()) }; checks nan; } + crate::Logb { domain _, [x, y] => inside([x, y], [Range::NONNEGATIVE; 2]); shared logb; native_complex |_, _, [b, x]| Some(crate::complex::divide(crate::complex::ln::(x), crate::complex::ln::(b))); native_real |_, _, [b, x]| Some(crate::ln::logb::(b, x)); slices NativeSlices::logarithm(); valid |_, _, o, x, r| logarithm_valid(o, x, r); } + crate::MathishLogb { domain _, [x, y] => inside([x, y], [Range::NONNEGATIVE; 2]); shared logb; native_complex |_, _, [b, x]| { let one = num_complex::Complex::new(F::one(), F::zero()); Some(if b == one && x == one { num_complex::Complex::new(F::zero(), F::zero()) } else { crate::complex::divide(crate::complex::ln::(x), crate::complex::ln::(b)) }) }; native_real |_, _, [b, x]| Some(crate::ln::logb_mode::(b, x)); slices NativeSlices::mathish_logarithm(); valid |_, _, o, x, r| logarithm_valid(o, x, r); } // `y` modulo `x`. Its slope in `x` is minus the quotient's floor, (y - r) ÷ x from the result, and 0 where `x` is 0, where the residue is `y` itself. Integers take a positive power of two that every item shares as a mask (`masked`). Integer kernels run its slice forms (`integer.rs`). Where a whole block shares one divisor, those take a mask or a reciprocal in place of a division. crate::tolerant::Residue { identity [Some(0.0), None]; @@ -623,12 +614,12 @@ objects!(2: o.negate(o.multiply(o.one_where(o.not(at_zero)), o.divide(o.subtract(y, r), o.add(x, o.one_where(at_zero))))) }, o.number(F::one())]; native_integer |_, _, [x, y]| Some(I::low(crate::integer::residue(x.into(), y.into()))); - native_complex |k, o, [x, y]| Some(crate::tolerant::TolerantComplex::complex_residue(native_tolerance(o.tolerance, k.0), x, y)); + native_complex |k, o, [x, y]| Some(crate::tolerant::TolerantComplex::complex_residue(k.0.overridden(o.tolerance), x, y)); native_real |k, o, [x, y]| { - let t = native_tolerance(o.tolerance, k.0); + let t = k.0.overridden(o.tolerance); Some(if t.is_exact() { crate::tolerant::exact_residue(x, y) } else { crate::tolerant::Tolerant::residue(t, x, y) }) }; - slices residue; + slices NativeSlices::residue(); exact |_, [x, y]| Ok(Some(crate::exact::residue(x, y))); } // C(n, k), with n first. @@ -669,15 +660,15 @@ objects!(2: crate::tolerant::Gcd { identity [Some(0.0), Some(0.0)]; table 14; native_integer |_, o, [x, y]| { let (x, y) = (x.into(), y.into()); Some(o.result(crate::integer::overflowed(o.variant, crate::integer::gcd(x, y), || crate::integer::gcd_wrapped(x, y)))) }; - native_complex |k, o, [x, y]| Some(native_tolerance(o.tolerance, k.0).complex_gcd(x, y)); - native_real |k, o, [x, y]| Some(native_tolerance(o.tolerance, k.0).gcd(x, y)); + native_complex |k, o, [x, y]| Some(k.0.overridden(o.tolerance).complex_gcd(x, y)); + native_real |k, o, [x, y]| Some(k.0.overridden(o.tolerance).gcd(x, y)); exact |_, [x, y]| Ok(Some(crate::exact::gcd(x, y))); checks nan; } crate::tolerant::Lcm { identity [Some(1.0), Some(1.0)]; table 8; native_integer |_, o, [x, y]| { let (x, y) = (x.into(), y.into()); Some(o.result(crate::integer::overflowed(o.variant, crate::integer::lcm(x, y), || crate::integer::lcm_wrapped(x, y)))) }; - native_complex |k, o, [x, y]| Some(lcm_from_gcd(x, y, native_tolerance(o.tolerance, k.0).complex_gcd(x, y))); - native_real |k, o, [x, y]| Some(lcm_from_gcd(x, y, native_tolerance(o.tolerance, k.0).gcd(x, y))); + native_complex |k, o, [x, y]| Some(lcm_from_gcd(x, y, k.0.overridden(o.tolerance).complex_gcd(x, y))); + native_real |k, o, [x, y]| Some(lcm_from_gcd(x, y, k.0.overridden(o.tolerance).gcd(x, y))); exact |_, [x, y]| Ok(Some(crate::exact::lcm(x, y))); checks nan; } crate::Nand { table 7; test |_, o, [x, y]| Some(o.not(o.and(boolean_one(o, x), boolean_one(o, y)))); integers formula; valid |_, v, o, x, _| truth_valid(v, o, x); } @@ -694,13 +685,15 @@ objects!(2: crate::tolerant::Equal { identity [Some(1.0), Some(1.0)]; table 9; characters Ordered; test |k, o, [x, y]| Some(tested(*k, k.0, o, x, y)); against |k| k.0; integers formula; - native_complex |k, o, [x, y]| Some(if native_tolerance(o.tolerance, k.0).complex_equal(x, y) { Complex::one() } else { Complex::zero() }); + native_complex |k, o, [x, y]| Some(if k.0.overridden(o.tolerance).complex_equal(x, y) { Complex::one() } else { Complex::zero() }); + complex_tests complex_equal_slice; exact |_, [x, y]| Ok(Some(crate::exact::boolean(x == y))); } crate::tolerant::NotEqual { identity [Some(0.0), Some(0.0)]; table 6; characters Ordered; test |k, o, [x, y]| Some(tested(*k, k.0, o, x, y)); against |k| k.0; integers formula; - native_complex |k, o, [x, y]| Some(if !native_tolerance(o.tolerance, k.0).complex_equal(x, y) { Complex::one() } else { Complex::zero() }); + native_complex |k, o, [x, y]| Some(if !k.0.overridden(o.tolerance).complex_equal(x, y) { Complex::one() } else { Complex::zero() }); + complex_tests complex_not_equal_slice; exact |_, [x, y]| Ok(Some(crate::exact::boolean(x != y))); } crate::tolerant::Less { identity [Some(0.0), None]; table 2; characters Ordered; test |k, o, [x, y]| Some(tested(*k, k.0, o, x, y)); against |k| k.0; integers formula; exact |_, [x, y]| Ok(Some(crate::exact::boolean(x < y))); } @@ -763,9 +756,6 @@ fn masked + ?Sized>(o: &O, x: O::V, y: Some(o.bit_and(y, o.number(F::narrow(mask as f64)))) } -#[inline] -fn native_tolerance(value: Option, fallback: Tolerance) -> Tolerance { value.and_then(Tolerance::new).unwrap_or(fallback) } - /// The tolerance at which `o` compares at tolerance `t`: the one that `o` sets, otherwise `t` (`Ops::tolerance`). `None` where it isn't known when the recipe is built. pub fn run_tolerance + ?Sized>(o: &O, t: crate::tolerant::Tolerance) -> Option { crate::tolerant::Tolerance::new(o.constant(o.tolerance(F::narrow(t.value())))?.wide()) diff --git a/core/src/float.rs b/core/src/float.rs index 67457d4..9465679 100644 --- a/core/src/float.rs +++ b/core/src/float.rs @@ -1,6 +1,5 @@ //! The float types that rustymath's functions take, and the operations that its algorithms build on. use crate::search::{Key, RadixKey}; -use crate::{parallel::Slots, slice::{Arg, Ternary}}; /// Whether `mla` fuses its multiply and add into one operation with one rounding. It does where the target has a fused multiply-add instruction, as on AArch64 and on x86-64 with FMA. Elsewhere, as in wasm and on baseline x86-64, `mul_add` calls a software `fma`, which is slow and stops the loop around it vectorizing. There `mla` multiplies and adds with two roundings. pub const FUSED: bool = cfg!(any(target_arch = "aarch64", target_feature = "fma")); @@ -31,8 +30,6 @@ pub trait Float: num_traits::Float + sealed::Sealed + crate::slice::KernelItem + fn total_key(self) -> Self::OrderKey; /// The float whose `order_key` or `total_key` is `k`. A total key reads back exactly, apart from a NaN's sign. An order key reads `¯0` back as `0`, and every NaN as one NaN. fn from_order_key(k: Self::OrderKey) -> Self; - /// Writes `x` where `mask`'s lanes aren't zero and `y` where they are into the next `n` slots of `out`, through [`Choose`](crate::Choose), and gives whether every result succeeded. Fused graphs pass a comparison's lanes of 1 and 0 as the mask. - fn choose(x: Arg<'_, Self>, y: Arg<'_, Self>, mask: Arg<'_, Self>, n: usize, out: &mut Slots<'_, Self>) -> bool; } /// A float type that rustymath's algorithms run in: `f32` or `f64`. @@ -149,8 +146,6 @@ macro_rules! native { #[inline] fn narrow(x: f64) -> $t { x as $t } #[inline] - fn choose(x: Arg<'_, $t>, y: Arg<'_, $t>, mask: Arg<'_, $t>, n: usize, out: &mut Slots<'_, $t>) -> bool { crate::Choose.many_into(x, y, mask, n, out) } - #[inline] fn next_up(self) -> $t { $t::next_up(self) } #[inline] fn next_down(self) -> $t { $t::next_down(self) } @@ -222,10 +217,6 @@ impl Float for half::f16 { fn wide(self) -> f64 { self.to_f64() } #[inline] fn narrow(x: f64) -> Self { half::f16::from_f64(x) } - #[inline] - fn choose(x: Arg<'_, Self>, y: Arg<'_, Self>, mask: Arg<'_, Self>, n: usize, out: &mut Slots<'_, Self>) -> bool { - crate::Choose.many_into(x, y, mask, n, out) - } // The float above a positive float has the next bit pattern. The float above a negative float has the previous one. Both zeros step to the least positive subnormal. #[inline] fn next_up(self) -> Self { diff --git a/core/src/integer.rs b/core/src/integer.rs index f8104b5..30cbacf 100644 --- a/core/src/integer.rs +++ b/core/src/integer.rs @@ -1,7 +1,7 @@ -//! Functions of `i64`s: greatest common divisor, least common multiple, binomial coefficients, powers, square roots and residues. Each gives `None` where its result isn't an `i64`. Their function objects are the ones that floats and complex numbers use: [`Gcd`], [`Lcm`] and [`Residue`] from [`tolerant`](crate::tolerant), whose integer forms ignore the tolerance they hold, and [`Binom`] and [`Pow`]. A slice form appends results to a vector and reports whether every item fitted. [`Divisor`] divides many integers by one, with multiplications in place of division instructions. Where every item has the same divisor, `Residue`'s slice form divides through `FloatDivisor`, which multiplies by the divisor's reciprocal in `f64` and vectorizes, for divisors and dividends below 2⁵¹ in magnitude. At types of 16 bits or fewer it multiplies 32-bit integers instead, and it takes `Divisor` for 64-bit dividends of 2⁵¹ or more. +//! Functions of `i64`s: greatest common divisor, least common multiple, binomial coefficients, powers, square roots and residues. Each gives `None` where its result isn't an `i64`. Their function objects are the ones that floats and complex numbers use: [`Gcd`](crate::tolerant::Gcd), [`Lcm`](crate::tolerant::Lcm) and [`Residue`](crate::tolerant::Residue) from [`tolerant`](crate::tolerant), whose integer formulas ignore the tolerance they hold, and [`Binom`](crate::Binom) and [`Pow`](crate::Pow). A native bulk routine writes results and reports whether every item fitted. [`Divisor`] divides many integers by one, with multiplications in place of division instructions. Where every item has the same divisor, `Residue`'s native bulk routine divides through `FloatDivisor`, which multiplies by the divisor's reciprocal in `f64` and vectorizes, for divisors and dividends below 2⁵¹ in magnitude. At types of 16 bits or fewer it multiplies 32-bit integers instead, and it takes `Divisor` for 64-bit dividends of 2⁵¹ or more. //! -//! `Gcd`, `Lcm`, `Binom`, `Pow` and `Residue` also take the narrower integer types of compact buffers, `u8`, `i16` and `i32`, which they compute through these functions. A result that doesn't fit its type fails. `Pow` floors a negative power. -use crate::{Arg, Binom, Pow, convert::{self, WIDE_LIMIT}, float::sealed::Core, tolerant::{Gcd, Lcm, Residue}}; +//! `Gcd`, `Lcm`, `Binom` and `Pow` compute the compact integer types `u8`, `i16` and `i32` through the `i64` helpers, with the variant's checked or wrapping result policy. `Residue` has a native bulk routine at each width. Integer powers floor negative exponents. +use crate::{Arg, convert::{self, WIDE_LIMIT}, float::sealed::Core, tolerant::Residue}; use strength_reduce::StrengthReducedU64; /// The greatest common divisor of x and y, which is never negative. `gcd(0, 0)` is 0. `None` only where it is 2^63, as for `gcd(i64::MIN, 0)`. @@ -211,15 +211,8 @@ fn decoded(bases: &[i64], digits: &[i64], step: impl Fn(i64, i64, i64) -> Option digits.chunks(bases.len()).map(|d| (1..d.len()).try_fold(d[0], |v, k| step(v, bases[k], d[k]))).collect() } -forms!(#[inline] Gcd at i64: |_, [x, y]| gcd(x, y)); -forms!(#[inline] Lcm at i64: |_, [x, y]| lcm(x, y)); -// On Booleans, the greatest common divisor is `or` and the least common multiple is `and`. -forms!(#[inline] Gcd at bool: |_, [x, y]| x | y); -forms!(#[inline] Lcm at bool: |_, [x, y]| x & y); -forms!(#[inline] Binom at i64: |_, [n, k]| binom(n, k)); -forms!(#[inline] Pow at i64: |_, [x, y]| wrapped_power(x, y, i64::MAX)); -forms!( - #[inline] Residue at i64: |_, [x, y]| residue(x, y), many |_, [x, y], n, out| if let Arg::One(d) = x && let Some(m) = mask(d) { +pub(crate) fn residue_long(_: &Residue, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x, y]: [Arg<'_, i64>; 2], n: usize, out: &mut [std::mem::MaybeUninit]) -> bool { + if let Arg::One(d) = x && let Some(m) = mask(d) { each!([y], n, out, |[y]| y & m) } else if let Arg::One(d) = x && let Some(r) = FloatDivisor::new(d) && fit(y, n) { each!([y], n, out, |[y]| r.residue(y)) @@ -228,18 +221,12 @@ forms!( } else { each!([x, y], n, out, |[x, y]| residue(x, y)) } -); +} -/// The narrower types' forms of `Gcd`, `Lcm`, `Binom`, `Pow` and `Residue`, which compute at `i64`. A result that doesn't fit the type fails, as one beyond `i64` does. A power that fails holds the power wrapped to the type. Every residue fits. -macro_rules! narrower { - ($($t:ty),+) => {$( - forms!(#[inline] Gcd at $t: |_, [x, y]| gcd(x.into(), y.into()).and_then(|r| <$t>::try_from(r).ok())); - forms!(#[inline] Lcm at $t: |_, [x, y]| lcm(x.into(), y.into()).and_then(|r| <$t>::try_from(r).ok())); - forms!(#[inline] Binom at $t: |_, [n, k]| binom(n.into(), k.into()).and_then(|r| <$t>::try_from(r).ok())); - forms!(#[inline] Pow at $t: |_, [x, y]| { let (p, exact) = wrapped_power(x.into(), y.into(), <$t>::MAX.into()); (p as $t, exact && <$t>::try_from(p).is_ok()) }); - forms!( - #[inline] Residue at $t: |_, [x, y]| residue(x.into(), y.into()) as $t, - many |_, [x, y], n, out| if let Arg::One(d) = x && let Some(m) = mask(d.into()) { +macro_rules! residues { + ($($name:ident: $t:ty),+) => {$( + pub(crate) fn $name(_: &Residue, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x, y]: [Arg<'_, $t>; 2], n: usize, out: &mut [std::mem::MaybeUninit<$t>]) -> bool { + if let Arg::One(d) = x && let Some(m) = mask(d.into()) { each!([y], n, out, |[y]| y & m as $t) } else if let Arg::One(d) = x && <$t>::BITS <= 16 && let Some(s) = Short::new(d as i16) { each!([y], n, out, |[y]| s.residue(y as i16) as $t) @@ -248,10 +235,10 @@ macro_rules! narrower { } else { each!([x, y], n, out, |[x, y]| residue(x.into(), y.into()) as $t) } - ); + } )+}; } -narrower!(u8, i16, i32); +residues!(residue_byte: u8, residue_short: i16, residue_int: i32); /// x^y at `variant`'s integer width, as [`wrapped_power`] gives it. 0 to a negative power gives the greatest value: i64's under dynamic kinds, which a narrower width reports as an overflow so the call widens, and otherwise the width's. Where `variant` fails on overflow, a power that doesn't fit the width reports an overflow. pub(crate) fn native_power(variant: crate::domain::Variant, x: i64, y: i64) -> i64 { diff --git a/core/src/lib.rs b/core/src/lib.rs index 8fe346f..55c1fe4 100644 --- a/core/src/lib.rs +++ b/core/src/lib.rs @@ -61,7 +61,7 @@ pub use hyperbolic::{RealAcosh, RealAtanh}; pub use ln::{Ln, Log1p, Log2, Log10, Logb, MathishLogb, RealLn}; pub use parts::{Abs, Angle, Choose, Clip, Conj, Imag, Isfinite, Isinf, Isnan, Real, Rint, Round, Scale, Sign, Trunc}; pub use pow::{Exponent, Pow, RealPow, Root}; -pub use slice::{Arg, Binary, CHUNK, Flipped, KernelItem, Results, Ternary, Unary, Wrapped, pairs}; +pub use slice::{Arg, CHUNK, Flipped, KernelItem, Wrapped, pairs}; pub use special::{BetaincAt, Erfc, Erfcinv, Erfinv}; pub use sqrts::{RealSqrt, RealSqrt1mx2, RealSqrtm1mx2, RealSqrtx2m1, Sqrt, Sqrtm1mx2}; pub use trig::{Cos, Sin, Tan}; @@ -71,7 +71,7 @@ pub use {half, num_complex}; /// How many items a function handles between checks for an interrupt ([`interrupt::check`]), which keeps the check out of the loop over items. pub const BETWEEN_CHECKS: usize = 1024; -/// Declares each function at the crate root, where it multiplies and adds as the target's `FUSED` says, and in `fused` and `unfused`, which choose for themselves. Each entry names the module whose function of the same name implements it, generic over `R: Native` and `const FUSED: bool`. Every wrapper is `#[inline(always)]`, because LLVM leaves a large function such as `gamma` out of line in a caller's loop, and the call stops that loop vectorizing. An entry marked `as` also declares a function object of that name, for the facts table's formula. An object with forms, such as [`Pow`] or [`Atan`], is written with `forms!` beside its function. `function_doc!(ln)` gives an entry's doc by its name. The Python package's functions take their docs from it. +/// Declares each scalar function at the crate root and in `fused`/`unfused`, with an optional Formula object. `function_doc!(ln)` gives an entry's docs for Python. macro_rules! functions { ($($(#[doc = $doc:literal])* $name:ident($($x:ident),+) in $module:ident $(as $object:ident)?;)+) => { $( diff --git a/core/src/ln.rs b/core/src/ln.rs index f86f3d8..479bd12 100644 --- a/core/src/ln.rs +++ b/core/src/ln.rs @@ -162,14 +162,31 @@ fn log1p_usual(x: R) -> R { let u = x + R::one(); #[inline(always)] fn usual1p(x: R) -> bool { usual(x + R::one()) & (x != R::zero()) } -usual_object!(/// [`ln`](crate::ln) as a function object. - Ln: ln, usual => #[inline(always)] |x| log_of::<_, { crate::FUSED }>(x.split(), natural())); -usual_object!(/// [`log2`](crate::log2) as a function object. - Log2: log2, usual => #[inline(always)] |x| log_of::<_, { crate::FUSED }>(x.split(), base2(&crate::recipe::Floats))); -usual_object!(/// [`log10`](crate::log10) as a function object. - Log10: log10, usual => #[inline(always)] |x| log_of::<_, { crate::FUSED }>(x.split(), base10(&crate::recipe::Floats))); -usual_object!(/// [`log1p`](crate::log1p) as a function object. - Log1p: log1p, usual1p => log1p_usual::<_, { crate::FUSED }>); +/// [`ln`](crate::ln) as a function object. +#[derive(Clone, Copy, Debug)] +pub struct Ln; +/// [`log2`](crate::log2) as a function object. +#[derive(Clone, Copy, Debug)] +pub struct Log2; +/// [`log10`](crate::log10) as a function object. +#[derive(Clone, Copy, Debug)] +pub struct Log10; +/// [`log1p`](crate::log1p) as a function object. +#[derive(Clone, Copy, Debug)] +pub struct Log1p; + +pub(crate) fn ln_slice(_: &Ln, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x]: [Arg<'_, R>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + each_usual(x, n, out, (usual, #[inline(always)] |x| log_of::<_, { crate::FUSED }>(x.split(), natural())), none(), ln::<_, { crate::FUSED }>) +} +pub(crate) fn log2_slice(_: &Log2, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x]: [Arg<'_, R>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + each_usual(x, n, out, (usual, #[inline(always)] |x| log_of::<_, { crate::FUSED }>(x.split(), base2(&crate::recipe::Floats))), none(), log2::<_, { crate::FUSED }>) +} +pub(crate) fn log10_slice(_: &Log10, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x]: [Arg<'_, R>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + each_usual(x, n, out, (usual, #[inline(always)] |x| log_of::<_, { crate::FUSED }>(x.split(), base10(&crate::recipe::Floats))), none(), log10::<_, { crate::FUSED }>) +} +pub(crate) fn log1p_slice(_: &Log1p, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x]: [Arg<'_, R>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + each_usual(x, n, out, (usual1p, log1p_usual::<_, { crate::FUSED }>), none(), log1p::<_, { crate::FUSED }>) +} /// ln|x|, the real part of ln x, as a function object. #[derive(Clone, Copy, Debug)] pub struct RealLn; diff --git a/core/src/parallel.rs b/core/src/parallel.rs index 05db062..461b3a2 100644 --- a/core/src/parallel.rs +++ b/core/src/parallel.rs @@ -1,5 +1,5 @@ //! Splitting work across threads. [`pieces`] cuts a call's items into pieces of the sizes its caller gives, and [`ranges`] into as many pieces as the time of its first unit justifies. The pieces run on a Rayon pool with Rayon's default number of threads: `RAYON_NUM_THREADS` where it's set, and otherwise one for each core. [`set_threads`] changes that number for the whole process. A call made on one of a pool's threads, such as a function that another call's piece runs, doesn't split again. In the browser build, which has no threads, every call runs its pieces in turn on the calling thread. -use crate::slice::{Arg, CHUNK}; +use crate::slice::Arg; use crate::{error::Error, interrupt::check}; /// The clock that `parallel`'s rounds and the speed tests' timings read: the browser's in a browser build, where `std::time::Instant` panics. #[cfg(not(all(target_arch = "wasm32", target_os = "unknown")))] @@ -258,13 +258,6 @@ impl Slots<'_, B> { } } -/// Runs `write` on each `CHUNK` of `slots` in turn, the results for the items from `start`, and gives whether every call said `true`. A slice form that makes several passes over buffers of its own keeps them in cache this way. -#[inline] -pub(crate) fn in_chunks(start: usize, slots: &mut [MaybeUninit], write: impl Fn(usize, &mut [MaybeUninit]) -> bool) -> bool { - let (n, mut ok) = (slots.len(), true); - for at in (0..n).step_by(CHUNK) { ok &= write(start + at, &mut slots[at..n.min(at + CHUNK)]) } - ok -} #[cfg(not(target_arch = "wasm32"))] mod pool { diff --git a/core/src/parts.rs b/core/src/parts.rs index 48b2b0d..1cea651 100644 --- a/core/src/parts.rs +++ b/core/src/parts.rs @@ -1,6 +1,4 @@ //! The magnitude, sign and angle of a number, as function objects. -use crate::slice::Total; -use num_complex::{Complex32, Complex64}; /// |x|, or the magnitude |z| of a complex z, which is real, as a function object. #[derive(Clone, Copy, Debug)] @@ -38,7 +36,7 @@ objects! { Rint; /// Clipping x between lower and upper bounds in one pass: `minimum(upper, maximum(lower, x))`. A lower bound above the upper one gives the upper. A NaN argument gives NaN, x or a bound, depending on the target. Clip; - /// `x` where `m` is nonzero and `y` where it's zero, item by item: a select, with the mask `m` as lanes of the items' type. + /// Choosing x where the mask is nonzero, otherwise y. Choose; } @@ -73,9 +71,3 @@ pub(crate) fn rounded_wide(x: i128, scale: f64) -> i128 { let q = q + if away > scale || away == scale && q & 1 != 0 { x.signum() } else { 0 }; q * scale } - -forms!(Choose at real: |_, [x, y, m]| Total(if m != 0.0 { x } else { y })); -forms!(Choose at Complex32, Complex64: |_, [x, y, m]| Total(if m.re != 0.0 { x } else { y })); -forms!(Choose at u8, i16, i32, i64: |_, [x, y, m]| if m != 0 { x } else { y }); -forms!(Choose at bool: |_, [x, y, m]| (m & x) | (!m & y)); -forms!(Choose at char: |_, [x, y, m]| if m == '\u{1}' { x } else { y }); diff --git a/core/src/pow.rs b/core/src/pow.rs index 992a3c6..73458ca 100644 --- a/core/src/pow.rs +++ b/core/src/pow.rs @@ -23,7 +23,7 @@ pub(crate) fn root_formula + ?Sized>(o: &O, n: O::V, rooted(o, Degree::of(k.wide()), k < F::zero(), y) } -/// The route [`Root`] takes for a degree n, chosen by n's magnitude. The slice form reads the route once for a constant n. +/// The route [`Root`] takes for a degree n, chosen by n's magnitude. The native bulk routine reads the route once for a constant n. #[derive(Clone, Copy, Debug, PartialEq, Eq)] enum Degree { One, Two, Three, Four, Eight, Half, General } @@ -142,7 +142,7 @@ pub(crate) const MULTIPLIED: u32 = 8; /// The greatest whole exponent that `pow` raises to by squaring and multiplying where the target fuses multiply-add (`FUSED`). Above `MULTIPLIED`, every step but those of bit position 0 runs in double words (`recipe::whole`). Larger exponents take the general route. From 64 on, squaring and multiplying is slower than the general route for some exponents. pub(crate) const COMPENSATED: u32 = 63; -/// How [`Pow`] raises to an exponent y. Its slice form reads it once for a constant y, and its scalar form for each item. The JIT's lowering of a constant power reads it too. +/// How [`Pow`] raises to an exponent y. Its native bulk routine resolves a constant y once; scalar evaluation resolves each call. The JIT's lowering of a constant power uses the same routes. #[derive(Clone, Copy, Debug, PartialEq, Eq)] pub enum Exponent { /// x^m, or x^-m where `negative`, through `recipe::power`. Positive m goes to `MULTIPLIED`, or `COMPENSATED` with fused multiply-add. Negative powers below -1 need fused multiply-add. No result fails. @@ -200,7 +200,7 @@ pub(crate) fn powers(e: Exponent, k: F, xs: &[F], out: &mut [st routes!(Whole 1, Whole 2, Whole 3, Whole 4, Whole 5, Whole 6, Whole 7, Whole 8, Halves 1, Halves 3); } -/// [`Pow`]'s slice form: [`powers`] for a shared exponent, otherwise the general block pass with scalar repair. Returns whether every result is valid: not NaN, or from a NaN argument. +/// [`Pow`]'s native bulk routine: [`powers`] for a shared exponent, otherwise the general block pass with scalar repair. Returns whether every result is valid: not NaN, or from a NaN argument. pub(crate) fn power_slice(_: &Pow, [x, y]: [crate::Arg<'_, F>; 2], n: usize, out: &mut [std::mem::MaybeUninit]) -> bool { if let (crate::Arg::Many(xs), crate::Arg::One(k)) = (x, y) { let e = Exponent::of(k.wide(), crate::FUSED); @@ -222,7 +222,7 @@ pub(crate) fn roots(n: F, ys: &[F], out: &mut [std::mem::MaybeU routes!(One, Two, Three, Four, Eight, Half); } -/// [`Root`]'s slice form: [`roots`] for a shared degree, otherwise the general block pass followed by per-item short-degree recipes. It gives final validity after repair. +/// [`Root`]'s native bulk routine: [`roots`] for a shared degree, otherwise the general block pass followed by per-item short-degree recipes. It gives final validity after repair. pub(crate) fn root_slice(_: &Root, [n, y]: [crate::Arg<'_, F>; 2], len: usize, out: &mut [std::mem::MaybeUninit]) -> bool { if let (crate::Arg::One(k), crate::Arg::Many(ys)) = (n, y) { roots(k, &ys[..len], &mut out[..len]); @@ -247,7 +247,7 @@ pub(crate) fn real_formula + ?Sized>(o: &O, x: O::V, o.select(o.less(x, o.number(F::zero())), real, r) } -/// [`root`](crate::root) as a function object, with a scalar form and a slice form that work as [`Pow`]'s do. A constant degree of 2 or 3 takes `sqrt` or `cbrt`. +/// [`root`](crate::root) as a function object, with a shared recipe and a native bulk routine as [`Pow`] has. A constant degree of 2 or 3 takes `sqrt` or `cbrt`. #[derive(Clone, Copy, Debug)] pub struct Root; diff --git a/core/src/reduce.rs b/core/src/reduce.rs index fe26e51..329a36e 100644 --- a/core/src/reduce.rs +++ b/core/src/reduce.rs @@ -81,16 +81,10 @@ pub fn minimum(x: F, y: F) -> F { F::from_working(Core::minimum(x.work /// [`maximum`] as a function object, which also takes integers and Booleans. The maximum of two Booleans is their `or`. #[derive(Clone, Copy, Debug)] pub struct Maximum; -forms!(#[inline] Maximum at real: |_, [x, y]| crate::slice::Total(maximum(x, y))); -forms!(#[inline] Maximum at u8, i16, i32, i64: |_, [x, y]| x.max(y)); -forms!(#[inline] Maximum at bool: |_, [x, y]| x | y); /// [`minimum`] as a function object, which also takes integers and Booleans. The minimum of two Booleans is their `and`. #[derive(Clone, Copy, Debug)] pub struct Minimum; -forms!(#[inline] Minimum at real: |_, [x, y]| crate::slice::Total(minimum(x, y))); -forms!(#[inline] Minimum at u8, i16, i32, i64: |_, [x, y]| x.min(y)); -forms!(#[inline] Minimum at bool: |_, [x, y]| x & y); /// A reduction by [`Sum`], [`Product`], [`ScaledProduct`], [`Max`], [`Min`], [`And`], [`Or`], [`Xor`], [`Xnor`], [`SumAs`], [`Squares`] or [`Alternated`], implemented inside rustymath once for each item type it takes. `slice` folds a whole slice of consecutive items, in any order, `merge` combines the folds of two neighbouring runs, and `scan` gives the running folds. [`along`], [`fold_into`] and [`scan_along`] fold along an axis whose items aren't consecutive. [`Sum`] and [`Product`] take floats, integers and complex numbers. [`Max`] and [`Min`] take floats, integers and Booleans, [`And`], [`Or`] and [`Xor`] take integers and Booleans, [`Xnor`] takes signed integers, bytes and Booleans, [`SumAs`] takes Booleans and integers, [`Squares`] takes floats and complex numbers, [`ScaledProduct`] takes floats, `Alternated` takes floats, integers and complex numbers, `Alternated` takes floats, and `Alternated` takes complex numbers. An `f16` reduction computes in `f32`. A sum or product of `f16`s gives an `f32` total, for its caller to round once. A total that doesn't fit its type fails, with its value wrapped to that type. pub trait Reduction { @@ -177,8 +171,6 @@ pub struct Xor; /// The fold of integers by bitwise XNOR, `!(a ^ b)`, and of Booleans by `==`: true where an even number of them are false. #[derive(Clone, Copy, Debug)] pub struct Xnor; -forms!(Xor at bool: |_, [x, y]| x != y); -forms!(Xnor at bool: |_, [x, y]| x == y); /// `x` folded by `add` from `unit`, through the running results that `running_results` gives, which `merge` combines. #[inline(always)] diff --git a/core/src/reduce/run.rs b/core/src/reduce/run.rs index 89d45f8..00ff9ba 100644 --- a/core/src/reduce/run.rs +++ b/core/src/reduce/run.rs @@ -95,7 +95,8 @@ fn inverse_result(result: Option>, items: usize) - fn native_read(shape: &[usize], axis: Option, checked: bool, input: Read<'_>) -> Result<(Buffer, bool, Option>), Error> where T: Cast, R: Reduction, R::Total: Totals { let dims = match axis { Some(axis) => Lanes::new(shape, axis)?.dims(), None => [1, crate::layout::count(shape), 1] }; - if let Some((view, mapping)) = input.positional(shape) && let Some(data) = T::of(view) { + let source = input.positional(shape).or_else(|| input.mapped_source().map(|(view, map)| (view, map.clone()))); + if let Some((view, mapping)) = source && let Some(data) = T::of(view) { if axis.is_none() && R::COMMUTATIVE && let Some(block) = mapping.block(shape) { return crate::error::checked(|| { let (dims, order) = block.fold(None); @@ -103,6 +104,13 @@ fn native_read(shape: &[usize], axis: Option, checked: bool, input: Ok((R::Total::buffer_of(totals), fits, order)) }); } + let flat = [crate::layout::count(shape)]; + if axis.is_none() && rereads_cheaply(&mapping, &flat) { + return crate::error::checked(|| { + let (totals, fits, order) = along_mapped::(data, &mapping, &flat, 0); + Ok((R::Total::buffer_of(totals), fits, order)) + }); + } if let Some(axis) = axis { return crate::error::checked(|| { let slide = if R::COMMUTATIVE { diff --git a/core/src/slice.rs b/core/src/slice.rs index 8f8e1cd..385f5b5 100644 --- a/core/src/slice.rs +++ b/core/src/slice.rs @@ -1,5 +1,4 @@ -//! Slice forms of rustymath's functions. Each function object, such as [`Pow`], has a scalar form `one` and a slice form `many`. Both run code that rustymath compiles once for each item type, inside its own crate. A caller's build settings can't change that code, and rustymath's speed test times what callers run. An edit to a function's body rebuilds only rustymath. -use crate::float::Float; +//! Shared slice loops, output storage and native bulk arguments. use num_complex::Complex; use std::mem::MaybeUninit; @@ -56,86 +55,14 @@ nan_items!(f32, f64, half::f16, Complex, Complex); macro_rules! exact_items { ($($t:ty),+) => {$( impl KernelItem for $t { const NAN: bool = false; #[inline(always)] fn has_nan(self) -> bool { false } } )+} } exact_items!(u8, i16, i32, i64, bool, char, std::cmp::Ordering); -/// The results of a function object of `N` items of type `T`: their type and whether any can fail. [`Unary`], [`Binary`] and [`Ternary`] each require it, so code that takes a function object of any arity reads these through one bound. -pub trait Results { - /// The type of each result. - type Out: KernelItem; - /// Whether no result fails: `checked` always gives `Some`, and `many` and `many_into` always give true. A total float function gives NaN from arguments that aren't NaN only where NaN is its value, as [`Subtract`](crate::Subtract) gives for ∞ - ∞ and [`Sin`](crate::Sin) for ±∞. - const TOTAL: bool; -} - -/// Native scalar and slice interfaces. The slice body initializes every output slot even when an item fails. -macro_rules! arities { - ($($(#[$doc:meta])* $name:ident<$n:literal>($($x:ident),+);)+) => {$( - $(#[$doc])* - /// - /// # Safety - /// The slice body initializes every slot in its output before returning. - pub unsafe trait $name: Results + Sync { - fn one(&self, $($x: T),+) -> Self::Out; - fn judged(&self, $($x: T),+) -> (Self::Out, bool); - fn slice(&self, args: [Arg<'_, T>; $n], n: usize, out: &mut [MaybeUninit]) -> bool; - #[inline] - fn checked(&self, $($x: T),+) -> Option { let (r, ok) = self.judged($($x),+); ok.then_some(r) } - fn many(&self, $($x: Arg<'_, T>,)+ n: usize, out: &mut Vec) -> bool where T: Copy + Sync { - unsafe { appended(out, n, |slots| crate::parallel::split(slots, CHUNK, &|start, slots| { - crate::parallel::in_chunks(start, slots, |at, part| self.slice([$($x.part(at, part.len())),+], part.len(), part)) - })) } - } - fn many_into(&self, $($x: Arg<'_, T>,)+ n: usize, out: &mut crate::parallel::Slots<'_, Self::Out>) -> bool where T: Copy { - // SAFETY: the unsafe object trait promises that `slice` initializes every output. - unsafe { out.fill(n, |slots| crate::parallel::in_chunks(0, slots, |at, part| self.slice([$($x.part(at, part.len())),+], part.len(), part))) } - } - } - )+}; -} - -arities! { - /// A function of one item as a value, with a scalar form and a slice form. rustymath's function objects, such as [`Exp`](crate::Exp), implement it at each item type they take: every [`Float`] type for a float function, and complex numbers, integers or Booleans for others. `Out` is the type of each result, which is the item type unless the function gives another, as [`Abs`](crate::Abs) gives the real magnitude of a complex number. Its methods take `&self`, so a function object can hold settings. - Unary<1>(x); - /// A function of two items as a value, as [`Unary`] is for one. - Binary<2>(x, y); - /// A function of three items as a value, as [`Unary`] is for one. - Ternary<3>(x, y, z); -} - -/// `K` with its two arguments in the other order, as a [`Binary`] at each item type where `K` is one: `Flipped(k).one(x, y)` is `k.one(y, x)`, and likewise for `many` and `checked`. +/// A binary formula with its arguments in the other order. #[derive(Clone, Copy, Debug)] pub struct Flipped(pub K); -impl> Results for Flipped { - type Out = K::Out; - const TOTAL: bool = K::TOTAL; -} -unsafe impl> Binary for Flipped { - fn one(&self, x: T, y: T) -> K::Out { self.0.one(y, x) } - fn judged(&self, x: T, y: T) -> (K::Out, bool) { self.0.judged(y, x) } - fn slice(&self, [x, y]: [Arg<'_, T>; 2], n: usize, out: &mut [MaybeUninit]) -> bool { self.0.slice([y, x], n, out) } -} - -/// `K` with no result failing. Each result is `K`'s, wrapped to its type, as static kinds and NumPy keep integer results. LLVM removes `K`'s test of whether a result fits, because the loop never reads it. The adapter suits an object whose scalar form is its whole computation, as integer arithmetic's is. +/// An integer reduction or product whose results wrap to its output type. #[derive(Clone, Copy, Debug)] pub struct Wrapped(pub K); -impl, const N: usize> Results for Wrapped { - type Out = K::Out; - const TOTAL: bool = true; -} -macro_rules! wrapped { - ($arity:ident, $n:literal, $($x:ident),+) => { - unsafe impl> $arity for Wrapped { - fn one(&self, $($x: T),+) -> K::Out { self.0.one($($x),+) } - fn judged(&self, $($x: T),+) -> (K::Out, bool) { (self.0.one($($x),+), true) } - fn slice(&self, x: [Arg<'_, T>; $n], n: usize, out: &mut [MaybeUninit]) -> bool { - each(x, n, out, #[inline(always)] |[$($x),+]| Total(self.0.one($($x),+))) - } - } - }; -} -wrapped!(Unary, 1, x); -wrapped!(Binary, 2, x, y); -wrapped!(Ternary, 3, x, y, z); - /// Whether `r`, the result for arguments `args`, keeps the rule that a NaN result needs a NaN argument. #[inline(always)] pub(crate) fn kept(r: O, args: [E; N]) -> bool { !r.has_nan() | args.iter().any(|a| a.has_nan()) } @@ -144,31 +71,18 @@ pub(crate) fn kept(r: O, args: [E; pub(crate) trait Outcome: Sized { /// Whether an output that `split` fails can still succeed, as a NaN output does where an argument is NaN. Only a plain output can. const RESCUE: bool; - /// Whether no output fails: a [`Total`] output, or a plain output of a type without NaN. - const TOTAL: bool; - /// The output, as a scalar form gives it. `None` panics, because a scalar form has no other way to report it. - fn value(self) -> O; /// The output, and whether it succeeded, judged by the output alone, so that a NaN output fails. That costs one test for each item, where testing each argument too would cost one more for each argument. fn split(self) -> (O, bool); - /// The output, and whether it succeeded for arguments `args`. - #[inline(always)] - fn judged(self, args: [E; N]) -> (O, bool) { let (r, ok) = self.split(); (r, ok || (Self::RESCUE && kept(r, args))) } } impl Outcome for O { const RESCUE: bool = true; - const TOTAL: bool = !O::NAN; - #[inline(always)] - fn value(self) -> O { self } #[inline(always)] fn split(self) -> (O, bool) { (self, !self.has_nan()) } } impl Outcome for Option { const RESCUE: bool = false; - const TOTAL: bool = false; - #[inline(always)] - fn value(self) -> O { self.expect("the result doesn't fit the output type") } #[inline(always)] fn split(self) -> (O, bool) { (self.unwrap_or_default(), self.is_some()) } } @@ -176,9 +90,6 @@ impl Outcome for Option { /// An output wrapped to its type, and whether it's the exact result, as an integer power gives it. A result that isn't exact fails, and its slot holds the wrapped output. That output is NumPy's result for fixed-width integers. impl Outcome for (O, bool) { const RESCUE: bool = false; - const TOTAL: bool = false; - #[inline(always)] - fn value(self) -> O { self.0 } #[inline(always)] fn split(self) -> (O, bool) { self } } @@ -189,15 +100,10 @@ pub(crate) struct Total(pub O); impl Outcome for Total { const RESCUE: bool = false; - const TOTAL: bool = true; - #[inline(always)] - fn value(self) -> O { self.0 } #[inline(always)] fn split(self) -> (O, bool) { (self.0, true) } } -/// Whether an item body never fails, as its outcome type `R` says. `forms!` passes the body as a closure, which names that type without running it. -pub(crate) const fn total, const N: usize>(_: fn(&K, [E; N]) -> R) -> bool { R::TOTAL } /// The items of `x` that result `i` takes. #[inline(always)] @@ -257,23 +163,8 @@ pub(crate) fn each_usual>( }) } -/// Declares a function object of one argument whose slice form runs `each_usual`. `$usual` says whether an item needs none of the function's special values, `$fast` computes such an item, and `$name` is the full function. An optional `$rare => $value` gives the form for a block whose items all take one rare value, as `exp` gives 0 below its underflow bound. A trailing `; total` wraps the full function's results in [`Total`]. -macro_rules! usual_object { - ($(#[doc = $doc:literal])* $object:ident: $name:ident, $usual:expr => $fast:expr $(, $rare:expr => $value:expr)? $(; $outcome:ident)?) => { - $(#[doc = $doc])* - #[derive(Clone, Copy, Debug)] - pub struct $object; - forms!($object at real: |_, [x]| usual_object!(@outcome [$($outcome)?], $name::<_, { crate::FUSED }>(x)), - many |_, [x], n, out| crate::slice::each_usual(x, n, out, ($usual, $fast), usual_object!(@rare $($rare, $value)?), - #[inline(always)] |x| usual_object!(@outcome [$($outcome)?], $name::<_, { crate::FUSED }>(x)))); - }; - (@rare) => { crate::slice::none() }; - (@rare $rare:expr, $value:expr) => { ($rare, $value) }; - (@outcome [], $e:expr) => { $e }; - (@outcome [total], $e:expr) => { crate::slice::Total($e) }; -} -/// [`each`] with the item body marked `#[inline(always)]`. Every slice form runs its loop this way: `forms!`'s default slice form, and a custom slice form once it has read its settings. LLVM leaves a large unmarked body out of line, and its call stops the loop vectorizing. +/// [`each`] with the item body marked `#[inline(always)]`, so native bulk loops inline their scalar body. macro_rules! each { ($x:expr, $n:expr, $out:expr, |$args:pat_param| $body:expr $(,)?) => { crate::slice::each( @@ -282,101 +173,6 @@ macro_rules! each { }; } -/// Writes a function object's forms, compiled inside rustymath. The list names the item types, such as `Complex` or `i64`. `real` stands for `f32` and `f64`, and also writes an `f16` form, which runs the `f32` form and narrows its outputs. `=> T` makes the output type `T` in place of the item type: for each listed type, or after `real` for every float type. The object is the first argument of the item body, and its items the bracketed ones. The body gives an output, an `Option` of one, or a [`Total`] output, as [`Outcome`] says. A `many` body replaces the default slice form, which runs the item body through [`each!`]. Its arguments are the object, the bracketed `Arg`s, the number of results, and the slots it writes them into. Attributes before the object, such as `#[inline]`, go on its scalar forms, `one` and `judged`, so that a caller in another crate can inline a small form into its own loop over items. -/// -/// A function object is a concrete type, never a generic one, so that its forms compile here. Settings that choose a loop, such as whether a tolerance is 0, are read once per call inside `many`, which then runs `each!`. -/// -/// ```text -/// forms!(Choose at real: |_, [x, y, m]| Total(if m != 0.0 { x } else { y })); -/// forms!(Choose at bool: |_, [x, y, m]| (m & x) | (!m & y)); -/// forms!(Gcd at i64: |_, [x, y]| gcd(x, y)); -/// forms!(Less at real => bool: |k, [x, y]| k.0.less(x, y), -/// many |k, [x, y], n, out| with_tolerance!(k.0, |t| each!([x, y], n, out, |[x, y]| t.less(x, y)))); -/// ``` -macro_rules! forms { - ($(#[$($m:tt)+])* $object:ident at real $(=> $o:ty)? : |$k:pat_param, $args:tt| $body:expr $(, many |$mk:pat_param, $margs:tt, $n:ident, $out:ident| $many:expr)?) => { - forms!(@types [$(#[$($m)+])*], $object, $k, $args, $body, [$((|$mk, $margs, $n, $out| $many))?], [f32 $(=> $o)?, f64 $(=> $o)?]); - forms!(@half [$(#[$($m)+])*], $object, forms!(@out half::f16 $(, $o)?), forms!(@narrow $($o)?), $args); - }; - ($(#[$($m:tt)+])* $object:ident at $($e:ty $(=> $o:ty)?),+ : |$k:pat_param, $args:tt| $body:expr $(, many |$mk:pat_param, $margs:tt, $n:ident, $out:ident| $many:expr)?) => { - forms!(@types [$(#[$($m)+])*], $object, $k, $args, $body, [$((|$mk, $margs, $n, $out| $many))?], [$($e $(=> $o)?),+]); - }; - // A generic type, such as `Indicator`, goes in brackets, because a type can't be followed by `at`. - ($(#[$($m:tt)+])* [$object:ty] at $($e:ty $(=> $o:ty)?),+ : |$k:pat_param, $args:tt| $body:expr $(, many |$mk:pat_param, $margs:tt, $n:ident, $out:ident| $many:expr)?) => { - forms!(@types [$(#[$($m)+])*], $object, $k, $args, $body, [$((|$mk, $margs, $n, $out| $many))?], [$($e $(=> $o)?),+]); - }; - // One type at a time, so that the optional `many` group, which repeats on its own, needn't sit inside the types' repetition. - (@types $attrs:tt, $object:ty, $k:pat_param, $args:tt, $body:expr, $many:tt, [$e:ty $(=> $o:ty)? $(, $($rest:tt)+)?]) => { - forms!(@concrete $attrs, $object, $e, forms!(@out $e $(, $o)?), $k, $args, $body, $many); - $(forms!(@types $attrs, $object, $k, $args, $body, $many, [$($rest)+]);)? - }; - (@concrete $attrs:tt, $object:ty, $e:ty, $o:ty, $k:pat_param, [$x:ident], $body:expr, $many:tt) => { - forms!(@impl $attrs, $object, Unary, 1, $e, $o, $k, [$x], $body, $many); - }; - (@concrete $attrs:tt, $object:ty, $e:ty, $o:ty, $k:pat_param, [$x:ident, $y:ident], $body:expr, $many:tt) => { - forms!(@impl $attrs, $object, Binary, 2, $e, $o, $k, [$x, $y], $body, $many); - }; - (@concrete $attrs:tt, $object:ty, $e:ty, $o:ty, $k:pat_param, [$x:ident, $y:ident, $z:ident], $body:expr, $many:tt) => { - forms!(@impl $attrs, $object, Ternary, 3, $e, $o, $k, [$x, $y, $z], $body, $many); - }; - (@impl [$(#[$($m:tt)+])*], $object:ty, $arity:ident, $arity_n:literal, $e:ty, $o:ty, $k:pat_param, [$($x:ident),+], $body:expr, $many:tt) => { - impl crate::slice::Results<$e, $arity_n> for $object { - type Out = $o; - const TOTAL: bool = crate::slice::total::(|$k, [$($x),+]| $body); - } - unsafe impl crate::slice::$arity<$e> for $object { - $(#[$($m)+])* - fn one(&self, $($x: $e),+) -> $o { - let $k = self; - crate::slice::Outcome::<$o>::value($body) - } - $(#[$($m)+])* - fn judged(&self, $($x: $e),+) -> ($o, bool) { - let $k = self; - crate::slice::Outcome::<$o>::judged($body, [$($x),+]) - } - forms!(@slice $e => $o, $k, [$($x),+], $body, $many); - } - }; - (@slice $e:ty => $o:ty, $k:pat_param, [$($x:ident),+], $body:expr, []) => { - fn slice(&self, args: [crate::Arg<'_, $e>; forms!(@n $($x),+)], n: usize, out: &mut [std::mem::MaybeUninit<$o>]) -> bool { - let $k = self; - each!(args, n, out, |[$($x),+]| $body) - } - }; - (@slice $e:ty => $o:ty, $k:pat_param, $args:tt, $body:expr, [(|$mk:pat_param, [$($mx:ident),+], $n:ident, $out:ident| $many:expr)]) => { - fn slice(&self, [$($mx),+]: [crate::Arg<'_, $e>; forms!(@n $($mx),+)], $n: usize, $out: &mut [std::mem::MaybeUninit<$o>]) -> bool { - let $mk = self; - $many - } - }; - (@half $attrs:tt, $object:ty, $o:ty, $narrow:expr, [$x:ident]) => { forms!(@half_impl $attrs, $object, Unary, 1, $o, $narrow, [$x]); }; - (@half $attrs:tt, $object:ty, $o:ty, $narrow:expr, [$x:ident, $y:ident]) => { forms!(@half_impl $attrs, $object, Binary, 2, $o, $narrow, [$x, $y]); }; - (@half $attrs:tt, $object:ty, $o:ty, $narrow:expr, [$x:ident, $y:ident, $z:ident]) => { forms!(@half_impl $attrs, $object, Ternary, 3, $o, $narrow, [$x, $y, $z]); }; - (@half_impl [$(#[$($m:tt)+])*], $object:ty, $arity:ident, $arity_n:literal, $o:ty, $narrow:expr, [$($x:ident),+]) => { - impl crate::slice::Results for $object { - type Out = $o; - const TOTAL: bool = >::TOTAL; - } - unsafe impl crate::slice::$arity for $object { - $(#[$($m)+])* - fn one(&self, $($x: half::f16),+) -> $o { ($narrow)(crate::slice::$arity::::one(self, $($x.to_f32()),+)) } - $(#[$($m)+])* - fn judged(&self, $($x: half::f16),+) -> ($o, bool) { - let (r, ok) = crate::slice::$arity::::judged(self, $($x.to_f32()),+); - (($narrow)(r), ok) - } - fn slice(&self, x: [crate::Arg<'_, half::f16>; $arity_n], n: usize, out: &mut [std::mem::MaybeUninit<$o>]) -> bool { - crate::slice::working(x, n, out, |x, m, o| crate::slice::$arity::::slice(self, x, m, o), $narrow) - } - } - }; - (@out $e:ty) => { $e }; - (@out $e:ty, $o:ty) => { $o }; - (@narrow) => { half::f16::from_f32 }; - (@narrow $o:ty) => { std::convert::identity }; - (@n $($x:ident),+) => { [$(stringify!($x)),+].len() }; -} /// Runs `$body` for each item `$i` of a result of `$n` items, with `$a` and `$b` the items of `$x` and `$y` that pair with it. Each pairing of `One` and `Many` has a loop of its own. The loops index slices of the result's length with no iterator adapters, and vectorize at any build settings. macro_rules! each_pair { @@ -515,14 +311,6 @@ pub(crate) fn two_pass( ok } -/// `f`, a slice form at `F`'s working type, for narrower floats, as `forms!`'s `f16` forms run it: the items convert to the working type and run, and `narrow` turns each output into the output at `F`. -pub(crate) fn working( - x: [Arg<'_, F>; N], n: usize, out: &mut [MaybeUninit

], f: impl Fn([Arg<'_, F::Working>; N], usize, &mut [MaybeUninit]) -> bool, - narrow: impl Fn(O) -> P) -> bool { - let mut bufs = [(); N].map(|()| Vec::new()); - let mut next = bufs.each_mut().into_iter(); - unsafe { buffered(n, out, |o| f(x.map(|a| a.converted(0, n, next.next().unwrap(), F::working)), n, o), narrow) } -} /// `f`, a slice form writing `n` outputs of type `O`, run into a buffer, with `each` turning each output into the `P` that `out` takes. It gives `f`'s result. /// diff --git a/core/src/special.rs b/core/src/special.rs index 7b2d571..8acfa08 100644 --- a/core/src/special.rs +++ b/core/src/special.rs @@ -627,14 +627,15 @@ fn erfcinv64(q: f64) -> f64 { maddock::((q - 1.0).abs(), q.min(2.0 - q)).copysign(1.0 - q) } -/// [`erfc`](crate::erfc) as a function object, with the slice form `many`, which runs `erfc_blocks` on a slice. +/// [`erfc`](crate::erfc) as a formula object, with a native bulk routine that groups inputs by region. #[derive(Clone, Copy, Debug)] pub struct Erfc; -forms!(Erfc at real: |_, [x]| crate::slice::Total(erfc::<_, { crate::FUSED }>(x)), -many |_, [x], n, out| match x { - Arg::Many(x) => erfc_blocks::<_, { crate::FUSED }>(&x[..n], out), - x => each!([x], n, out, |[x]| crate::slice::Total(erfc::<_, { crate::FUSED }>(x))), -}); +pub(crate) fn erfc_slice(_: &Erfc, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [x]: [Arg<'_, R>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + match x { + Arg::Many(x) => erfc_blocks::<_, { crate::FUSED }>(&x[..n], out), + x => each!([x], n, out, |[x]| crate::slice::Total(erfc::<_, { crate::FUSED }>(x))), + } +} /// `erfc` on each item of `x`, written into `out`, through `slice::blocks` of 512 items. A block whose magnitudes all lie below `ERFC_SPLIT` computes only `erfc_near`, and one whose magnitudes all reach it computes only `erfc_far`. Where `branches` holds, a mixed block's items are partitioned: `placed` puts the positions of those below `ERFC_SPLIT` at the front and the others' at the back, and each formula runs on its own items, gathered into one buffer, as `regions` does, because a branch on each item was mispredicted on mixed arguments. Elsewhere a mixed block runs `erfc`, which computes both formulas in every lane. Partitioning every block added about 12 instructions an item, more than computing in place where a block lies in one region. Blocks of 64, the size `each_usual` uses, cost the normal CDF about 7% more instructions on shuffled and on sorted items than blocks of 512. fn erfc_blocks(x: &[R], out: &mut [MaybeUninit]) -> bool { @@ -662,17 +663,19 @@ fn erfc_blocks(x: &[R], out: &mut [MaybeUninit] true }) } -/// [`erfinv`](crate::erfinv) as a function object, with the slice form `many`, which runs `regions` on a slice of f64s. +/// [`erfinv`](crate::erfinv) as a formula object, with a native bulk routine that groups f64 inputs by region. #[derive(Clone, Copy, Debug)] pub struct Erfinv; -forms!(Erfinv at real: |_, [p]| erfinv::<_, { crate::FUSED }>(p), - many |_, [p], n, out| inverses(p, n, out, #[inline(always)] |p| erfinv::<_, { crate::FUSED }>(p), #[inline(always)] |p| (p.abs(), 1.0 - p.abs(), p))); +pub(crate) fn erfinv_slice(_: &Erfinv, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [p]: [Arg<'_, R>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + inverses(p, n, out, #[inline(always)] |p| erfinv::<_, { crate::FUSED }>(p), #[inline(always)] |p| (p.abs(), 1.0 - p.abs(), p)) +} -/// [`erfcinv`](crate::erfcinv) as a function object, with the slice form `many`, which runs `regions` on a slice of f64s. +/// [`erfcinv`](crate::erfcinv) as a formula object, with a native bulk routine that groups f64 inputs by region. #[derive(Clone, Copy, Debug)] pub struct Erfcinv; -forms!(Erfcinv at real: |_, [q]| erfcinv::<_, { crate::FUSED }>(q), - many |_, [q], n, out| inverses(q, n, out, #[inline(always)] |q| erfcinv::<_, { crate::FUSED }>(q), #[inline(always)] |q| ((q - 1.0).abs(), q.min(2.0 - q), 1.0 - q))); +pub(crate) fn erfcinv_slice(_: &Erfcinv, _: &crate::recipe::NativeOps<{ crate::FUSED }>, [q]: [Arg<'_, R>; 1], n: usize, out: &mut [MaybeUninit]) -> bool { + inverses(q, n, out, #[inline(always)] |q| erfcinv::<_, { crate::FUSED }>(q), #[inline(always)] |q| ((q - 1.0).abs(), q.min(2.0 - q), 1.0 - q)) +} /// The slice form of [`Erfinv`] and [`Erfcinv`]: `regions` on a slice of f64s, and `f` on each item otherwise. `parts` gives an f64 item's `maddock` arguments and the value whose sign its result takes. fn inverses(x: Arg<'_, R>, n: usize, out: &mut [MaybeUninit], f: impl Fn(R) -> R, parts: impl Fn(f64) -> (f64, f64, f64)) -> bool { diff --git a/core/src/structure.rs b/core/src/structure.rs index 8d9d0ad..7957164 100644 --- a/core/src/structure.rs +++ b/core/src/structure.rs @@ -19,29 +19,74 @@ impl Count { pub fn of(n: i64) -> Self { Count { negative: n < 0, magnitude: usize::try_from(n.unsigned_abs()).unwrap_or(usize::MAX) } } } +/// A slice bound: a positional offset or a key on the selected axis. +#[derive(Clone, Debug, PartialEq)] +pub enum SliceBound { Position(isize), Key(Arc) } + +/// A slice before the selected axis's length and keys are known. A keyed stop includes its position. +#[derive(Clone, Debug, PartialEq)] +pub struct SliceSpec { pub start: Option, pub stop: Option, pub step: isize } + +impl std::fmt::Display for SliceSpec { + /// Diagnostic spelling, not source syntax for a language's selectors. + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let bound = |v: &Option| match v { + Some(SliceBound::Position(i)) => i.to_string(), Some(SliceBound::Key(k)) => format!("{k:?}"), None => "None".into() + }; + write!(f, "slice({}, {}, {})", bound(&self.start), bound(&self.stop), self.step) + } +} + +impl SliceSpec { + /// Resolve this slice against the actual axis, clipping positional bounds to its length. + pub fn selector(&self, layout: &Layout, axis: usize) -> Result { + if self.step == 0 { return Err(Error::new(Kind::Domain, "slice step cannot be zero")); } + let n = *layout.shape().get(axis).ok_or_else(|| Error::new(Kind::Index, "too many axis selectors"))? as i128; + let direction = self.step.signum() as i128; + let (low, high) = if direction > 0 { (0, n) } else { (-1, n - 1) }; + let bound = |value: Option<&SliceBound>, default: i128, stop: bool| -> Result { + let at = match value { + None => return Ok(default), + Some(SliceBound::Position(i)) => { let i = *i as i128; if i < 0 { i + n } else { i } }, + Some(SliceBound::Key(name)) => { + let i = layout.keys(axis).and_then(|keys| keys.position(name)) + .ok_or_else(|| Error::new(Kind::Index, format!("missing key: {name}")))?; + i as i128 + if stop { direction } else { 0 } + } + }; + Ok(at.clamp(low, high)) + }; + let start = bound(self.start.as_ref(), if direction > 0 { 0 } else { n - 1 }, false)?; + let stop = bound(self.stop.as_ref(), if direction > 0 { n } else { -1 }, true)?; + let distance = (stop - start) * direction; + let len = if distance <= 0 { 0 } else { ((distance - 1) / (self.step as i128).abs() + 1) as usize }; + Ok(Selector::Range { start: start.max(0) as usize, len, step: self.step }) + } +} + /// Written index parts, before the selected axes and their lengths are known. #[derive(Clone, Debug, PartialEq)] -pub enum IndexPart { Axis(T), All, NewAxis, Rest } +pub enum IndexPart { Axis(T), Slice(Arc), All, NewAxis, Rest } impl IndexPart { /// The part with its axis payload borrowed. pub fn as_ref(&self) -> IndexPart<&T> { match self { IndexPart::Axis(t) => IndexPart::Axis(t), IndexPart::All => IndexPart::All, IndexPart::NewAxis => IndexPart::NewAxis, - IndexPart::Rest => IndexPart::Rest + IndexPart::Rest => IndexPart::Rest, IndexPart::Slice(slice) => IndexPart::Slice(slice.clone()) } } /// The part with `f` applied to its axis payload. The structural markers stay as they are. pub fn map(self, f: impl FnOnce(T) -> U) -> IndexPart { match self { IndexPart::Axis(t) => IndexPart::Axis(f(t)), IndexPart::All => IndexPart::All, IndexPart::NewAxis => IndexPart::NewAxis, - IndexPart::Rest => IndexPart::Rest + IndexPart::Rest => IndexPart::Rest, IndexPart::Slice(slice) => IndexPart::Slice(slice) } } /// The part with fallible `f` applied to its axis payload, or `f`'s error. pub fn try_map(self, f: impl FnOnce(T) -> Result) -> Result, E> { Ok(match self { IndexPart::Axis(t) => IndexPart::Axis(f(t)?), IndexPart::All => IndexPart::All, IndexPart::NewAxis => IndexPart::NewAxis, - IndexPart::Rest => IndexPart::Rest + IndexPart::Rest => IndexPart::Rest, IndexPart::Slice(slice) => IndexPart::Slice(slice) }) } /// The axis payload, or `None` for a structural marker. @@ -50,7 +95,7 @@ impl IndexPart { /// Expands Rest once geometry is known, with one source claim per consuming part. Axis payloads identify the original written part. pub fn index_claims<'a>(rank: usize, parts: &[IndexPart<(Option<&'a str>, usize)>]) -> Result<(Vec>, Vec<(Option<&'a str>, usize)>), Error> { - let used: usize = parts.iter().map(|p| match p { IndexPart::Axis((_, n)) => *n, IndexPart::All => 1, _ => 0 }).sum(); + let used: usize = parts.iter().map(|p| match p { IndexPart::Axis((_, n)) => *n, IndexPart::All | IndexPart::Slice(_) => 1, _ => 0 }).sum(); if used > rank { return Err(Error::new(Kind::Index, "too many axis selectors")); } if parts.iter().filter(|p| matches!(p, IndexPart::Rest)).count() > 1 { return Err(Error::new(Kind::Index, "only one rest selector is allowed")); } let mut written = Vec::new(); @@ -65,6 +110,7 @@ pub fn index_claims<'a>(rank: usize, parts: &[IndexPart<(Option<&'a str>, usize) written.push(IndexPart::All); claims.push((None, 1)); } + IndexPart::Slice(slice) => { written.push(IndexPart::Slice(slice.clone())); claims.push((None, 1)); } IndexPart::NewAxis => written.push(IndexPart::NewAxis), IndexPart::Rest => { for _ in used..rank { @@ -136,6 +182,7 @@ pub fn routes( let (selector, single) = decode(i, axis)?; (selector, width, single) } + IndexPart::Slice(slice) => (slice.selector(layout, axis).map_err(&err)?, 1, false), _ => (Selector::All, 1, false), }; one &= single; diff --git a/core/src/tolerant.rs b/core/src/tolerant.rs index 3acb5f8..429d26c 100644 --- a/core/src/tolerant.rs +++ b/core/src/tolerant.rs @@ -1,10 +1,10 @@ //! Tolerant comparison of floats: `x` equals `y` within a relative tolerance `t` when `|x-y|` is less than `t` times the larger of `|x|` and `|y|`. APL's definition says at most, which differs only where `|x-y|` lands exactly on that bound. Dyalog APL calls `t` the comparison tolerance, `⎕CT`, and its default is 1e-14. A tolerance has no effect when either number is zero. //! -//! A tolerance is a value of a type that implements `Tolerant`: a `Tolerance`, or `Exactly`, whose methods compile to plain IEEE comparisons and `floor`. A function written once over `Tolerant` compiles to a version for each. `with_tolerance!` picks one when the program runs, outside any loop. The methods are generic over the float type and compute at its width. `TolerantComplex` gives the complex versions. The function objects `Equal`, `NotEqual`, `Less`, `LessEqual`, `Greater`, `GreaterEqual`, `Matches`, `Compare`, `Floor`, `Ceil`, `Residue`, `FloorDivide`, `Gcd` and `Lcm` each hold a `Tolerance` and give the same results as the methods, one item at a time or over slices. A slice form picks its loop with `with_tolerance!` once. A comparison's slice form compares many items with one number through the range of floats equal to it (`Tolerant::range`), so each item takes one or two plain comparisons. All but `Matches` and `Compare` also have integer forms, from the `integer` module, which ignore the tolerance. `Floor` and `Ceil` also take Booleans. +//! `Tolerant` supplies real comparison and rounding; `TolerantComplex` supplies complex operations. Formula objects hold a Tolerance. Native bulk loops choose Exactly or Tolerance once through `with_tolerance!`, and comparison kernels use `against` for one constant argument. //! //! Tolerant equality isn't transitive: `x` can equal `y` and `y` equal `z` while `x` doesn't equal `z`. Hash maps and sorting need a transitive `PartialEq`. This module therefore has no float type with a tolerant `==`. Each comparison is a method of the tolerance. use crate::{ - Arg, BETWEEN_CHECKS, buffer::Element, convert::{float, int}, error::Error, float::{Float, Native}, interrupt::check, recipe::{self, Floats}, + Arg, BETWEEN_CHECKS, convert::{float, int}, error::Error, float::{Float, Native}, interrupt::check, recipe::{self, Floats}, search::{self, Chains, Key, Links}, slice::{Total, pairs_into} }; use num_bigint::BigInt; @@ -25,6 +25,9 @@ impl Tolerance { /// The tolerance `t`, when it lies from 0 up to but not including 1. `None` for any other `t`, including NaN. #[inline] pub fn new(t: f64) -> Option { (0.0..1.0).contains(&t).then_some(Self(t)) } + /// A valid runtime tolerance in place of this formula's declared tolerance. + #[inline] + pub(crate) fn overridden(self, value: Option) -> Self { value.and_then(Self::new).unwrap_or(self) } /// Whether this tolerance is 0. #[inline] pub fn is_exact(self) -> bool { self.0 == 0.0 } @@ -290,45 +293,6 @@ macro_rules! with_tolerance { }}; } -/// Writes a tolerant function object's forms from one item body over a tolerance `$t`, at `real` or at a list of concrete types, as `forms!` takes them, with its attributes, such as `#[inline]`. The object holds a `Tolerance`. Its slice form reads it once and runs the body's loop at `Exactly` where the tolerance is 0. Its scalar forms run the body at the `Tolerance` itself, which gives the exact result at 0 with no test for each item, so a caller's own loop over items has no branch to hoist. -macro_rules! tolerant { - ($(#[$($m:tt)+])* $object:ident at real $(=> $o:ty)?: $($rest:tt)+) => { tolerant!(@forms [$(#[$($m)+])*], $object, (real $(=> $o)?), $($rest)+); }; - ($(#[$($m:tt)+])* $object:ident at $($e:ty $(=> $o:ty)?),+: $($rest:tt)+) => { tolerant!(@forms [$(#[$($m)+])*], $object, ($($e $(=> $o)?),+), $($rest)+); }; - (@forms [$(#[$($m:tt)+])*], $object:ident, ($($at:tt)+), |$t:ident, [$($x:ident),+]| $body:expr) => { - forms!($(#[$($m)+])* $object at $($at)+: |k, [$($x),+]| { let $t = k.0; $body }, - many |k, [$($x),+], n, out| with_tolerance!(k.0, |$t| each!([$($x),+], n, out, |[$($x),+]| $body))); - }; -} - -/// `b` as 1 or 0 for an [`Indicator`]. Choosing between the two constants turns a comparison's lane mask into the result with one `and`. [`Element::truth`] converts through an integer instead, which suits the stored Booleans that `slice::masked` reads but would narrow and widen a lane mask here. -#[inline(always)] -fn indicated(b: bool) -> E where bool: crate::buffer::Source { - if b { crate::buffer::Source::::read(true) } else { crate::buffer::Source::::read(false) } -} - -/// Writes an [`Indicator`]'s forms for the tolerant comparison `$object` at each float type `$e`, through `forms!`: its item body over a tolerance `$t` as 1 or 0, and a slice form that compares many items with one number through `against`, writing each float in the same loop. -macro_rules! indicator { - ($object:ident: |$t:ident, [$x:ident, $y:ident]| $body:expr, $left:ident, $right:ident, $($e:ty),+) => {$( - forms!([Indicator<$object>] at $e: |k, [$x, $y]| crate::slice::Total(indicated::<$e>({ let $t = k.0.0; $body })), - many |k, [$x, $y], n, out| match against(k.0.0, $x, $y, n, out, Band::$left, Band::$right, |xs, lo, hi, band, out| in_band(xs, lo, hi, band, out, indicated::<$e>)) { - Some(ok) => ok, - None => with_tolerance!(k.0.0, |$t| each!([$x, $y], n, out, |[$x, $y]| crate::slice::Total(indicated::<$e>($body)))), - }); - )+}; -} - -/// Writes a tolerant comparison's forms at `real`: its item body over a tolerance `$t`, at the `Tolerance` itself, as `tolerant!` writes them, and a slice form that compares many items with one number through `against`. `$left` is the band an item must lie in when it is on the left of the one number, and `$right` the band when it is on the right. It also writes the comparison's [`Indicator`] forms at `f32` and `f64`. -macro_rules! compared { - ($object:ident: |$t:ident, [$x:ident, $y:ident]| $body:expr, $left:ident, $right:ident) => { - forms!($object at real => bool: |k, [$x, $y]| { let $t = k.0; $body }, - many |k, [$x, $y], n, out| match against(k.0, $x, $y, n, out, Band::$left, Band::$right, Banded::banded) { - Some(ok) => ok, - None => with_tolerance!(k.0, |$t| each!([$x, $y], n, out, |[$x, $y]| $body)), - }); - indicator!($object: |$t, [$x, $y]| $body, $left, $right, f32, f64); - impl Banding for $object { const BANDS: (Band, Band) = (Band::$left, Band::$right); } - }; -} /// Where an item must lie against the range of floats equal to one number (`Tolerant::range`) for a comparison with that number to hold. `NotAbove` and `NotBelow` include the range itself. #[derive(Clone, Copy)] @@ -396,9 +360,6 @@ pub fn against( Some(true) } -/// A comparison's results as floats, 1 where it holds and 0 where it doesn't, as a [`Binary`](crate::Binary) at `f32` and `f64` for each comparison here. A caller that computes in floats reads a comparison's results this way, rather than storing Booleans and converting them. Its slice form compares each item and writes its float in one loop. -#[derive(Clone, Copy, Debug)] -pub struct Indicator(pub K); /// Wasm has no instruction that narrows a 64-bit lane to a byte. LLVM writes a float comparison's results two items at a time, in about ten instructions for each item. Where it can prove that lanes hold masks, it also replaces wasm's narrowing instructions with slower code. Here each block of 16 items compares with vector instructions, then shuffles the first byte of each mask into place: 7 shuffles for `f64`'s 8 vectors, and 3 for `f32`'s 4. #[cfg(all(target_arch = "wasm32", target_feature = "simd128"))] @@ -463,16 +424,20 @@ mod narrowed { /// [`Tolerant::equal`], or [`TolerantComplex::complex_equal`] for complex items, at the tolerance it holds, as a function object. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Equal(pub Tolerance); -compared!(Equal: |t, [x, y]| t.equal(x, y), Inside, Inside); -forms!(Equal at Complex => bool, Complex => bool: |k, [x, y]| k.0.complex_equal(x, y), - many |k, [x, y], n, out| with_tolerance!(k.0, |t| complex_equals(t, [x, y], n, out, false))); +impl Banding for Equal { const BANDS: (Band, Band) = (Band::Inside, Band::Inside); } +pub(crate) fn complex_equal_slice(k: &Equal, o: &recipe::NativeOps<{ crate::FUSED }>, args: [Arg<'_, Complex>; 2], n: usize, out: &mut [std::mem::MaybeUninit]) { + let tolerance = k.0.overridden(o.tolerance); + with_tolerance!(tolerance, |t| { complex_equals(t, args, n, out, false); }); +} /// The negation of [`Equal`] at the tolerance it holds, as a function object, for real and complex items. #[derive(Clone, Copy, Debug, PartialEq)] pub struct NotEqual(pub Tolerance); -compared!(NotEqual: |t, [x, y]| !t.equal(x, y), Outside, Outside); -forms!(NotEqual at Complex => bool, Complex => bool: |k, [x, y]| !k.0.complex_equal(x, y), - many |k, [x, y], n, out| with_tolerance!(k.0, |t| complex_equals(t, [x, y], n, out, true))); +impl Banding for NotEqual { const BANDS: (Band, Band) = (Band::Outside, Band::Outside); } +pub(crate) fn complex_not_equal_slice(k: &NotEqual, o: &recipe::NativeOps<{ crate::FUSED }>, args: [Arg<'_, Complex>; 2], n: usize, out: &mut [std::mem::MaybeUninit]) { + let tolerance = k.0.overridden(o.tolerance); + with_tolerance!(tolerance, |t| { complex_equals(t, args, n, out, true); }); +} /// `Equal`'s slice form on complex items, or with `negate` `NotEqual`'s. Each `CHUNK` of items splits into slices of real and imaginary parts, as in `slice::split`, so `squared_equal` runs over plain slices with no branch and vectorizes. `scaled_equal` then takes the pairs it leaves undecided, one at a time. No result fails. #[inline(always)] @@ -490,22 +455,22 @@ fn complex_equals( /// [`Tolerant::less`] at the tolerance it holds, as a function object. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Less(pub Tolerance); -compared!(Less: |t, [x, y]| t.less(x, y), Below, Above); +impl Banding for Less { const BANDS: (Band, Band) = (Band::Below, Band::Above); } /// [`Tolerant::less_equal`] at the tolerance it holds, as a function object. #[derive(Clone, Copy, Debug, PartialEq)] pub struct LessEqual(pub Tolerance); -compared!(LessEqual: |t, [x, y]| t.less_equal(x, y), NotAbove, NotBelow); +impl Banding for LessEqual { const BANDS: (Band, Band) = (Band::NotAbove, Band::NotBelow); } /// [`Tolerant::greater`] at the tolerance it holds, as a function object. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Greater(pub Tolerance); -compared!(Greater: |t, [x, y]| t.greater(x, y), Above, Below); +impl Banding for Greater { const BANDS: (Band, Band) = (Band::Above, Band::Below); } /// [`Tolerant::greater_equal`] at the tolerance it holds, as a function object. #[derive(Clone, Copy, Debug, PartialEq)] pub struct GreaterEqual(pub Tolerance); -compared!(GreaterEqual: |t, [x, y]| t.greater_equal(x, y), NotBelow, NotAbove); +impl Banding for GreaterEqual { const BANDS: (Band, Band) = (Band::NotBelow, Band::NotAbove); } /// A tolerant comparison as steps of `recipe`'s formulas, which a compiler records through its `recipe::Ops`: the composition that `Tolerant`'s method of the same name uses, with `recipe::close` at the object's tolerance and the arguments in the same order. pub trait Stepped { @@ -564,38 +529,26 @@ impl Stepped for GreaterEqual { fn steps + ?Sized>(self, o: &O, x: O::V, y: O::V) -> O::B { LessEqual(self.0).steps(o, y, x) } } -/// [`Tolerant::matches`] at the tolerance it holds, as a function object: `equal`, except that NaN matches NaN. -#[derive(Clone, Copy, Debug, PartialEq)] -pub struct Matches(pub Tolerance); -compared!(Matches: |t, [x, y]| t.matches(x, y), Inside, Inside); - -/// [`Tolerant::compare`] at the tolerance it holds, as a function object. NaN orders after every number and equal to NaN, as [`Tolerant::matches`] has NaN match NaN. -#[derive(Clone, Copy, Debug, PartialEq)] -pub struct Compare(pub Tolerance); -tolerant!(Compare at real => Ordering: |t, [x, y]| t.compare(x, y).unwrap_or(x.is_nan().cmp(&y.is_nan()))); /// [`Tolerant::floor`], or [`TolerantComplex::complex_floor`] for complex items, at the tolerance it holds, as a function object. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Floor(pub Tolerance); -tolerant!(#[inline] Floor at real: |t, [y]| Total(t.floor(y))); -tolerant!(Floor at Complex, Complex: |t, [y]| Total(t.complex_floor(y))); /// [`Tolerant::ceil`], or [`TolerantComplex::complex_ceil`] for complex items, at the tolerance it holds, as a function object. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Ceil(pub Tolerance); -tolerant!(#[inline] Ceil at real: |t, [y]| Total(t.ceil(y))); -tolerant!(Ceil at Complex, Complex: |t, [y]| Total(t.complex_ceil(y))); /// [`Tolerant::residue`], or [`TolerantComplex::complex_residue`] for complex items, at the tolerance it holds, as a function object: `y` modulo `x`, for arguments `x` and `y` in that order. At no tolerance, where `y/x` passes the float's precision ([`far`]), the rounded quotient has lost its low digits, and [`exact_residue`] gives the residue. The slice form recomputes those items in a second pass. Its `i64` form is [`integer::residue`](crate::integer::residue), which needs no tolerance. #[derive(Clone, Copy, Debug, PartialEq)] pub struct Residue(pub Tolerance); -forms!(#[inline] Residue at real: |k, [x, y]| Total(if k.0.is_exact() { exact_residue(x, y) } else { k.0.residue(x, y) }), -many |k, [x, y], n, out| with_tolerance!(k.0, |t| { - if t.value() != 0.0 { return each!([x, y], n, out, |[x, y]| Total(t.residue(x, y))); } - if !each!([x, y], n, out, |[x, y]| (t.residue(x, y), near(x, y))) { replace_far(x, y, n, out); } - true -})); -tolerant!(Residue at Complex, Complex: |t, [x, y]| t.complex_residue(x, y)); +pub(crate) fn residue_slice(k: &Residue, o: &recipe::NativeOps<{ crate::FUSED }>, [x, y]: [Arg<'_, R>; 2], n: usize, out: &mut [std::mem::MaybeUninit]) -> bool { + let tolerance = k.0.overridden(o.tolerance); + with_tolerance!(tolerance, |t| { + if t.value() != 0.0 { return each!([x, y], n, out, |[x, y]| Total(t.residue(x, y))); } + if !each!([x, y], n, out, |[x, y]| (t.residue(x, y), near(x, y))) { replace_far(x, y, n, out); } + true + }) +} /// Recomputes through [`exact_residue`] the residue in `out` of each item that isn't [`near`]. Of those, only the items that [`far`] picks change. The slice form runs this only where its loop finds an item that isn't near. fn replace_far(x: Arg<'_, F>, y: Arg<'_, F>, n: usize, out: &mut [std::mem::MaybeUninit]) { @@ -632,14 +585,10 @@ pub struct MathishFloorDivide(pub Tolerance); /// [`Tolerant::gcd`], or [`TolerantComplex::complex_gcd`] for complex items, at the tolerance it holds, as a function object. A non-finite argument gives NaN, which marks the item failed. Its `i64` form is [`integer::gcd`](crate::integer::gcd). #[derive(Clone, Copy, Debug, PartialEq)] pub struct Gcd(pub Tolerance); -tolerant!(Gcd at real: |t, [x, y]| t.gcd(x, y)); -tolerant!(Gcd at Complex, Complex: |t, [x, y]| t.complex_gcd(x, y)); /// The least common multiple x y / gcd(x, y), with the sign of x y, at the tolerance it holds, as a function object. Where the gcd is 0, so is the result. Its `i64` form is [`integer::lcm`](crate::integer::lcm). #[derive(Clone, Copy, Debug, PartialEq)] pub struct Lcm(pub Tolerance); -tolerant!(Lcm at real: |t, [x, y]| lcm_from_gcd(x, y, t.gcd(x, y))); -tolerant!(Lcm at Complex, Complex: |t, [x, y]| lcm_from_gcd(x, y, t.complex_gcd(x, y))); #[inline] pub(crate) fn lcm_from_gcd(x: T, y: T, g: T) -> T { if g.is_zero() { g } else { x / g * y } } diff --git a/kernels/src/operation.rs b/kernels/src/operation.rs index 7fd0312..11e8193 100644 --- a/kernels/src/operation.rs +++ b/kernels/src/operation.rs @@ -80,7 +80,7 @@ impl Map { if inputs.len() != self.arity() { output.fill(range.len()); return Err(Error::new(Kind::Length, "map arity does not match its operands")); } match self { Self::Kernel(kernel, settings) => kernel.map_in(inputs, range, *settings, scratch, output), - Self::Convert(conversion, settings, checked) => conversion.map_in(inputs[0], range, *checked, settings.tolerance, scratch, output) + Self::Convert(conversion, settings, checked) => conversion.map_in(inputs[0], range, *checked, settings.tolerance, Some(scratch), output) } } } @@ -656,7 +656,7 @@ fn with_projected_view(view: View<'_>, projection: Projection, f: impl FnOnce fn with_projected( inputs: &[crate::agreement::Read<'_>], range: Option>, projection: Projection, call: impl FnOnce(&[crate::agreement::Read<'_>]) -> Result) -> Result { - use crate::{agreement::{Binding, Read, ReadPlan}, buffer::{FloatBuf, IntBuf}}; + use crate::{agreement::{Binding, Read, ReadPlan}, buffer::{Floats, IntBuf}}; let Some(t) = projection.wholes() else { return call(inputs) }; let codes = matches!(projection, Projection::Codes(_)); let sources: Vec> = if codes { inputs.iter().map(|r| r.sources.iter().map(|s| s.as_codes()).collect()).collect() } else { Vec::new() }; @@ -677,9 +677,8 @@ fn with_projected( let converted = projected.iter().map(|read| { if !read.sources.iter().any(|s| matches!(s.view.dtype(), DType::Float(_))) { return Ok(None); } let range = range.clone().unwrap_or(0..read.plan.len()); - let mut floats = Buffer::with_capacity(DType::Float(FloatWidth::F64), range.len()); - read.read_into(range.clone(), floats.dtype(), true, &mut floats)?; - let Buffer::Floats(FloatBuf::F64(values)) = &floats else { unreachable!("a read into 64-bit floats gives them") }; + let mut floats = Buffer::with_capacity(DType::Float(FloatWidth::F64), 0); + let View::Floats(Floats::F64(values)) = read.part_view(range.clone(), floats.dtype(), true, &mut floats)? else { unreachable!("a read into 64-bit floats gives them") }; Ok(Some((range.start, Buffer::Integers(IntBuf::I64(whole_numbers(t, values)?))))) }).collect::, Error>>()?; let bindings = converted.iter().zip(inputs) diff --git a/kernels/src/operation/map.rs b/kernels/src/operation/map.rs index bc15fbf..35fa3ec 100644 --- a/kernels/src/operation/map.rs +++ b/kernels/src/operation/map.rs @@ -21,17 +21,17 @@ pub(super) fn map_formula, const N: usize>( if let Some(slices) = &slices { let valid = |ok: bool| if judge.checked && !ok { Err(Kind::Domain) } else { Ok(()) }; let blocks = match variant.dtype { - DType::Float(FloatWidth::F64) => sliced(object, &ops, slices.double(), inputs, range.clone(), variant.checked, &valid, output), - DType::Float(FloatWidth::F32) => sliced(object, &ops, slices.single(), inputs, range.clone(), variant.checked, &valid, output), - DType::Float(FloatWidth::F16) => slices.single().map(|single| { + DType::Float(FloatWidth::F64) => sliced(object, &ops, slices.double, inputs, range.clone(), variant.checked, &valid, output), + DType::Float(FloatWidth::F32) => sliced(object, &ops, slices.single, inputs, range.clone(), variant.checked, &valid, output), + DType::Float(FloatWidth::F16) => slices.single.map(|single| { map_blocks::(inputs, range.clone(), variant.checked, output, &|args, n, out| { valid(unsafe { crate::slice::buffered(n, out, |out| single(object, &ops, args, n, out), half::f16::from_f32) }) }) }), - DType::Integer(Width::U8) => sliced(object, &ops, slices.byte(), inputs, range.clone(), variant.checked, &valid, output), - DType::Integer(Width::I16) => sliced(object, &ops, slices.short(), inputs, range.clone(), variant.checked, &valid, output), - DType::Integer(Width::I32) => sliced(object, &ops, slices.int(), inputs, range.clone(), variant.checked, &valid, output), - DType::Integer(Width::I64) => sliced(object, &ops, slices.long(), inputs, range.clone(), variant.checked, &valid, output), + DType::Integer(Width::U8) => sliced(object, &ops, slices.byte, inputs, range.clone(), variant.checked, &valid, output), + DType::Integer(Width::I16) => sliced(object, &ops, slices.short, inputs, range.clone(), variant.checked, &valid, output), + DType::Integer(Width::I32) => sliced(object, &ops, slices.int, inputs, range.clone(), variant.checked, &valid, output), + DType::Integer(Width::I64) => sliced(object, &ops, slices.long, inputs, range.clone(), variant.checked, &valid, output), _ => None, }; if let Some(blocks) = blocks { return blocks; } @@ -45,16 +45,16 @@ pub(super) fn map_formula, const N: usize>( }, DType::Integer(width) => with_width!(width, I => map_integer::(object, variant, inputs, range, judge, output)), DType::Complex(crate::buffer::ComplexWidth::F32) if K::REAL_COMPLEX && inputs.real() => { - map_real_complex::(object, variant, inputs, range, &ops, slices.and_then(|s| s.real_complex_single()), output) + map_real_complex::(object, variant, inputs, range, &ops, slices.and_then(|s| s.real_complex_single), output) } DType::Complex(crate::buffer::ComplexWidth::F64) if K::REAL_COMPLEX && inputs.real() => { - map_real_complex::(object, variant, inputs, range, &ops, slices.and_then(|s| s.real_complex_double()), output) + map_real_complex::(object, variant, inputs, range, &ops, slices.and_then(|s| s.real_complex_double), output) } DType::Complex(crate::buffer::ComplexWidth::F32) => { - map_complex::(object, variant, inputs, range, &ops, slices.and_then(|s| s.complex_single()), judge, output) + map_complex::(object, variant, inputs, range, &ops, slices.as_ref().and_then(|s| s.complex_single), slices.and_then(|s| s.complex_test_single), judge, output) } DType::Complex(crate::buffer::ComplexWidth::F64) => { - map_complex::(object, variant, inputs, range, &ops, slices.and_then(|s| s.complex_double()), judge, output) + map_complex::(object, variant, inputs, range, &ops, slices.as_ref().and_then(|s| s.complex_double), slices.and_then(|s| s.complex_test_double), judge, output) } DType::Character => Err(Error::new(Kind::Domain, "characters are not numerical kernel operands")), } diff --git a/kernels/src/operation/map/complex.rs b/kernels/src/operation/map/complex.rs index 9a32742..8786472 100644 --- a/kernels/src/operation/map/complex.rs +++ b/kernels/src/operation/map/complex.rs @@ -2,11 +2,17 @@ use super::*; pub(super) fn map_complex( object: &K, variant: crate::domain::Variant, inputs: Inputs<'_, '_>, range: std::ops::Range, ops: &NativeOps<{ FUSED }>, - slice: Option, N, ()>>, judge: crate::domain::Variant, output: &mut Output<'_, '_>) -> Result<(), Error> + slice: Option, N, ()>>, test: Option, N, (), bool>>, + judge: crate::domain::Variant, output: &mut Output<'_, '_>) -> Result<(), Error> where K: crate::facts::Formula, F: crate::Native + crate::buffer::Cast + OutputItem, num_complex::Complex: crate::buffer::Cast + OutputItem { let value = |x| object.complex(ops, x).unwrap_or(num_complex::Complex::new(F::zero(), F::zero())); let invalid = |x, r| (!complex_valid::(ops, x, r)).then_some(Kind::Domain); if const { K::BOOLEAN } { + if let Some(test) = test { + return map_blocks::, bool, N>(inputs, range, variant.checked, output, &|args, n, out| { + test(object, ops, args, n, out); Ok(()) + }); + } map_values::, bool, N>( inputs, range, variant.checked, judge.checked, output, |x| value(x).re != F::zero(), |x, r| invalid(x, num_complex::Complex::new(if r { F::one() } else { F::zero() }, F::zero()))) diff --git a/nbs/index.ipynb b/nbs/index.ipynb index 6af2498..ada0840 100644 --- a/nbs/index.ipynb +++ b/nbs/index.ipynb @@ -40,7 +40,7 @@ "pip install rustymath\n", "```\n", "\n", - "Rust users can add the crate to `Cargo.toml`. See [Using rustymath from Rust](https://answerdotai.github.io/rustymath/rust.html) for scalar and slice forms, build settings, accuracy and speed." + "Rust users can add the crate to `Cargo.toml`. See [Using rustymath from Rust](https://answerdotai.github.io/rustymath/rust.html) for function objects and native bulk routines, build settings, accuracy and speed." ] }, { diff --git a/src/array.rs b/src/array.rs index 5514ade..361e6bc 100644 --- a/src/array.rs +++ b/src/array.rs @@ -2,7 +2,7 @@ use crate::agreement::{Agreement, Binding, Mapping, Missing, Read, ReadPlan, ReadRoute}; use crate::buffer::{Buffer, DType, Element, Items, View}; use crate::error::{Error, Kind}; -use crate::gather::{Steps, Writes}; +use crate::gather::Writes; use crate::layout::{Keys, Lanes, Layout, element_count}; use crate::operation::{Kernel, Settings, WriteRoot}; pub use rustymath_core::array::ArrayView; @@ -94,19 +94,13 @@ impl Array { if view.dtype() != source.dtype() { return Err(Error::new(Kind::Domain, "concatenate items must have one dtype")); } layout = crate::structure::catenate(&layout, view.layout(), axis)?; } - let lanes = arrays.iter().map(|a| Lanes::new(a.layout().shape(), axis)).collect::, _>>()?; - let lens: Vec = lanes.iter().map(|l| l.len * l.inner).collect(); - let (outer, stride) = (lanes[0].outer, lens.iter().sum::()); - let buffer = with_dtype!(source.dtype(), T => { - let parts = arrays.iter().map(|a| a.items::()).collect::, _>>()?; - let mut out = parts.iter().find_map(|p| p.first()).map_or_else(Vec::new, |&x| vec![x; outer * stride]); - let mut at = 0; - for (part, &n) in parts.iter().zip(&lens) { - crate::gather::scatter_walk(&mut out, part, at, &[Steps::along(outer, stride), Steps::along(n, 1)], crate::gather::Same); - at += n; - } - T::buffer(out) - }); + let dtype = source.dtype(); + let outer = Lanes::new(source.layout().shape(), axis)?.outer; + let parts = arrays.iter().map(|a| a.cast(dtype)).collect::, _>>()?; + let sources: Vec<_> = parts.iter().map(|part| Binding::full(part.view())).collect(); + let plan = ReadPlan::joined(outer, sources.iter().enumerate().map(|(i, source)| ReadPlan::input(i, source.items)))?; + let mut buffer = Buffer::with_capacity(dtype, layout.count()); + Read { sources: &sources, plan: &plan }.read_into(0..layout.count(), dtype, true, &mut buffer)?; Self::new(layout, buffer) } pub fn choose(condition: ArrayView<'_>, x: ArrayView<'_>, y: ArrayView<'_>) -> Result { diff --git a/src/lib.md b/src/lib.md index a93270b..de1f7e8 100644 --- a/src/lib.md +++ b/src/lib.md @@ -16,7 +16,7 @@ rustflags = ["-C", "target-cpu=x86-64-v3"] Rust's default x86-64 target has no fused multiply-add and no vector instruction for `floor` or `trunc`. Without the setting, rustymath multiplies and adds separately. Some of its loops also run one item at a time. -Each function comes in two forms. The free function, such as `exp`, is generic over `f32`, `f64` and `half::f16`. The compiler inlines it into the loop that calls it and vectorizes that loop: +Free functions such as `exp` are generic over `f32`, `f64` and `half::f16`. The compiler inlines it into the loop that calls it and vectorizes that loop: ```rust use half::f16; @@ -27,16 +27,27 @@ let e32 = exp(1.0f32); // 2.718282 let e16 = exp(f16::from_f32(1.0)); // 2.719 ``` -The function object, such as `Exp` or `Pow`, has a scalar form `one` and a slice form `many`. Both run code that rustymath compiles itself. Your build settings don't change their speed. Each argument of `many` is one item, which every result pairs with, or a slice with one item for each result. `many` appends the results to a vector and returns whether every item succeeded. A NaN result from arguments that aren't NaN counts as a failure. An object implements `Unary`, `Binary` or `Ternary` at each item type it takes: `f16`, `f32` and `f64`, complex numbers where the function has complex values, and integers for some, such as `Gcd` at `i64`. Code that is generic over the float type states the bound it needs, such as `where Exp: Unary`. The free functions need no bound. +For whole arrays, resolve a formula object such as `Pow` to an `operation::Kernel`. `Kernel::map` reads each argument through a `ReadPlan` and replaces the contents of an output `Buffer`. A plan can read a source in order or repeat one item, as the exponent does here: ```rust -use rustymath::{Arg, Binary, Pow}; - -let mut squares = Vec::new(); -Pow.many(Arg::Many(&[1.5, 2.0, 3.0]), Arg::One(2.0), 3, &mut squares); // squares: [2.25, 4.0, 9.0] +use rustymath::{ + Pow, agreement::{Binding, Mapping, Missing, Read, ReadPlan}, buffer::{Buffer, Element}, + domain::Variant, operation::{Kernel, Settings}, tolerant::Tolerance, +}; + +let x = [1.5, 2.0, 3.0]; +let exponent = [2.0]; +let sources = [Binding::full(f64::view(&x)), Binding::full(f64::view(&exponent))]; +let plans = [ReadPlan::input(0, 3), ReadPlan::copy(1, 1, Mapping::Single, 3, Missing::Default)]; +let reads = plans.each_ref().map(|plan| Read { sources: &sources, plan }); +let pow = Kernel::resolve(Pow, Variant::checked(f64::DTYPE), [f64::DTYPE; 2]).unwrap(); +let settings = Settings { tolerance: Tolerance::EXACT }; +let mut squares = Buffer::with_capacity(pow.result_dtype(), 3); +pow.map(&reads, 0..3, settings, &mut squares).unwrap(); +assert_eq!(f64::vec(squares).unwrap(), [2.25, 4.0, 9.0]); ``` -`facts::Formula` holds an object's numerical and kind rules. `Kernel::resolve` selects a working dtype and returns an `operation::Kernel` for callers choosing dtypes at run time. Its `map`, `fold_read`, `fold_pieces` and `scan` methods consume shared `agreement::Read` plans. `Kernel::real`, `complex` and `exact` evaluate scalar arguments through the same object. +`facts::Formula` owns each object's numerical rules, supported kinds, validity, ranges and derivatives. `Kernel::resolve` selects its computing and result dtypes. `Kernel::of` selects a formula from a runtime `facts::Op`. Kernels use native bulk routines where a formula supplies them, and otherwise run its item formula. Their `map`, `fold_read`, `fold_pieces` and `scan` methods consume shared `agreement::Read` plans. `real`, `integers`, `complex` and `exact` evaluate scalar arguments with the kernel's numerical policy and return a `Result`. ```rust use rustymath::{Exp, domain::Variant, operation::Kernel, buffer::{DType, FloatWidth, Width}}; @@ -48,19 +59,23 @@ assert_eq!(exp.variant().dtype, double); `Sin`, `Cos` and `Tan` use short reduction and polynomials for their common range. An object's repair predicate marks arguments beyond 2^23 at f64 or 2^20 at f32. Its native slice driver applies the object's scalar repair routine to those items before returning the results. `Kernel::ranged` omits repair when the same predicate is proved false over the supplied ranges. -`atan2(y, x)` and its object `Atan2` give the angle of the point (x, y), from -π to π. `y` comes first, as in `f64::atan2`. +`atan2(y, x)` and a kernel of `Atan2` give the angle of the point (x, y), from -π to π. `y` comes first, as in `f64::atan2`. `divide(x, y)` divides complex numbers, inline in the caller's loop, as the free functions do. Two real arguments divide as reals, so a real zero divisor gives ±∞ or NaN with an imaginary part of 0. Its formula's |y|² overflows beyond about 1e154 at f64 and 1e19 at f32. ## Tolerant comparison -The `tolerant` module has APL's tolerant comparison: x equals y when |x - y| is less than the tolerance times the larger of |x| and |y|. Its objects, such as `Equal`, `Less`, `Floor`, `Residue` and `Gcd`, each hold a `Tolerance`. With `Tolerance::EXACT`, a slice form runs a loop of plain IEEE comparisons, chosen once per call. +The `tolerant` module has APL's tolerant comparison: x equals y when |x - y| is less than the tolerance times the larger of |x| and |y|. Its formula objects, such as `Equal`, `Less`, `Floor`, `Residue` and `Gcd`, hold a `Tolerance`. `operation::Settings` supplies the tolerance for a kernel execution. `Tolerance::EXACT` selects ordinary comparison once per call. A comparison kernel produces Boolean array results; its scalar `real` method represents true and false as 1 and 0: ```rust -use rustymath::{tolerant::{Equal, Tolerance}, Arg, Binary}; - -let mut same = Vec::new(); -Equal(Tolerance::DEFAULT).many(Arg::Many(&[0.3, 0.1 + 0.2, 0.31]), Arg::One(0.3), 3, &mut same); // same: [true, true, false] +use rustymath::{ + buffer::Element, domain::Variant, operation::{Kernel, Settings}, tolerant::{Equal, Tolerance}, +}; + +let equal = Kernel::resolve(Equal(Tolerance::DEFAULT), Variant::checked(f64::DTYPE), [f64::DTYPE; 2]).unwrap(); +let settings = Settings { tolerance: Tolerance::DEFAULT }; +let same = [0.3, 0.1 + 0.2, 0.31].map(|x| equal.real(&[x, 0.3], settings).unwrap()); +assert_eq!(same, [1.0, 1.0, 0.0]); ``` ## Other modules @@ -90,11 +105,11 @@ Equal(Tolerance::DEFAULT).many(Arg::Many(&[0.3, 0.1 + 0.2, 0.31]), Arg::One(0.3) The special functions `gamma`, `lgamma`, `digamma`, `erf`, `erfc`, their inverses, `binom` with its slopes `binom_dn` and `binom_dk`, the incomplete gamma and beta functions and Lambert W have free functions and objects at the crate root, as the elementary functions do. -`Add`, `Subtract`, `Multiply` and the other arithmetic objects at the crate root have no scalar or slice forms. `Kernel` runs them at every float, complex and integer type, with its explicit `Variant` checking and wrapping policy. Unbounded rational callers use `Formula::exact` through `Kernel::exact`. +`Add`, `Subtract`, `Multiply` and the other arithmetic objects at the crate root implement `Formula`. `Kernel` runs them at every float, complex and integer type, with its explicit `Variant` checking and wrapping policy. Unbounded rational callers use `Formula::exact` through `Kernel::exact`. ## Accuracy and speed -The table gives nanoseconds per item for a loop over 65,536 inputs on an Apple M5 Max, in October 2026, beside the standard library's version of each function. A row named after an object times its slice form. `Pow.one` and `Root.one` time the scalar forms. To print the table on your own machine, run `RUSTYMATH_COSTS=1 cargo test --test speed`. +The table gives nanoseconds per item for a loop over 65,536 inputs on an Apple M5 Max, in October 2026, beside the standard library's version of each function. A row named after an object times its bulk kernel. The labels `Pow.one` and `Root.one` time scalar formulas pair by pair. To print the table on your own machine, run `RUSTYMATH_COSTS=1 cargo test --test speed`. | Function | f64 | std f64 | f32 | std f32 | |---|---|---|---|---| diff --git a/src/lib.rs b/src/lib.rs index 2c16d2d..e2dfc4e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -10,7 +10,7 @@ pub use rustymath_core::{ fused, unfused, half, num_complex, }; pub use rustymath_core::{ - FUSED, Float, Native, Arg, Unary, Binary, Ternary, + FUSED, Float, Native, Arg, Abs, Acos, Acosh, Add, Angle, Asin, Asinh, Atan, Atan2, Atanh, Betainc, BetaincAt, Betaincinv, Binom, BinomDk, BinomDn, BitAnd, BitOr, Bitwise, Cbrt, Choose, Cis, Clip, Conj, Cos, Cosh, Digamma, Divide, Erf, Erfc, Erfcinv, Erfinv, Exp, Expm1, Factorial, Gamma, Gammainc, Gammaincc, Gammaincinv, Hypot, Imag, Isfinite, Isinf, Isnan, LambertW, LeftShift, Lgamma, Ln, Log1p, Log2, diff --git a/tests/accuracy.rs b/tests/accuracy.rs index c1ad28e..b0a9be1 100644 --- a/tests/accuracy.rs +++ b/tests/accuracy.rs @@ -462,7 +462,7 @@ fn error_functions_accuracy() { /// `Erfc`'s slice form gives `erfc`'s bits at both widths. Its blocks of 256 items lie wholly below 0.84375 (with NaN), wholly from it up (with -∞ and magnitudes past 28), on both sides, and in a short block of zeros and ∞. #[test] fn erfc_blocks() { - use rustymath::{Arg, Erfc, Unary}; + use rustymath::{Erfc, buffer::Element, domain::Variant, operation::Kernel}; let mut xs = vec![f64::NAN]; xs.extend((1..256).map(|i| -0.84 + 1.68 * f64::from(i) / 256.0)); xs.push(f64::NEG_INFINITY); @@ -470,8 +470,8 @@ fn erfc_blocks() { xs.extend((0..256).map(|i| if i % 2 == 0 { f64::from(i) / 1024.0 } else { -0.9 - f64::from(i) / 16.0 })); xs.extend([0.0, -0.0, f64::INFINITY]); let xs32: Vec = xs.iter().map(|&x| x as f32).collect(); - let (mut out, mut out32) = (Vec::new(), Vec::new()); - assert!(Erfc.many(Arg::Many(&xs), xs.len(), &mut out) & Erfc.many(Arg::Many(&xs32), xs32.len(), &mut out32)); + let out: Vec = crate::support::mapped(&Kernel::new(Erfc, Variant::checked(f64::DTYPE)), [&xs[..]], xs.len()).unwrap(); + let out32: Vec = crate::support::mapped(&Kernel::new(Erfc, Variant::checked(f32::DTYPE)), [&xs32[..]], xs32.len()).unwrap(); let same = |a: f64, b: f64| a.to_bits() == b.to_bits() || (a.is_nan() && b.is_nan()); assert!(out.iter().zip(&xs).all(|(&r, &x)| same(r, rustymath::erfc(x))), "{out:?}"); assert!(out32.iter().zip(&xs32).all(|(&r, &x)| same(r.into(), rustymath::erfc(x).into())), "{out32:?}"); diff --git a/tests/linalg.rs b/tests/linalg.rs index 2de7919..c2aa6b0 100644 --- a/tests/linalg.rs +++ b/tests/linalg.rs @@ -336,10 +336,7 @@ mod factor { pub(super) fn wide(x: &[T]) -> Vec { x.iter().map(|v| v.wide()).collect() } - /// The product of the `m`×`k` matrix `a` and the `k`×`n` matrix `b`, in f64, item by item. - pub(super) fn product(a: &[f64], b: &[f64], [m, k, n]: [usize; 3]) -> Vec { - (0..m * n).map(|ij| (0..k).map(|p| a[ij / n * k + p] * b[p * n + ij % n]).sum()).collect() - } + pub(super) fn product(a: &[f64], b: &[f64], size: [usize; 3]) -> Vec { PlusTimes.inner(a, b, size) } pub(super) fn norm(x: &[f64]) -> f64 { x.iter().map(|v| v * v).sum::().sqrt() } diff --git a/tests/numeric.rs b/tests/numeric.rs index 4c94417..e1469a4 100644 --- a/tests/numeric.rs +++ b/tests/numeric.rs @@ -118,7 +118,14 @@ mod convert { mod distribution { use rustymath::distribution::*; - use rustymath::{Arg, Unary}; + use rustymath::{buffer::{Buffer, Element}, domain::Variant, facts::Formula, operation::Kernel}; + + fn many>(f: K, items: &[F]) -> (bool, Vec) { + let kernel = Kernel::new(f, Variant::checked(F::DTYPE)); + let mut out = Buffer::with_capacity(kernel.result_dtype(), items.len()); + let ok = rustymath::error::checked(|| crate::support::map(&kernel, [F::view(items)], items.len(), &mut out)).is_ok(); + (ok, F::vec(out).unwrap()) + } /// One of each distribution, with each sampling method's case: binomials that draw by inversion, by BTRD and reflected, and Poisson distributions that draw by multiplying and by PTRS. A second Weibull distribution has an odd whole shape, whose root of a negative number is finite, and a density whose power overflows in its upper tail. fn cases() -> Vec { @@ -177,10 +184,6 @@ mod distribution { /// Each distribution's density and cdf succeed at every x, and its quantile at every p from 0 to 1. A quantile fails at a p outside 0 to 1, where it gives NaN. #[test] fn domains() { - fn many(f: impl Unary, items: &[f64]) -> (bool, Vec) { - let mut out = Vec::new(); - (f.many(Arg::Many(items), items.len(), &mut out), out) - } let (xs, ps) = ([-1.0, 0.0, 0.25, 0.5, 1.0, 3.0, 100.0, f64::NAN], [0.001, 0.25, 0.5, 0.999]); for d in cases() { assert!(many(Pdf(d), &xs).0 && many(Cdf(d), &xs).0 && many(Quantile(d), &ps).0, "{d:?}"); @@ -190,9 +193,10 @@ mod distribution { } /// The worst relative error of `f`'s f32 slice form against `wide`, the same law at f64, over `xs`. The error is relative where the f64 result rounds to a normal f32 number. Elsewhere the f32 result must be a NaN beside a NaN, the same infinity, or at most f32's least normal number in magnitude where the f64 result rounds to zero or a subnormal. A result outside its class, or a success verdict other than the one the f64 results imply, counts as an infinite error. - fn worst_at_f32(f: impl Unary, wide: impl Fn(f64) -> f64, xs: &[f32]) -> f64 { - let (mut out, wants): (Vec, Vec) = (Vec::new(), xs.iter().map(|&x| wide(x.into())).collect()); - if f.many(Arg::Many(xs), xs.len(), &mut out) != xs.iter().zip(&wants).all(|(x, w)| !w.is_nan() || x.is_nan()) { return f64::INFINITY; } + fn worst_at_f32(f: impl Formula<1>, wide: impl Fn(f64) -> f64, xs: &[f32]) -> f64 { + let wants: Vec = xs.iter().map(|&x| wide(x.into())).collect(); + let (ok, out) = many(f, xs); + if ok != xs.iter().zip(&wants).all(|(x, w)| !w.is_nan() || x.is_nan()) { return f64::INFINITY; } wants.iter().zip(&out).map(|(&want, &r)| { let (near, r) = (want as f32, f64::from(r)); let same = r == f64::from(near) || (r.is_nan() && near.is_nan()) || (near.abs() < f32::MIN_POSITIVE && r.abs() <= f64::from(f32::MIN_POSITIVE)); @@ -246,7 +250,7 @@ mod integer { use num_bigint::BigInt; use rustymath::integer::*; use rustymath::{ - Abs, Arg, Binary, Binom, Pow, buffer::{DType, Element}, domain::Variant, operation::Kernel, + Abs, Binom, Pow, buffer::{DType, Element}, domain::Variant, operation::Kernel, tolerant::{Ceil, Equal, Floor, FloorDivide, Gcd, Greater, GreaterEqual, Lcm, LessEqual, NotEqual, Residue} }; @@ -348,14 +352,12 @@ mod integer { // The integer forms ignore the tolerance. let t = rustymath::tolerant::Tolerance::DEFAULT; for d in [7, -12, 1, -1, 0, i64::MIN] { - let mut out = Vec::new(); - assert!(Residue(t).many(Arg::One(d), Arg::Many(&ys), ys.len(), &mut out)); + let out: Vec = mapped(&Kernel::new(Residue(t), Variant::checked(i64::DTYPE)), [&[d], &ys[..]], ys.len()).unwrap(); assert_eq!(out, ys.iter().map(|&y| residue(d, y)).collect::>()); } let ys16: Vec = (i16::MIN..=i16::MAX).collect(); for d in [7i16, -12, 1, -1, 0, 3, i16::MAX, i16::MIN] { - let mut out = Vec::new(); - assert!(Residue(t).many(Arg::One(d), Arg::Many(&ys16), ys16.len(), &mut out)); + let out: Vec = mapped(&Kernel::new(Residue(t), Variant::checked(i16::DTYPE)), [&[d], &ys16[..]], ys16.len()).unwrap(); assert_eq!(out, ys16.iter().map(|&y| residue(d.into(), y.into()) as i16).collect::>()); } // Every `u8` and `i16` divisor of magnitude 2 or more divides through a 16-bit high product. A quotient too large would show first one below a multiple of the divisor's magnitude, or at 2¹⁵, the largest dividend's magnitude. A quotient too small would show first at a multiple. These dividends therefore test every dividend. @@ -363,14 +365,12 @@ mod integer { let m = i32::from(d.unsigned_abs()); let ys: Vec = (1..).flat_map(|q| [q * m - 1, q * m]).take_while(|&a| a < 1 << 15).flat_map(|a| [a as i16, -a as i16]).chain([i16::MIN]).collect(); - let mut out = Vec::new(); - assert!(Residue(t).many(Arg::One(d), Arg::Many(&ys), ys.len(), &mut out)); + let out: Vec = mapped(&Kernel::new(Residue(t), Variant::checked(i16::DTYPE)), [&[d], &ys[..]], ys.len()).unwrap(); assert_eq!(out, ys.iter().map(|&y| residue(d.into(), y.into()) as i16).collect::>(), "divisor {d}"); } let ys8: Vec = (0..=u8::MAX).collect(); for d in 2..=u8::MAX { - let mut out = Vec::new(); - assert!(Residue(t).many(Arg::One(d), Arg::Many(&ys8), ys8.len(), &mut out)); + let out: Vec = mapped(&Kernel::new(Residue(t), Variant::checked(u8::DTYPE)), [&[d], &ys8[..]], ys8.len()).unwrap(); assert_eq!(out, ys8.iter().map(|&y| residue(d.into(), y.into()) as u8).collect::>(), "divisor {d}"); } // Below 2⁵¹ in magnitude, 64-bit dividends take the residue by one divisor through its reciprocal, and every `i32` dividend does. Dividends beside multiples of each divisor, and beside 2⁵¹, test the step that corrects the quotient. @@ -381,38 +381,40 @@ mod integer { if let Some(m) = k.checked_mul(d) { ys.extend([-1, 0, 1].into_iter().flat_map(|e| [m + e, e - m])); } } ys.retain(|y| y.abs() <= limit); - let mut out = Vec::new(); - assert!(Residue(t).many(Arg::One(d), Arg::Many(&ys), ys.len(), &mut out)); + let out: Vec = mapped(&Kernel::new(Residue(t), Variant::checked(i64::DTYPE)), [&[d], &ys[..]], ys.len()).unwrap(); assert_eq!(out, ys.iter().map(|&y| residue(d, y)).collect::>(), "divisor {d}"); } let ys32: Vec = samples().into_iter().map(|v| v as i32).chain([i32::MIN, i32::MIN + 1, i32::MAX]).collect(); for d in [3i32, -3, 7, -12, 1, -1, i32::MAX, i32::MIN] { - let mut out = Vec::new(); - assert!(Residue(t).many(Arg::One(d), Arg::Many(&ys32), ys32.len(), &mut out)); + let out: Vec = mapped(&Kernel::new(Residue(t), Variant::checked(i32::DTYPE)), [&[d], &ys32[..]], ys32.len()).unwrap(); assert_eq!(out, ys32.iter().map(|&y| residue(d.into(), y.into()) as i32).collect::>(), "divisor {d}"); } - let mut out = vec![5]; - let ok = Gcd(t).many(Arg::Many(&xs), Arg::Many(&ys), xs.len(), &mut out); + let mut out = i64::buffer(vec![5]); + let gcd_kernel = Kernel::new(Gcd(t), Variant::checked(i64::DTYPE)); + let ok = rustymath::error::checked(|| crate::support::map(&gcd_kernel, [i64::view(&xs), i64::view(&ys)], xs.len(), &mut out)).is_ok(); assert_eq!(ok, xs.iter().zip(&ys).all(|(&x, &y)| gcd(x, y).is_some())); - assert_eq!(out[0], 5); - assert_eq!(Binary::::checked(&Pow, 3, 4), Some(81)); - assert_eq!(Binary::::checked(&Pow, 3, 40), None); - assert_eq!(Binom.one(5i64, 2), 10); - assert_eq!(Lcm(t).one(-4i64, 6), -12); + assert_eq!(i64::vec(out).unwrap().len(), xs.len()); + let power = Kernel::new(Pow, Variant::checked(i64::DTYPE)); + assert_eq!(power.integers(&[3, 4]).ok(), Some(81)); + assert!(power.integers(&[3, 40]).is_err()); + assert_eq!(Kernel::new(Binom, Variant::checked(i64::DTYPE)).integers(&[5, 2]).unwrap(), 10); + assert_eq!(Kernel::new(Lcm(t), Variant::checked(i64::DTYPE)).integers(&[-4, 6]).unwrap(), -12); } #[test] fn other_widths() { let t = rustymath::tolerant::Tolerance::DEFAULT; // Every pair of `u8`s agrees with the `i64` functions where the result fits, and fails where it doesn't. + let v = Variant::checked(u8::DTYPE); + let (gcd_kernel, lcm_kernel, pow_kernel, residue_kernel) = (Kernel::new(Gcd(t), v), Kernel::new(Lcm(t), v), Kernel::new(Pow, v), Kernel::new(Residue(t), v)); for x in 0..=u8::MAX { for y in 0..=u8::MAX { let (a, b) = (i64::from(x), i64::from(y)); let fit = |r: Option| r.and_then(|r| u8::try_from(r).ok()); - assert_eq!(Binary::::checked(&Gcd(t), x, y), fit(gcd(a, b)), "gcd({x}, {y})"); - assert_eq!(Binary::::checked(&Lcm(t), x, y), fit(lcm(a, b)), "lcm({x}, {y})"); - assert_eq!(Binary::::checked(&Pow, x, y), fit(power(a, b)), "power({x}, {y})"); - assert_eq!(i64::from(Residue(t).one(x, y)), residue(a, b), "residue({x}, {y})"); + assert_eq!(gcd_kernel.integers(&[a, b]).ok().map(|r| r as u8), fit(gcd(a, b)), "gcd({x}, {y})"); + assert_eq!(lcm_kernel.integers(&[a, b]).ok().map(|r| r as u8), fit(lcm(a, b)), "lcm({x}, {y})"); + assert_eq!(pow_kernel.integers(&[a, b]).ok().map(|r| r as u8), fit(power(a, b)), "power({x}, {y})"); + assert_eq!(residue_kernel.integers(&[a, b]).unwrap(), residue(a, b), "residue({x}, {y})"); } } // Floor division of every `u8` dividend by every nonzero `u8` divisor, through each width's checked kernel. @@ -557,7 +559,7 @@ mod number_theory { mod operation { use rustymath::{ agreement::{Agreement, Binding, Mapping, Missing, Read, ReadPlan}, array::{Array, ArrayView}, buffer::{Buffer, DType, Element, Width}, - domain::{Range, Variant}, error::{self, Error, ErrorSlot}, facts::Formula, gather::Steps, layout::Layout, operation::{Kernel, Settings}, + domain::{Range, Variant}, error::{self, Error, ErrorSlot}, facts::{self, Formula}, gather::Steps, layout::Layout, operation::{Kernel, Settings}, reduce::{Direction, Seed} }; @@ -584,6 +586,14 @@ mod operation { let map = rustymath::operation::Map::Kernel(kernel.clone(), Settings { tolerance: Tolerance::EXACT }); Slots::append(values, n, |slots| map.run(&reads, 0..n, &mut Scratch::default(), MapOutput::F64(slots))) } + fn convert(data: rustymath::buffer::View<'_>, dtype: DType, values: &mut Buffer) -> Result<(), Error> { + let sources = [Binding::full(data)]; + let plan = ReadPlan::input(0, data.len()); + let map = rustymath::operation::Map::Convert(rustymath::convert::Conversion { + input: data.dtype(), output: dtype, kind: rustymath::convert::ConversionKind::Cast + }, Settings { tolerance: Tolerance::EXACT }, true); + map.run(&[Read { sources: &sources, plan: &plan }], 0..data.len(), &mut Scratch::default(), MapOutput::Buffer(values)) + } let float = f64::DTYPE; let add = Kernel::resolve(Add, Variant::checked(float), [float; 2]).unwrap(); let row_shape = Layout::from([2, 32]); @@ -609,6 +619,13 @@ mod operation { assert_eq!(i16::vec(Buffer::floored(data, i16::DTYPE, Tolerance::EXACT).unwrap()).unwrap(), [-2, 2]); assert_eq!(i16::vec(Buffer::cast_results(data, i16::DTYPE, true).unwrap()).unwrap(), [-1, 2]); assert_eq!(u8::vec(Buffer::floored(i64::view(&[300]), u8::DTYPE, Tolerance::EXACT).unwrap()).unwrap(), [44]); + let precise = f64::view(&[1. + f64::EPSILON, 2.5]); + let mut narrow = f32::buffer(vec![9.]); + convert(precise, f32::DTYPE, &mut narrow).unwrap(); + assert_eq!(f32::of(narrow.view()).unwrap(), [9., 1., 2.5]); + let mut wide = f64::buffer(vec![9.]); + convert(precise, f32::DTYPE, &mut wide).unwrap(); + assert_eq!(f64::of(wide.view()).unwrap(), [9., 1., 2.5]); let errors = rustymath::error::ErrorSlot::default(); errors.within(|| { assert!(Buffer::kept(data, i16::DTYPE).is_ok()); @@ -617,12 +634,17 @@ mod operation { assert_eq!(errors.take().unwrap().kind, rustymath::error::Kind::Overflow); assert!(Buffer::cast_results(i64::view(&[0xd800]), char::DTYPE, true).is_ok()); assert_eq!(errors.take().unwrap().kind, rustymath::error::Kind::Domain); + let mut partial = f64::buffer(vec![9.]); + let complex = [num_complex::Complex64::new(1., 0.), num_complex::Complex64::new(2., 1.)]; + assert!(convert(num_complex::Complex64::view(&complex), float, &mut partial).is_ok()); + assert_eq!((partial.len(), f64::of(partial.view()).unwrap()[0]), (3, 9.)); + assert_eq!(errors.take().unwrap().kind, rustymath::error::Kind::Domain); }); } #[test] fn native_operations_use_input_kinds_and_run_settings() { - use rustymath::{Add, Binary, Binom, Bitwise, Exp, Subtract, tolerant::{Residue, Tolerance}}; + use rustymath::{Add, Binom, Bitwise, Exp, Subtract, tolerant::{Residue, Tolerance}}; let (small, wide) = (DType::Integer(Width::I16), DType::Integer(Width::I64)); let exp = Kernel::resolve(Exp, Variant::checked(wide), [wide]).unwrap(); let mut out = Buffer::with_capacity(exp.result_dtype(), 2); @@ -677,11 +699,11 @@ mod operation { let values = [[3.], [3. - 1e-14]]; let bindings = values.each_ref().map(|x| Binding::full(f64::view(x))); let reads = bindings.each_ref().map(|x| Read { sources: std::slice::from_ref(x), plan: &scalar }); - assert_ne!(Residue(Tolerance::EXACT).one(3., values[1][0]), Residue(Tolerance::DEFAULT).one(3., values[1][0])); + assert_ne!(facts::formula_at(&Residue(Tolerance::EXACT), [3., values[1][0]]), facts::formula_at(&Residue(Tolerance::DEFAULT), [3., values[1][0]])); for tolerance in [Tolerance::EXACT, Tolerance::DEFAULT] { let mut out = Buffer::with_capacity(dtype, 1); residue.map(&reads, 0..1, Settings { tolerance }, &mut out).unwrap(); - assert_eq!(f64::vec(out).unwrap(), [Residue(tolerance).one(3., values[1][0])]); + assert_eq!(f64::vec(out).unwrap(), [facts::formula_at(&Residue(tolerance), [3., values[1][0]])]); } for table in 0..16 { let kernel = Kernel::resolve(Bitwise(table), Variant::checked(DType::Boolean), [DType::Boolean; 2]).unwrap(); @@ -759,7 +781,7 @@ mod operation { #[test] fn folds_keep_uncalled_values_and_scans_keep_their_carry() { - use rustymath::{Add, Binary, tolerant::{Less, Tolerance}}; + use rustymath::{Add, tolerant::{Less, Tolerance, Tolerant}}; let dtype = DType::Integer(Width::I64); let less = Kernel::resolve(Less(Tolerance::EXACT), Variant::checked(dtype), [dtype; 2]).unwrap(); let input = [Binding::full(i64::view(&[3, 2, 1]))]; @@ -812,7 +834,7 @@ mod operation { let (mut carry, mut expected) = (initial, Vec::new()); for x in ordered { let value = carry.map_or(x, |c| { - u8::from(if direction == Direction::Left { Less(tolerance).one(c, x) } else { Less(tolerance).one(x, c) }) as f64 + u8::from(if direction == Direction::Left { tolerance.less(c, x) } else { tolerance.less(x, c) }) as f64 }); expected.push(value); carry = Some(value); @@ -842,6 +864,7 @@ mod operation { let extended = view(&extended, f64::view(&data), &Mapping::Linear(1)); assert_eq!(extended.read().plan.contiguous(1..4), Some((0, 1..4))); assert_eq!(extended.read().plan.contiguous(1..1), None); + assert!(view(&shape, f64::view(&[]), &Mapping::Single).read().mapped_source().is_none()); let fallback = ReadPlan::copy(1, 0, Mapping::Linear(1), 2, Missing::Input(ReadPlan::input(0, 1))); assert!(fallback.contiguous(0..2).is_none()); assert_eq!(extended.items::().unwrap().as_ref(), &[0., 1., 2., 3., 4., 5., 0., 0.]); @@ -973,6 +996,8 @@ mod operation { kernel.map(&[chars.read(), offsets.read()], 1..2, Settings { tolerance: rustymath::tolerant::Tolerance::DEFAULT }, &mut output).unwrap(); assert_eq!(char::of(output.view()), Some(&['c'][..])); assert!(kernel.map(&[chars.read(), offsets.read()], 0..1, Settings { tolerance: rustymath::tolerant::Tolerance::DEFAULT }, &mut output).is_err()); + let repeated = view(&layout, f64::view(&[1.]), &Mapping::Single); + assert_eq!(char::vec(apply(&kernel, [ArrayView::compact(&layout, ab.view()), repeated]).unwrap()), Some(vec!['b', 'c'])); let seed = f64::buffer(vec![0.5]); let seed = ArrayView::compact(&Layout::EMPTY, seed.view()); assert!( diff --git a/tests/slices.rs b/tests/slices.rs index 15c443c..7156d06 100644 --- a/tests/slices.rs +++ b/tests/slices.rs @@ -183,14 +183,14 @@ mod forms { #[test] fn outputs_can_differ_from_items() { use rustymath::tolerant::{Greater, Tolerance}; - fn positive(x: &[F]) -> Vec where Greater: Binary { - let mut out = Vec::new(); - assert!(Greater(Tolerance::EXACT).many(Arg::Many(x), Arg::One(F::zero()), x.len(), &mut out)); - out + fn positive(x: &[F]) -> Vec { + let k = Kernel::new(Greater(Tolerance::EXACT), domain::Variant::checked(F::DTYPE)); + crate::support::mapped(&k, [x, &[F::zero()]], x.len()).unwrap() } assert_eq!(positive(&[1.0f64, -1.0, f64::NAN]), [true, false, false]); assert_eq!(positive(&[f16::from_f32(2.0), f16::NAN]), [true, false]); - assert_eq!(Greater(Tolerance::EXACT).checked(f32::NAN, 0.0), Some(false)); + let k = Kernel::new(Greater(Tolerance::EXACT), domain::Variant::checked(::DTYPE)); + assert_eq!(k.real(&[f64::NAN, 0.0], operation::Settings { tolerance: Tolerance::EXACT }).unwrap(), 0.0); } } @@ -796,12 +796,19 @@ mod tolerant { use half::f16; use num_complex::Complex64; use num_traits::Float as _; - use rustymath::{Arg, Binary, Float, Unary, search, tolerant::*, with_tolerance}; + use rustymath::{Float, buffer::{Buffer, Element}, domain::Variant, facts::Formula, operation::{Kernel, Settings}, search, tolerant::*, with_tolerance}; use std::cmp::Ordering; const T: Tolerance = Tolerance::DEFAULT; fn c(re: f64, im: f64) -> Complex64 { Complex64::new(re, im) } + fn many, const N: usize>(k: K, t: Tolerance, args: [&[T]; N], n: usize) -> Vec { + let kernel = Kernel::new(k, Variant { checked: false, ..Variant::checked(T::DTYPE) }); + let mut out = Buffer::with_capacity(kernel.result_dtype(), n); + crate::support::map_with(&kernel, args.map(T::view), n, &mut out, Settings { tolerance: t }).unwrap(); + O::vec(out).unwrap() + } + #[test] fn equal_and_match() { // The rounded pair 0.3 and 0.1+0.2 lies inside the default tolerance, and 1 and 1+3e-14 outside it. @@ -862,7 +869,7 @@ mod tolerant { let f64s = DType::Float(FloatWidth::F64); let quotient = Kernel::resolve(FloorDivide(Tolerance::EXACT), Variant::checked(f64s), [f64s; 2]).unwrap(); assert_eq!( - (quotient.real(&[x, y], rustymath::operation::Settings { tolerance: Tolerance::EXACT }).unwrap(), Residue(Tolerance::EXACT).one(y, x)), + (quotient.real(&[x, y], rustymath::operation::Settings { tolerance: Tolerance::EXACT }).unwrap(), Exactly.residue(y, x)), (2.0, x - 2.0 * y) ); if rustymath::FUSED { assert_ne!(Exactly.residue(0.1, 1.0), 0.0); } @@ -899,10 +906,11 @@ mod tolerant { // The slice forms decide most pairs from squared magnitudes, in a loop with no branch, then redo the rest. let zs = [c(2.0, 1.0), huge, tiny, c(0.0, 0.0), c(f64::INFINITY, 1.0), c(f64::NAN, 0.0), c(-1e300, 0.0)]; let (xs, ys): (Vec<_>, Vec<_>) = zs.iter().flat_map(|&x| zs.iter().map(move |&y| (x, y * (1.0 + 1e-15)))).unzip(); - let (mut equal, mut unequal) = (Vec::new(), Vec::new()); - assert!(Equal(T).many(Arg::Many(&xs), Arg::Many(&ys), xs.len(), &mut equal)); - assert!(NotEqual(T).many(Arg::Many(&xs), Arg::Many(&ys), xs.len(), &mut unequal)); - for i in 0..xs.len() { assert!(equal[i] == T.complex_equal(xs[i], ys[i]) && unequal[i] != equal[i], "{:?} {:?}", xs[i], ys[i]) } + for t in [T, Tolerance::EXACT] { + let equal: Vec = many(Equal(T), t, [&xs[..], &ys[..]], xs.len()); + let unequal: Vec = many(NotEqual(T), t, [&xs[..], &ys[..]], xs.len()); + for i in 0..xs.len() { assert!(equal[i] == t.complex_equal(xs[i], ys[i]) && unequal[i] != equal[i], "{:?} {:?}", xs[i], ys[i]) } + } } /// `range(c)` gives the last floats equal to `c` on each side. @@ -929,104 +937,71 @@ mod tolerant { assert_eq!(T.range(f16::from_f32(0.3)), Some((f16::from_f32(0.3), f16::from_f32(0.3)))); } - /// Checks each real comparison's slice form against its scalar form for many items against one number, in either order, with each item of `items` as the number. - fn against_one(t: Tolerance, items: &[f64]) where Equal: Binary, NotEqual: Binary, - Less: Binary, LessEqual: Binary, Greater: Binary, GreaterEqual: Binary, - Matches: Binary { - fn check + std::fmt::Debug>(k: K, x: &[F]) { + /// Native comparison kernels keep broadcast comparisons and each scalar tolerance method in agreement. + fn against_one(t: Tolerance, items: &[f64]) { + fn check>(k: K, t: Tolerance, x: &[F], one: impl Fn(F, F) -> bool) { for &c in x { - let (mut left, mut right) = (Vec::new(), Vec::new()); - assert!(k.many(Arg::Many(x), Arg::One(c), x.len(), &mut left) & k.many(Arg::One(c), Arg::Many(x), x.len(), &mut right)); - assert_eq!(left, x.iter().map(|&v| k.one(v, c)).collect::>(), "{k:?} with {c:?} on the right"); - assert_eq!(right, x.iter().map(|&v| k.one(c, v)).collect::>(), "{k:?} with {c:?} on the left"); + let left: Vec = many(k.clone(), t, [x, &[c]], x.len()); + let right: Vec = many(k.clone(), t, [&[c], x], x.len()); + assert_eq!(left, x.iter().map(|&v| one(v, c)).collect::>()); + assert_eq!(right, x.iter().map(|&v| one(c, v)).collect::>()); } } let x: Vec = items.iter().map(|&v| F::narrow(v)).collect(); - check(Equal(t), &x); - check(NotEqual(t), &x); - check(Less(t), &x); - check(LessEqual(t), &x); - check(Greater(t), &x); - check(GreaterEqual(t), &x); - check(Matches(t), &x); + check(Equal(t), t, &x, |a, b| t.equal(a, b)); + check(NotEqual(t), t, &x, |a, b| !t.equal(a, b)); + check(Less(t), t, &x, |a, b| t.less(a, b)); + check(LessEqual(t), t, &x, |a, b| t.less_equal(a, b)); + check(Greater(t), t, &x, |a, b| t.greater(a, b)); + check(GreaterEqual(t), t, &x, |a, b| t.greater_equal(a, b)); } - /// Each tolerant object gives its method's result for each item, at the default tolerance and at 0, for slices beside slices and one item beside a slice, at f32, f64 and complex items. Only `Gcd`'s infinite argument fails. + + /// Tolerant kernels give the tolerance methods' results for native arrays, broadcasts and complex items. #[test] fn objects() { let (xs, ys) = (&[3.0, -3.0, 0.3, 0.0, f64::NAN, 2.0][..], &[-5.0, 5.0, 0.1 + 0.2, 7.0, 1.0, 2.9999999999999996][..]); let agree = |a: &[f64], b: &[f64]| a.len() == b.len() && a.iter().zip(b).all(|(x, y)| x == y || x.is_nan() && y.is_nan()); for t in [T, Tolerance::EXACT] { - let (mut floors, mut residues, mut equal, mut order) = (Vec::new(), Vec::new(), Vec::new(), Vec::new()); - assert!(Floor(t).many(Arg::Many(ys), 6, &mut floors) & Ceil(t).many(Arg::Many(ys), 6, &mut floors)); - assert!(Residue(t).many(Arg::Many(xs), Arg::Many(ys), 6, &mut residues) & Residue(t).many(Arg::One(3.0), Arg::Many(ys), 6, &mut residues)); - assert!(Equal(t).many(Arg::Many(xs), Arg::Many(ys), 6, &mut equal) & Matches(t).many(Arg::Many(xs), Arg::One(f64::NAN), 6, &mut equal)); - assert!(Compare(t).many(Arg::Many(xs), Arg::Many(ys), 6, &mut order)); + let mut floors: Vec = many(Floor(t), t, [ys], 6); + floors.extend(many::<_, f64, _, 1>(Ceil(t), t, [ys], 6)); + let mut residues: Vec = many(Residue(t), t, [xs, ys], 6); + residues.extend(many::<_, f64, _, 2>(Residue(t), t, [&[3.0], ys], 6)); let floor: Vec = ys.iter().map(|&y| t.floor(y)).chain(ys.iter().map(|&y| t.ceil(y))).collect(); let residue: Vec = (0..6).map(|i| t.residue(xs[i], ys[i])).chain(ys.iter().map(|&y| t.residue(3.0, y))).collect(); - let eq: Vec = (0..6).map(|i| t.equal(xs[i], ys[i])).chain(xs.iter().map(|&x| t.matches(x, f64::NAN))).collect(); assert!(agree(&floors, &floor) && agree(&residues, &residue), "{floors:?} {residues:?}"); - assert_eq!(equal, eq); - assert_eq!(order, (0..6).map(|i| t.compare(xs[i], ys[i]).unwrap_or(Ordering::Greater)).collect::>()); - let check = |many: &dyn Fn(&mut Vec) -> bool, one: &dyn Fn(f64, f64) -> bool| { - let mut got = Vec::new(); - assert!(many(&mut got)); + let check = |got: Vec, one: &dyn Fn(f64, f64) -> bool| { assert_eq!(got, (0..6).map(|i| one(xs[i], ys[i])).collect::>()); }; - check(&|o| NotEqual(t).many(Arg::Many(xs), Arg::Many(ys), 6, o), &|x, y| !t.equal(x, y)); - check(&|o| Less(t).many(Arg::Many(xs), Arg::Many(ys), 6, o), &|x, y| t.less(x, y)); - check(&|o| LessEqual(t).many(Arg::Many(xs), Arg::Many(ys), 6, o), &|x, y| t.less_equal(x, y)); - check(&|o| Greater(t).many(Arg::Many(xs), Arg::Many(ys), 6, o), &|x, y| t.greater(x, y)); - check(&|o| GreaterEqual(t).many(Arg::Many(xs), Arg::Many(ys), 6, o), &|x, y| t.greater_equal(x, y)); - let (ys32, mut floors32) = (ys.iter().map(|&y| y as f32).collect::>(), Vec::new()); - assert!(Floor(t).many(Arg::Many(&ys32), 6, &mut floors32)); + check(many(Equal(t), t, [xs, ys], 6), &|x, y| t.equal(x, y)); + check(many(NotEqual(t), t, [xs, ys], 6), &|x, y| !t.equal(x, y)); + check(many(Less(t), t, [xs, ys], 6), &|x, y| t.less(x, y)); + check(many(LessEqual(t), t, [xs, ys], 6), &|x, y| t.less_equal(x, y)); + check(many(Greater(t), t, [xs, ys], 6), &|x, y| t.greater(x, y)); + check(many(GreaterEqual(t), t, [xs, ys], 6), &|x, y| t.greater_equal(x, y)); + let ys32 = ys.iter().map(|&y| y as f32).collect::>(); + let floors32: Vec = many(Floor(t), t, [&ys32], 6); assert_eq!(floors32, ys32.iter().map(|&y| t.floor(y)).collect::>()); - let (zs, mut complex) = (&[c(2.3, -1.7), c(0.5, 0.5), c(-1.0, 2.0)][..], Vec::new()); - assert!(Floor(t).many(Arg::Many(zs), 3, &mut complex) & Residue(t).many(Arg::One(c(1.0, 1.0)), Arg::Many(zs), 3, &mut complex)); + let zs = &[c(2.3, -1.7), c(0.5, 0.5), c(-1.0, 2.0)][..]; + let mut complex: Vec = many(Floor(t), t, [zs], 3); + complex.extend(many::<_, Complex64, _, 2>(Residue(t), t, [&[c(1.0, 1.0)], zs], 3)); let want: Vec = zs.iter().map(|&z| t.complex_floor(z)).chain(zs.iter().map(|&z| t.complex_residue(c(1.0, 1.0), z))).collect(); assert_eq!(complex, want); } - // Many items against one number take the range of floats equal to it, in either order. Each test item serves as the number. let items = [3.0, -3.0, 0.3, 0.1 + 0.2, 0.0, -0.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY, 2.0, 2.9999999999999996, 1e-310, 1.0 + 3e-14]; for t in [T, Tolerance::EXACT] { against_one::(t, &items); against_one::(t, &items); against_one::(t, &items); } - assert_eq!(Gcd(T).one(12.0, 18.0), 6.0); - assert!(!Gcd(T).many(Arg::One(f64::INFINITY), Arg::Many(&[2.0][..]), 1, &mut Vec::new())); - // Lcm has the sign of x y, gives 0 where the gcd is 0, and divides at the tolerance's gcd. - assert_eq!([Lcm(T).one(-4.0, 6.0), Lcm(T).one(0.0, 5.0), Lcm(T).one(0.6, 0.4)], [-12.0, 0.0, 1.2000000000000002]); - assert_eq!(Lcm(T).one(c(2.0, 0.0), c(0.0, 3.0)), c(0.0, 6.0)); + let real = |k| Kernel::new(k, Variant::checked(f64::DTYPE)); + assert_eq!(real(Gcd(T)).real(&[12.0, 18.0], Settings { tolerance: T }).unwrap(), 6.0); + assert!(real(Gcd(T)).real(&[f64::INFINITY, 2.0], Settings { tolerance: T }).is_err()); + let lcm = Kernel::new(Lcm(T), Variant::checked(f64::DTYPE)); + assert_eq!([(-4.0, 6.0), (0.0, 5.0), (0.6, 0.4)].map(|(x, y)| lcm.real(&[x, y], Settings { tolerance: T }).unwrap()), [-12.0, 0.0, 1.2000000000000002]); + let lcm = Kernel::new(Lcm(T), Variant::checked(Complex64::DTYPE)); + assert_eq!(lcm.complex(&[c(2.0, 0.0), c(0.0, 3.0)], Settings { tolerance: T }).unwrap(), c(0.0, 6.0)); } - /// An `Indicator` gives each comparison's Booleans as 1 and 0 at each float width, over more items than its buffer holds, with one number on either side and with slices on both. - #[test] - fn indicator() { - fn check + Copy>(k: K) where Indicator: Binary { - let x: Vec = (0..700).map(|i| F::narrow(f64::from(i % 7) - 3.0)).collect(); - let zero = F::narrow(0.0); - for (a, b) in [(Arg::Many(&x[..]), Arg::One(zero)), (Arg::One(zero), Arg::Many(&x[..])), (Arg::Many(&x[..]), Arg::Many(&x[1..]))] { - let (mut bools, mut lanes) = (Vec::new(), Vec::new()); - assert!(k.many(a, b, 699, &mut bools) & Indicator(k).many(a, b, 699, &mut lanes)); - assert_eq!(lanes, bools.iter().map(|&b| F::narrow(if b { 1.0 } else { 0.0 })).collect::>()); - } - } - macro_rules! comparisons { - ($f:ty, $t:expr) => { - check::<$f, _>(Equal($t)); - check::<$f, _>(NotEqual($t)); - check::<$f, _>(Less($t)); - check::<$f, _>(LessEqual($t)); - check::<$f, _>(Greater($t)); - check::<$f, _>(GreaterEqual($t)); - check::<$f, _>(Matches($t)); - }; - } - for t in [T, Tolerance::EXACT] { - comparisons!(f64, t); - comparisons!(f32, t); - } - } /// Tolerant search finds the first match that a brute-force search finds, at every width and at small and large tolerances, among floats crowded near a few values, and among zeros, infinities and NaN, whether it scans for a few items or indexes many. #[test] diff --git a/tests/speed_python.txt b/tests/speed_python.txt index 52c13ea..830c3ff 100644 --- a/tests/speed_python.txt +++ b/tests/speed_python.txt @@ -79,6 +79,8 @@ np.asarray(x) 0.5 0.5 instructions linux 0.29 xf+yf 3.4 3.2 instructions linux 0.735 # 5.7 copies, for Exp's 5.5 copies at f32, plus the call. rm.exp(xf) 12 6.5 instructions linux 4.193 +# Degree-five roots of nonnegative f32 values, matching the native Root row; public bulk cost is uncalibrated. +rm.root(root_nf, root_yf) - - # Real-input cis with complex output at f32; bulk cost is uncalibrated. rm.cis(xf) - - # Error-function slice forms at f32; bulk costs are uncalibrated. diff --git a/tests/support/mod.rs b/tests/support/mod.rs index 9ed7e34..cdae223 100644 --- a/tests/support/mod.rs +++ b/tests/support/mod.rs @@ -32,6 +32,13 @@ pub fn threads_lock() -> std::sync::MutexGuard<'static, ()> { pub fn map( kernel: &rustymath::operation::Kernel, args: [rustymath::buffer::View<'_>; N], n: usize, out: &mut rustymath::buffer::Buffer +) -> Result<(), rustymath::error::Error> { + map_with(kernel, args, n, out, rustymath::operation::Settings { tolerance: rustymath::tolerant::Tolerance::DEFAULT }) +} + +pub fn map_with( + kernel: &rustymath::operation::Kernel, args: [rustymath::buffer::View<'_>; N], n: usize, out: &mut rustymath::buffer::Buffer, + settings: rustymath::operation::Settings ) -> Result<(), rustymath::error::Error> { use rustymath::agreement::{Binding, Mapping, Missing, Read, ReadPlan}; let sources = args.map(Binding::full); @@ -39,7 +46,7 @@ pub fn map( ReadPlan::input(i, args[i].len()).remapped(if args[i].len() == 1 { Mapping::Single } else { Mapping::Linear(1) }, n, Missing::Default) }); let reads = plans.each_ref().map(|plan| Read { sources: &sources, plan }); - kernel.map(&reads, 0..n, rustymath::operation::Settings { tolerance: rustymath::tolerant::Tolerance::DEFAULT }, out) + kernel.map(&reads, 0..n, settings, out) } /// `kernel`'s `n` results on `args`, each of one item or of `n`, or `None` where any result fails. diff --git a/tests/test_speed.py b/tests/test_speed.py index 231c8ca..add88cd 100644 --- a/tests/test_speed.py +++ b/tests/test_speed.py @@ -13,8 +13,10 @@ def arrays(): u = lambda *shape: rng.uniform(-1, 1, shape) x, y = u(N), u(N) xs = np.sort(x) + root_nf = rm.array(np.full(N, 5, dtype=np.float32)) + root_yf = rm.array((1000 * ((np.arange(N) * 0.414_213_562_373_095_1) % 1)).astype(np.float32)) return dict(np=np, rm=rm, xn=x, x=rm.array(x), y=rm.array(y), xf=rm.array(x.astype(np.float32)), yf=rm.array(y.astype(np.float32)), - xs=rm.array(xs), xsf=rm.array(xs.astype(np.float32)), + xs=rm.array(xs), xsf=rm.array(xs.astype(np.float32)), root_nf=root_nf, root_yf=root_yf, m=rm.array(u(256, 256)), a=rm.array(u(256, 256)), b=rm.array(u(256, 256)), row=rm.array(u(256)), idx=rm.array(rng.permutation(N)), img=rm.array(u(256, 256)), k=rm.array(u(3, 3)), kern=rm.array(u(16)), frames=rm.array(u(1024, 64)), r=rm.array(rng.integers(0, 256, N)), c=rm.array(rng.integers(0, 256, N)), mask=rm.array(u(N) > 0)) @@ -93,3 +95,24 @@ def test_speed(capsys): failures += [line for _, line in over] if timing.checks_budgets: assert not failures, '\n'.join(failures) for line in failures: print(f'{line}: over', file=sys.stderr) + + +def test_conversion_controls(capsys): + "Public conversion workloads, including Python calls and result release; fill converts one scalar, and float-width astype uses ArrayView's reader." + with capsys.disabled(), held(rm.timing()): + for n in (16, 65_536): + x = rm.array((np.arange(n) % 31) / 4) + whole, outside = rm.array(np.arange(n) % 31, dtype='int16'), rm.array([n]) + cases = [('floor', lambda: x.astype('int16'), (n,)), + ('float-width', lambda: x.astype('float32'), (n,)), + ('seeded-bool-scan', lambda: rm.accumulate(rm.equal, x, start=0), (n,)), + ('scalar-fill-keep', lambda: whole.take(outside, mode='fill', fill=3.).compact(), (1,))] + for name, run, shape in cases: + result = run() + assert result.shape == shape + del result + steps = max(1, 2048 // n) + def repeated(): + for _ in range(steps): run() + ns = fastest(repeated) * 1e9 / steps + print(f'conversion {name} n={n} output={shape[0]} ns/call={ns:.3f} ns/item={ns/shape[0]:.6f}', file=sys.stderr) diff --git a/tests/test_structure.py b/tests/test_structure.py index 0f0349d..49c3389 100644 --- a/tests/test_structure.py +++ b/tests/test_structure.py @@ -19,6 +19,10 @@ def test_selection_and_metadata(): assert a[[2]].tolist() == a[[0.5]].tolist() == [[0, 0, 0]] repeated = keyed[['a', 'a']] assert repeated.tolist() == [[1, 2], [1, 2]] and repeated.keys == (None, None) + for s in (slice(None, -1, -1), slice(-1, None, -1), slice(-10**100, 10**100, 2), slice(10**100, -10**100, -2), slice(None, None, -10**100)): + assert a[:, s].tolist() == np.asarray(a)[:, s].tolist() + assert keyed['b':'a':-1].tolist() == [[3, 4], [1, 2]] + with pytest.raises(KeyError): keyed['missing':] def test_named_index_arrays_pair_and_address_by_name(): a = rm.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]) @@ -100,6 +104,9 @@ def test_join_split_repeat_and_padding(a): assert rm.concatenate([a, a], axis=1).tolist() == [[1, 2, 1, 2], [3, 4, 3, 4]] and rm.concatenate([a, a]).shape == (4, 2) assert rm.concatenate([a, [5]], axis=None).tolist() == [1, 2, 3, 4, 5] assert rm.stack([a, a], axis=-1).tolist() == [[[1, 1], [2, 2]], [[3, 3], [4, 4]]] + pieces = [a, rm.array([[5], [6]]), rm.array([[7, 8], [9, 10]])] + assert rm.concatenate(pieces, axis=1).tolist() == [[1, 2, 5, 7, 8], [3, 4, 6, 9, 10]] + assert rm.concatenate([a[:, :0], a, a[:, :0]], axis=1).tolist() == [[1, 2], [3, 4]] assert [x.tolist() for x in rm.split(a, 2, axis=1)] == [[[1], [3]], [[2], [4]]] assert [x.tolist() for x in rm.split(a, 2, axis=None)] == [[1, 2], [3, 4]] assert [x.tolist() for x in rm.split([1, 2, 3], [1, 2])] == [[1], [2], [3]] diff --git a/tests/test_ufuncs.py b/tests/test_ufuncs.py index a136e42..7aee9ad 100644 --- a/tests/test_ufuncs.py +++ b/tests/test_ufuncs.py @@ -340,6 +340,8 @@ def test_view_folds(threads): for cols, in_place in [(slice(None, 512), True), (slice(None, None, 2), False)]: with rm.explain() as report: assert np.isclose(rm.sum(r[:, cols]), m[:, cols].sum()) + assert ('in place' in report.routes) == in_place + with rm.explain() as report: assert all(np.allclose(rm.sum(r[:, cols], axis=a), m[:, cols].sum(a)) for a in [0, 1]) assert 'mapped fold' in report.routes and ('in place' in report.routes) == in_place diff --git a/timing/src/lib.rs b/timing/src/lib.rs index 097f732..8179d17 100644 --- a/timing/src/lib.rs +++ b/timing/src/lib.rs @@ -13,16 +13,25 @@ pub const RUNS: usize = 5; const ENOUGH: Duration = Duration::from_millis(20); /// The time of the fastest of `run`'s calls, in seconds. It stops after the second call, or any later one, once the calls have taken `ENOUGH` since the first began. Otherwise it makes `RUNS` calls, and goes on calling until the calls have taken `least`: a short case's fastest call then comes from many calls, and varies less between test runs. Call `k` is passed `k`, so it can read its own inputs. A caller whose inputs hold one set for each of `RUNS` calls passes `Duration::ZERO`. -pub fn fastest(least: Duration, mut run: impl FnMut(usize)) -> f64 { - let (begun, mut best) = (Instant::now(), Duration::MAX); - for k in 0.. { - let start = Instant::now(); - run(k); - best = best.min(start.elapsed()); - let spent = begun.elapsed(); - if k > 0 && spent >= ENOUGH || k + 1 >= RUNS && spent >= least { break; } +pub fn fastest(least: Duration, run: impl FnMut(usize)) -> f64 { Sampling::Measure.fastest(least, run) } + +/// Routine budget checks need two completed samples; measurements retain the full sampling window. +#[derive(Clone, Copy)] +pub enum Sampling { Check, Measure } + +impl Sampling { + pub fn fastest(self, least: Duration, mut run: impl FnMut(usize)) -> f64 { + let (least, runs) = match self { Self::Check => (Duration::ZERO, 2), Self::Measure => (least, RUNS) }; + let (begun, mut best) = (Instant::now(), Duration::MAX); + for k in 0.. { + let start = Instant::now(); + run(k); + best = best.min(start.elapsed()); + let spent = begun.elapsed(); + if k > 0 && spent >= ENOUGH || k + 1 >= runs && spent >= least { break; } + } + best.as_secs_f64() } - best.as_secs_f64() } /// What `f` gives, and the instructions retired in user space while it ran. On Linux a `perf_event_open` counter counts the calling thread, without any threads that `f` starts. On macOS `proc_pid_rusage` counts the whole process. The speed tests run each case on one thread, where the two agree. `None` where the platform has no counter or refuses one, as Linux does where `/proc/sys/kernel/perf_event_paranoid` is above 2, and in the browser build.