easel3d-wasm/fractional/src/fractional.rs

//
// Some code to support fractional numbers for full precision rational number
// 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>
//
// Copyright © 2019 Georg Hopp
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// 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::cmp::Ordering;
use std::ops::{Add,Sub,Neg,Mul,Div};
use std::fmt;
use std::convert::{TryFrom, TryInto};
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 {
        0 => x,
        _ => hcf(y, x % y),
    }
}

// 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 {
        let Fractional(_, d1) = self;
        let Fractional(_, d2) = other;
        (d1 * d2) / hcf(d1, d2)
    }

    #[inline]
    pub fn reduce(self) -> Self {
        let Fractional(n, d) = self;
        let (_n, _d) = if n > d { (n, d) } else { (d, n) };

        // if the difference from _n % _d to _n is very big we are close to
        // a whole number and can ignore the fractional part... this reduces
        // the precision but ensures smaller numbers for numerator and
        // denominator.
        if _d > 1 && (_n % _d) * 10000000 < _n {
            Self(_n / _d, 1)
        } else {
            Self(n / hcf(n, d), d / hcf(n, d))
        }
    }

    #[inline]
    pub fn numerator(self) -> i64 {
        self.0
    }

    #[inline]
    pub fn denominator(self) -> i64 {
        self.1
    }

    // 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).
        // 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 {
                    (&sqrt_cfrac(x, l - 1)).into()
                } else {
                    let m = (l + h) / 2;
                    match x.cmp(&(m * m)) {
                        Ordering::Equal   => m.into(),
                        Ordering::Less    => inner(l, m - 1, x),
                        Ordering::Greater => inner(m + 1, h, x),
                    }
                }
            }

            match x {
                0 => 0.into(),
                1 => 1.into(),
                _ => inner(1, x / 2, x),
            }
        }

        let Fractional(n, d) = self;

        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   => 0.into(),
            Ordering::Less    => return Err("sqrt on negative undefined"),
            Ordering::Greater => floor_sqrt(d),
        };

        Ok(n / d)
    }
}

impl fmt::Display for Fractional {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "({}/{})", self.0, self.1)
    }
}

impl From<i64> for Fractional {
    fn from(x: i64) -> Self {
        Self(x, 1)
    }
}

impl TryFrom<usize> for Fractional {
    type Error = &'static str;

    fn try_from(x: usize) -> Result<Self, Self::Error> {
        let v = i64::try_from(x);
        match v {
            Err(_) => Err("Conversion from usize to i32 failed"),
            Ok(_v) => Ok(Self(_v, 1)),
        }
    }
}

pub fn from_vector(xs: &Vec<i64>) -> Vec<Fractional> {
    xs.iter().map(|x| Fractional(*x, 1)).collect()
}

impl TryInto<f64> for Fractional {
    type Error = TryFromIntError;

    fn try_into(self) -> Result<f64, Self::Error> {
        let n :i32 = self.0.try_into()?;
        let d :i32 = self.1.try_into()?;
        Ok(f64::from(n) / f64::from(d))
    }
}

impl TryInto<(i32, i32)> for Fractional {
    type Error = TryFromIntError;

    fn try_into(self) -> Result<(i32, i32), Self::Error> {
        let a :i32 = (self.0 / self.1).try_into()?;
        let b :i32 = (self.0 % self.1).try_into()?;
        Ok((a, b))
    }
}

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;
        let Fractional(n2, d2) = other;
        n1 * (self.gcd(*other) / d1) == n2 * (self.gcd(*other) / d2)
    }
}

impl PartialOrd for Fractional {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for Fractional {
    fn cmp(&self, other: &Self) -> Ordering {
        let Fractional(n1, d1) = self;
        let Fractional(n2, d2) = other;
        let x = n1 * (self.gcd(*other) / d1);
        let y = n2 * (self.gcd(*other) / d2);
        x.cmp(&y)
    }
}

impl Add for Fractional {
    type Output = Self;

    fn add(self, other: Self) -> Self {
        let Fractional(n1, d1) = self;
        let Fractional(n2, d2) = other;
        let n = n1 * (self.gcd(other) / d1) + n2 * (self.gcd(other) / d2);
        Self(n, self.gcd(other)).reduce()
    }
}

impl Sub for Fractional {
    type Output = Self;

    fn sub(self, other: Self) -> Self {
        self + -other
    }
}

impl Neg for Fractional {
    type Output = Self;

    fn neg(self) -> Self {
        let Fractional(n, d) = self;
        Self(-n, d).reduce()
    }
}

impl Mul for Fractional {
    type Output = Self;

    fn mul(self, other :Self) -> Self {
        let Fractional(n1, d1) = self;
        let Fractional(n2, d2) = other;
        Self(n1 * n2, d1 * d2).reduce()
    }
}

impl Div for Fractional {
    type Output = Self;

    fn div(self, other: Self) -> Self {
        let Fractional(n, d) = other;
        self * Fractional(d, n)
    }
}

    /* some stuff that could be tested...
    let x = Fractional(1, 3);
    let y = Fractional(1, 6);

    println!(
        "Greatest common denominator of {} and {}: {}", x, y, x.gcd(y));
    println!("Numerator of {}: {}", x, x.numerator());
    println!("Denominator of {}: {}", x, x.denominator());
    assert_eq!(Fractional(1, 3), Fractional(2, 6));
    assert_eq!(Fractional(1, 3), Fractional(1, 3));
    assert_eq!(y < x, true);
    assert_eq!(y > x, false);
    assert_eq!(x == y, false);
    assert_eq!(x == x, true);
    assert_eq!(x + y, Fractional(1, 2));
    println!("{} + {} = {}", x, y, x + y);
    assert_eq!(x - y, Fractional(1, 6));
    println!("{} - {} = {}", x, y, x - y);
    assert_eq!(y - x, Fractional(-1, 6));
    println!("{} - {} = {}", y, x, y - x);
    assert_eq!(-x, Fractional(-1, 3));
    println!("-{} = {}", x, -x);
    assert_eq!(x * y, Fractional(1, 18));
    println!("{} * {} = {}", x, y, x * y);
    assert_eq!(x / y, Fractional(2, 1));
    println!("{} / {} = {}", x, y, x / y);
    assert_eq!(y / x, Fractional(1, 2));
    println!("{} / {} = {}", y, x, y / x);

    println!("Fractional from 3: {}", Fractional::from(3));
    let z :f64 = Fractional::into(x);
    println!("Floating point of {}: {}", x, z);
    let (d, r) = Fractional::into(x);
    println!("(div, rest) of {}: ({}, {})", x, d, r);
    */