Add continued fraction and use them to get a better sqrt approximation

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 <georg@steffers.org>
 //
@@ -31,6 +31,10 @@ use std::num::TryFromIntError;
 #[derive(Debug, Eq, Clone, Copy)]
 pub struct Fractional (pub i64, pub i64);
 
+pub type Continuous = Vec<i64>;
+
+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<Self, Error> {
         // 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<Continuous> 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;

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 <http://www.gnu.org/licenses/>.
 //
-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<i64>) -> Result<Fractional, TryFromIntError> {
     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::<f64>::try_into(c).unwrap());
-    println!("        sqrt f : {}"     , f.sqrt());
-    println!(" sqrt f as f64 : {}"     , TryInto::<f64>::try_into(f.sqrt()).unwrap());
-    println!("        sqrt g : {}"     , g.sqrt());
-    println!(" sqrt g as f64 : {}"     , TryInto::<f64>::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::<f64>::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::<Fractional>::into(&dc));
+    println!("   reverted ec : {:?} {}", ec, Into::<Fractional>::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::<f64>::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::<f64>::try_into(PI * PI).unwrap(), PI * PI);
-    println!("         sin 9 : {}"     , sin(9));
-    println!("         sin 9 : {}"     , TryInto::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64>::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();
 }

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 <georg@steffers.org>
 //