Share native math execution, selector geometry and conversions

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"] }

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

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<usize> {
-    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<SliceSpec> {
+    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::<PyString>() { py.None().into_bound(py) } else { v.clone() };
+    let unpack = py.get_type::<PySlice>().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<Option<SliceBound>> {
+        Ok(if v.is_none() { None }
+        else if let Ok(s) = v.cast::<PyString>() { 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::<PySlice>() || obj.is_instance_of::<PyInt>() || 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<Vec<String>> {
@@ -20,21 +35,8 @@ pub fn key_names(obj: &Bound<'_, PyAny>) -> Option<Vec<String>> {
 /// 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::<PySlice>() {
-        let py = obj.py();
-        let bound = |name| -> PyResult<Bound<'_, PyAny>> {
-            let value = slice.getattr(name)?;
-            if let Ok(s) = value.cast::<PyString>() { 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::<PyString>() {
-            let reverse = !step.is_none() && step.extract::<isize>()? < 0;
-            let at = stop.extract::<isize>()?;
-            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::<PySlice>().call1((start, stop, step))?.cast_into::<PySlice>()?.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::<PyString>();
@@ -45,7 +47,7 @@ pub fn selector(layout: &Layout, axis: usize, obj: &Bound<'_, PyAny>, array: imp
         let i = obj.extract::<i64>()?;
         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<Vec<IndexPart<Bound<'py
         else if p.is(&ellipsis) {
             IndexPart::Rest
         } else if p.is_instance_of::<PySlice>() && ["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::<PySlice>() { 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`.

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<T: crate::buffer::Cast>(&self, range: Range<usize>, checked: bool, output: &mut Vec<T>) -> 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<usize>) -> 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::<usize>();
         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<A: Copy>(mapping: &Mapping, data: &[A], start: usize, n: usize) -> Option<A> {
-    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 [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<Arg<'a, A>> {
+fn in_place<'a, A: Copy>(mapping: &Walked, data: &'a [A], start: usize, n: usize) -> Option<Arg<'a, A>> {
     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<A>) -> 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<A: Element>(data: &[A], mapping: &Mapping, run: usize, buf: &mut Vec<A>) -> Option<Mapping> {
@@ -1272,7 +1266,7 @@ pub fn tiled_cycle<A: Element>(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 [A] {
     buf.clear();
     let reason = copy_walked(mapping, data, start, n, A::FILL, buf);

core/src/arc.rs

@@ -70,8 +70,9 @@ pub(crate) fn atan<R: Native, const FUSED: bool>(x: R) -> R {
     (hi + ((u * z).mla::<FUSED>(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<R: Native>(_: &Atan, _: &crate::recipe::NativeOps<{ crate::FUSED }>, x: [crate::Arg<'_, R>; 1], n: usize, out: &mut [std::mem::MaybeUninit<R>]) -> 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]