|
| 1 | +use std::fmt::{self, Display, Formatter}; |
| 2 | +use std::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign}; |
| 3 | + |
| 4 | +use kurbo::PathSeg; |
| 5 | + |
| 6 | +/// A struct that represents a polynomial with a maximum degree of `N-1`. |
| 7 | +/// |
| 8 | +/// It provides basic mathematical operations for polynomials like addition, multiplication, differentiation, integration, etc. |
| 9 | +#[derive(Copy, Clone, Debug, PartialEq)] |
| 10 | +pub struct Polynomial<const N: usize> { |
| 11 | + coefficients: [f64; N], |
| 12 | +} |
| 13 | + |
| 14 | +impl<const N: usize> Polynomial<N> { |
| 15 | + /// Create a new polynomial from the coefficients given in the array. |
| 16 | + /// |
| 17 | + /// The coefficient for nth degree is at the nth index in array. Therefore the order of coefficients are reversed than the usual order for writing polynomials mathematically. |
| 18 | + pub fn new(coefficients: [f64; N]) -> Polynomial<N> { |
| 19 | + Polynomial { coefficients } |
| 20 | + } |
| 21 | + |
| 22 | + /// Create a polynomial where all its coefficients are zero. |
| 23 | + pub fn zero() -> Polynomial<N> { |
| 24 | + Polynomial { coefficients: [0.; N] } |
| 25 | + } |
| 26 | + |
| 27 | + /// Return an immutable reference to the coefficients. |
| 28 | + /// |
| 29 | + /// The coefficient for nth degree is at the nth index in array. Therefore the order of coefficients are reversed than the usual order for writing polynomials mathematically. |
| 30 | + pub fn coefficients(&self) -> &[f64; N] { |
| 31 | + &self.coefficients |
| 32 | + } |
| 33 | + |
| 34 | + /// Return a mutable reference to the coefficients. |
| 35 | + /// |
| 36 | + /// The coefficient for nth degree is at the nth index in array. Therefore the order of coefficients are reversed than the usual order for writing polynomials mathematically. |
| 37 | + pub fn coefficients_mut(&mut self) -> &mut [f64; N] { |
| 38 | + &mut self.coefficients |
| 39 | + } |
| 40 | + |
| 41 | + /// Evaluate the polynomial at `value`. |
| 42 | + pub fn eval(&self, value: f64) -> f64 { |
| 43 | + self.coefficients.iter().rev().copied().reduce(|acc, x| acc * value + x).unwrap() |
| 44 | + } |
| 45 | + |
| 46 | + /// Return the same polynomial but with a different maximum degree of `M-1`.\ |
| 47 | + /// |
| 48 | + /// Returns `None` if the polynomial cannot fit in the specified size. |
| 49 | + pub fn as_size<const M: usize>(&self) -> Option<Polynomial<M>> { |
| 50 | + let mut coefficients = [0.; M]; |
| 51 | + |
| 52 | + if M >= N { |
| 53 | + coefficients[..N].copy_from_slice(&self.coefficients); |
| 54 | + } else if self.coefficients.iter().rev().take(N - M).all(|&x| x == 0.) { |
| 55 | + coefficients.copy_from_slice(&self.coefficients[..M]) |
| 56 | + } else { |
| 57 | + return None; |
| 58 | + } |
| 59 | + |
| 60 | + Some(Polynomial { coefficients }) |
| 61 | + } |
| 62 | + |
| 63 | + /// Computes the derivative in place. |
| 64 | + pub fn derivative_mut(&mut self) { |
| 65 | + self.coefficients.iter_mut().enumerate().for_each(|(index, x)| *x *= index as f64); |
| 66 | + self.coefficients.rotate_left(1); |
| 67 | + } |
| 68 | + |
| 69 | + /// Computes the antiderivative at `C = 0` in place. |
| 70 | + /// |
| 71 | + /// Returns `None` if the polynomial is not big enough to accommodate the extra degree. |
| 72 | + pub fn antiderivative_mut(&mut self) -> Option<()> { |
| 73 | + if self.coefficients[N - 1] != 0. { |
| 74 | + return None; |
| 75 | + } |
| 76 | + self.coefficients.rotate_right(1); |
| 77 | + self.coefficients.iter_mut().enumerate().skip(1).for_each(|(index, x)| *x /= index as f64); |
| 78 | + Some(()) |
| 79 | + } |
| 80 | + |
| 81 | + /// Computes the polynomial's derivative. |
| 82 | + pub fn derivative(&self) -> Polynomial<N> { |
| 83 | + let mut ans = *self; |
| 84 | + ans.derivative_mut(); |
| 85 | + ans |
| 86 | + } |
| 87 | + |
| 88 | + /// Computes the antiderivative at `C = 0`. |
| 89 | + /// |
| 90 | + /// Returns `None` if the polynomial is not big enough to accommodate the extra degree. |
| 91 | + pub fn antiderivative(&self) -> Option<Polynomial<N>> { |
| 92 | + let mut ans = *self; |
| 93 | + ans.antiderivative_mut()?; |
| 94 | + Some(ans) |
| 95 | + } |
| 96 | +} |
| 97 | + |
| 98 | +impl<const N: usize> Default for Polynomial<N> { |
| 99 | + fn default() -> Self { |
| 100 | + Self::zero() |
| 101 | + } |
| 102 | +} |
| 103 | + |
| 104 | +impl<const N: usize> Display for Polynomial<N> { |
| 105 | + fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result { |
| 106 | + let mut first = true; |
| 107 | + for (index, coefficient) in self.coefficients.iter().enumerate().rev().filter(|&(_, &coefficient)| coefficient != 0.) { |
| 108 | + if first { |
| 109 | + first = false; |
| 110 | + } else { |
| 111 | + f.write_str(" + ")? |
| 112 | + } |
| 113 | + |
| 114 | + coefficient.fmt(f)?; |
| 115 | + if index == 0 { |
| 116 | + continue; |
| 117 | + } |
| 118 | + f.write_str("x")?; |
| 119 | + if index == 1 { |
| 120 | + continue; |
| 121 | + } |
| 122 | + f.write_str("^")?; |
| 123 | + index.fmt(f)?; |
| 124 | + } |
| 125 | + |
| 126 | + Ok(()) |
| 127 | + } |
| 128 | +} |
| 129 | + |
| 130 | +impl<const N: usize> AddAssign<&Polynomial<N>> for Polynomial<N> { |
| 131 | + fn add_assign(&mut self, rhs: &Polynomial<N>) { |
| 132 | + self.coefficients.iter_mut().zip(rhs.coefficients.iter()).for_each(|(a, b)| *a += b); |
| 133 | + } |
| 134 | +} |
| 135 | + |
| 136 | +impl<const N: usize> Add for &Polynomial<N> { |
| 137 | + type Output = Polynomial<N>; |
| 138 | + |
| 139 | + fn add(self, other: &Polynomial<N>) -> Polynomial<N> { |
| 140 | + let mut output = *self; |
| 141 | + output += other; |
| 142 | + output |
| 143 | + } |
| 144 | +} |
| 145 | + |
| 146 | +impl<const N: usize> Neg for &Polynomial<N> { |
| 147 | + type Output = Polynomial<N>; |
| 148 | + |
| 149 | + fn neg(self) -> Polynomial<N> { |
| 150 | + let mut output = *self; |
| 151 | + output.coefficients.iter_mut().for_each(|x| *x = -*x); |
| 152 | + output |
| 153 | + } |
| 154 | +} |
| 155 | + |
| 156 | +impl<const N: usize> Neg for Polynomial<N> { |
| 157 | + type Output = Polynomial<N>; |
| 158 | + |
| 159 | + fn neg(mut self) -> Polynomial<N> { |
| 160 | + self.coefficients.iter_mut().for_each(|x| *x = -*x); |
| 161 | + self |
| 162 | + } |
| 163 | +} |
| 164 | + |
| 165 | +impl<const N: usize> SubAssign<&Polynomial<N>> for Polynomial<N> { |
| 166 | + fn sub_assign(&mut self, rhs: &Polynomial<N>) { |
| 167 | + self.coefficients.iter_mut().zip(rhs.coefficients.iter()).for_each(|(a, b)| *a -= b); |
| 168 | + } |
| 169 | +} |
| 170 | + |
| 171 | +impl<const N: usize> Sub for &Polynomial<N> { |
| 172 | + type Output = Polynomial<N>; |
| 173 | + |
| 174 | + fn sub(self, other: &Polynomial<N>) -> Polynomial<N> { |
| 175 | + let mut output = *self; |
| 176 | + output -= other; |
| 177 | + output |
| 178 | + } |
| 179 | +} |
| 180 | + |
| 181 | +impl<const N: usize> MulAssign<&Polynomial<N>> for Polynomial<N> { |
| 182 | + fn mul_assign(&mut self, rhs: &Polynomial<N>) { |
| 183 | + for i in (0..N).rev() { |
| 184 | + self.coefficients[i] = self.coefficients[i] * rhs.coefficients[0]; |
| 185 | + for j in 0..i { |
| 186 | + self.coefficients[i] += self.coefficients[j] * rhs.coefficients[i - j]; |
| 187 | + } |
| 188 | + } |
| 189 | + } |
| 190 | +} |
| 191 | + |
| 192 | +impl<const N: usize> Mul for &Polynomial<N> { |
| 193 | + type Output = Polynomial<N>; |
| 194 | + |
| 195 | + fn mul(self, other: &Polynomial<N>) -> Polynomial<N> { |
| 196 | + let mut output = *self; |
| 197 | + output *= other; |
| 198 | + output |
| 199 | + } |
| 200 | +} |
| 201 | + |
| 202 | +/// Returns two [`Polynomial`]s representing the parametric equations for x and y coordinates of the bezier curve respectively. |
| 203 | +/// The domain of both the equations are from t=0.0 representing the start and t=1.0 representing the end of the bezier curve. |
| 204 | +pub fn pathseg_to_parametric_polynomial(segment: PathSeg) -> (Polynomial<4>, Polynomial<4>) { |
| 205 | + match segment { |
| 206 | + PathSeg::Line(line) => { |
| 207 | + let term1 = line.p0 - line.p1; |
| 208 | + (Polynomial::new([line.p0.x, term1.x, 0., 0.]), Polynomial::new([line.p0.y, term1.y, 0., 0.])) |
| 209 | + } |
| 210 | + PathSeg::Quad(quad_bez) => { |
| 211 | + let term1 = 2. * (quad_bez.p1 - quad_bez.p0); |
| 212 | + let term2 = quad_bez.p0 - 2. * quad_bez.p1.to_vec2() + quad_bez.p2.to_vec2(); |
| 213 | + |
| 214 | + (Polynomial::new([quad_bez.p0.x, term1.x, term2.x, 0.]), Polynomial::new([quad_bez.p0.y, term1.y, term2.y, 0.])) |
| 215 | + } |
| 216 | + PathSeg::Cubic(cubic_bez) => { |
| 217 | + let term1 = 3. * (cubic_bez.p1 - cubic_bez.p0); |
| 218 | + let term2 = 3. * (cubic_bez.p2 - cubic_bez.p1) - term1; |
| 219 | + let term3 = cubic_bez.p3 - cubic_bez.p0 - term2 - term1; |
| 220 | + |
| 221 | + ( |
| 222 | + Polynomial::new([cubic_bez.p0.x, term1.x, term2.x, term3.x]), |
| 223 | + Polynomial::new([cubic_bez.p0.y, term1.y, term2.y, term3.y]), |
| 224 | + ) |
| 225 | + } |
| 226 | + } |
| 227 | +} |
| 228 | + |
| 229 | +#[cfg(test)] |
| 230 | +mod test { |
| 231 | + use super::*; |
| 232 | + |
| 233 | + #[test] |
| 234 | + fn evaluation() { |
| 235 | + let p = Polynomial::new([1., 2., 3.]); |
| 236 | + |
| 237 | + assert_eq!(p.eval(1.), 6.); |
| 238 | + assert_eq!(p.eval(2.), 17.); |
| 239 | + } |
| 240 | + |
| 241 | + #[test] |
| 242 | + fn size_change() { |
| 243 | + let p1 = Polynomial::new([1., 2., 3.]); |
| 244 | + let p2 = Polynomial::new([1., 2., 3., 0.]); |
| 245 | + |
| 246 | + assert_eq!(p1.as_size(), Some(p2)); |
| 247 | + assert_eq!(p2.as_size(), Some(p1)); |
| 248 | + |
| 249 | + assert_eq!(p2.as_size::<2>(), None); |
| 250 | + } |
| 251 | + |
| 252 | + #[test] |
| 253 | + fn addition_and_subtaction() { |
| 254 | + let p1 = Polynomial::new([1., 2., 3.]); |
| 255 | + let p2 = Polynomial::new([4., 5., 6.]); |
| 256 | + |
| 257 | + let addition = Polynomial::new([5., 7., 9.]); |
| 258 | + let subtraction = Polynomial::new([-3., -3., -3.]); |
| 259 | + |
| 260 | + assert_eq!(&p1 + &p2, addition); |
| 261 | + assert_eq!(&p1 - &p2, subtraction); |
| 262 | + } |
| 263 | + |
| 264 | + #[test] |
| 265 | + fn multiplication() { |
| 266 | + let p1 = Polynomial::new([1., 2., 3.]).as_size().unwrap(); |
| 267 | + let p2 = Polynomial::new([4., 5., 6.]).as_size().unwrap(); |
| 268 | + |
| 269 | + let multiplication = Polynomial::new([4., 13., 28., 27., 18.]); |
| 270 | + |
| 271 | + assert_eq!(&p1 * &p2, multiplication); |
| 272 | + } |
| 273 | + |
| 274 | + #[test] |
| 275 | + fn derivative_and_antiderivative() { |
| 276 | + let mut p = Polynomial::new([1., 2., 3.]); |
| 277 | + let p_deriv = Polynomial::new([2., 6., 0.]); |
| 278 | + |
| 279 | + assert_eq!(p.derivative(), p_deriv); |
| 280 | + |
| 281 | + p.coefficients_mut()[0] = 0.; |
| 282 | + assert_eq!(p_deriv.antiderivative().unwrap(), p); |
| 283 | + |
| 284 | + assert_eq!(p.antiderivative(), None); |
| 285 | + } |
| 286 | + |
| 287 | + #[test] |
| 288 | + fn display() { |
| 289 | + let p = Polynomial::new([1., 2., 0., 3.]); |
| 290 | + |
| 291 | + assert_eq!(format!("{:.2}", p), "3.00x^3 + 2.00x + 1.00"); |
| 292 | + } |
| 293 | +} |
0 commit comments