Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2)

This task performs the basic mathematical functions on 2 continued fractions. This requires the full version of matrix NG:

{\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}

I may perform perform the following operations:

Input the next term of continued fraction N₁

Input the next term of continued fraction N₂

Output a term of the continued fraction resulting from the operation.

I output a term if the integer parts of ${\frac {a}{b}}$ and ${\frac {a_{1}}{b_{1}}}$ and ${\frac {a_{2}}{b_{2}}}$ and ${\frac {a_{12}}{b_{12}}}$ are equal. Otherwise I input a term from continued fraction N₁ or continued fraction N₂. If I need a term from N but N has no more terms I inject $\infty$ .

When I input a term t from continued fraction N₁ I change my internal state:

{\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}

is transposed thus

{\begin{bmatrix}a_{2}+a_{12}*t&a+a_{1}*t&a_{12}&a_{1}\\b_{2}+b_{12}*t&b+b_{1}*t&b_{12}&b_{1}\end{bmatrix}}

When I need a term from exhausted continued fraction N₁ I change my internal state:

{\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}

is transposed thus

{\begin{bmatrix}a_{12}&a_{1}&a_{12}&a_{1}\\b_{12}&b_{1}&b_{12}&b_{1}\end{bmatrix}}

When I input a term t from continued fraction N₂ I change my internal state:

{\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}

is transposed thus

{\begin{bmatrix}a_{1}+a_{12}*t&a_{12}&a+a_{2}*t&a_{2}\\b_{1}+b_{12}*t&b_{12}&b+b_{2}*t&b_{2}\end{bmatrix}}

When I need a term from exhausted continued fraction N₂ I change my internal state:

{\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}

is transposed thus

{\begin{bmatrix}a_{12}&a_{12}&a_{2}&a_{2}\\b_{12}&b_{12}&b_{2}&b_{2}\end{bmatrix}}

When I output a term t I change my internal state:

{\begin{bmatrix}a_{12}&a_{1}&a_{2}&a\\b_{12}&b_{1}&b_{2}&b\end{bmatrix}}

is transposed thus

{\begin{bmatrix}b_{12}&b_{1}&b_{2}&b\\a_{12}-b_{12}*t&a_{1}-b_{1}*t&a_{2}-b_{2}*t&a-b*t\end{bmatrix}}

When I need to choose to input from N₁ or N₂ I act:

if b and b2 are zero I choose N₁

if b is zero I choose N₂

if b2 is zero I choose N₂

if abs(

{\frac {a_{1}}{b_{1}}}-{\frac {a}{b}})

is greater than abs(

{\frac {a_{2}}{b_{2}}}-{\frac {a}{b}})

I choose N₁

otherwise I choose N₂

When performing arithmetic operation on two potentially infinite continued fractions it is possible to generate a rational number. eg ${\sqrt {2}}$ * ${\sqrt {2}}$ should produce 2. This will require either that I determine that my internal state is approaching infinity, or limiting the number of terms I am willing to input without producing any output.

C++

<lang cpp>/* Implement matrix NG

  Nigel Galloway, February 12., 2013

/

class NG_8 : public matrixNG {

 private: int a12, a1, a2, a, b12, b1, b2, b, t;
          double ab, a1b1, a2b2, a12b12;
 const int chooseCFN(){return fabs(a1b1-ab) > fabs(a2b2-ab)? 0 : 1;}
 const bool needTerm() {
   if (b1==0 and b==0 and b2==0 and b12==0) return false;
   if (b==0){cfn = b2==0? 0:1; return true;} else ab = ((double)a)/b;
   if (b2==0){cfn = 1; return true;} else a2b2 = ((double)a2)/b2;
   if (b1==0){cfn = 0; return true;} else a1b1 = ((double)a1)/b1;
   if (b12==0){cfn = chooseCFN(); return true;} else a12b12 = ((double)a12)/b12;
   thisTerm = (int)ab;
   if (thisTerm==(int)a1b1 and thisTerm==(int)a2b2 and thisTerm==(int)a12b12){
     t=a; a=b; b=t-b*thisTerm; t=a1; a1=b1; b1=t-b1*thisTerm; t=a2; a2=b2; b2=t-b2*thisTerm; t=a12; a12=b12; b12=t-b12*thisTerm;
     haveTerm = true; return false;
   }
   cfn = chooseCFN();
   return true;
 }
 void consumeTerm(){if(cfn==0){a=a1; a2=a12; b=b1; b2=b12;} else{a=a2; a1=a12; b=b2; b1=b12;}}
 void consumeTerm(int n){
   if(cfn==0){t=a; a=a1; a1=t+a1*n; t=a2; a2=a12; a12=t+a12*n; t=b; b=b1; b1=t+b1*n; t=b2; b2=b12; b12=t+b12*n;}
   else{t=a; a=a2; a2=t+a2*n; t=a1; a1=a12; a12=t+a12*n; t=b; b=b2; b2=t+b2*n; t=b1; b1=b12; b12=t+b12*n;}
 }
 public:
 NG_8(int a12, int a1, int a2, int a, int b12, int b1, int b2, int b): a12(a12), a1(a1), a2(a2), a(a), b12(b12), b1(b1), b2(b2), b(b){

}};</lang>

Testing

[3;7] + [0;2] <lang cpp>int main() {

 NG_8 a(0,1,1,0,0,0,0,1);
 r2cf n2(22,7);
 r2cf n1(1,2);
 for(NG n(&a, &n1, &n2); n.moreTerms(); std::cout << n.nextTerm() << " ");
 std::cout << std::endl;

 NG_4 a3(2,1,0,2);
 r2cf n3(22,7);
 for(NG n(&a3, &n3); n.moreTerms(); std::cout << n.nextTerm() << " ");
 std::cout << std::endl;
 return 0;

}</lang>

Output:

3 1 1 1 4
3 1 1 1 4

[1:5,2] * [3;7] <lang cpp>int main() {

 NG_8 b(1,0,0,0,0,0,0,1);
 r2cf b1(13,11);
 r2cf b2(22,7);
 for(NG n(&b, &b1, &b2); n.moreTerms(); std::cout << n.nextTerm() << " ");
 std::cout << std::endl;
 for(NG n(&a, &b2, &b1); n.moreTerms(); std::cout << n.nextTerm() << " ");
 std::cout << std::endl;
 for(r2cf cf(286,77); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
 std::cout << std::endl;
 return 0;

}</lang>

Output:

3 1 2 2
3 1 2 2

[1:5,2] - [3;7] <lang cpp>int main() {

 NG_8 c(0,1,-1,0,0,0,0,1);
 r2cf c1(13,11);
 r2cf c2(22,7);
 for(NG n(&c, &c1, &c2); n.moreTerms(); std::cout << n.nextTerm() << " ");
 std::cout << std::endl;
 for(r2cf cf(-151,77); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
 std::cout << std::endl;
 return 0;

}</lang>

Output:

-1 -1 -24 -1 -2
-1 -1 -24 -1 -2

Divide [] by [3;7] <lang cpp>int main() {

 NG_8 d(0,1,0,0,0,0,1,0);
 r2cf d1(22*22,7*7);
 r2cf d2(22,7);
 for(NG n(&d, &d1, &d2); n.moreTerms(); std::cout << n.nextTerm() << " ");
 std::cout << std::endl;
 return 0;

}</lang>

Output:

3 7

([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2]) <lang cpp>int main() {

 r2cf a1(2,7);
 r2cf a2(13,11);
 NG_8 na(0,1,1,0,0,0,0,1);
 NG A(&na, &a1, &a2); //[0;3,2] + [1;5,2]
 r2cf b1(2,7);
 r2cf b2(13,11);
 NG_8 nb(0,1,-1,0,0,0,0,1);
 NG B(&nb, &b1, &b2); //[0;3,2] - [1;5,2]
 NG_8 nc(1,0,0,0,0,0,0,1); //A*B
 for(NG n(&nc, &A, &B); n.moreTerms(); std::cout << n.nextTerm() << " ");
 std::cout << std::endl;
 for(r2cf cf(2,7); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
 std::cout << std::endl;
 for(r2cf cf(13,11); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
 std::cout << std::endl;
 for(r2cf cf(-7797,5929); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
 std::cout << std::endl;
 return 0;

}</lang>

Perl 6

Works with: Rakudo version 2013.05

<lang perl6>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 = 100) {
       my @cfs = [@cf1], [@cf2];
       gather {
           while @cfs[0].elems or @cfs[1].elems {
               take self!extract unless self!needterm;
               my $from = self!from;
               $from = @cfs[$from].elems ?? $from !! $from +^ 1; 
               self!inject($from, @cfs[$from].shift);
           }
           take self!drain while $!b;
       }[ ^ $limit ];
   }

   # 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 ];</lang>

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] = [-1;-1,-24,-1,-2] => -151/77

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

The NG2 object can work with infinitely long continued fractions, it does lazy evaluation. By default it is limited to returning the first 100 terms. Pass in a limit value if you want something other than default. For example, lets do the square of √2. √2 as a cf is [1;2,2,2,2,2,...] repeating indefinitely. First we'll construct a lazy infinite continued fraction, then multiply it by itself and limit the result to 50 terms. We'll then convert that continued fraction back to an arbitrary precision 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.) We'll the convert that number to a normal precision Rat, which is accurate to the nearest 1 / 2^64. (Note that the resulting continued fraction has been manually wrapped to fit on this page.)

<lang perl6>say "First 50 terms of √2 expressed as a continued fraction: "; my @root2 = 1, 2 xx *; my @result = NG2.new.operator(|%ops{'*'}).apply( @root2, @root2, limit => 50 ); say @root2.&ppcf, "² = \n"; say @result.&ppcf; say "\nConverted back to an arbitrary precision Rational: ", @result.&cf2r.&pprat; say "\nCoerced to a standard precision Rational: ", @result.&cf2r.Num.&pprat;</lang>

First 50 terms of √2 expressed as a continued fraction: 
[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,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,...]² =
 
[1;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-58451683124983302025,0,-1927184886226364356176,0,
-65467555105469489418600,0,-2223969688699736275876224,0,-75549501860685563890373016,0,-2566459093574609435996806320,
0,-87184059679676035260001041864,0,-2961691570015410589404038617056,0,-100610329320844284004477311938040,0,
-3417789505338690245562824567276304,0,-116104232852194624065131557975456296,0,-3944126127469278527968910146598237760]

Converted back to an arbitrary precision Rational: 8127503623352124799931806454985399361/4063751811676062399965903227492699681

Coerced to a standard precision Rational: 2

Tcl

This uses the Generator class, R2CF class and printcf procedure from the r2cf task.

Works with: Tcl version 8.6

<lang tcl>oo::class create NG2 {

   variable a b a1 b1 a2 b2 a12 b12 cf1 cf2
   superclass Generator
   constructor {args} {

lassign $args a12 a1 a2 a b12 b1 b2 b next

   }
   method operands {N1 N2} {

set cf1 $N1 set cf2 $N2 return [self]

   method Ingress1 t {

lassign [list [expr {$a2+$a12*$t}] [expr {$a+$a1*$t}] $a12 $a1 \ [expr {$b2+$b12*$t}] [expr {$b+$b1*$t}] $b12 $b1] \ a12 a1 a2 a b12 b1 b2 b

   }
   method Exhaust1 {} {

lassign [list $a12 $a1 $a12 $a1 $b12 $b1 $b12 $b1] \ a12 a1 a2 a b12 b1 b2 b

   }
   method Ingress2 t {

lassign [list [expr {$a1+$a12*$t}] $a12 [expr {$a+$a2*$t}] $a2 \ [expr {$b1+$b12*$t}] $b12 [expr {$b+$b2*$t}] $b2] \ a12 a1 a2 a b12 b1 b2 b

   }
   method Exhaust2 {} {

lassign [list $a12 $a12 $a2 $a2 $b12 $b12 $b2 $b2] \ a12 a1 a2 a b12 b1 b2 b

   }
   method Egress {} {

set t [expr {$a/$b}] lassign [list $b12 $b1 $b2 $b \ [expr {$a12 - $b12*$t}] [expr {$a1 - $b1*$t}] \ [expr {$a2 - $b2*$t}] [expr {$a - $b*$t}]] \ a12 a1 a2 a b12 b1 b2 b return $t

   method DoIngress1 {} {

try {tailcall my Ingress1 [$cf1]} on break {} {} oo::objdefine [self] forward DoIngress1 my Exhaust1 set cf1 "" tailcall my Exhaust1

   }
   method DoIngress2 {} {

try {tailcall my Ingress2 [$cf2]} on break {} {} oo::objdefine [self] forward DoIngress2 my Exhaust2 set cf2 "" tailcall my Exhaust2

   }
   method Ingress {} {

if {$b==0} { if {$b2 == 0} { tailcall my DoIngress1 } else { tailcall my DoIngress2 } } if {!$b2} { tailcall my DoIngress2 } if {!$b1} { tailcall my DoIngress1 } if {[my FirstSource?]} { tailcall my DoIngress1 } else { tailcall my DoIngress2 }

   method FirstSource? {} {

expr {abs($a1*$b*$b2 - $a*$b1*$b2) > abs($a2*$b*$b1 - $a*$b1*$b2)}

   }
   method NeedTerm? {} {

expr { ($b*$b1*$b2*$b12==0) || !($a/$b == $a1/$b1 && $a/$b == $a2/$b2 && $a/$b == $a12/$b12) }

   }
   method Done? {} {

expr {$b==0 && $b1==0 && $b2==0 && $b12==0}

   method Produce {} {

# Until we've drained both continued fractions... while {$cf1 ne "" || $cf2 ne ""} { if {[my NeedTerm?]} { my Ingress } else { yield [my Egress] } } # Drain our internal state while {![my Done?]} { yield [my Egress] }

}</lang> Demonstrating: <lang tcl>set op [[NG2 new 0 1 1 0 0 0 0 1] operands [R2CF new 1/2] [R2CF new 22/7]] printcf "\[3;7\] + \[0;2\]" $op

set op [[NG2 new 1 0 0 0 0 0 0 1] operands [R2CF new 13/11] [R2CF new 22/7]] printcf "\[1:5,2\] * \[3;7\]" $op

set op [[NG2 new 0 1 -1 0 0 0 0 1] operands [R2CF new 13/11] [R2CF new 22/7]] printcf "\[1:5,2\] - \[3;7\]" $op

set op [[NG2 new 0 1 0 0 0 0 1 0] operands [R2CF new 484/49] [R2CF new 22/7]] printcf "div test" $op

set op1 [[NG2 new 0 1 1 0 0 0 0 1] operands [R2CF new 2/7] [R2CF new 13/11]] set op2 [[NG2 new 0 1 -1 0 0 0 0 1] operands [R2CF new 2/7] [R2CF new 13/11]] set op3 [[NG2 new 1 0 0 0 0 0 0 1] operands $op1 $op2] printcf "layered test" $op3</lang>

Output:

[3;7] + [0;2]  -> 3,1,1,1,4
[1:5,2] * [3;7]-> 3,1,2,2
[1:5,2] - [3;7]-> -2,25,1,2
div test       -> 3,7
layered test   -> -2,1,2,5,1,2,1,26,3