From 6d356f0c681cb2bb0f26d848ed467f65c49769a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Manuel=20Tom=C3=A1s?= Date: Sat, 9 Jan 2021 19:57:43 -0300 Subject: Implement polynomial gcd This was done to find the common roots between two polynomials. --- src/lib.rs | 124 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- src/main.rs | 4 ++ 2 files changed, 126 insertions(+), 2 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 31f230e..c70a40f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,3 +1,5 @@ +use std::cmp::Ordering; + pub type Number = f32; pub type Poly = Vec; @@ -29,13 +31,13 @@ pub fn lp(l: Box) -> Poly { Lerp::Node(a, b) => { let a = lp(a); let b = lp(b); - let c = poly_sub(&b, &a); + let c = sub(&b, &a); skewed_sum(a, c) } } } -fn poly_sub(a: &Poly, b: &Poly) -> Poly { +fn sub(a: &Poly, b: &Poly) -> Poly { let mut r = a.clone(); for i in 0..r.len() { r[i] -= b[i]; @@ -51,3 +53,121 @@ fn skewed_sum(a: Poly, b: Poly) -> Poly { r.push(b[b.len() - 1]); r } + +fn is_zero(p: &Poly) -> bool { + for i in 0..p.len() { + if p[i] != 0.0 { + return false + } + } + true +} + +fn degree(p: &Poly) -> usize { + let mut i = p.len() - 1; + while p[i] == 0.0 && i > 0 { + i -= 1; + } + i +} + +fn prod(p: &Poly, n: Number) -> Poly { + let mut r = p.clone(); + for i in 0..r.len() { + r[i] *= n; + } + r +} + +fn shift(p: &Poly, amount: usize) -> Poly { + let mut r = vec![0.0; p.len()]; + for i in 0..p.len() { + if i + amount < r.len() { + r[i + amount] = p[i]; + } + } + r +} + +fn rem(a: &Poly, b: &Poly) -> Poly { + let an = a[degree(a)]; + let bn = b[degree(b)]; + match degree(a).cmp(°ree(b)) { + Ordering::Equal => sub(a, &prod(b, an / bn)), + Ordering::Greater => sub(a, &prod(&shift(b, (degree(a) - degree(b)) as usize), an / bn)), + Ordering::Less => vec![0.0; b.len()] + } +} + +pub fn gcd(a: &Poly, b: &Poly) -> Poly { + let c = rem(a, b); + if is_zero(&c) { + b.clone() + } else { + gcd(b, &c) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn poly_is_zero() { + let p = vec![0.0; 5]; + assert_eq!(is_zero(&p), true); + } + + #[test] + fn poly_is_not_zero() { + let mut p = vec![0.0; 5]; + p[2] = 8.0; + assert_eq!(is_zero(&p), false); + } + + #[test] + fn degree_is_five() { + let mut p = vec![0.0; 6]; + p[5] = 2.0; + assert_eq!(degree(&p), 5); + } + + #[test] + fn degree_is_zero() { + let p = vec![0.0; 6]; + assert_eq!(degree(&p), 0); + } + + #[test] + fn prod_test() { + let p = vec![1.0, 2.0, 3.0, 4.0, 5.0]; + assert_eq!(prod(&p, 2.0), vec![2.0, 4.0, 6.0, 8.0, 10.0]); + } + + #[test] + fn shift_test() { + let p = vec![1.0, 2.0, 3.0, 0.0, 0.0]; + assert_eq!(shift(&p, 2), vec![0.0, 0.0, 1.0, 2.0, 3.0]); + } + + #[test] + fn rem_equal() { + let a = vec![6.0, 7.0, 1.0]; + let b = vec![-6.0, -5.0, 1.0]; + assert_eq!(rem(&a, &b), vec![12.0, 12.0, 0.0]); + } + + #[test] + fn rem_greater() { + let a = vec![-6.0, -5.0, 1.0]; + let b = vec![12.0, 12.0, 0.0]; + assert_eq!(rem(&a, &b), vec![-6.0, -6.0, 0.0]); + } + + #[test] + fn gcd_test() { + let a = vec![6.0, 7.0, 1.0]; + let b = vec![-6.0, -5.0, 1.0]; + assert_eq!(gcd(&a, &b), vec![-6.0, -6.0, 0.0]); + } +} diff --git a/src/main.rs b/src/main.rs index 13f13da..8831976 100644 --- a/src/main.rs +++ b/src/main.rs @@ -23,6 +23,7 @@ fn main() { let b = Lerp::new(vec![400.0, 140.0, 500.0, 300.0]); let pa = poly::lp(a); let pb = poly::lp(b); + let p = poly::gcd(&pa, &pb); let sdl_context = sdl2::init().unwrap(); let video_subsystem = sdl_context.video().unwrap(); @@ -51,10 +52,13 @@ fn main() { for t in 0..800 { let x = eval_poly(&pa, t as Number / 800.0); let y = eval_poly(&pb, t as Number / 800.0); + let z = eval_poly(&p, t as Number / 800.0); canvas.set_draw_color(Color::RGB(180, 20, 20)); canvas.fill_rect(Rect::new(t, x as i32, 5, 5)).unwrap(); canvas.set_draw_color(Color::RGB(20, 180, 20)); canvas.fill_rect(Rect::new(t, y as i32, 5, 5)).unwrap(); + canvas.set_draw_color(Color::RGB(20, 20, 180)); + canvas.fill_rect(Rect::new(t, z as i32, 5, 5)).unwrap(); } canvas.present(); -- cgit v1.2.3