rustymath_core/
polynomial.rs

1//! Polynomials of one variable. A polynomial is the slice of its coefficients in ascending order, so `c` stands for `c[0] + c[1]x + c[2]x² + …`. Trailing zero coefficients don't count towards the degree. `evaluate` therefore never multiplies an infinite x by one of them. `from_roots`, `derivative`, `value` and `factored` check for an interrupt between roots, or between groups of [`BETWEEN_CHECKS`](crate::BETWEEN_CHECKS) coefficients or roots, and return its error as soon as they get one, as [`linalg::exact`](crate::linalg::exact)'s functions do. A call too short for a second group makes no check, so evaluating a short polynomial at each of many points costs no checks.
2use crate::{BETWEEN_CHECKS, error::Error, field::{minus, times}, float::{FUSED, Float, sealed::Core}, interrupt::check};
3use num_complex::Complex;
4use num_traits::{NumAssignRef, Zero};
5
6/// `c` without its trailing zeros.
7pub fn trimmed<T: Zero>(c: &[T]) -> &[T] { &c[..c.iter().rposition(|k| !k.is_zero()).map_or(0, |i| i + 1)] }
8
9/// Appends to `out` the value of the polynomial `c` at each item of `x`, by Horner's rule. An f16 polynomial computes in f32 and rounds each result once.
10pub fn evaluate<F: Float>(c: &[F], x: &[F], out: &mut Vec<F>) {
11    let Some((&top, rest)) = trimmed(c).split_last() else { return out.resize(out.len() + x.len(), F::zero()) };
12    unsafe {
13        crate::slice::appended(out, x.len(), |out| {
14            if crate::slice::running(size_of::<F::Working>()) == 16 { blocks::<F, 16>(top, rest, x, out) } else { blocks::<F, 32>(top, rest, x, out) }
15        })
16    }
17}
18
19/// `evaluate` through `horner` on blocks of `L` items, as many as `slice::running` gives.
20#[inline]
21fn blocks<F: Float, const L: usize>(top: F, rest: &[F], x: &[F], out: &mut [std::mem::MaybeUninit<F>]) {
22    let ((blocks, tail), (outs, last)) = (x.as_chunks::<L>(), out.as_chunks_mut::<L>());
23    for (xs, ys) in blocks.iter().zip(outs) { horner(top, rest, xs, ys) }
24    let mut padded = [F::zero(); L];
25    padded[..tail.len()].copy_from_slice(tail);
26    let mut values = [std::mem::MaybeUninit::uninit(); L];
27    horner(top, rest, &padded, &mut values);
28    last.copy_from_slice(&values[..tail.len()]);
29}
30
31/// Writes into `out` the value at each of `xs` of the polynomial whose last coefficient is `top` and whose others are `rest`. The items' running values stay in registers, and each coefficient takes one multiply-add per item. The loop over the items vectorizes, and the items' operations don't wait for each other. The block's fixed size keeps its loops free of calls to `memcpy`, which a variable length brings in.
32#[inline]
33fn horner<F: Float, const L: usize>(top: F, rest: &[F], xs: &[F; L], out: &mut [std::mem::MaybeUninit<F>; L]) {
34    let mut y = [top.working(); L];
35    for &k in rest.iter().rev() { let k = k.working(); for i in 0..L { y[i] = y[i].mla::<FUSED>(xs[i].working(), k); } }
36    for i in 0..L { out[i].write(F::from_working(y[i])); }
37}
38
39/// The roots of the polynomial `c`, as the eigenvalues of its companion matrix, computed in f64 by [`eigenvalues`](crate::linalg::eigenvalues), sorted by real part, then imaginary part. Real coefficients take the real solver, so a real root has an imaginary part of exactly 0. `None` when the eigenvalue solver doesn't converge.
40pub fn roots<F: Float>(c: &[Complex<F>]) -> Option<Vec<Complex<F>>> {
41    let c: Vec<Complex<f64>> = trimmed(c).iter().map(|z| Complex::new(z.re.wide(), z.im.wide())).collect();
42    let n = c.len().saturating_sub(1);
43    if n == 0 { return Some(Vec::new()); }
44    // The companion matrix has ones below its diagonal and -c[i]/c[n] down its last column.
45    let a: Vec<Complex<f64>> = c[..n].iter().map(|&k| -crate::complex::divide(k, c[n])).collect();
46    let companion: Vec<Complex<f64>> = (0..n * n).map(|k| if k % n == n - 1 { a[k / n] } else { Complex::new(f64::from(k / n == k % n + 1), 0.0) }).collect();
47    let mut values = if a.iter().all(|z| z.im == 0.0) { crate::linalg::eigenvalues(&companion.iter().map(|z| z.re).collect::<Vec<f64>>(), n)? }
48    else { crate::linalg::eigenvalues(&companion, n)? };
49    values.sort_by(|x, y| x.re.total_cmp(&y.re).then(x.im.total_cmp(&y.im)));
50    Some(values.into_iter().map(|z| Complex::new(F::narrow(z.re), F::narrow(z.im))).collect())
51}
52
53/// The polynomial `m` times the product of `x - r` over the roots `r`. It checks for an interrupt between roots.
54pub fn from_roots<T: Clone + NumAssignRef>(m: T, roots: &[T]) -> Result<Vec<T>, Error> {
55    let mut c = vec![m];
56    for (i, r) in roots.iter().enumerate() {
57        if i > 0 { check()?; }
58        // `a - b r`
59        let less = |a: &T, b: &T| minus(a, &times(b, r));
60        c.push(T::zero());
61        for i in (1..c.len()).rev() { c[i] = less(&c[i - 1], &c[i]); }
62        c[0] = less(&T::zero(), &c[0]);
63    }
64    Ok(c)
65}
66
67/// The derivative of the polynomial `c`. It checks for an interrupt between groups of `BETWEEN_CHECKS` coefficients.
68pub fn derivative<T: Clone + NumAssignRef>(c: &[T]) -> Result<Vec<T>, Error> {
69    let (mut n, mut d) = (T::zero(), Vec::with_capacity(c.len().saturating_sub(1)));
70    for (i, group) in c.get(1..).unwrap_or_default().chunks(BETWEEN_CHECKS).enumerate() {
71        if i > 0 { check()?; }
72        d.extend(group.iter().map(|k| { n += T::one(); times(&n, k) }));
73    }
74    Ok(d)
75}
76
77/// The value of the polynomial `c` at `x`, by Horner's rule, for any `T` with the field operations, such as an exact rational. [`evaluate`] computes the values of a float polynomial at many points at once. It checks for an interrupt between groups of `BETWEEN_CHECKS` coefficients.
78pub fn value<T: Clone + NumAssignRef>(c: &[T], x: &T) -> Result<T, Error> {
79    let Some((top, rest)) = trimmed(c).split_last() else { return Ok(T::zero()) };
80    let mut y = top.clone();
81    for (i, group) in rest.rchunks(BETWEEN_CHECKS).enumerate() {
82        if i > 0 { check()?; }
83        for k in group.iter().rev() {
84            y *= x;
85            y += k;
86        }
87    }
88    Ok(y)
89}
90
91/// The value at `x` of `m` times the product of `x - r` over the roots `r`. It checks for an interrupt between groups of `BETWEEN_CHECKS` roots.
92pub fn factored<T: Clone + NumAssignRef>(m: &T, roots: &[T], x: &T) -> Result<T, Error> {
93    let mut y = m.clone();
94    for (i, group) in roots.chunks(BETWEEN_CHECKS).enumerate() { if i > 0 { check()?; } for r in group { y *= &minus(x, r); } }
95    Ok(y)
96}
97
98/// The value of each cell of `bases.len()` digits in `digits`, in the mixed radix of `bases`, as APL's `⊥` gives it: Horner's rule, with each step's own base. The first base is never used. `integer::decode` gives `i64`s, checked for overflow.
99pub fn decode<T: Clone + NumAssignRef>(bases: &[T], digits: &[T]) -> Vec<T> {
100    digits.chunks(bases.len()).map(|d| {
101        d[1..].iter().zip(&bases[1..]).fold(d[0].clone(), |mut v, (k, b)| {
102            v *= b;
103            v += k;
104            v
105        })
106    }).collect()
107}
108
109/// The digits of each of `values` in the mixed radix of `bases`, as APL's `⊤` gives them: each value's digits together, the most significant first. `residue(b, v)` is `v`'s residue by the base `b`, which for a base of 0 is `v`, so a base of 0 takes all of what remains as its digit. `integer::encode` gives `i64`s, checked for overflow.
110pub fn encode<T: Clone + NumAssignRef>(bases: &[T], values: &[T], residue: impl Fn(&T, &T) -> T) -> Vec<T> {
111    encoded(bases, values, |b, mut v| {
112        let digit = residue(b, &v);
113        if b.is_zero() { v.set_zero() } else {
114            v -= &digit;
115            v /= b;
116        }
117        Some((v, digit))
118    }).expect("a step in a field gives a value")
119}
120
121/// `encode`'s rule, where `div_mod(base, value)` gives the value that the next base divides and the digit, and `None` where it gives `None`. `integer::encode` divides through its `Divisor`s this way.
122pub(crate) fn encoded<B, T: Clone + Zero>(bases: &[B], values: &[T], div_mod: impl Fn(&B, T) -> Option<(T, T)>) -> Option<Vec<T>> {
123    let mut digits = vec![T::zero(); bases.len() * values.len()];
124    for (out, value) in digits.chunks_mut(bases.len().max(1)).zip(values) {
125        let mut value = value.clone();
126        for (digit, base) in out.iter_mut().zip(bases).rev() { (value, *digit) = div_mod(base, value)?; }
127    }
128    Some(digits)
129}
130
131/// The polynomial of degree `degree` that fits the points `(x[i], y[i])` by least squares, which [`lstsq`](crate::linalg::lstsq) solves in f64. Where fewer than `degree + 1` distinct points fix the polynomial, it's the fit whose coefficients have the least norm. Panics if `x` and `y` differ in length.
132pub fn fit<F: Float>(x: &[F], y: &[F], degree: usize) -> Vec<F> {
133    assert_eq!(x.len(), y.len(), "fit needs one y for each x");
134    let powers: Vec<f64> = x.iter().flat_map(|x| {
135        let x = x.wide();
136        std::iter::successors(Some(1.0), move |p| Some(p * x)).take(degree + 1)
137    }).collect();
138    let y: Vec<f64> = y.iter().map(|y| y.wide()).collect();
139    crate::linalg::lstsq(&powers, &y, [x.len(), degree + 1], 1).0.into_iter().map(F::narrow).collect()
140}
rustymath_core/polynomial.rs

rustymath_core/
polynomial.rs
