num_irrational/cont_frac/
simple.rs

1//! Implementation of simple continued fractions
2
3use super::block::Block;
4use super::infinite::InfiniteContinuedFraction;
5use crate::quadratic::{QuadraticBase, QuadraticSurd};
6use crate::traits::{Approximation, Computable, WithSigned, WithUnsigned};
7use core::convert::TryFrom;
8use core::str::FromStr;
9use num_integer::Integer;
10use num_rational::Ratio;
11use num_traits::{CheckedAdd, CheckedMul, Num, NumRef, One, RefNum, Signed, Zero};
12use std::fmt;
13use std::ops::{Add, AddAssign, Div, Mul, Neg, Sub};
14
15/// This struct represents a simple continued fraction `a0 + 1/(a1 + 1/ (a2 + ...))`
16/// Where a0 is an signed integer, a1, a2, .. are positive integers
17/// It's capable of representing rational numbers and quadratic surds
18#[derive(Clone, Debug, PartialEq)]
19pub struct ContinuedFraction<T> {
20    /// Coefficients of aperiodic part
21    a_coeffs: Vec<T>,
22
23    /// Coefficients of periodic part
24    p_coeffs: Vec<T>,
25
26    /// Sign of the fraction
27    negative: bool,
28}
29
30impl<T> ContinuedFraction<T> {
31    #[inline]
32    pub fn aperiodic_coeffs(&self) -> &[T] {
33        &self.a_coeffs[..]
34    }
35
36    #[inline]
37    pub fn periodic_coeffs(&self) -> &[T] {
38        &self.p_coeffs[..]
39    }
40
41    #[inline]
42    pub fn is_negative(&self) -> bool {
43        self.negative
44    }
45
46    #[inline]
47    pub fn is_rational(&self) -> bool {
48        self.p_coeffs.len() == 0
49    }
50
51    #[inline]
52    pub fn is_integer(&self) -> bool {
53        self.a_coeffs.len() == 1 && self.p_coeffs.len() == 0
54    }
55}
56
57impl<U> ContinuedFraction<U> {
58    /// Create a continued fraction from
59    ///
60    /// This function will make sure that if two numbers are equal,
61    /// their internal representation in continued fraction will be the same
62    pub fn new<T: Num + PartialOrd + AddAssign + WithUnsigned<Unsigned = U>>(
63        a_coeffs: Vec<T>,
64        p_coeffs: Vec<T>,
65        negate: bool,
66    ) -> Self {
67        if a_coeffs.len() == 0 && p_coeffs.len() == 0 {
68            panic!("at least one coefficient is required!");
69        }
70
71        let mut negate = negate; // carry the flag when we process through the coefficients
72        let mut negative = None; // final sign for the number
73        let mut dedup_a: Vec<T> = Vec::with_capacity(a_coeffs.len());
74        let mut it = a_coeffs.into_iter();
75
76        // iteratively consume aperiodic coefficients
77        loop {
78            match it.next() {
79                Some(a) => {
80                    // apply negate flag
81                    let a = if negate { T::zero() - a } else { a };
82
83                    // special cases
84                    if dedup_a.len() == 0 {
85                        // if first element
86                        // make sure the first element is non-negative
87                        if a < T::zero() {
88                            dedup_a.push(T::zero() - a);
89                            negative = Some(true);
90                            negate = !negate;
91                        } else {
92                            if !a.is_zero() {
93                                negative = Some(false);
94                            }
95                            dedup_a.push(a);
96                        }
97                        continue;
98                    } else if a.is_zero() {
99                        // if a is zero
100                        if dedup_a.last().unwrap().is_zero() {
101                            dedup_a.pop(); // remove consecutive 2 zeros
102                        } else {
103                            dedup_a.push(a); // dealt by the next value
104                        }
105                        continue;
106                    }
107
108                    // make sure all coefficients are non-negative
109                    if a < T::zero() {
110                        // In this case we use the following formula
111                        // [.., a_{n-1}, -a_n, a_{n+1}, ..] = [.., a_{n-1}-1, 1, a_n-1, -[a_{n+1}, ..]]
112                        let mut target = a;
113                        loop {
114                            let am1 = dedup_a.pop().unwrap(); // a_{n-1}
115                            let has_am2 = dedup_a.len() != 0; // has a_{n-2}
116                            debug_assert!(am1 >= T::zero());
117
118                            if am1.is_one() {
119                                if has_am2 {
120                                    // [.., a_{n-2}, 1-1, 1, a_n-1, -[a_{n+1}, ..]]
121                                    // = [.., a_{n-2}+1, a_n-1, -[a_{n+1}, ..]]
122                                    *dedup_a.last_mut().unwrap() += T::one();
123                                    dedup_a.push(T::zero() - target - T::one());
124                                } else {
125                                    // a_{n-1}-1 = 0 will be the first term
126                                    dedup_a.push(T::zero());
127                                    dedup_a.push(T::one());
128                                    dedup_a.push(T::zero() - target - T::one());
129                                    negative = Some(false);
130                                }
131                                negate = !negate;
132                                break;
133                            } else if am1.is_zero() {
134                                if has_am2 {
135                                    // [.., a_{n-2}, 0-1, 1, a_n-1, -[a_{n+1}, ..]]
136                                    // = [.., a_{n-2}-1, 1, 0, -[1, a_n-1, -[a_{n+1}, ..]]]
137                                    // = [.., a_{n-2}-a_n, a_{n+1}, ..]
138                                    let am2 = dedup_a.pop().unwrap();
139                                    let new_target = am2 + target;
140                                    if new_target < T::zero() {
141                                        // propagate to previous terms
142                                        target = new_target;
143                                    } else {
144                                        dedup_a.push(new_target);
145                                        break;
146                                    }
147                                } else {
148                                    // a_n is the first non-zero coefficient, flip the sign
149                                    dedup_a.push(T::zero());
150                                    dedup_a.push(T::zero() - target);
151                                    negative = Some(true);
152                                    negate = !negate;
153                                    break;
154                                }
155                            } else {
156                                // normal case
157                                debug_assert!(am1 > T::one());
158                                dedup_a.push(am1 - T::one());
159                                dedup_a.push(T::one());
160                                dedup_a.push(T::zero() - target - T::one());
161                                negate = !negate;
162                                break;
163                            }
164                        }
165                    } else {
166                        if dedup_a.last().unwrap().is_zero() && dedup_a.len() > 1 {
167                            // combine current value with the value before zero
168                            dedup_a.pop();
169                            *dedup_a.last_mut().unwrap() += a;
170                        } else {
171                            debug_assert!(!a.is_zero());
172                            if matches!(negative, None) {
173                                negative = Some(a < T::zero());
174                            }
175                            dedup_a.push(a);
176                        }
177                    }
178                }
179                None => {
180                    break;
181                }
182            }
183        }
184
185        if dedup_a.len() == 0 && p_coeffs.len() == 0 {
186            panic!("no effective coefficient!")
187        }
188
189        // when the last term is one or still zero
190        if p_coeffs.len() == 0 {
191            while dedup_a.len() >= 2 {
192                let last_a = dedup_a.last().unwrap();
193                if last_a.is_zero() {
194                    dedup_a.pop();
195                    dedup_a.pop();
196                } else if last_a.is_one() {
197                    let one = dedup_a.pop().unwrap();
198                    *dedup_a.last_mut().unwrap() += one;
199                } else {
200                    break;
201                }
202            }
203
204            if dedup_a.len() == 1 && dedup_a[0].is_zero() {
205                negative = Some(false);
206            }
207        }
208
209        if matches!(negative, None) {
210            if matches!(p_coeffs.first(), Some(p) if p <= &T::zero()) {
211                // TODO: how to handle negative coefficients in the periodic parts
212                unimplemented!()
213            } else {
214                negative = Some(negate);
215            }
216        }
217
218        // collect the results
219        let a_coeffs: Vec<U> = dedup_a.into_iter().map(|v| v.to_unsigned()).collect();
220        let p_coeffs: Vec<U> = p_coeffs.into_iter().map(|v| v.to_unsigned()).collect();
221        ContinuedFraction {
222            a_coeffs,
223            p_coeffs,
224            negative: negative.unwrap(),
225        }
226    }
227}
228
229/// Iterator of coeffcients in a [ContinuedFraction]
230#[derive(Debug, Clone)]
231pub struct Coefficients<'a, T> {
232    a_iter: Option<std::slice::Iter<'a, T>>, // None if aperiodic part has been consumed
233    p_ref: &'a Vec<T>,
234    p_iter: Option<std::slice::Iter<'a, T>>, // None before aperiodic part is consumed, or when periodic part is empty
235}
236
237impl<'a, T> Iterator for Coefficients<'a, T> {
238    type Item = &'a T;
239
240    fn next(&mut self) -> Option<Self::Item> {
241        if let Some(it) = self.a_iter.as_mut() {
242            // in aperiodic part
243            match it.next() {
244                Some(v) => Some(v),
245                None => {
246                    self.a_iter = None;
247                    if self.p_ref.len() > 0 {
248                        let mut new_iter = self.p_ref.iter();
249                        let result = new_iter.next();
250                        self.p_iter = Some(new_iter);
251                        result
252                    } else {
253                        None
254                    }
255                }
256            }
257        } else {
258            if let Some(it) = self.p_iter.as_mut() {
259                // in periodic part
260                match it.next() {
261                    Some(v) => Some(v),
262                    None => {
263                        let mut new_iter = self.p_ref.iter();
264                        let result = new_iter.next();
265                        self.p_iter = Some(new_iter);
266                        result
267                    }
268                }
269            } else {
270                None
271            }
272        }
273    }
274}
275
276/// Iterator that converts coefficients of simple continued fraction to
277/// coefficients of a general continued fraction. Then it can be consumed by
278/// [GeneralContinuedFraction][crate::GeneralContinuedFraction]
279pub struct GeneralCoefficients<T> {
280    coeffs: T,
281    negative: bool, // store the sign
282}
283
284impl<'r, I: Iterator<Item = &'r T>, T: 'r + WithSigned<Signed = U> + Clone, U: Signed> Iterator
285    for GeneralCoefficients<I>
286{
287    type Item = (U, U);
288
289    fn next(&mut self) -> Option<Self::Item> {
290        if self.negative {
291            self.coeffs
292                .next()
293                .map(|v| (U::one(), -v.clone().to_signed()))
294        } else {
295            self.coeffs
296                .next()
297                .map(|v| (U::one(), v.clone().to_signed()))
298        }
299    }
300}
301
302/// Iterator of convergents of a [ContinuedFraction]
303pub struct Convergents<'a, T> {
304    coeffs: Coefficients<'a, T>,
305    block: Block<T>,
306    neg: bool, // store the sign
307}
308
309impl<
310        'a,
311        T: Integer + Clone + CheckedAdd + CheckedMul + WithSigned<Signed = U>,
312        U: Integer + Clone + Signed,
313    > Iterator for Convergents<'a, T>
314{
315    type Item = Ratio<U>;
316
317    fn next(&mut self) -> Option<Self::Item> {
318        let a = self.coeffs.next()?;
319        let (p, q) = self.block.checked_rmove(a.clone())?;
320        self.block.update(p.clone(), q.clone());
321
322        let r = Ratio::new(p.to_signed(), q.to_signed());
323        if self.neg {
324            Some(-r)
325        } else {
326            Some(r)
327        }
328    }
329}
330
331/// Iterator of coefficients in a [ContinuedFraction] with sign applied. This iterator
332/// can be used to construct a [InfiniteContinuedFraction]
333pub struct SignedCoefficients<T> {
334    coeffs: T,
335    negative: bool, // store the sign
336}
337
338impl<'r, I: Iterator<Item = &'r T>, T: 'r + WithSigned<Signed = U> + Clone, U: Signed> Iterator
339    for SignedCoefficients<I>
340{
341    type Item = U;
342
343    fn next(&mut self) -> Option<Self::Item> {
344        if self.negative {
345            self.coeffs.next().map(|v| -v.clone().to_signed())
346        } else {
347            self.coeffs.next().map(|v| v.clone().to_signed())
348        }
349    }
350}
351
352impl<T> ContinuedFraction<T> {
353    /// Returns an iterator of the coefficients in the continued fraction
354    /// Note that for a negative number, the coefficients of it's absolute value is returned
355    pub fn coeffs(&self) -> Coefficients<T> {
356        Coefficients {
357            a_iter: Some(self.a_coeffs.iter()),
358            p_ref: &self.p_coeffs,
359            p_iter: None,
360        }
361    }
362
363    /// Returns an iterator of generalized coefficients, that can be consumed
364    /// by GeneralContinuedFraction
365    pub fn generalized(&self) -> GeneralCoefficients<Coefficients<T>> {
366        GeneralCoefficients {
367            coeffs: self.coeffs(),
368            negative: self.negative,
369        }
370    }
371
372    /// Returns an iterator of the coefficients in the continued fraction.
373    /// The coefficients will be negative if the number is negative
374    pub fn coeffs_signed(&self) -> SignedCoefficients<Coefficients<T>> {
375        SignedCoefficients {
376            coeffs: self.coeffs(),
377            negative: self.negative,
378        }
379    }
380}
381
382impl<T: WithSigned<Signed = U> + Clone, U: Signed> ContinuedFraction<T> {
383    /// Wrap the continued fraction object as `InfiniteContinuedFraction`
384    pub fn expanded(&self) -> InfiniteContinuedFraction<SignedCoefficients<Coefficients<T>>> {
385        self.coeffs_signed().into()
386    }
387}
388
389impl<
390        T: Integer + Clone + CheckedAdd + CheckedMul + WithSigned<Signed = U>,
391        U: Integer + Clone + Signed,
392    > ContinuedFraction<T>
393{
394    /// Returns an iterator of the convergents. The iterator will stop
395    /// if all coefficients are consumed, or numeric overflow happened.
396    pub fn convergents(&self) -> Convergents<T> {
397        Convergents {
398            coeffs: self.coeffs(),
399            block: Block::identity(),
400            neg: self.negative,
401        }
402    }
403
404    /// Converts to an integer by using the first coefficient
405    #[inline]
406    pub fn to_integer(&self) -> Approximation<U> {
407        let a = self.a_coeffs.first().unwrap().clone().to_signed();
408        let a = if self.negative { -a } else { a };
409        if self.is_integer() {
410            Approximation::Exact(a)
411        } else {
412            Approximation::Approximated(a)
413        }
414    }
415
416    #[inline]
417    /// This method returns the corresponding rational number if it's rational,
418    /// returns the expansion until the first repeating occurence
419    pub fn to_rational(&self) -> Approximation<Ratio<U>> {
420        if self.is_rational() {
421            Approximation::Exact(self.convergents().last().unwrap())
422        } else {
423            Approximation::Approximated(
424                self.convergents()
425                    .nth(self.a_coeffs.len() + self.p_coeffs.len())
426                    .unwrap(),
427            )
428        }
429    }
430}
431
432impl<
433        T: Integer + Clone + CheckedAdd + CheckedMul + WithSigned<Signed = U>,
434        U: Integer + Clone + Signed + CheckedAdd,
435    > Computable<U> for ContinuedFraction<T>
436{
437    fn approximated(&self, limit: &U) -> Approximation<Ratio<U>> {
438        let mut convergents = self.convergents();
439        let mut last_conv = convergents.next().unwrap();
440        if last_conv.denom() > limit {
441            let i = self.a_coeffs.first().unwrap().clone();
442            return Approximation::Approximated(Ratio::from(i.to_signed()));
443        }
444        loop {
445            last_conv = match convergents.next() {
446                Some(v) => {
447                    if v.denom() < limit {
448                        v
449                    } else {
450                        return Approximation::Approximated(last_conv);
451                    }
452                }
453                None => return Approximation::Exact(last_conv),
454            }
455        }
456    }
457}
458
459impl<T: fmt::Display> fmt::Display for ContinuedFraction<T> {
460    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
461        if self.negative {
462            write!(f, "-")?;
463        }
464
465        write!(f, "[{}", self.a_coeffs.first().unwrap())?;
466        if self.a_coeffs.len() == 1 {
467            if self.p_coeffs.len() == 0 {
468                return write!(f, "]");
469            } else {
470                write!(f, "; ")?;
471            }
472        } else {
473            let mut aiter = self.a_coeffs.iter().skip(1);
474            write!(f, "; {}", aiter.next().unwrap())?;
475            while let Some(v) = aiter.next() {
476                write!(f, ", {}", v)?;
477            }
478            if self.p_coeffs.len() > 0 {
479                write!(f, ", ")?;
480            }
481        }
482
483        if self.p_coeffs.len() > 0 {
484            let mut piter = self.p_coeffs.iter();
485            write!(f, "({}", piter.next().unwrap())?;
486            while let Some(v) = piter.next() {
487                write!(f, ", {}", v)?;
488            }
489            write!(f, ")]")
490        } else {
491            write!(f, "]")
492        }
493    }
494}
495
496// trait bound free implementations for several simple functions
497impl<T: Zero> ContinuedFraction<T> {
498    #[inline]
499    fn _zero() -> Self {
500        ContinuedFraction {
501            a_coeffs: vec![T::zero()],
502            p_coeffs: Vec::new(),
503            negative: false,
504        }
505    }
506
507    #[inline]
508    fn _is_zero(&self) -> bool {
509        self.a_coeffs.len() == 1 && self.a_coeffs[0].is_zero() && self.p_coeffs.len() == 0
510    }
511}
512
513impl<
514        T: QuadraticBase + AddAssign + WithUnsigned<Unsigned = U>,
515        U: Integer + Clone + NumRef + CheckedAdd + CheckedMul + WithSigned<Signed = T>,
516    > Zero for ContinuedFraction<U>
517where
518    for<'r> &'r T: RefNum<T>,
519    for<'r> &'r U: RefNum<U>,
520{
521    #[inline]
522    fn zero() -> Self {
523        Self::_zero()
524    }
525
526    #[inline]
527    fn is_zero(&self) -> bool {
528        self._is_zero()
529    }
530}
531
532impl<
533        T: QuadraticBase + AddAssign + WithUnsigned<Unsigned = U>,
534        U: Integer + Clone + NumRef + CheckedAdd + CheckedMul + WithSigned<Signed = T>,
535    > One for ContinuedFraction<U>
536where
537    for<'r> &'r T: RefNum<T>,
538    for<'r> &'r U: RefNum<U>,
539{
540    fn one() -> Self {
541        ContinuedFraction {
542            a_coeffs: vec![U::one()],
543            p_coeffs: Vec::new(),
544            negative: false,
545        }
546    }
547
548    fn is_one(&self) -> bool {
549        self.a_coeffs.len() == 1 && self.a_coeffs[0].is_one() && self.p_coeffs.len() == 0
550    }
551}
552
553// TODO: implement from_float, using FloatCore trait?
554
555impl<T: Num + WithUnsigned<Unsigned = U> + PartialOrd, U> From<T> for ContinuedFraction<U> {
556    fn from(t: T) -> Self {
557        if t < T::zero() {
558            ContinuedFraction {
559                a_coeffs: vec![(T::zero() - t).to_unsigned()],
560                p_coeffs: Vec::new(),
561                negative: true,
562            }
563        } else {
564            ContinuedFraction {
565                a_coeffs: vec![t.to_unsigned()],
566                p_coeffs: Vec::new(),
567                negative: false,
568            }
569        }
570    }
571}
572
573impl<T: Integer + Clone + WithUnsigned<Unsigned = U>, U: Zero> From<Ratio<T>>
574    for ContinuedFraction<U>
575{
576    fn from(r: Ratio<T>) -> Self {
577        if r.is_zero() {
578            return Self::_zero();
579        }
580
581        let mut coeffs = Vec::new();
582        let (mut n, mut d) = r.into();
583
584        let negative = if n < T::zero() {
585            n = T::zero() - n;
586            true
587        } else {
588            false
589        };
590
591        if n < d {
592            std::mem::swap(&mut n, &mut d);
593            coeffs.push(U::zero());
594        }
595
596        while !d.is_zero() {
597            let (quo, rem) = n.div_rem(&d);
598            coeffs.push(quo.to_unsigned());
599            n = d;
600            d = rem;
601        }
602
603        ContinuedFraction {
604            a_coeffs: coeffs,
605            p_coeffs: Vec::new(),
606            negative,
607        }
608    }
609}
610
611// TODO: expose this after implemented
612mod parse {
613    use super::{ContinuedFraction, FromStr};
614
615    pub enum ParseContFracError {}
616
617    impl<T> FromStr for ContinuedFraction<T> {
618        type Err = ParseContFracError;
619
620        /// Parse from standard format (like 355/113 = "[3; 7, 16]", (1+sqrt(5))/2 = "[1; (1)]")
621        fn from_str(s: &str) -> Result<Self, Self::Err> {
622            unimplemented!()
623        }
624    }
625}
626
627impl<
628        T: Num + PartialOrd + AddAssign + Signed + WithUnsigned<Unsigned = U>,
629        U: NumRef + PartialOrd + WithSigned<Signed = T>,
630    > Add<T> for ContinuedFraction<U>
631where
632    for<'r> &'r U: RefNum<U>,
633{
634    type Output = Self;
635
636    fn add(self, rhs: T) -> Self {
637        let mut new_a = self.a_coeffs;
638        let i = new_a.first().unwrap();
639
640        // find if the signed is flipped after adding
641        let (new_i, flipped) = match (self.negative, rhs.is_negative()) {
642            (true, true) => {
643                // neg + neg
644                let rhs = (-rhs).to_unsigned();
645                (rhs + i, false)
646            }
647            (true, false) => {
648                // neg + pos
649                let rhs = rhs.to_unsigned();
650                if i < &rhs {
651                    (rhs - i, true)
652                } else {
653                    (i - rhs, false)
654                }
655            }
656            (false, true) => {
657                // pos + neg
658                let rhs = (-rhs).to_unsigned();
659                if i < &rhs {
660                    (rhs - i, true)
661                } else {
662                    (i - rhs, false)
663                }
664            }
665            (false, false) => {
666                // pos + pos
667                (i + rhs.to_unsigned(), false)
668            }
669        };
670
671        if flipped {
672            // delegate to the constructor to handle the negative sign
673            let mut new_a: Vec<T> = new_a.into_iter().map(|u| u.to_signed()).collect();
674            let new_p: Vec<T> = self.p_coeffs.into_iter().map(|u| u.to_signed()).collect();
675            *new_a.first_mut().unwrap() = -new_i.to_signed();
676            Self::new(new_a, new_p, self.negative)
677        } else {
678            *new_a.first_mut().unwrap() = new_i;
679
680            let mut result = ContinuedFraction {
681                a_coeffs: new_a,
682                p_coeffs: self.p_coeffs,
683                negative: self.negative,
684            };
685            if result._is_zero() {
686                // clear negative flag for 0
687                result.negative = false;
688            }
689            result
690        }
691    }
692}
693
694impl<
695        T: Num + PartialOrd + AddAssign + Signed + WithUnsigned<Unsigned = U>,
696        U: NumRef + PartialOrd + WithSigned<Signed = T>,
697    > Sub<T> for ContinuedFraction<U>
698where
699    for<'r> &'r U: RefNum<U>,
700{
701    type Output = Self;
702
703    fn sub(self, rhs: T) -> Self::Output {
704        self.add(-rhs)
705    }
706}
707
708impl<T: Zero> Neg for ContinuedFraction<T> {
709    type Output = ContinuedFraction<T>;
710
711    fn neg(self) -> Self::Output {
712        let mut result = self;
713        if !result._is_zero() {
714            // don't negate when the number is zero
715            result.negative = !result.negative;
716        }
717        result
718    }
719}
720
721impl<
722        T: QuadraticBase + AddAssign + WithUnsigned<Unsigned = U>,
723        U: Integer + Clone + NumRef + CheckedAdd + CheckedMul + WithSigned<Signed = T>,
724    > Mul<T> for ContinuedFraction<U>
725where
726    for<'r> &'r T: RefNum<T>,
727    for<'r> &'r U: RefNum<U>,
728{
729    type Output = Self;
730
731    fn mul(self, rhs: T) -> Self {
732        if self.is_rational() {
733            Self::from(self.to_rational().value() * rhs)
734        } else {
735            Self::try_from(QuadraticSurd::<T>::from(self) * rhs).unwrap()
736        }
737    }
738}
739
740impl<
741        T: QuadraticBase + AddAssign + WithUnsigned<Unsigned = U>,
742        U: Integer + Clone + NumRef + CheckedAdd + CheckedMul + WithSigned<Signed = T>,
743    > Div<T> for ContinuedFraction<U>
744where
745    for<'r> &'r T: RefNum<T>,
746    for<'r> &'r U: RefNum<U>,
747{
748    type Output = Self;
749
750    fn div(self, rhs: T) -> Self {
751        if self.is_rational() {
752            Self::from(self.to_rational().value() / rhs)
753        } else {
754            Self::try_from(QuadraticSurd::<T>::from(self) / rhs).unwrap()
755        }
756    }
757}
758
759macro_rules! impl_binop_for_ratio_surd {
760    (impl $imp:ident, $method:ident) => {
761        impl<
762                T: QuadraticBase + AddAssign + WithUnsigned<Unsigned = U>,
763                U: Integer + Clone + NumRef + CheckedAdd + CheckedMul + WithSigned<Signed = T>,
764            > $imp<Ratio<T>> for ContinuedFraction<U>
765        where
766            for<'r> &'r T: RefNum<T>,
767            for<'r> &'r U: RefNum<U>,
768        {
769            type Output = Self;
770
771            fn $method(self, rhs: Ratio<T>) -> Self {
772                if self.is_rational() {
773                    Self::from(self.to_rational().value().$method(rhs))
774                } else {
775                    Self::try_from(QuadraticSurd::<T>::from(self).$method(rhs)).unwrap()
776                }
777            }
778        }
779
780        impl<
781                T: QuadraticBase + AddAssign + WithUnsigned<Unsigned = U>,
782                U: Integer + Clone + NumRef + CheckedAdd + CheckedMul + WithSigned<Signed = T>,
783            > $imp<ContinuedFraction<U>> for ContinuedFraction<U>
784        where
785            for<'r> &'r T: RefNum<T>,
786            for<'r> &'r U: RefNum<U>,
787        {
788            type Output = Self;
789
790            fn $method(self, rhs: ContinuedFraction<U>) -> Self {
791                if rhs.is_rational() {
792                    self.$method(rhs.to_rational().value())
793                } else {
794                    Self::try_from(
795                        QuadraticSurd::<T>::from(self).$method(QuadraticSurd::<T>::from(rhs)),
796                    )
797                    .unwrap()
798                }
799            }
800        }
801    };
802}
803
804impl_binop_for_ratio_surd!(impl Add, add);
805impl_binop_for_ratio_surd!(impl Sub, sub);
806impl_binop_for_ratio_surd!(impl Mul, mul);
807impl_binop_for_ratio_surd!(impl Div, div);
808
809// XXX: implement sqrt for InfiniteContinuedFraction
810// REF: https://crypto.stanford.edu/pbc/notes/contfrac/algebraic.html
811//      https://github.com/blynn/frac/blob/master/newton.c#L74
812
813/// This trait provide conversion from iterator of ASCII chars to
814/// continued fraction. This can be used for accepting high-precision
815/// decimal, or infinite continued fraction representation
816trait ParseContinuedFraction {
817    // TODO: implement followings and make the trait public
818    // fn parse_as_decimals() -> InfiniteContinuedFraction;
819    // fn parse_as_cfrac() -> InfiniteContinuedFraction;
820}
821
822impl<I: Iterator<Item = u8>> ParseContinuedFraction for I {}
823
824#[cfg(test)]
825mod tests {
826    use super::*;
827
828    #[test]
829    fn creation_test() {
830        let cf = ContinuedFraction::new(vec![0, 1, 2], vec![3, 4], false);
831        assert_eq!(
832            cf,
833            ContinuedFraction::new(vec![0, 1, 0, 0, 2], vec![3, 4], false)
834        );
835        assert_eq!(
836            cf,
837            ContinuedFraction::new(vec![0, 1, 1, 0, 1], vec![3, 4], false)
838        );
839
840        let cf = ContinuedFraction::new(vec![2, 3, 4], vec![], true);
841        assert_eq!(
842            cf,
843            ContinuedFraction::new(vec![2, 1, 0, 2, 4], vec![], true)
844        );
845        assert_eq!(
846            cf,
847            ContinuedFraction::new(vec![3, -1, -2, -4], vec![], true)
848        );
849        assert_eq!(cf, ContinuedFraction::new(vec![-3, 1, 2, 4], vec![], false));
850    }
851
852    #[test]
853    fn iter_test() {
854        let one = ContinuedFraction::<u32>::new(vec![1], vec![], false);
855        assert_eq!(one.coeffs().cloned().collect::<Vec<_>>(), vec![1]);
856        assert_eq!(one.convergents().collect::<Vec<_>>(), vec![Ratio::from(1)]);
857
858        let n_one = ContinuedFraction::<u32>::new(vec![1], vec![], true);
859        assert_eq!(n_one.coeffs().cloned().collect::<Vec<_>>(), vec![1]);
860        assert_eq!(
861            n_one.convergents().collect::<Vec<_>>(),
862            vec![Ratio::from(-1)]
863        );
864
865        let sq2 = ContinuedFraction::<u32>::new(vec![1], vec![2], false);
866        assert_eq!(
867            sq2.coeffs().take(5).cloned().collect::<Vec<_>>(),
868            vec![1, 2, 2, 2, 2]
869        );
870        assert_eq!(
871            sq2.convergents().take(5).collect::<Vec<_>>(),
872            vec![
873                Ratio::from(1),
874                Ratio::new(3, 2),
875                Ratio::new(7, 5),
876                Ratio::new(17, 12),
877                Ratio::new(41, 29)
878            ]
879        );
880
881        let n_sq2 = ContinuedFraction::<u32>::new(vec![1], vec![2], true);
882        assert_eq!(
883            n_sq2.coeffs().take(5).cloned().collect::<Vec<_>>(),
884            vec![1, 2, 2, 2, 2]
885        );
886        assert_eq!(
887            n_sq2.convergents().take(5).collect::<Vec<_>>(),
888            vec![
889                Ratio::from(-1),
890                Ratio::new(-3, 2),
891                Ratio::new(-7, 5),
892                Ratio::new(-17, 12),
893                Ratio::new(-41, 29)
894            ]
895        );
896    }
897
898    #[test]
899    fn conversion_test() {
900        assert_eq!(
901            ContinuedFraction::from(Ratio::from(3)),
902            ContinuedFraction::new(vec![3u32], vec![], false)
903        );
904        assert_eq!(
905            ContinuedFraction::from(Ratio::new(22, 7)),
906            ContinuedFraction::new(vec![3u32, 7], vec![], false)
907        );
908        assert_eq!(
909            ContinuedFraction::from(Ratio::new(-22, 7)),
910            ContinuedFraction::new(vec![3i32, 7], vec![], true)
911        );
912        assert_eq!(
913            ContinuedFraction::from(Ratio::new(7, 22)),
914            ContinuedFraction::new(vec![0u32, 3, 7], vec![], false)
915        );
916        assert_eq!(
917            ContinuedFraction::from(Ratio::new(-7, 22)),
918            ContinuedFraction::new(vec![0i32, 3, 7], vec![], true)
919        );
920        assert_eq!(
921            ContinuedFraction::from(Ratio::new(355, 113)),
922            ContinuedFraction::new(vec![3u32, 7, 16], vec![], false)
923        );
924    }
925
926    #[test]
927    fn fmt_test() {
928        assert_eq!(
929            format!("{}", ContinuedFraction::new(vec![1], vec![], false)),
930            "[1]"
931        );
932        assert_eq!(
933            format!("{}", ContinuedFraction::new(vec![1, 2, 3], vec![], false)),
934            "[1; 2, 3]"
935        );
936        assert_eq!(
937            format!("{}", ContinuedFraction::new(vec![1], vec![2], false)),
938            "[1; (2)]"
939        );
940        assert_eq!(
941            format!(
942                "{}",
943                ContinuedFraction::new(vec![1, 2, 3], vec![3, 2], false)
944            ),
945            "[1; 2, 3, (3, 2)]"
946        );
947    }
948
949    #[test]
950    fn arithmetic_test() {
951        // TODO: add more tests to cover all cases
952
953        let one = ContinuedFraction::one();
954        let n_one = ContinuedFraction::<u32>::new(vec![1], vec![], true);
955        let sq2 = ContinuedFraction::new(vec![1], vec![2], false);
956        let n_sq2 = ContinuedFraction::new(vec![1], vec![2], true);
957
958        // Add
959        assert_eq!(one.clone() + 1, ContinuedFraction::from(2));
960        assert_eq!(one.clone() + (-1), ContinuedFraction::zero());
961        assert_eq!(n_one.clone() + 1, ContinuedFraction::zero());
962        assert_eq!(n_one.clone() + (-1), ContinuedFraction::from(-2));
963
964        assert_eq!(
965            sq2.clone() + 1,
966            ContinuedFraction::new(vec![2u32], vec![2], false)
967        );
968        assert_eq!(
969            sq2.clone() + (-1),
970            ContinuedFraction::new(vec![0u32], vec![2], false)
971        );
972        assert_eq!(
973            n_sq2.clone() + 1,
974            ContinuedFraction::new(vec![0i32], vec![2], true)
975        );
976        assert_eq!(
977            n_sq2.clone() + (-1),
978            ContinuedFraction::new(vec![2i32], vec![2], true)
979        );
980
981        // Mul
982        assert_eq!(one.clone() * 2, ContinuedFraction::from(2));
983        assert_eq!(one.clone() * -2, ContinuedFraction::from(-2));
984        assert_eq!(n_one.clone() * 2, ContinuedFraction::from(-2));
985        assert_eq!(n_one.clone() * -2, ContinuedFraction::from(2));
986
987        assert_eq!(
988            sq2.clone() * 2,
989            ContinuedFraction::new(vec![2u32], vec![1, 4], false)
990        );
991        assert_eq!(
992            sq2.clone() * -2,
993            ContinuedFraction::new(vec![2i32], vec![1, 4], true)
994        );
995        assert_eq!(
996            n_sq2.clone() * 2,
997            ContinuedFraction::new(vec![2i32], vec![1, 4], true)
998        );
999        assert_eq!(
1000            n_sq2.clone() * -2,
1001            ContinuedFraction::new(vec![2u32], vec![1, 4], false)
1002        );
1003    }
1004}
num_irrational/cont_frac/simple.rs

num_irrational/cont_frac/
simple.rs
