easel3d/fractional/src/continuous.rs

//
// A «continued fraction» is a representation of a fraction as a vector
// of integrals… Irrational fractions will result in infinite most of the
// time repetitive vectors. They can be used to get a resonable approximation
// for sqrt on fractionals.
//
// 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 crate::Fractional;

#[derive(Debug)]
pub struct Continuous (Vec<i64>);

impl Continuous {
    // calculate a sqrt as continued fraction sequence. Taken from:
    // https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#
    //  Continued_fraction_expansion
    pub fn sqrt(x :i64, a0 :i64) -> Self {
        fn inner(v :&mut [i64], x :i64, a0 :i64, mn :i64, dn :i64, an :i64) {
            let mn_1 = dn * an - mn;
            let dn_1 = (x - mn_1 * mn_1) / dn;
            let an_1 = (a0 + mn_1) / dn_1;

            v[0] = an;
            // The convergence criteria „an_1 == 2 * a0“ is not good for
            // very small x thus I decided to break the iteration at constant
            // time. Which is the 5 below.
            if v.len() > 1 {
                inner(&mut v[1..], x, a0, mn_1, dn_1, an_1);
            }
        }

        let mut v :Vec<i64> = vec!(0; 5);
        inner(&mut v, x, a0, 0, 1, a0);
        Continuous(v)
    }

    // general continous fraction form of a fractional...
    pub fn from_prec(f :&Fractional, prec :Option<usize>) -> Self {
        fn inner(v :&mut Vec<i64>, f :Fractional, prec :Option<usize>) {
            let mut process = |prec :Option<usize>| {
                let Fractional(n, d)   = f;
                let a                  = n / d;
                let Fractional(_n, _d) = f.noreduce_sub(a.into());

                v.push(a);
                match _n {
                    1 => v.push(_d),
                    0 => {},
                    _ => inner(v, Fractional(_d, _n), prec),
                }
            };

            match prec {
                Some(0) => {},
                None    => process(None),
                Some(p) => process(Some(p - 1)),
            }
        }

        let mut v = match prec {
            None    => Vec::with_capacity(100),
            Some(p) => Vec::with_capacity(p + 1),
        };

        inner(&mut v, *f, prec);
        Continuous(v)
    }

    pub fn into_prec(&self, prec :Option<usize>) -> Fractional {
        let Continuous(c) = self;
        let             p = match prec {
            Some(p) => if p <= c.len() { p } else { c.len() },
            None    => c.len(),
        };

        let to_frac = |acc :Fractional, x :&i64| {
            let Fractional(an, ad) = acc.noreduce_add((*x).into());
            Fractional(ad, an)
        };

        let Fractional(n, d) = c[..p]
                             . into_iter()
                             . rev()
                             . fold(Fractional(0, 1), to_frac);
        Fractional(d, n)
    }
}

impl From<&Fractional> for Continuous {
    fn from(x :&Fractional) -> Self {
        Self::from_prec(x, None)
    }
}

impl Into<Fractional> for &Continuous {
    fn into(self) -> Fractional {
        (&self).into_prec(None)
    }
}