Continued fraction/Arithmetic/G(matrix ng, continued fraction n)
This task investigates mathmatical operations that can be performed on a single continued fraction. This requires only a baby version of NG:
I may perform perform the following operations:
- Input the next term of N1
- Output a term of the continued fraction resulting from the operation.
I output a term if the integer parts of and are equal. Otherwise I input a term from N. If I need a term from N but N has no more terms I inject .
When I input a term t my internal state: is transposed thus
When I output a term t my internal state: is transposed thus
When I need a term t but there are no more my internal state: is transposed thus
I am done when b1 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 (e.g., from Continued fraction) calculate . 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 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
First I generate 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
First I generate 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
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";
}
- 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
Python
Python: NG
<lang python>class NG:
def __init__(self, a1, a, b1, b): self.a1, self.a, self.b1, self.b = a1, a, b1, b
def ingress(self, n): self.a, self.a1 = self.a1, self.a + self.a1 * n self.b, self.b1 = self.b1, self.b + self.b1 * n
@property def needterm(self): return (self.b == 0 or self.b1 == 0) or not self.a//self.b == self.a1//self.b1
@property def egress(self): n = self.a // self.b self.a, self.b = self.b, self.a - self.b * n self.a1, self.b1 = self.b1, self.a1 - self.b1 * n return n
@property def egress_done(self): if self.needterm: self.a, self.b = self.a1, self.b1 return self.egress
@property def done(self): return self.b == 0 and self.b1 == 0</lang>
Python: Testing
Uses r2cf method from here. <lang python>data = [["[1;5,2] + 1/2", [2,1,0,2], [13,11]],
["[3;7] + 1/2", [2,1,0,2], [22, 7]],
["[3;7] divided by 4", [1,0,0,4], [22, 7]]]
for string, ng, r in data:
print( "%-20s->" % string, end= ) op = NG(*ng) for n in r2cf(*r): if not op.needterm: print( " %r" % op.egress, end= ) op.ingress(n) while True: print( " %r" % op.egress_done, end= ) if op.done: break print()</lang>
- Output:
[1;5,2] + 1/2 -> 1 1 2 7 [3;7] + 1/2 -> 3 1 1 1 4 [3;7] divided by 4 -> 0 1 3 1 2
Racket
Main part of the NG-baby matrices. They are implemented as mutable structs. <lang Racket>#lang racket/base
(struct ng (a1 a b1 b) #:transparent #:mutable)
(define (ng-ingress! v t)
(define a (ng-a v)) (define a1 (ng-a1 v)) (define b (ng-b v)) (define b1 (ng-b1 v)) (set-ng-a! v a1) (set-ng-a1! v (+ a (* a1 t))) (set-ng-b! v b1) (set-ng-b1! v (+ b (* b1 t))))
(define (ng-needterm? v)
(or (zero? (ng-b v))
(zero? (ng-b1 v))
(not (= (quotient (ng-a v) (ng-b v)) (quotient (ng-a1 v) (ng-b1 v))))))
(define (ng-egress! v)
(define t (quotient (ng-a v) (ng-b v))) (define a (ng-a v)) (define a1 (ng-a1 v)) (define b (ng-b v)) (define b1 (ng-b1 v)) (set-ng-a! v b) (set-ng-a1! v b1) (set-ng-b! v (- a (* b t))) (set-ng-b1! v (- a1 (* b1 t))) t)
(define (ng-infty! v)
(when (ng-needterm? v) (set-ng-a! v (ng-a1 v)) (set-ng-b! v (ng-b1 v))))
(define (ng-done? v)
(and (zero? (ng-b v)) (zero? (ng-b1 v))))</lang>
Auxiliary functions to create producers of well known continued fractions. The function rational->cf is copied from r2cf task. <lang Racket>(define ((rational->cf n d))
(and (not (zero? d))
(let-values ([(q r) (quotient/remainder n d)])
(set! n d)
(set! d r)
q)))
(define (sqrt2->cf)
(define first? #t)
(lambda ()
(if first?
(begin
(set! first? #f)
1)
2)))</lang>
The function combine-ng-cf->cf combines a ng-matrix and a cf- producer and creates a cf-producer. The cf-producers can represent infinitely long continued fractions. The function cf-showln shows the first coefficients of a continued fraction represented in a cf-producer. <lang Racket>(define (combine-ng-cf->cf ng cf)
(define empty-producer? #f)
(lambda ()
(let loop ()
(cond
[(not empty-producer?) (define t (cf))
(cond
[t (ng-ingress! ng t)
(if (ng-needterm? ng)
(loop)
(ng-egress! ng))]
[else (set! empty-producer? #t)
(loop)])]
[(ng-done? ng) #f]
[(ng-needterm? ng) (ng-infty! ng)
(loop)]
[else (ng-egress! ng)]))))
(define (cf-showln cf n)
(for ([i (in-range n)])
(define val (cf))
(when val
(printf " ~a" val)))
(when (cf)
(printf " ..."))
(printf "~n"))</lang>
Some test <lang Racket>(display "[1;5,2] + 1/2 ->") (cf-showln (combine-ng-cf->cf (ng 2 1 0 2) (rational->cf 13 11)) 20)
(display "[3;7] + 1/2 ->") (cf-showln (combine-ng-cf->cf (ng 2 1 0 2) (rational->cf 22 7)) 20)
(display "[3;7] / 4 ->") (cf-showln (combine-ng-cf->cf (ng 1 0 0 4) (rational->cf 22 7)) 20)
(display "sqrt(2)/2 ->") (cf-showln (combine-ng-cf->cf (ng 1 0 0 2) (sqrt2->cf)) 20)
(display "1/sqrt(2) ->") (cf-showln (combine-ng-cf->cf (ng 0 1 1 0) (sqrt2->cf)) 20)
(display "(1+sqrt(2))/2 ->") (cf-showln (combine-ng-cf->cf (ng 1 1 0 2) (sqrt2->cf)) 20)</lang>
Sample output:
[1;5,2] + 1/2 -> 1 1 2 7 [3;7] + 1/2 -> 3 1 1 1 4 [3;7] / 4 -> 0 1 3 1 2 sqrt(2)/2 -> 0 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ... 1/sqrt(2) -> 0 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ... (1+sqrt(2))/2 -> 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 ...
Ruby
NG
<lang ruby># I define a class to implement baby NG class NG
def initialize(a1, a, b1, b) @a1, @a, @b1, @b = a1, a, b1, b end def ingress(n) @a, @a1 = @a1, @a + @a1 * n @b, @b1 = @b1, @b + @b1 * n end def needterm? return true if @b == 0 or @b1 == 0 return true unless @a/@b == @a1/@b1 false end def egress n = @a / @b @a, @b = @b, @a - @b * n @a1, @b1 = @b1, @a1 - @b1 * n n end def egress_done @a, @b = @a1, @b1 if needterm? egress end def done? @b == 0 and @b1 == 0 end
end</lang>
Testing
Uses r2cf method from here. <lang ruby>data = [["[1;5,2] + 1/2", [2,1,0,2], [13,11]],
["[3;7] + 1/2", [2,1,0,2], [22, 7]],
["[3;7] divided by 4", [1,0,0,4], [22, 7]]]
data.each do |str, ng, r|
printf "%-20s->", str
op = NG.new(*ng)
r2cf(*r) do |n|
print " #{op.egress}" unless op.needterm?
op.ingress(n)
end
print " #{op.egress_done}" until op.done?
puts
end</lang>
- Output:
[1;5,2] + 1/2 -> 1 1 2 7 [3;7] + 1/2 -> 3 1 1 1 4 [3;7] divided by 4 -> 0 1 3 1 2
Tcl
This uses the Generator class, R2CF class and printcf procedure from the r2cf task.
<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