1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127

/* This file is part of bacon.
 * Copyright (c) Wyatt Campbell.
 *
 * See repository LICENSE for information.
 */

use super::polynomial::Polynomial;
use nalgebra::ComplexField;
use num_traits::FromPrimitive;

mod spline;
pub use spline::*;

/// Create a Lagrange interpolating polynomial.
///
/// Create an nth degree polynomial matching the n points (xs[i], ys[i])
/// using Neville's iterated method for Lagrange polynomials. The result will
/// match no derivatives.
///
/// # Examples
/// ```
/// use bacon_sci::interp::lagrange;
/// use bacon_sci::polynomial::Polynomial;
/// fn example() {
///     let xs: Vec<_> = (0..10).map(|i| i as f64).collect();
///     let ys: Vec<_> = xs.iter().map(|x| x.cos()).collect();
///     let poly = lagrange(&xs, &ys, 1e-6).unwrap();
///     for x in xs {
///         assert!((x.cos() - poly.evaluate(x)).abs() < 0.00001);
///     }
/// }
/// ```
pub fn lagrange<N: ComplexField + FromPrimitive + Copy>(
    xs: &[N],
    ys: &[N],
    tol: N::RealField,
) -> Result<Polynomial<N>, String>
where
    <N as ComplexField>::RealField: FromPrimitive + Copy,
{
    if xs.len() != ys.len() {
        return Err("lagrange: slices have mismatched dimension".to_owned());
    }

    let mut qs = vec![Polynomial::with_tolerance(tol)?; xs.len() * xs.len()];
    for (ind, y) in ys.iter().enumerate() {
        qs[ind] = polynomial![*y];
    }

    for i in 1..xs.len() {
        let mut poly_2 = polynomial![N::one(), -xs[i]];
        poly_2.set_tolerance(tol)?;
        for j in 1..=i {
            let mut poly_1 = polynomial![N::one(), -xs[i - j]];
            poly_1.set_tolerance(tol)?;
            let idenom = N::one() / (xs[i] - xs[i - j]);
            let numer =
                &poly_1 * &qs[i + xs.len() * (j - 1)] - &poly_2 * &qs[(i - 1) + xs.len() * (j - 1)];
            qs[i + xs.len() * j] = numer * idenom;
        }
    }

    for i in 0..=qs[xs.len() * xs.len() - 1].order() {
        if qs[xs.len() * xs.len() - 1].get_coefficient(i).abs() < tol {
            qs[xs.len() * xs.len() - 1].purge_coefficient(i);
        }
    }

    qs[xs.len() * xs.len() - 1].purge_leading();
    Ok(qs[xs.len() * xs.len() - 1].clone())
}

pub fn hermite<N: ComplexField + FromPrimitive + Copy>(
    xs: &[N],
    ys: &[N],
    derivs: &[N],
    tol: N::RealField,
) -> Result<Polynomial<N>, String>
where
    <N as ComplexField>::RealField: FromPrimitive + Copy,
{
    if xs.len() != ys.len() {
        return Err("hermite: slices have mismatched dimension".to_owned());
    }
    if xs.len() != derivs.len() {
        return Err("hermite: derivatives have mismatched dimension".to_owned());
    }

    let mut zs = vec![N::zero(); 2 * xs.len()];
    let mut qs = vec![N::zero(); 4 * xs.len() * xs.len()];

    for i in 0..xs.len() {
        zs[2 * i] = xs[i];
        zs[2 * i + 1] = xs[i];
        qs[2 * i] = ys[i];
        qs[2 * i + 1] = ys[i];
        qs[2 * i + 1 + (2 * xs.len())] = derivs[i];

        if i != 0 {
            qs[2 * i + (2 * xs.len())] = (qs[2 * i] - qs[2 * i - 1]) / (zs[2 * i] - zs[2 * i - 1]);
        }
    }

    for i in 2..2 * xs.len() {
        for j in 2..=i {
            qs[i + j * (2 * xs.len())] = (qs[i + (j - 1) * (2 * xs.len())]
                - qs[i - 1 + (j - 1) * (2 * xs.len())])
                / (zs[i] - zs[i - j]);
        }
    }

    let mut hermite = polynomial![N::zero()];
    for i in (1..2 * xs.len()).rev() {
        hermite += qs[i + i * (2 * xs.len())];
        hermite *= polynomial![N::one(), -xs[(i - 1) / 2]];
    }
    hermite += qs[0];

    for i in 0..=hermite.order() {
        if hermite.get_coefficient(i).abs() < tol {
            hermite.purge_coefficient(i);
        }
    }

    hermite.purge_leading();
    Ok(hermite)
}