Continued fraction/Arithmetic/G(matrix ng, continued fraction n)

Continued fraction/Arithmetic/G(matrix ng, continued fraction n)
You are encouraged to solve this task according to the task description, using any language you may know.

This task investigates mathmatical operations that can be performed on a single continued fraction. This requires only a baby version of NG:

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

I may perform perform the following operations:

Input the next term of 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}}}$ are equal. Otherwise I input a term from N. If I need a term from N but N has no more terms I inject $\infty$.

When I input a term t my internal state: ${\begin{bmatrix}a_{1}&a\\b_{1}&b\end{bmatrix}}$ is transposed thus ${\begin{bmatrix}a+a_{1}*t&a_{1}\\b+b_{1}*t&b_{1}\end{bmatrix}}$

When I output a term t my internal state: ${\begin{bmatrix}a_{1}&a\\b_{1}&b\end{bmatrix}}$ is transposed thus ${\begin{bmatrix}b_{1}&b\\a_{1}-b_{1}*t&a-b*t\end{bmatrix}}$

When I need a term t but there are no more my internal state: ${\begin{bmatrix}a_{1}&a\\b_{1}&b\end{bmatrix}}$ is transposed thus ${\begin{bmatrix}a_{1}&a_{1}\\b_{1}&b_{1}\end{bmatrix}}$

I am done when b₁ and b are zero.

Demonstrate your solution by calculating:

[1;5,2] + 1/2

[3;7] + 1/2

[3;7] divided by 4

Using a generator for ${\sqrt {2}}$ (e.g., from Continued fraction) calculate ${\frac {1}{\sqrt {2}}}$. You are now at the starting line for using Continued Fractions to implement Arithmetic-geometric mean without ulps and epsilons.

The first step in implementing Arithmetic-geometric mean is to calculate ${\frac {1+{\frac {1}{\sqrt {2}}}}{2}}$ do this now to cross the starting line and begin the race.

C++

<lang cpp>/* Interface for all matrixNG classes

  Nigel Galloway, February 10th., 2013.

• /

class matrixNG {

 private:
 virtual void consumeTerm(){}
 virtual void consumeTerm(int n){}
 virtual const bool needTerm(){}
 protected: int cfn = 0, thisTerm;
            bool haveTerm = false;
 friend class NG;

}; /* Implement the babyNG matrix

  Nigel Galloway, February 10th., 2013.

• /

class NG_4 : public matrixNG {

 private: int a1, a, b1, b, t;
 const bool needTerm() {
   if (b1==0 and b==0) return false;
   if (b1==0 or b==0) return true; else thisTerm = a/b;
   if (thisTerm==(int)(a1/b1)){
     t=a; a=b; b=t-b*thisTerm; t=a1; a1=b1; b1=t-b1*thisTerm;
     haveTerm=true; return false;
   }
   return true;
 }
 void consumeTerm(){a=a1; b=b1;}
 void consumeTerm(int n){t=a; a=a1; a1=t+a1*n; t=b; b=b1; b1=t+b1*n;}
 public:
 NG_4(int a1, int a, int b1, int b): a1(a1), a(a), b1(b1), b(b){}

}; /* Implement a Continued Fraction which returns the result of an arithmetic operation on

  1 or more Continued Fractions (Currently 1 or 2).
  Nigel Galloway, February 10th., 2013.

• /

class NG : public ContinuedFraction {

 private:
  matrixNG* ng;
  ContinuedFraction* n[2];
 public:
 NG(NG_4* ng, ContinuedFraction* n1): ng(ng){n[0] = n1;}
 NG(NG_8* ng, ContinuedFraction* n1, ContinuedFraction* n2): ng(ng){n[0] = n1; n[1] = n2;}
 const int nextTerm() {ng->haveTerm = false; return ng->thisTerm;}
 const bool moreTerms(){
   while(ng->needTerm()) if(n[ng->cfn]->moreTerms()) ng->consumeTerm(n[ng->cfn]->nextTerm()); else ng->consumeTerm();
   return ng->haveTerm;
 }

};</lang>

Testing

[1;5,2] + 1/2

<lang cpp>int main() {

 NG_4 a1(2,1,0,2);
 r2cf n1(13,11);
 for(NG n(&a1, &n1); n.moreTerms(); std::cout << n.nextTerm() << " ");
 std::cout << std::endl;
 return 0;

}</lang>

Output:

1 1 2 7

[3;7] * 7/22

<lang cpp>int main() {

 NG_4 a2(7,0,0,22);
 r2cf n2(22,7);
 for(NG n(&a2, &n2); n.moreTerms(); std::cout << n.nextTerm() << " ");
 std::cout << std::endl;
 return 0;

}</lang>

Output:

1

[3;7] + 1/22

<lang cpp>int main() {

 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;7] divided by 4

<lang cpp>int main() {

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

}</lang>

Output:

0 1 3 1 2

${\frac {1}{\sqrt {2}}}$

First I generate ${\frac {1}{\sqrt {2}}}$ as a continued fraction, then I obtain an approximate value using r2cf for comparison. <lang cpp>int main() {

 NG_4 a5(0,1,1,0);
 SQRT2 n5;
 int i = 0;
 for(NG n(&a5, &n5); n.moreTerms() and i++ < 20; std::cout << n.nextTerm() << " ");
 std::cout << "..." << std::endl;
 for(r2cf cf(10000000, 14142136); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
 std::cout << std::endl;
 return 0;

}</lang>

Output:

0 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ...
0 1 2 2 2 2 2 2 2 2 2 6 1 2 4 1 1 2

${\frac {1+{\sqrt {2}}}{2}}$

First I generate ${\frac {1+{\sqrt {2}}}{2}}$ as a continued fraction, then I obtain an approximate value using r2cf for comparison. <lang cpp>int main() {

 int i = 0;
 NG_4 a6(1,1,0,2);
 SQRT2 n6;
 for(NG n(&a6, &n6); n.moreTerms() and i++ < 20; std::cout << n.nextTerm() << " ");
 std::cout << "..." << std::endl;
 for(r2cf cf(24142136, 20000000); cf.moreTerms(); std::cout << cf.nextTerm() << " ");
 std::cout << std::endl;
 return 0;

}</lang>

Output:

1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 ...
1 4 1 4 1 4 1 4 1 4 3 2 1 9 5

Perl 6

Works with: Rakudo version 2013.05

All the important stuff takes place in the NG object. Everything else is helper subs for testing and display. <lang perl6>class NG {

   has ( $!a1, $!a, $!b1, $!b );
   submethod BUILD ( :$!a1, :$!a, :$!b1, :$!b ) { }

   # Public methods
   method new( $a1, $a, $b1, $b ) { self.bless( *, :$a1, :$a, :$b1, :$b ) }
   method apply(@cf, :$limit = Inf) {
       (gather {
           map { take self!extract unless self!needterm; self!inject($_) }, @cf;
           take self!drain until self!done;
       })[ ^ $limit ]
   }

   # Private methods
   method !inject ($n) {
       sub xform($n, $x, $y) { $x, $n * $x + $y }
       ( $!a, $!a1 ) = xform( $n, $!a1, $!a );
       ( $!b, $!b1 ) = xform( $n, $!b1, $!b );
   }
   method !extract {
       sub xform($n, $x, $y) { $y, $x - $y * $n }
       my $n = $!a div $!b;
       ($!a,  $!b ) = xform( $n, $!a,  $!b  );
       ($!a1, $!b1) = xform( $n, $!a1, $!b1 );
       $n
   }
   method !drain { $!a = $!a1, $!b = $!b1 if self!needterm; self!extract }
   method !needterm { so [||] !$!b, !$!b1, $!a/$!b != $!a1/$!b1 }
   method !done { not [||] $!b, $!b1 }

}

sub r2cf(Rat $x is copy) { # Rational to continued fraction

   gather loop {

$x -= take $x.floor; last if !$x; $x = 1 / $x;

   }

}

sub cf2r(@a) { # continued fraction to Rational

   my $x = @a[* - 1]; # Use FatRats for arbitrary precision
   $x = ( @a[$_- 1] + 1 / $x ).FatRat for reverse 1 ..^ @a;
   $x

}

sub ppcf(@cf) { # format continued fraction for pretty printing

   "[{ @cf.join(',').subst(',',';') }]"

}

sub pprat($a) { # format Rational for pretty printing

  # Use FatRats for arbitrary precision
  $a.FatRat.denominator == 1 ?? $a !! $a.FatRat.nude.join('/')

}

sub test_NG ($rat, @ng, $op) {

   my @cf = $rat.Rat(1e-18).&r2cf;
   my @op = NG.new( |@ng ).apply( @cf );
   say $rat.perl, ' as a cf: ', @cf.&ppcf, " $op = ",
       @op.&ppcf, "\tor ", @op.&cf2r.&pprat, "\n";

}

1. Testing

test_NG(|$_) for (

   [ 13/11, [<2 1 0 2>], '+ 1/2 '    ],
   [ 22/7,  [<2 1 0 2>], '+ 1/2    ' ],
   [ 22/7,  [<1 0 0 4>], '/ 4      ' ],
   [ 22/7,  [<7 0 0 22>], '* 7/22   ' ],
   [ 2**.5, [<1 1 0 2>], "\n(1+√2)/2 (approximately)" ]

); </lang>

Output

<13/11> as a cf: [1;5,2] + 1/2  = [1;1,2,7]	or 37/22

<22/7> as a cf: [3;7] + 1/2     = [3;1,1,1,4]	or 51/14

<22/7> as a cf: [3;7] / 4       = [0;1,3,1,2]	or 11/14

<22/7> as a cf: [3;7] * 7/22    = [1]	or 1

1.4142135623731e0 as a cf: [1;2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2] 
(1+√2)/2 (approximately) = [1;4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4]	or 225058681/186444716

The cf for (1+√2)/2 in the testing routine is an approximation. The NG object is capable of working with infinitely long continued fractions, but displaying them can be problematic. You can pass in a limit to the apply method to get a fixed maximum number of terms though. Here are the first 100 terms from the infinite cf (1+√2)/2 and its Rational representation. <lang perl6>my @continued-fraction = NG.new( 1,1,0,2 ).apply( ( 1, 2 xx * ), limit => 100 ); say @continued-fraction.&ppcf.comb(/ . ** 1..80/).join("\n"); say @continued-fraction.&cf2r.&pprat;</lang>

[1;4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4
,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4
,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4]
161733217200188571081311986634082331709/133984184101103275326877813426364627544

Ruby

NG

<lang ruby>=begin

 I define a class to implement baby NG
 Nigel Galloway February 6th., 2013

=end class NG

 def initialize(a1, a, b1, b)
   @a1 = a1; @a = a; @b1 = b1; @b = b;
 end
 def ingress(n)
   t=@a; @a=@a1; @a1=t + @a1 * n; t=@b; @b=@b1; @b1=t + @b1 * n;
 end
 def needterm?
   return true if @b1 == 0 or @b == 0
   return true unless @a/@b == @a1/@b1
   return false
 end
 def egress
   n = @a/@b
   t=@a; @a=@b; @b=t - @b * n; t=@a1; @a1=@b1; @b1=t - @b1 * n;
   return n
 end
 def egress_done
   if needterm? then @a=@a1; @b=@b1 end
   return egress
 end
 def done?
   if @b1 == 0 and @b == 0 then return true else return false end
 end

end</lang>

Testing

[1;5,2] + 1/2

<lang ruby>op = NG.new(2,1,0,2) r2cf(13,11) {|n| if op.needterm? then op.ingress(n) else print "#{op.egress} "; op.ingress(n) end} while not op.done? do print "#{op.egress_done} " end</lang>

Output:

1 1 2 7

[3;7] + 1/2

<lang ruby>op = NG.new(2,1,0,2) r2cf(22,7) {|n| if op.needterm? then op.ingress(n) else print "#{op.egress} "; op.ingress(n) end} while not op.done? do print "#{op.egress_done} " end</lang>

Output:

3 1 1 1 4

[3;7] divided by 4

<lang ruby>op = NG.new(1,0,0,4) r2cf(22,7) {|n| if op.needterm? then op.ingress(n) else print "#{op.egress} "; op.ingress(n) end} while not op.done? do print "#{op.egress_done} " end</lang>

Output:

0 1 3 1 2

Tcl

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

Works with: Tcl version 8.6

Translation of: Ruby

<lang tcl># The single-operand version of the NG operator, using our little generator framework oo::class create NG1 {

   superclass Generator

   variable a1 a b1 b cf
   constructor args {

next lassign $args a1 a b1 b

   }
   method Ingress n {

lassign [list [expr {$a + $a1*$n}] $a1 [expr {$b + $b1*$n}] $b1] \ a1 a b1 b

   }
   method NeedTerm? {} {

expr {$b1 == 0 || $b == 0 || $a/$b != $a1/$b1}

   }
   method Egress {} {

set n [expr {$a/$b}] lassign [list $b1 $b [expr {$a1 - $b1*$n}] [expr {$a - $b*$n}]] \ a1 a b1 b return $n

   }
   method EgressDone {} {

if {[my NeedTerm?]} { set a $a1 set b $b1 } tailcall my Egress

   }
   method Done? {} {

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

   }

   method operand {N} {

set cf $N return [self]

   }
   method Produce {} {

while 1 { set n [$cf] if {![my NeedTerm?]} { yield [my Egress] } my Ingress $n } while {![my Done?]} { yield [my EgressDone] }

   }

}</lang> Demonstrating: <lang tcl># The square root of 2 as a continued fraction in the framework oo::class create Root2 {

   superclass Generator
   method apply {} {

yield 1 while {[self] ne ""} { yield 2 }

   }

}

set op [[NG1 new 2 1 0 2] operand [R2CF new 13/11]] printcf "\[1;5,2\] + 1/2" $op

set op [[NG1 new 7 0 0 22] operand [R2CF new 22/7]] printcf "\[3;7\] * 7/22" $op

set op [[NG1 new 2 1 0 2] operand [R2CF new 22/7]] printcf "\[3;7\] + 1/2" $op

set op [[NG1 new 1 0 0 4] operand [R2CF new 22/7]] printcf "\[3;7\] / 4" $op

set op [[NG1 new 0 1 1 0] operand [Root2 new]] printcf "1/\u221a2" $op 20

set op [[NG1 new 1 1 0 2] operand [Root2 new]] printcf "(1+\u221a2)/2" $op 20 printcf "approx val" [R2CF new 24142136 20000000]</lang>

Output:

[1;5,2] + 1/2  -> 1,1,2,7
[3;7] * 7/22   -> 1
[3;7] + 1/2    -> 3,1,1,1,4
[3;7] / 4      -> 0,1,3,1,2
1/√2           -> 0,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,…
(1+√2)/2       -> 1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,…
approx val     -> 1,4,1,4,1,4,1,4,1,4,3,2,1,9,5