#![warn(missing_docs)]
#![warn(missing_crate_level_docs)]
#![warn(macro_use_extern_crate)]
#![warn(invalid_html_tags)]
#![warn(missing_copy_implementations)]
#![warn(missing_debug_implementations)]
#![warn(unreachable_pub)]
#![warn(unused_extern_crates)]
#![warn(unused_lifetimes)]
use core::{
ops::{
Add,
Mul,
AddAssign,
SubAssign,
MulAssign,
},
borrow::{
Borrow,
BorrowMut,
},
iter::{
self,
FromIterator,
IntoIterator,
ExactSizeIterator,
Extend,
},
marker::{
PhantomData,
},
convert::{
self,
},
};
use num_traits::{
Zero,
One,
Inv,
};
use num_complex::{
Complex,
Complex32,
Complex64,
};
use custom_error::custom_error;
custom_error!{
#[derive(PartialEq)]
#[allow(missing_docs)]
pub Error
CoeffConv{cause: String}
= "Could not convert coefficient: {cause}",
Singular = "Encountered singular matrix when expecting nonsingular; perhaps sparsity was set too low",
SparsityTooLow = "Sparsity bound set too low in sparse interpolation",
MissingExponents = "The exponent of a non-zero term was missing from the list",
}
pub type Result<T> = core::result::Result<T, Error>;
pub trait OneWay {
type Source;
type Dest;
fn one_way(src: Self::Source) -> Result<Self::Dest>;
}
pub trait TwoWay: OneWay {
fn other_way(src: <Self as OneWay>::Dest) -> Result<<Self as OneWay>::Source>;
}
#[derive(Debug)]
pub struct DefConv<S,D>(PhantomData<(S,D)>);
impl<S,D> TwoWay for DefConv<S,D>
where Self: OneWay<Source=S, Dest=D>,
DefConv<D,S>: OneWay<Source=D, Dest=S>,
{
fn other_way(src: D) -> Result<S> {
<DefConv<D,S> as OneWay>::one_way(src)
}
}
impl<S> OneWay for DefConv<S,S>
{
type Source = S;
type Dest = S;
#[inline(always)]
fn one_way(src: S) -> Result<S> {
Ok(src)
}
}
macro_rules! one_way_as {
($s:ty, $($ds:ty),+) => {
$(
impl OneWay for DefConv<$s,$ds> {
type Source = $s;
type Dest = $ds;
fn one_way(src: $s) -> Result<$ds> {
Ok(src as $ds)
}
}
)+
};
}
macro_rules! one_way_as_check {
($s:ty, $($ds:ty),+) => {
$(
impl OneWay for DefConv<$s,$ds> {
type Source = $s;
type Dest = $ds;
fn one_way(src: $s) -> Result<$ds> {
let dest = src as $ds;
if (dest as $s) == src {
Ok(dest)
} else {
Err(Error::CoeffConv{cause: format!("{} not representable as {}", src, stringify!($ds))})
}
}
}
)+
}
}
macro_rules! one_way_round {
($s:ty, $($ds:ty),+) => {
$(
impl OneWay for DefConv<$s,$ds> {
type Source = $s;
type Dest = $ds;
fn one_way(src: $s) -> Result<$ds> {
let rounded = src.round();
let dest = rounded as $ds;
if (dest as $s) == rounded {
Ok(dest)
} else {
Err(Error::CoeffConv{cause: format!("{} not representable as {}", src, stringify!($ds))})
}
}
}
)+
}
}
macro_rules! one_way_try {
($s:ty, $($ds:ty),+) => {
$(
impl OneWay for DefConv<$s,$ds> {
type Source = $s;
type Dest = $ds;
fn one_way(src: $s) -> Result<$ds> {
convert::TryInto::<$ds>::try_into(src).map_err(|e| Error::CoeffConv{cause: e.to_string()})
}
}
)+
}
}
one_way_try!(i8, u8);
one_way_as!(i8, i16, u16, i32, u32, i64, u64, i128, u128, f32, f64);
one_way_try!(u8, i8);
one_way_as!(u8, i16, u16, i32, u32, i64, u64, i128, u128, f32, f64);
one_way_try!(i16, i8, u8, u16);
one_way_as!(i16, i32, u32, i64, u64, i128, u128, f32, f64);
one_way_try!(u16, i8, u8, i16);
one_way_as!(u16, i32, u32, i64, u64, i128, u128, f32, f64);
one_way_try!(i32, i8, u8, i16, u16, u32);
one_way_as!(i32, i64, u64, i128, u128, f64);
one_way_as_check!(i32, f32);
one_way_try!(u32, i8, u8, i16, u16, i32);
one_way_as!(u32, i64, u64, i128, u128, f64);
one_way_as_check!(u32, f32);
one_way_try!(i64, i8, u8, i16, u16, i32, u32, u64);
one_way_as!(i64, i128, u128);
one_way_as_check!(i64, f32, f64);
one_way_try!(u64, i8, u8, i16, u16, i32, u32, i64);
one_way_as!(u64, i128, u128);
one_way_as_check!(u64, f32, f64);
one_way_try!(i128, i8, u8, i16, u16, i32, u32, i64, u64, u128);
one_way_as_check!(i128, f32, f64);
one_way_try!(u128, i8, u8, i16, u16, i32, u32, i64, u64, i128);
one_way_as_check!(u128, f32, f64);
one_way_round!(f32, i8, u8, i16, u16, i32, u32, i64, u64, i128, u128);
one_way_as!(f32, f64);
one_way_round!(f64, i8, u8, i16, u16, i32, u32, i64, u64, i128, u128);
impl OneWay for DefConv<f64, f32> {
type Source = f64;
type Dest = f32;
fn one_way(src: f64) -> Result<f32> {
let dest = src as f32;
if dest.log2().round() == (src.log2().round() as f32) {
Ok(dest)
} else {
Err(Error::CoeffConv{cause: format!("f64->f32 exponent overflow on {}", src)})
}
}
}
macro_rules! two_way_complex {
($t:ty, $($us:ty),+) => {
$(
impl OneWay for DefConv<$us, Complex<$t>> {
type Source = $us;
type Dest = Complex<$t>;
#[inline(always)]
fn one_way(src: Self::Source) -> Result<Self::Dest> {
Ok(Complex::from(DefConv::<$us,$t>::one_way(src)?))
}
}
impl OneWay for DefConv<Complex<$t>, $us> {
type Source = Complex<$t>;
type Dest = $us;
#[inline(always)]
fn one_way(src: Self::Source) -> Result<Self::Dest> {
DefConv::<$t,$us>::one_way(src.re)
}
}
)+
};
}
two_way_complex!(f32, i8, u8, i16, u16, i32, u32, i64, u64, i128, u128, f32, f64);
two_way_complex!(f64, i8, u8, i16, u16, i32, u32, i64, u64, i128, u128, f32, f64);
pub trait CloseTo {
type Item;
fn close_to(&self, x: &Self::Item, y: &Self::Item) -> bool;
fn close_to_zero(&self, x: &Self::Item) -> bool;
fn close_to_iter<'a, Iter1, Iter2>(&'a self, x: Iter1, y: Iter2) -> bool
where Iter1: Iterator<Item=&'a Self::Item>,
Iter2: Iterator<Item=&'a Self::Item>,
{
x.zip(y).all(|(xi, yi)| self.close_to(xi, yi))
}
}
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
pub struct CloseToEq<T>{
pub zero: T,
}
impl<T: Zero> CloseToEq<T> {
pub fn new() -> Self {
Self {
zero: T::zero(),
}
}
}
impl<T: PartialEq> CloseTo for CloseToEq<T> {
type Item = T;
#[inline(always)]
fn close_to(&self, x: &Self::Item, y: &Self::Item) -> bool {
x.eq(y)
}
#[inline(always)]
fn close_to_zero(&self, x: &Self::Item) -> bool {
x.eq(&self.zero)
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct RelativeParams<T, E=T>
{
pub zero_thresh: E,
pub max_relative: E,
phantom: PhantomData<T>,
}
macro_rules! impl_relative_params {
($t:ident, $e:ident, $norm:ident) => {
impl RelativeParams<$t,$e> {
#[inline(always)]
pub fn new(zero_thresh: Option<$e>, max_relative: Option<$e>) -> Self {
Self {
zero_thresh: zero_thresh.unwrap_or($e::EPSILON),
max_relative: max_relative.unwrap_or($e::EPSILON),
phantom: PhantomData,
}
}
}
impl Default for RelativeParams<$t,$e> {
#[inline(always)]
fn default() -> Self {
Self::new(None, None)
}
}
impl CloseTo for RelativeParams<$t,$e> {
type Item = $t;
#[inline]
fn close_to(&self, x: &Self::Item, y: &Self::Item) -> bool {
if x == y {
return true
};
if $t::is_infinite(*x) || $t::is_infinite(*y) {
return false
};
let absx = $t::$norm(*x);
let absy = $t::$norm(*y);
let largest = if absx <= absy {
absy
} else {
absx
};
if largest <= self.zero_thresh {
true
} else {
$t::$norm(x - y) / largest <= self.max_relative
}
}
#[inline(always)]
fn close_to_zero(&self, x: &Self::Item) -> bool {
$t::$norm(*x) <= self.zero_thresh
}
}
};
}
impl_relative_params!(f32, f32, abs);
impl_relative_params!(f64, f64, abs);
impl_relative_params!(Complex32, f32, norm);
impl_relative_params!(Complex64, f64, norm);
pub trait PolyTraits {
type Coeff;
type SparseInterpEval: EvalTypes<Coeff=Self::Coeff>;
type SparseInterpInfo;
fn slice_mul(out: &mut [Self::Coeff], lhs: &[Self::Coeff], rhs: &[Self::Coeff]);
#[inline(always)]
fn mp_eval_prep<U>(pts: impl Iterator<Item=U>) -> <EvalTrait<Self,U> as EvalTypes>::EvalInfo
where EvalTrait<Self,U>: EvalTypes<Coeff=Self::Coeff, Eval=U>
{
<EvalTrait<Self,U> as EvalTypes>::prep(pts)
}
#[inline(always)]
fn mp_eval_slice<U>(out: &mut impl Extend<U>,
coeffs: &[Self::Coeff],
info: &<EvalTrait<Self,U> as EvalTypes>::EvalInfo) -> Result<()>
where EvalTrait<Self,U>: EvalTypes<Coeff=Self::Coeff, Eval=U>
{
<EvalTrait<Self,U> as EvalTypes>::post(out, coeffs, info)
}
fn sparse_interp_prep(sparsity: usize, expons: impl Iterator<Item=usize>, max_coeff: &Self::Coeff)
-> (<Self::SparseInterpEval as EvalTypes>::EvalInfo, Self::SparseInterpInfo)
;
fn sparse_interp_slice(evals: &[<Self::SparseInterpEval as EvalTypes>::Eval],
info: &Self::SparseInterpInfo)
-> Result<Vec<(usize, Self::Coeff)>>;
}
#[derive(Debug)]
pub struct EvalTrait<T: ?Sized, U>(PhantomData<T>,PhantomData<U>);
pub trait EvalTypes {
type Coeff;
type Eval;
type EvalInfo;
fn prep(pts: impl Iterator<Item=Self::Eval>) -> Self::EvalInfo;
fn post(out: &mut impl Extend<Self::Eval>, coeffs: &[Self::Coeff], info: &Self::EvalInfo) -> Result<()>;
}
impl<C,U> EvalTypes for EvalTrait<ClassicalTraits<C>, U>
where C: Clone,
U: Clone + Zero + MulAssign + AddAssign,
DefConv<C,U>: OneWay<Source=C, Dest=U>,
{
type Coeff = C;
type Eval = U;
type EvalInfo = Vec<U>;
#[inline(always)]
fn prep(pts: impl Iterator<Item=Self::Eval>) -> Self::EvalInfo {
pts.collect()
}
fn post(out: &mut impl Extend<Self::Eval>, coeffs: &[Self::Coeff], info: &Self::EvalInfo) -> Result<()> {
ErrorIter::new(info.iter().map(|x| horner_eval(coeffs.iter(), x)))
.try_use(|evals| out.extend(evals))
}
}
#[derive(Debug,Default,Clone,Copy,PartialEq,Eq)]
pub struct ClassicalTraits<C>(PhantomData<C>);
impl<C> PolyTraits for ClassicalTraits<C>
where C: Clone + Mul<Output=C> + AddAssign,
DefConv<C, Complex64>: TwoWay<Source=C, Dest=Complex64>,
{
type Coeff = C;
type SparseInterpEval = EvalTrait<Self, Complex64>;
type SparseInterpInfo = (usize, Vec<(usize, Complex64)>, RelativeParams<Complex64, f64>);
#[inline(always)]
fn slice_mul(out: &mut [Self::Coeff], lhs: &[Self::Coeff], rhs: &[Self::Coeff]) {
classical_slice_mul(out, lhs, rhs);
}
fn sparse_interp_prep(sparsity: usize, expons: impl Iterator<Item=usize>, _max_coeff: &Self::Coeff)
-> (<Self::SparseInterpEval as EvalTypes>::EvalInfo, Self::SparseInterpInfo)
{
let mut theta_pows: Vec<_> = expons.map(|d| (d, Complex64::default())).collect();
let theta = match theta_pows.iter().map(|pair| pair.0).max() {
Some(max_pow) => Complex64::from_polar(1., 2.*core::f64::consts::PI / (max_pow as f64 + 1.)),
None => Complex64::default(),
};
for (ref expon, ref mut power) in theta_pows.iter_mut() {
*power = theta.powu(*expon as u32);
}
(EvalTrait::<Self, Complex64>::prep((0..2*sparsity).map(|e| theta.powu(e as u32))),
(sparsity, theta_pows, RelativeParams::<Complex64,f64>::new(Some(1e-8), Some(1e-8)))
)
}
fn sparse_interp_slice(
evals: &[Complex64],
info: &Self::SparseInterpInfo,
) -> Result<Vec<(usize, Self::Coeff)>>
{
let (sparsity, pows, close) = info;
assert_eq!(evals.len(), 2*sparsity);
let mut lambda = bad_berlekamp_massey(evals, close)?;
lambda.push(Complex64::from(-1.));
let (degs, roots): (Vec<usize>, Vec<Complex64>) = pows.iter().filter_map(
|(deg, rootpow)| {
classical_sparse_interp_helper(&lambda, rootpow).ok().and_then(
|eval| match close.close_to_zero(&eval) {
true => Some((deg, rootpow)),
false => None,
}
)}
).unzip();
if degs.len() != lambda.len() - 1 {
Err(Error::MissingExponents)
} else {
let evslice = &evals[..degs.len()];
ErrorIter::new(
bad_trans_vand_solve(&roots, evslice, close)?.into_iter()
.map(|c| DefConv::<C,Complex64>::other_way(c))
).try_use(|it| degs.into_iter().zip(it).collect())
}
}
}
fn classical_sparse_interp_helper(x: &Vec<Complex64>, y: &Complex64) -> Result<Complex64> {
horner_eval(x.iter(), y)
}
#[derive(Debug)]
#[repr(transparent)]
pub struct Poly<T,U> {
rep: T,
traits: PhantomData<U>,
}
pub type ClassicalPoly<C> = Poly<Vec<C>, ClassicalTraits<C>>;
impl<T,U,V,W> PartialEq<Poly<V,W>> for Poly<T,U>
where U: PolyTraits,
W: PolyTraits,
T: Borrow<[U::Coeff]>,
V: Borrow<[W::Coeff]>,
U::Coeff: PartialEq<W::Coeff>,
{
fn eq(&self, other: &Poly<V,W>) -> bool {
self.rep.borrow() == other.rep.borrow()
}
}
impl<T,U> Eq for Poly<T,U>
where U: PolyTraits,
T: Borrow<[U::Coeff]>,
U::Coeff: Eq,
{ }
impl<U> Poly<Vec<U::Coeff>,U>
where U: PolyTraits,
U::Coeff: Zero,
{
pub fn new(rep: Vec<U::Coeff>) -> Self {
let mut out = Self {
rep,
traits: PhantomData,
};
out.normalize();
out
}
fn normalize(&mut self) {
while self.rep.last().map_or(false, Zero::is_zero) {
self.rep.pop();
}
}
}
impl<T,U> Poly<T,U>
where U: PolyTraits,
T: Borrow<[U::Coeff]>,
U::Coeff: Zero,
{
pub fn new_norm(rep: T) -> Self {
assert!(rep.borrow().last().map_or(true, |x| !x.is_zero()));
Self {
rep,
traits: PhantomData,
}
}
}
impl<T,U> Default for Poly<T,U>
where T: Default,
{
#[inline(always)]
fn default() -> Self {
Self {
rep: T::default(),
traits: PhantomData,
}
}
}
impl<T,U> FromIterator<U::Coeff> for Poly<T,U>
where U: PolyTraits,
T: FromIterator<U::Coeff> + Borrow<[U::Coeff]>,
U::Coeff: Zero,
{
fn from_iter<V>(iter: V) -> Self
where V: IntoIterator<Item=U::Coeff>,
{
struct Trimmed<I: Iterator>(I, Option<(usize, I::Item)>);
impl<I: Iterator> Iterator for Trimmed<I>
where I::Item: Zero,
{
type Item = I::Item;
fn next(&mut self) -> Option<Self::Item> {
match self.1.take() {
Some((count, nzcoeff)) => Some(
if count == 0 {
nzcoeff
} else {
self.1 = Some((count - 1, nzcoeff));
Zero::zero()
}
),
None => self.0.next().and_then(|coeff| {
if coeff.is_zero() {
let mut count = 1;
while let Some(coeff) = self.0.next() {
if ! coeff.is_zero() {
self.1 = Some((count-1, coeff));
return Some(Zero::zero());
}
count += 1;
}
None
} else {
Some(coeff)
}
}),
}
}
}
Self::new_norm(Trimmed(iter.into_iter(), None).collect())
}
}
impl<T,U> Poly<T,U>
where U: PolyTraits,
T: Borrow<[U::Coeff]>,
U::Coeff: Clone + Zero + AddAssign + MulAssign,
{
#[inline(always)]
pub fn eval<V>(&self, x: &V) -> Result<V>
where DefConv<U::Coeff, V>: OneWay<Source=U::Coeff, Dest=V>,
V: Clone + Zero + AddAssign + MulAssign,
{
horner_eval(self.rep.borrow().iter(), x)
}
#[inline(always)]
pub fn mp_eval_prep<V>(pts: impl Iterator<Item=V>) -> <EvalTrait<U,V> as EvalTypes>::EvalInfo
where EvalTrait<U, V>: EvalTypes<Coeff=U::Coeff, Eval=V>,
{
U::mp_eval_prep(pts)
}
#[inline]
pub fn mp_eval<V>(&self, info: &<EvalTrait<U,V> as EvalTypes>::EvalInfo) -> Result<Vec<V>>
where EvalTrait<U, V>: EvalTypes<Coeff=U::Coeff, Eval=V>,
{
let mut out = Vec::new();
U::mp_eval_slice(&mut out, self.rep.borrow(), info)?;
Ok(out)
}
}
impl<T,U> Poly<T,U>
where U: PolyTraits,
{
#[inline(always)]
pub fn sparse_interp_prep(sparsity: usize, expons: impl Iterator<Item=usize>, max_coeff: &U::Coeff)
-> (<U::SparseInterpEval as EvalTypes>::EvalInfo, U::SparseInterpInfo)
{
U::sparse_interp_prep(sparsity, expons, max_coeff)
}
#[inline(always)]
pub fn sparse_interp(evals: &[<U::SparseInterpEval as EvalTypes>::Eval], info: &U::SparseInterpInfo)
-> Result<Vec<(usize, U::Coeff)>>
{
U::sparse_interp_slice(evals, info)
}
}
#[must_use]
struct LongZip<A,B,F,G,H>(A,B,F,G,H);
impl<A,B,F,G,H,O> Iterator for LongZip<A,B,F,G,H>
where A: Iterator,
B: Iterator,
F: FnMut(A::Item, B::Item) -> O,
G: FnMut(A::Item) -> O,
H: FnMut(B::Item) -> O,
{
type Item = O;
fn next(&mut self) -> Option<Self::Item> {
match self.0.next() {
Some(x) => Some(match self.1.next() {
Some(y) => self.2(x,y),
None => self.3(x),
}),
None => self.1.next().map(|y| self.4(y)),
}
}
}
impl<'a,'b,T,U,V,W> Add<&'b Poly<V,W>> for &'a Poly<T,U>
where U: PolyTraits,
W: PolyTraits,
&'a T: IntoIterator<Item=&'a U::Coeff>,
&'b V: IntoIterator<Item=&'b W::Coeff>,
&'a U::Coeff: Add<&'b W::Coeff, Output=U::Coeff>,
U::Coeff: AddAssign<&'b W::Coeff>,
U::Coeff: Clone + Zero
{
type Output = Poly<Vec<U::Coeff>, U>;
fn add(self, rhs: &'b Poly<V,W>) -> Self::Output {
Poly::from_iter(LongZip(
self.rep.into_iter(),
rhs.rep.into_iter(),
Add::add,
Clone::clone,
|x| {
let mut sum = U::Coeff::zero();
sum += x;
sum
},
))
}
}
impl<T,U,V,W> Add<Poly<V,W>> for Poly<T,U>
where U: PolyTraits,
W: PolyTraits,
T: IntoIterator<Item=U::Coeff>,
V: IntoIterator<Item=W::Coeff>,
U::Coeff: Zero + From<W::Coeff> + Add<W::Coeff, Output=U::Coeff>,
{
type Output = Poly<Vec<U::Coeff>, U>;
fn add(self, rhs: Poly<V,W>) -> Self::Output {
Poly::from_iter(LongZip(
self.rep.into_iter(),
rhs.rep.into_iter(),
Add::add,
convert::identity,
From::from,
))
}
}
impl<'a,'b,T,U,V> Mul<&'b Poly<V,U>> for &'a Poly<T,U>
where U: PolyTraits,
T: Borrow<[U::Coeff]>,
V: Borrow<[U::Coeff]>,
U::Coeff: Zero,
{
type Output = Poly<Box<[U::Coeff]>, U>;
fn mul(self, rhs: &'b Poly<V,U>) -> Self::Output {
let aslice = self.rep.borrow();
let bslice = rhs.rep.borrow();
if aslice.is_empty() || bslice.is_empty() {
return Poly::new_norm(Box::new([]));
}
let clen = aslice.len() + bslice.len() - 1;
let mut outbox = {
let mut vec = Vec::with_capacity(clen);
vec.resize_with(clen, U::Coeff::zero);
vec.into_boxed_slice()
};
U::slice_mul(outbox.borrow_mut(), aslice, bslice);
Poly::new_norm(outbox)
}
}
fn dot_product<'a,'b,T,U>(lhs: impl Iterator<Item=&'a T>, rhs: impl Iterator<Item=&'b U>) -> T
where T: 'a + Clone + Mul<Output=T> + AddAssign,
U: 'b + Clone + Into<T>,
{
let mut zipit = lhs.cloned().zip(rhs.cloned());
let (a,b) = zipit.next().expect("dot product must be with non-empty iterators");
let mut out = a * b.into();
for (a,b) in zipit {
out += a * b.into();
}
out
}
fn classical_slice_mul<T,U>(output: &mut [T], lhs: &[T], rhs: &[U])
where T: Clone + Mul<Output=T> + AddAssign,
U: Clone + Into<T>,
{
if lhs.len() == 0 || rhs.len() == 0 {
assert_eq!(output.len(), 0);
return;
}
assert_eq!(lhs.len() + rhs.len(), output.len() + 1);
let mut i = 0;
while i < lhs.len() && i < rhs.len() {
output[i] = dot_product(lhs[..i+1].iter(), rhs[..i+1].iter().rev());
i = i + 1;
}
let mut j = 1;
while i < lhs.len() {
output[i] = dot_product(lhs[j..i+1].iter(), rhs.iter().rev());
i = i + 1;
j = j + 1;
}
let mut k = 1;
while i < rhs.len() {
output[i] = dot_product(lhs.iter(), rhs[k..i+1].iter().rev());
i = i + 1;
k = k + 1;
}
while i < output.len() {
output[i] = dot_product(lhs[j..].iter(), rhs[k..].iter().rev());
i = i + 1;
j = j + 1;
k = k + 1;
}
}
fn horner_eval<'a,'b,T,U,I>(mut coeffs: I, x: &'b U)
-> Result<U>
where T: 'a + Clone,
U: Clone + Zero + MulAssign + AddAssign,
I: DoubleEndedIterator<Item=&'a T>,
DefConv<T,U>: OneWay<Source=T, Dest=U>,
{
if let Some(leading) = coeffs.next_back() {
let mut out = DefConv::<T,U>::one_way(leading.clone())?;
for coeff in coeffs.rev() {
out *= x.clone();
out += DefConv::<T,U>::one_way(coeff.clone())?;
}
Ok(out)
} else {
Ok(U::zero())
}
}
#[inline]
pub fn refpow<T>(base: &T, exp: usize) -> T
where T: Clone + One + Mul<Output=T> + MulAssign,
{
match exp {
0 => T::one(),
1 => base.clone(),
_ => {
let mut acc = base.clone() * base.clone();
let mut curexp = exp >> 1;
while curexp & 1 == 0 {
acc *= acc.clone();
curexp >>= 1;
}
if curexp > 1 {
let mut basepow = acc.clone() * acc.clone();
curexp >>= 1;
while curexp > 1 {
if curexp & 1 == 1 {
acc *= basepow.clone();
}
basepow *= basepow.clone();
curexp >>= 1;
}
acc *= basepow.clone();
}
if exp & 1 == 1 {
acc *= base.clone();
}
acc
}
}
}
fn bad_linear_solve<'a,'b,M,T>(
matrix: M,
rhs: impl IntoIterator<Item=&'b T>,
close: &impl CloseTo<Item=T>,
) -> Result<Vec<T>>
where M: IntoIterator,
M::Item: IntoIterator<Item=&'a T>,
T: 'a + 'b + Clone + One + Mul<Output=T> + SubAssign + MulAssign + Inv<Output=T>,
{
let mut workmat: Vec<Vec<_>> =
matrix .into_iter()
.map(|row| row.into_iter().map(Clone::clone).collect())
.collect();
let mut sol: Vec<_> = rhs.into_iter().map(Clone::clone).collect();
let n = workmat.len();
assert!(workmat.iter().all(|row| row.len() == n));
assert_eq!(sol.len(), n);
for i in 0..n {
{
let mut j = i;
while close.close_to_zero(&workmat[j][i]) {
j += 1;
if j == n {
return Err(Error::Singular);
}
}
if i != j {
workmat.swap(i, j);
sol.swap(i, j);
}
}
{
let pivot_inverse = workmat[i][i].clone().inv();
for x in workmat[i].split_at_mut(i+1).1 {
*x *= pivot_inverse.clone();
}
sol[i] *= pivot_inverse;
workmat[i][i].set_one();
}
{
let (top, midbot) = workmat.split_at_mut(i);
let (pivrow, bottom) = midbot.split_first_mut().expect("row index i must exist");
let (soltop, solmidbot) = sol.split_at_mut(i);
let (solpiv, solbot) = solmidbot.split_first_mut().expect("row index imust exist");
for (row, solx) in top.iter_mut().chain(bottom.iter_mut())
.zip(soltop.iter_mut().chain(solbot.iter_mut()))
{
if !close.close_to_zero(&row[i]) {
let (left, right) = row.split_at_mut(i+1);
for (j, x) in (i+1..n).zip(right) {
*x -= left[i].clone() * pivrow[j].clone();
}
*solx -= left[i].clone() * solpiv.clone();
}
}
}
}
Ok(sol)
}
fn bad_berlekamp_massey<T>(seq: &[T], close: &impl CloseTo<Item=T>) -> Result<Vec<T>>
where T: Clone + One + Mul<Output=T> + SubAssign + MulAssign + Inv<Output=T>,
{
assert_eq!(seq.len() % 2, 0);
let n = seq.len() / 2;
for k in (1..=n).rev() {
if let Ok(res) = bad_linear_solve(
(0..k).map(|i| &seq[i..i+k]),
&seq[k..2*k],
close)
{
return Ok(res);
}
}
Err(Error::SparsityTooLow)
}
fn bad_trans_vand_solve<'a,'b,T,U>(
roots: U,
rhs: impl IntoIterator<Item=&'b T>,
close: &impl CloseTo<Item=T>)
-> Result<Vec<T>>
where T: 'a + 'b + Clone + One + Mul<Output=T> + SubAssign + MulAssign + Inv<Output=T>,
U: Copy + IntoIterator<Item=&'a T>,
U::IntoIter: ExactSizeIterator,
{
let n = roots.into_iter().len();
let mat: Vec<_> = iter::successors(
Some(iter::repeat_with(One::one).take(n).collect::<Vec<T>>()),
|row| Some(row.iter().cloned().zip(roots.into_iter().cloned()).map(|(x,y)| x * y).collect())
).take(n).collect();
bad_linear_solve(&mat, rhs, close)
}
struct ErrorIter<I,E> {
iter: I,
err: Option<E>,
}
impl<'a,I,E,T> Iterator for &'a mut ErrorIter<I,E>
where
I: 'a + Iterator<Item=core::result::Result<T,E>>,
T: 'a,
E: 'a,
{
type Item = T;
fn next(&mut self) -> Option<Self::Item> {
match self.err {
Some(_) => None,
None => {
match self.iter.next() {
Some(Ok(x)) => Some(x),
Some(Err(e)) => {
self.err = Some(e);
None
}
None => None,
}
}
}
}
}
impl<I,E> ErrorIter<I,E> {
fn new(iter: I) -> Self {
Self {
iter,
err: None,
}
}
fn try_use<F,T>(mut self, f: F) -> core::result::Result<T,E>
where F: FnOnce(&mut Self) -> T
{
let x = f(&mut self);
match self.err {
Some(e) => Err(e),
None => Ok(x),
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn bad_linear_algebra() {
let very_close = RelativeParams::<f64>::default();
{
let mat :Vec<Vec<f64>> = vec![
vec![2., 2., 2.],
vec![-3., -3., -1.],
vec![5., 3., -3.],
];
let rhs = vec![12., -12., 10.];
assert!(very_close.close_to_iter(
bad_linear_solve(&mat, &rhs, &very_close).unwrap().iter(),
[5., -2., 3.].iter()));
}
{
let mat :Vec<Vec<f64>> = vec![
vec![2., 8., 10., 52.],
vec![1., 4., 5., 24.],
vec![-2., -6., -4., -40.],
vec![3., 11., 12., 77.],
];
let rhs = vec![1., 2., 3., 4.];
assert_eq!(bad_linear_solve(&mat, &rhs, &very_close), Err(Error::Singular));
}
{
let seq = vec![1., 0., 5., -2., 12., -1.];
assert!(very_close.close_to_iter(
bad_berlekamp_massey(&seq, &very_close).unwrap().iter(),
[3., 2., -1.].iter()));
}
{
let roots = vec![2., -1., 3.];
let rhs = vec![8.5, 29., 70.];
assert!(very_close.close_to_iter(
bad_trans_vand_solve(&roots, &rhs, &very_close).unwrap().iter(),
[4.5,-2.,6.].iter()));
}
}
#[test]
fn add() {
let a = vec![10, 20, 30, 40];
let b = vec![3, 4, 5];
let c = vec![13, 24, 35, 40];
let ap: ClassicalPoly<_> = a.iter().copied().collect();
let bp: ClassicalPoly<_> = b.iter().copied().collect();
let cp: ClassicalPoly<_> = c.iter().copied().collect();
assert_eq!(ap + bp, cp);
}
#[test]
fn mul() {
let a = vec![1, 2, 3, 4, 5];
let b = vec![300, 4, 50000];
let c = vec![300, 604, 50908, 101212, 151516, 200020, 250000];
let ap: ClassicalPoly<_> = a.iter().copied().collect();
let bp: ClassicalPoly<_> = b.iter().copied().collect();
let cp: ClassicalPoly<_> = c.iter().copied().collect();
assert_eq!(&ap * &bp, cp);
assert_eq!(&bp * &ap, cp);
}
#[test]
fn eval() {
type CP = ClassicalPoly<i16>;
let f = CP::new(vec![-5,3,-1,2]);
assert_eq!(f.eval(&-3.0f32),
Ok((-5 + 3*-3 + -1*-3*-3 + 2*-3*-3*-3) as f32));
{
let pts = vec![-2,0,7];
let prep = CP::mp_eval_prep(pts.iter().copied());
assert!(f.mp_eval(&prep).unwrap().into_iter().eq(
pts.iter().map(|x| f.eval(x).unwrap())));
}
}
#[test]
fn pow() {
assert_eq!(refpow(&3u64, 0), 1);
assert_eq!(refpow(&25u64, 1), 25);
assert_eq!(refpow(&4u64, 2), 16);
assert_eq!(refpow(&-5i32, 3), -125);
assert_eq!(refpow(&2u16, 8), 256);
let mut pow3 = 1u64;
for d in 1..41 {
pow3 *= 3u64;
assert_eq!(refpow(&3u64, d), pow3);
}
}
#[test]
fn classical_sparse_interp_exact() {
type CP = ClassicalPoly<i32>;
let f = CP::new(vec![3, 0, -2, 0, 0, 0, -1]);
let t = 3;
let (eval_info, interp_info) = CP::sparse_interp_prep(t, 0..10, &100);
let evals = f.mp_eval(&eval_info).unwrap();
let expected_sparse: Vec<_> = f.rep.iter().enumerate().filter(|(_,c)| **c != 0).map(|(e,c)| (e,*c)).collect();
let sparse_f = CP::sparse_interp(&evals, &interp_info).unwrap();
assert_eq!(expected_sparse, sparse_f);
}
#[test]
fn classical_sparse_interp_overshoot() {
type CP = ClassicalPoly<i32>;
let f = CP::new(vec![3, 0, -2, 0, 0, 0, -1]);
let t = 5;
let (eval_info, interp_info) = CP::sparse_interp_prep(t, 0..10, &100);
let evals = f.mp_eval(&eval_info).unwrap();
let expected_sparse: Vec<_> = f.rep.iter().enumerate().filter(|(_,c)| **c != 0).map(|(e,c)| (e,*c)).collect();
let sparse_f = CP::sparse_interp(&evals, &interp_info).unwrap();
assert_eq!(expected_sparse, sparse_f);
}
#[test]
fn classical_spinterp_big() {
type CP = ClassicalPoly<f64>;
let mut coeffs = vec![0.0f64; 200];
let expected_sparse: Vec<(usize, f64)> = vec![
(3, 1.9),
(4, -3.),
(31, 18.),
(101, 19.4),
(108, 0.5),
(109, 0.6),
(110, -1.),
(151, -12.6),
(199, 2.5),
];
for (e, c) in expected_sparse.iter() {
coeffs[*e] = *c;
}
let f = CP::new(coeffs.clone());
let (eval_info, interp_info) = CP::sparse_interp_prep(15, 0..200, &20000.);
let evals = f.mp_eval(&eval_info).unwrap();
let (computed_expons, computed_coeffs): (Vec<_>, Vec<_>)
= CP::sparse_interp(&evals, &interp_info).unwrap().into_iter().unzip();
let (expected_expons, expected_coeffs): (Vec<_>, Vec<_>)
= expected_sparse.into_iter().unzip();
assert_eq!(expected_expons, computed_expons);
let close = RelativeParams::<f64>::new(Some(0.001), Some(0.001));
assert!(close.close_to_iter(expected_coeffs.iter(), computed_coeffs.iter()));
}
}