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

From Rosetta Code
Continued fraction/Arithmetic/G(matrix NG, Contined Fraction N1, Contined Fraction N2) is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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

I may perform perform the following operations:

Input the next term of continued fraction N1
Input the next term of continued fraction N2
Output a term of the continued fraction resulting from the operation.

I output a term if the integer parts of and and and are equal. Otherwise I input a term from continued fraction N1 or continued fraction N2. If I need a term from N but N has no more terms I inject .

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

is transposed thus

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

is transposed thus

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

is transposed thus

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

is transposed thus

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

is transposed thus

When I need to choose to input from N1 or N2 I act:

if b and b2 are zero I choose N1
if b is zero I choose N2
if b2 is zero I choose N2
if abs( is greater than abs( I choose N1
otherwise I choose N2

When performing arithmetic operation on two potentially infinite continued fractions it is possible to generate a rational number. eg * 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++[edit]

/* 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){
}};

Testing[edit]

[3;7] + [0;2]

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;
}
Output:
3 1 1 1 4
3 1 1 1 4

[1:5,2] * [3;7]

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;
}
Output:
3 1 2 2
3 1 2 2

[1:5,2] - [3;7]

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;
}
Output:
-1 -1 -24 -1 -2
-1 -1 -24 -1 -2

Divide [] by [3;7]

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;
}
Output:
3 7

([0;3,2] + [1;5,2]) * ([0;3,2] - [1;5,2])

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;
}

Perl 6[edit]

Works with: Rakudo version 2016.01
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 ];

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 30 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 6 terms for brevitys' 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: ";
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;
√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,...]² = 

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

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

Coerced to a standard precision Rational: 2

Tcl[edit]

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

Works with: Tcl version 8.6
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]
}
}
}

Demonstrating:

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
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