Continued fraction/Arithmetic/G(matrix NG, Contined Fraction N1, Contined Fraction N2)

The NG2 object can work with infinitely long continued fractions, it does lazy evaluation. By default, it is limited to returning the first 30 terms. Pass in a limit value if you want something other than default.

class NG2 {
    has ( $!a12, $!a1, $!a2, $!a, $!b12, $!b1, $!b2, $!b );

    # Public methods
    method operator($!a12, $!a1, $!a2, $!a, $!b12, $!b1, $!b2, $!b ) { self }

    method apply(@cf1, @cf2, :$limit = 30) {
        my @cfs = [@cf1], [@cf2];
        gather {
            while @cfs[0] or @cfs[1] {
                my $term;
                (take $term if $term = self!extract) unless self!needterm;
                my $from = self!from;
                $from = @cfs[$from] ?? $from !! $from +^ 1;
                self!inject($from, @cfs[$from].shift);
            }
            take self!drain while $!b;
        }[ ^$limit ].grep: *.defined;
    }

    # Private methods
    method !inject ($n, $t) {
        multi sub xform(0, $t, $x12, $x1, $x2, $x) { $x2 + $x12 * $t, $x + $x1 * $t, $x12, $x1 }
        multi sub xform(1, $t, $x12, $x1, $x2, $x) { $x1 + $x12 * $t, $x12, $x + $x2 * $t, $x2 }
        ( $!a12, $!a1, $!a2, $!a ) = xform($n, $t, $!a12, $!a1, $!a2, $!a );
        ( $!b12, $!b1, $!b2, $!b ) = xform($n, $t, $!b12, $!b1, $!b2, $!b );
    }
    method !extract {
        my $t = $!a div $!b;
        ( $!a12, $!a1, $!a2, $!a, $!b12, $!b1, $!b2, $!b ) =
          $!b12, $!b1, $!b2, $!b,
                                  $!a12 - $!b12 * $t,
                                         $!a1 - $!b1 * $t,
                                               $!a2 - $!b2 * $t,
                                                     $!a - $!b * $t;
        $t;
    }
    method !from {
        return $!b == $!b2 == 0 ?? 0 !!
           $!b == 0 || $!b2 == 0 ?? 1 !!
           abs($!a1*$!b*$!b2 - $!a*$!b1*$!b2) > abs($!a2*$!b*$!b1 - $!a*$!b1*$!b2) ?? 0 !! 1;
    }
    method !needterm {
        so !([&&] $!b12, $!b1, $!b2, $!b) or $!a/$!b != $!a1/$!b1 != $!a2/$!b2 != $!a12/$!b1;
    }
    method !noterms($which) {
        $which ?? (($!a1, $!a, $!b1, $!b ) = $!a12, $!a2, $!b12, $!b2)
               !! (($!a2, $!a, $!b2, $!b ) = $!a12, $!a1, $!b12, $!b1);
    }
    method !drain {
    self!noterms(self!from) if self!needterm;
    self!extract;
    }
}

sub r2cf(Rat $x is copy) { # Rational to continued fraction
    gather loop {
    $x -= take $x.floor;
    last unless $x;
    $x = 1 / $x;
    }
}

sub cf2r(@a) { # continued fraction to Rational
    my $x = @a[* - 1].FatRat; # Use FatRats for arbitrary precision
    $x = @a[$_- 1] + 1 / $x for reverse 1 ..^ @a;
    $x
}

# format continued fraction for pretty printing
sub ppcf(@cf) { "[{ @cf.join(',').subst(',',';') }]" }

# format Rational for pretty printing. Use FatRats for arbitrary precision
sub pprat($a) { $a.FatRat.denominator == 1 ?? $a !! $a.FatRat.nude.join('/') }

my %ops = ( # convenience hash of NG matrix operators
    '+' => (0,1,1,0,0,0,0,1),
    '-' => (0,1,-1,0,0,0,0,1),
    '*' => (1,0,0,0,0,0,0,1),
    '/' => (0,1,0,0,0,0,1,0)
);

sub test_NG2 ($rat1, $op, $rat2) {
    my @cf1 = $rat1.&r2cf;
    my @cf2 = $rat2.&r2cf;
    my @result = NG2.new.operator(|%ops{$op}).apply( @cf1, @cf2 );
    say "{$rat1.&pprat} $op {$rat2.&pprat} => {@cf1.&ppcf} $op ",
        "{@cf2.&ppcf} = {@result.&ppcf} => {@result.&cf2r.&pprat}\n";
}

# Testing
test_NG2(|$_) for
   [   22/7, '+',  1/2 ],
   [  23/11, '*', 22/7 ],
   [  13/11, '-', 22/7 ],
   [ 484/49, '/', 22/7 ];


# Sometimes you may want to limit the terms in the continued fraction to something other than default.
# Here a lazy infinite continued fraction for  √2, then multiply it by itself. We'll limit the result
# to 6 terms for brevity’s' sake. We'll then convert that continued fraction back to an arbitrary precision
# FatRat Rational number. (Perl 6 stores FatRats internally as a ratio of two arbitrarily long integers.
# We need to exercise a little caution because they can eat up all of your memory if allowed to grow unchecked,
# hence the limit of 6 terms in continued fraction.) We'll then convert that number to a normal precision
# Rat, which is accurate to the nearest 1 / 2^64, 

say "√2 expressed as a continued fraction, then squared: ";
my @root2 = lazy flat 1, 2 xx *;
my @result = NG2.new.operator(|%ops{'*'}).apply( @root2, @root2, limit => 6 );
say @root2.&ppcf, "² = \n";
say @result.&ppcf;
say "\nConverted back to an arbitrary (ludicrous) precision Rational: ";
say @result.&cf2r.nude.join(" /\n");
say "\nCoerced to a standard precision Rational: ", @result.&cf2r.Num.Rat;

Output:

22/7 + 1/2 => [3;7] + [0;2] = [3;1,1,1,4] => 51/14

23/11 * 22/7 => [2;11] * [3;7] = [6;1,1,3] => 46/7

13/11 - 22/7 => [1;5,2] - [3;7] = [-2;25,1,2] => -151/77

484/49 / 22/7 => [9;1,7,6] / [3;7] = [3;7] => 22/7

√2 expressed as a continued fraction, then squared: 
[1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]² = 

[1;1,-58451683124983302025,-1927184886226364356176,-65467555105469489418600,-2223969688699736275876224]

Converted back to an arbitrary (ludicrous) precision Rational: 
32802382178012409621354320392819425499699206367450594986122623570838188983519955166754002 /
16401191089006204810536863200564985394427741343927508600629139291039556821665755787817601

Coerced to a standard precision Rational: 2