diff --git a/fractional/src/fractional.rs b/fractional/src/fractional.rs index 387f907..f72e4f7 100644 --- a/fractional/src/fractional.rs +++ b/fractional/src/fractional.rs @@ -1,9 +1,9 @@ // // Some code to support fractional numbers for full precision rational number -// calculations. -// TODO -// - maybe this could be build as a generic for all integral numbers. -// (Question, how can I assure that it is build from integral numbers? +// calculations. (At least for the standard operations.) +// This also implements a sqrt on fractional numbers, which can not be precise +// because of the irrational nature of most sqare roots. +// Fractions can only represent rational numbers precise. // // Georg Hopp // @@ -31,6 +31,10 @@ use std::num::TryFromIntError; #[derive(Debug, Eq, Clone, Copy)] pub struct Fractional (pub i64, pub i64); +pub type Continuous = Vec; + +pub type Error = &'static str; + #[inline] fn hcf(x :i64, y :i64) -> i64 { match y { @@ -39,6 +43,38 @@ fn hcf(x :i64, y :i64) -> i64 { } } +// calculate a sqrt as continued fraction sequence. Taken from: +// https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Continued_fraction_expansion +fn sqrt_cfrac(x :i64, a0 :i64) -> Continuous { + let v :Continuous = Vec::new(); + + fn inner(mut v :Continuous, + x :i64, + a0 :i64, + mn :i64, + dn :i64, + an :i64) -> Continuous { + let mn_1 = dn * an - mn; + let dn_1 = (x - mn_1 * mn_1) / dn; + let an_1 = (a0 + mn_1) / dn_1; + + v.push(an); + match v.len() { + 10 => v, + _ => inner(v, x, a0, mn_1, dn_1, an_1), + } + + // This convergence criterium is not good for very small x + // thus I decided to break the iteration at constant time. +// match an_1 == 2 * a0 { +// true => v, +// _ => inner(v, x, a0, mn_1, dn_1, an_1), +// } + } + + inner(v, x, a0, 0, 1, a0) +} + impl Fractional { #[inline] pub fn gcd(self, other: Self) -> i64 { @@ -63,31 +99,54 @@ impl Fractional { self.1 } - pub fn sqrt(self) -> Self { + // This is a really bad approximation of sqrt for a fractional... + // for (9/3) it will result 3 which if way to far from the truth, + // which is ~1.7320508075 + // BUT we can use this value as starting guess for creating a + // continous fraction for the sqrt... and create a much better + // fractional representation of the sqrt. + // So, if inner converges, but is not a perfect square (does not + // end up in an Ordering::Equal - which is the l > h case) + // we use the l - 1 as starting guess for sqrt_cfrac. + // taken from: + // https://www.geeksforgeeks.org/square-root-of-an-integer/ + pub fn sqrt(self) -> Result { // find the sqrt of x in O(log x/2). - fn floor_sqrt(x :i64) -> i64 { - fn inner(l :i64, h :i64, x :i64) -> i64 { + // This stops if a perfect sqare was found. Else it passes + // the found value as starting guess to the continous fraction + // sqrt function. + fn floor_sqrt(x :i64) -> Fractional { + fn inner(l :i64, h :i64, x :i64) -> Fractional { if l > h { - l - 1 + (&sqrt_cfrac(x, l - 1)).into() } else { let m = (l + h) / 2; match x.cmp(&(m * m)) { - Ordering::Equal => m, - Ordering::Less => inner(l, m, x), + Ordering::Equal => m.into(), + Ordering::Less => inner(l, m - 1, x), Ordering::Greater => inner(m + 1, h, x), } } } - match x.cmp(&0) { - Ordering::Equal => 0, - Ordering::Less => -inner(1, -x / 2, -x), - Ordering::Greater => inner(1, x / 2, x), - } + inner(1, x / 2, x) } let Fractional(n, d) = self; - Fractional(floor_sqrt(n), floor_sqrt(d)).reduce() + + let n = match n.cmp(&0) { + Ordering::Equal => 0.into(), + Ordering::Less => return Err("sqrt on negative undefined"), + Ordering::Greater => floor_sqrt(n), + }; + + let d = match d.cmp(&0) { + Ordering::Equal => return Err("division by zero"), + Ordering::Less => return Err("sqrt on negative undefined"), + Ordering::Greater => floor_sqrt(d), + }; + + Ok(n / d) } } @@ -139,6 +198,37 @@ impl TryInto<(i32, i32)> for Fractional { } } +impl Into for Fractional { + // general continous fraction form of a fractional... + fn into(self) -> Continuous { + let v :Continuous = Vec::new(); + + fn inner(mut v :Continuous, f :Fractional) -> Continuous { + let Fractional(n, d) = f; + let a = n / d; + let Fractional(_n, _d) = f - a.into(); + + v.push(a); + match _n { + 1 => { v.push(_d); v }, + _ => inner(v, Fractional(_d, _n)), + } + } + + inner(v, self) + } +} + +impl From<&Continuous> for Fractional { + fn from(x :&Continuous) -> Self { + let Self(n, d) = x.iter().rev().fold(Fractional(0, 1), |acc, x| { + let Self(an, ad) = acc + (*x).into(); + Self(ad, an) + }); + Self(d, n) + } +} + impl PartialEq for Fractional { fn eq(&self, other: &Self) -> bool { let Fractional(n1, d1) = self; diff --git a/fractional/src/main.rs b/fractional/src/main.rs index 4539b93..8fa1284 100644 --- a/fractional/src/main.rs +++ b/fractional/src/main.rs @@ -18,11 +18,11 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . // -use std::convert::{TryFrom, TryInto}; +use std::convert::{TryFrom, TryInto, Into}; use std::num::TryFromIntError; use std::f64::consts::PI as FPI; -use fractional::fractional::{Fractional, from_vector}; +use fractional::fractional::{Fractional, from_vector, Continuous}; use fractional::trigonometry::{sin, cos, PI}; struct Vector { @@ -37,84 +37,96 @@ fn mean(v: &Vec) -> Result { Ok(Fractional(r, l)) } -fn main() { - let a = vec![3, 6, 1, 9]; - let b = from_vector(&a); - let c = mean(&a).unwrap(); // This might fail if the len of the - // vector (usize) does not fit into i32. - let d :f64 = c.try_into().unwrap(); - let e :f64 = Fractional::try_into(c).unwrap(); - - let f = Fractional(9, 4); - let g = Fractional(-9, 4); - - println!(" [i32] : {:?}" , a); - println!(" [Fractional] : {:?}" , b); - println!(" mean of [i32] : {}" , c); - println!(" as f64 : {}" , d); - println!(" and as f64 : {}" , e); - println!(" again as f64 : {}" , TryInto::::try_into(c).unwrap()); - println!(" sqrt f : {}" , f.sqrt()); - println!(" sqrt f as f64 : {}" , TryInto::::try_into(f.sqrt()).unwrap()); - println!(" sqrt g : {}" , g.sqrt()); - println!(" sqrt g as f64 : {}" , TryInto::::try_into(g.sqrt()).unwrap()); +fn common_fractional() { + let a = vec![3, 6, 1, 9]; + let b = from_vector(&a); + let c = mean(&a).unwrap(); // This might fail if the len of the + // vector (usize) does not fit into i32. + let cr :f64 = c.try_into().unwrap(); + + println!(" [i32] : {:?}", a); + println!(" [Fractional] : {:?}", b); + println!(" mean of [i32] : {}" , c); + println!(" as f64 : {}" , cr); + println!(" again as f64 : {}" , TryInto::::try_into(c).unwrap()); +} + +fn continuous() { + let d = Fractional(45, 16); + let e = Fractional(16, 45); + + let dc :Continuous = d.into(); + let ec :Continuous = e.into(); + + println!("cont frac of d : {} => {:?}", d, dc); + println!("cont frac of e : {} => {:?}", e, ec); + println!(" reverted dc : {:?} {}", dc, Into::::into(&dc)); + println!(" reverted ec : {:?} {}", ec, Into::::into(&ec)); +} + +fn sqrt() { + let f = Fractional(-9, 4); + let fr :f64 = f.try_into().unwrap(); + let sq = f.sqrt(); + let _sq = f64::sqrt(fr); + + println!("{:>14} : {:?} / {}", format!("sqrt {}", f), sq, _sq); + + for f in [ Fractional(9, 4) + , Fractional(45, 16) + , Fractional(16, 45) + , Fractional(9, 3) ].iter() { + let fr :f64 = (*f).try_into().unwrap(); + let sq = f.sqrt().unwrap(); + let sqr :f64 = sq.try_into().unwrap(); + let _sqr = f64::sqrt(fr); + + println!("{:>14} : {} {} / {}", format!("sqrt {}", f), sq, sqr, _sqr); + } +} + +fn pi() { + let pir :f64 = PI.try_into().unwrap(); + let pit :(i32, i32) = PI.try_into().unwrap(); + let pi2r :f64 = (PI * PI).try_into().unwrap(); + println!(" Rust π : {}" , FPI); - println!(" π : {} {}" , TryInto::::try_into(PI).unwrap(), PI); - println!(" π as tuple : {:?}" , TryInto::<(i32, i32)>::try_into(PI).unwrap()); + println!(" π : {} {}" , PI, pir); + println!(" π as tuple : {:?}" , pit); println!(" Rust π² : {}" , FPI * FPI); - println!(" π² : {} {}" , TryInto::::try_into(PI * PI).unwrap(), PI * PI); - println!(" sin 9 : {}" , sin(9)); - println!(" sin 9 : {}" , TryInto::::try_into(sin(9)).unwrap()); - println!(" Rust sin 9 : {}" , f64::sin(9.0 * FPI / 180.0)); - println!(" sin 17 : {}" , sin(17)); - println!(" sin 17 : {}" , TryInto::::try_into(sin(17)).unwrap()); - println!(" Rust sin 17 : {}" , f64::sin(17.0 * FPI / 180.0)); - println!(" sin 31 : {}" , sin(31)); - println!(" sin 31 : {}" , TryInto::::try_into(sin(31)).unwrap()); - println!(" Rust sin 31 : {}" , f64::sin(31.0 * FPI / 180.0)); - println!(" sin 45 : {}" , sin(45)); - println!(" sin 45 : {}" , TryInto::::try_into(sin(45)).unwrap()); - println!(" Rust sin 45 : {}" , f64::sin(45.0 * FPI / 180.0)); - println!(" sin 73 : {}" , sin(73)); - println!(" sin 73 : {}" , TryInto::::try_into(sin(73)).unwrap()); - println!(" Rust sin 73 : {}" , f64::sin(73.0 * FPI / 180.0)); - println!(" sin 123 : {}" , sin(123)); - println!(" sin 123 : {}" , TryInto::::try_into(sin(123)).unwrap()); - println!(" Rust sin 123 : {}" , f64::sin(123.0 * FPI / 180.0)); - println!(" sin 213 : {}" , sin(213)); - println!(" sin 213 : {}" , TryInto::::try_into(sin(213)).unwrap()); - println!(" Rust sin 213 : {}" , f64::sin(213.0 * FPI / 180.0)); - println!(" sin 312 : {}" , sin(312)); - println!(" sin 312 : {}" , TryInto::::try_into(sin(312)).unwrap()); - println!(" Rust sin 312 : {}" , f64::sin(312.0 * FPI / 180.0)); - println!(" sin 876 : {}" , sin(876)); - println!(" sin 876 : {}" , TryInto::::try_into(sin(876)).unwrap()); - println!(" Rust sin 876 : {}" , f64::sin(876.0 * FPI / 180.0)); - println!(" cos 9 : {}" , cos(9)); - println!(" cos 9 : {}" , TryInto::::try_into(cos(9)).unwrap()); - println!(" Rust cos 9 : {}" , f64::cos(9.0 * FPI / 180.0)); - println!(" cos 17 : {}" , cos(17)); - println!(" cos 17 : {}" , TryInto::::try_into(cos(17)).unwrap()); - println!(" Rust cos 17 : {}" , f64::cos(17.0 * FPI / 180.0)); - println!(" cos 31 : {}" , cos(31)); - println!(" cos 31 : {}" , TryInto::::try_into(cos(31)).unwrap()); - println!(" Rust cos 31 : {}" , f64::cos(31.0 * FPI / 180.0)); - println!(" cos 45 : {}" , cos(45)); - println!(" cos 45 : {}" , TryInto::::try_into(cos(45)).unwrap()); - println!(" Rust cos 45 : {}" , f64::cos(45.0 * FPI / 180.0)); - println!(" cos 73 : {}" , cos(73)); - println!(" cos 73 : {}" , TryInto::::try_into(cos(73)).unwrap()); - println!(" Rust cos 73 : {}" , f64::cos(73.0 * FPI / 180.0)); - println!(" cos 123 : {}" , cos(123)); - println!(" cos 123 : {}" , TryInto::::try_into(cos(123)).unwrap()); - println!(" Rust cos 123 : {}" , f64::cos(123.0 * FPI / 180.0)); - println!(" cos 213 : {}" , cos(213)); - println!(" cos 213 : {}" , TryInto::::try_into(cos(213)).unwrap()); - println!(" Rust cos 213 : {}" , f64::cos(213.0 * FPI / 180.0)); - println!(" cos 312 : {}" , cos(312)); - println!(" cos 312 : {}" , TryInto::::try_into(cos(312)).unwrap()); - println!(" Rust cos 312 : {}" , f64::cos(312.0 * FPI / 180.0)); - println!(" cos 876 : {}" , cos(876)); - println!(" cos 876 : {}" , TryInto::::try_into(cos(876)).unwrap()); - println!(" Rust cos 876 : {}" , f64::cos(876.0 * FPI / 180.0)); + println!(" π² : {} {}" , PI * PI, pi2r); +} + +fn _sin() { + for d in [9, 17, 31, 45, 73, 123, 213, 312, 876].iter() { + let s = sin(*d as i32); + let sr :f64 = s.try_into().unwrap(); + let _s = f64::sin(*d as f64 * FPI / 180.0); + + println!("{:>14} : {} {} / {}", format!("sin {}", d), s, sr, _s); + } +} + +fn _cos() { + for d in [9, 17, 31, 45, 73, 123, 213, 312, 876].iter() { + let s = cos(*d as i32); + let sr :f64 = s.try_into().unwrap(); + let _s = f64::cos(*d as f64 * FPI / 180.0); + + println!("{:>14} : {} {} / {}", format!("cos {}", d), s, sr, _s); + } +} + +fn main() { + common_fractional(); + println!(); + continuous(); + println!(); + sqrt(); + println!(); + pi(); + println!(); + _sin(); + println!(); + _cos(); } diff --git a/fractional/src/trigonometry.rs b/fractional/src/trigonometry.rs index 9f137dd..cc809a1 100644 --- a/fractional/src/trigonometry.rs +++ b/fractional/src/trigonometry.rs @@ -1,5 +1,10 @@ // -// Test our fractional crate / module... +// Some trigonometic functions with Fractions results. +// Currently only sin and cos are implemented. As I was unable +// to find a really good integral approximation for them I +// implement them as tables which are predefined using the +// floating point function f64::sin and then transformed into +// a fraction of a given PRECISION. // // Georg Hopp //