Struct find_real_roots_of_polynomial::SturmChain

pub struct SturmChain<K> { /* fields omitted */ }

calculate Sturm chain and find all real roots by Sturm method.

use find_real_roots_of_polynomial::SturmChain;
use polynomial_ring::{polynomial, Polynomial};
use num::{BigRational, BigInt};
let v = [-2, 0, 1].iter().map(|x| BigRational::from(BigInt::from(*x))).collect::<Vec<_>>();
let f = Polynomial::new(v); // f=x^2-2=(x-√2)(x+√2)
let limit = BigRational::new(BigInt::from(1), BigInt::from(16)); // 1/16
let sc = SturmChain::<BigRational>::new(f);
let roots = sc.find_all_real_roots(&limit);
assert_eq!(roots.len(), 2);
assert!(&roots[0].1 - &roots[0].0 < limit);
assert!(&roots[1].1 - &roots[1].0 < limit);
// -√2 ∈ [-93/64, -45/32] = [-1.453125, -1.40625]
//  √2 ∈ [45/32, 93/64] = [1.40625, 1.453125]
assert_eq!(roots[0].0, BigRational::new(BigInt::from(-93), BigInt::from(64)));
assert_eq!(roots[0].1, BigRational::new(BigInt::from(-45), BigInt::from(32)));
assert_eq!(roots[1].0, BigRational::new(BigInt::from(45), BigInt::from(32)));
assert_eq!(roots[1].1, BigRational::new(BigInt::from(93), BigInt::from(64)));

Implementations

impl<K: Sized> SturmChain<K>

pub fn new(poly: Polynomial<K>) -> Self where
    K: Clone + Ord + Zero + One + Mul<Output = K> + Neg<Output = K> + for<'x> AddAssign<&'x K> + for<'x> SubAssign<&'x K> + for<'x> DivAssign<&'x K>,
    for<'x> &'x K: Mul<Output = K> + Div<Output = K>, 

make Sturm chain

Input polynomial is converted to square-free polynomial. Then Sturm chain calcurated.

pub fn strum_method(&self, left: &K, right: &K) -> Vec<(K, K)> where
    K: Clone + Ord + Zero + One + Mul<Output = K> + Div<Output = K> + for<'x> AddAssign<&'x K> + for<'x> MulAssign<&'x K>,
    for<'x> &'x K: Add<Output = K> + Mul<Output = K>, 

This function returns half-open intervals. Each of them has exact one root.

Let polynomial has $n$ real roots in (left, right]. Let the roots be $r_0, r_1, \dots, r_{n-1}$. Where, $r_0 < r_1 < \dots < r_{n-1}$. (Note. input ploynomial was converted to square-free polynomial. So, all root are simple root.) This function returns $n$ half-open intervals. Let this intervals $(l_0, u_0], (l_1, u_1], \dots, (l_{n-1}, u_{n-1}]$. $\forall i, r_i \in (l_i, u_i]$.

use find_real_roots_of_polynomial::SturmChain;
use polynomial_ring::{polynomial, Polynomial};
use num::{BigRational, BigInt, Zero, One};
let v = [0, -1, 0, 1].iter().map(|x| BigRational::from(BigInt::from(*x))).collect::<Vec<_>>();
let f = Polynomial::new(v); // f=x^3-1=(x-1)x(x+1)
let m_two = BigRational::from(BigInt::from(-2));
let p_two = BigRational::from(BigInt::from(2));
let sc = SturmChain::<BigRational>::new(f);
let intervals = sc.strum_method(&m_two, &p_two);
assert_eq!(intervals.len(), 3);
// -1 ∈ (-2, -1]
assert_eq!(intervals[0].0, m_two);
assert_eq!(intervals[0].1, -BigRational::one());
// 0 ∈ (-1, 0]
assert_eq!(intervals[1].0, -BigRational::one());
assert_eq!(intervals[1].1, BigRational::zero());
// 1 ∈ (0, 2]
assert_eq!(intervals[2].0, BigRational::zero());
assert_eq!(intervals[2].1, p_two);

pub fn find_all_real_roots(&self, limit: &K) -> Vec<(K, K)> where
    K: Clone + Ord + Zero + One + Mul<Output = K> + Div<Output = K> + Neg<Output = K> + for<'x> AddAssign<&'x K> + for<'x> MulAssign<&'x K> + for<'x> DivAssign<&'x K>,
    for<'x> &'x K: Add<Output = K> + Sub<Output = K> + Neg<Output = K> + Mul<Output = K> + Div<Output = K>, 

Find all real roots.

size of interval is less than limit.

use find_real_roots_of_polynomial::SturmChain;
use polynomial_ring::{polynomial, Polynomial};
use num::{BigRational, BigInt};
let v = [-1, -1, 1].iter().map(|x| BigRational::from(BigInt::from(*x))).collect::<Vec<_>>();
let f = Polynomial::new(v); // f=x^2-x-1=(x-(1+√5)/2)(x-(1-√5)/2)
let limit = BigRational::new(BigInt::from(1), BigInt::from(8)); // 1/8
let sc = SturmChain::<BigRational>::new(f);
let roots = sc.find_all_real_roots(&limit);
assert_eq!(roots.len(), 2);
assert!(&roots[0].1 - &roots[0].0 < limit);
assert!(&roots[1].1 - &roots[1].0 < limit);
// (1-√5)/2(≒-0.61803) ∈ [-5/8, -9/16] = [-0.625, -0.5625]
// (1+√5)/2(≒ 1.61803) ∈ [45/32, 3/2] = [1.5625, 1.625]
assert_eq!(roots[0].0, BigRational::new(BigInt::from(-5), BigInt::from(8)));
assert_eq!(roots[0].1, BigRational::new(BigInt::from(-9), BigInt::from(16)));
assert_eq!(roots[1].0, BigRational::new(BigInt::from(25), BigInt::from(16)));
assert_eq!(roots[1].1, BigRational::new(BigInt::from(13), BigInt::from(8)));