DGtal  0.9.4beta
Irreducible fraction, continued fractions

Table of Contents

Author(s) of this documentation:
Jacques-Olivier Lachaud

Part of the Arithmetic package.

This part of the manual describes how to create and use irreducible fractions. More precisely, several representations of irreducible fractions are provided. They are based on the Stern-Brocot tree structure. With these fractions, amortized constant time operations are provided for computing reduced fractions.

Part of the code is a backport from ImaGene. [39]

Introduction to positive irreducible fraction.

A positive irreducible fraction is a fraction p / q, such that p and q are positive integers and gcd(p, q)=1. The fractions \(\frac{0}{1}\) (zero) and \(\frac{1}{0}\) (infinite) are added for convenience. Standard operations like '+', '-', '*', or '/' can be defined over these fractions: it is a matter of playing with numerators and denominators.

However, these fractions have other properties that are useful in the context of digital geometry, especially in digital straightness. These other properties are related to Euclid algorithm, the Stern-Brocot tree, and simple continued fractions.

Euclid algorithm

The well-known Euclid algorithm is useful to compute the greatest common divisor of two integers p and q. However, its intermediate computations are also useful. Let us play with 8/3 (see IntegerComputer::gcd).

p =u *q +r
8 =2 *3 +2
3 =1 *2 +1
2 =2 *1 +0
1 *0

At each iteration, p takes the former value of q, q takes the former value of r. The algorithm stops when q is null. Then p is the gcd (1 here).

The values of u are called the quotients of \(\frac{p}{q}\), and are numbered from 0. A unique way to describe \(\frac{p}{q}\) is the sequence \(u_0,u_1,\ldots,u_k\), which is in fact denoted as \([u_0;u_1,\ldots,u_k]\). This sequence characterizes the fraction.

Continued fraction

In fact, one can prove that \(\frac{p}{q}=u_0+\frac{1}{u_1+\frac{1}{\ldots+\frac{1}{u_k}}}\), which is the simple continued fraction of p / q.

The successive fractions \(\frac{p_i}{q_i}=[u_0;u_1,\ldots,u_i]\), for \(k < i\) are called the principal convergents of \(\frac{p}{q}\), and are the best approximations of this fraction for the size of their denominator.

A i-reduced of a continued fraction is its k - i-th convergent. The depth of the fraction is k.

A nice property of convergents with respect to their fraction is that even convergents are smaller than the fraction while odd convergents are greater.

Stern-Brocot tree

Another (but one-to-one) way of seeing irreducible fractions is the Stern-Brocot tree. Place the fractions 0/1 and 1/0 at the top, then for each pair of consecutive fractions \(\frac{p}{q}\) and \(\frac{p'}{q'}\) compute its mediant \(\frac{p+p'}{q+q'}\). The fraction 1/1 appears first, then 1/2 and 2/1, then 1/3, 2/3, 3/2, 3/1, etc.

Sternbrocot.png
Stern-Brocot tree of irreducible fractions.

Now the sequence of quotients \(u_0, u_1,\ldots, u_k \) is exactly the sequence of right-then-left moves in the Stern-Brocot tree from node 1/1 when \( u_0 \ge 1 \), except for the last \( u_k \) where there is one less movement. Otherwise, when \( u_0 = 0 \), the sequence \( u_1,\ldots, u_k \) is exactly the sequence of left-then-right moves in the Stern-Brocot tree from node 1/1, except for the last \( u_k \) where there is one less movement.

Looking back at \( 8/3=[2;1,2] \) gives 2 right moves, 1 left move, then 2-1 right move.

Implementation of irreducible fractions.

When playing with irreducible fractions, computing reduced fractions or partial quotients is one of the main task. For instance, splitting formula and Berstel splitting formula is intensively used by algorithms related to digital straightness. Therefore, we must have representations of irreducible fractions which are efficient for that kind of requests.

Several representations have been implemented in the Arithmetic package. Their interface are common (see Irreducible fraction concept) so that a user can choose which representation suits him best.

Irreducible fraction concept

Irreducible fractions are first described abstractly with the concept CPositiveIrreducibleFraction. Note that the following operations are defined for each irreducible fraction.

Name Expression Type requirements Return type Precondition Semantics Post condition Complexity
Constructor Fraction( p, q )X creates the fraction p'/q', where p'=p/g, q'=q/g, g=gcd(p,q) o(p+q)
numerator x.p() Integer ! x.null() returns the numeratorO(1)
denominator x.q() Integer ! x.null() returns the denominatorO(1)
quotient x.u() Size ! x.null() returns the quotient \(u_k\) O(1)
depth x.k() Size ! x.null() returns the depth k O(1)
null test x.null() bool returns 'true' if the fraction is null 0/0 (default fraction) O(1)
even parity x.even() bool ! x.null() returns 'true' iff the fraction is even, i.e. k is even O(1)
odd parity x.odd() bool ! x.null() returns 'true' iff the fraction is odd, i.e. k is odd O(1)
left descendantx.left() X ! x.null() returns the left descendant of p/q in the Stern-Brocot tree O(1)
right descendantx.right() X ! x.null() returns the right descendant of p/q in the Stern-Brocot tree O(1)
father x.father() X ! x.null() returns the father of this fraction, ie \([u_0,...,u_k - 1]\) O(1)
m-father x.father(m)X ! x.null(), m>=0 returns the m-father of this fraction, ie \([u_0,...,u_{k-1}, m]\) O( m)
previousPartialx.previousPartial() X ! x.null() returns the previous partial of this fraction, ie \([u_0,...,u_{k-1}]\) O(1)
inverse x.inverse() X ! x.null() returns the inverse of this fraction, ie \([0,u_0,...,u_k]\) if \(u_0 \neq 0 \) or \([u_1,...,u_k]\) otherwise O(1)
m-th partial x.partial(m) X ! x.null() returns the m-th partial of this fraction, ie \([u_0,...,u_m]\) O(1)
m-th reduced x.reduced(m) X ! x.null() returns the m-th reduced of this fraction, equivalently the \(k-m\) partial, ie \([u_0,...,u_{k-m}]\) O(1)
splitting formula x.getSplit(x1, x2)void ! x.null() modifies fractions x1 and x2 such that \( x1 \oplus x2 = x \)O(1)
Berstel splitting formula x.getSplitBerstel(x1, n1, x2, n2)void ! x.null() modifies fractions x1 and x2 and integers n1 and n2 such that \( (x1)^{n1} \oplus (x2)^{n2} = x \)O(1)
Continued fraction coefficients x.getCFrac(quots)void modifies the vector quots such that it contains the quotients \(u_0,u_1,...,u_k \)O(k)
equality x.equals(p, q)bool returns 'true' iff the fraction is equal to \( p / q \). O(1)
less than x.lessThan(p, q)bool returns 'true' iff the fraction is inferior to \( p / q \). O(1)
more than x.moreThan(p, q)bool returns 'true' iff the fraction is superior to \( p / q \). O(1)
equality == x == y bool returns 'true' iff the fraction is equal to y. O(1)
inequality != x != y bool returns 'true' iff the fraction is different from y. O(1)
less than < x < y bool returns 'true' iff the fraction is inferior to y. O(1)
more than > x > y bool returns 'true' iff the fraction is superior to y. O(1)
Next continued fraction x.pushBack( pair ) transforms this fraction \([0,u_0,...,u_k]\) into \([0,u_0,...,u_k,m]\), where pair is \((m,k+1)\) O(m)
Next continued fraction x.push_back( pair ) transforms this fraction \([0,u_0,...,u_k]\) into \([0,u_0,...,u_k,m]\), where pair is \((m,k+1)\) O(m)
Begin visiting quotients x.begin() ConstIterator returns a forward iterator on the beginning of the sequence of quotients \([u_0,...,u_k]\)
End visiting quotients x.end() ConstIterator returns a forward iterator after the end of the sequence of quotients \([u_0,...,u_k]\)

Inner types are:

Notations are:

Naive approach

A naive approach is to represent an irreducible fraction with three arrays of integers: the numerators (p_i), the denominators (q_i), the quotients (u_i).

However, each time we duplicate the fraction, or each time we compute a reduced fraction, the complexity of the operation is proportional to the depth of the fraction, i.e. in O(k).

We wish to have amortized O(1) complexity for these operations. Therefore, we take another path.

First approach with Stern-Brocot tree

This approach is implemented in the class SternBrocot::Fraction (see also class SternBrocot). We construct on demand parts of the Stern-Brocot tree structure. For instance, if one wish to manipule fraction 8/3, the following nodes of the tree are computed once and for all: 0/1, 1/0, 1/1, 2/1, 1/2, 3/1, 1/3, 5/2, 2/5, 8/3, 3/8.

Each node in the Stern-Brocot is created once. The node stores information on the irreducible fraction itself (p/q, the partial quotient u, the depth k), but also pointers to ascendants, descendants and inverse in the Stern-Brocot tree. Nodes are constructed on demand, when the user ask for descendant or for a specific fraction.

A fraction SternBrocot::Fraction is just a pointer to the corresponding node in the Stern-Brocot tree. Looking for a reduced partial is just moving from node to node.

#include "DGtal/arithmetic/SternBrocot.h"
...
typedef SternBrocot<DGtal::int64_t,DGtal::int32_t> SB;
typedef SB::Fraction Fraction;
typedef Fraction::ConstIterator ConstIterator;
Fraction f1( 117, 45 ); // 117/45
Fraction f2; // null 0/0
f2.push_back( std::make_pair( 2, 0 ) ); // u_0 = 2
f2.push_back( std::make_pair( 1, 1 ) ); // u_1 = 1
f2.push_back( std::make_pair( 2, 2 ) ); // u_2 = 2
// f2 = 8/3
// Displays quotients of f1.
for ( ConstIterator it = f1.begin(), itend = f1.end();
it != itend; ++it )
std::cout << "u_" << (*it).second
<< " = " << (*it).first << std::endl;

You may have a look at testSternBrocot.cpp to see how to use these fractions.

Also inverses and reduced fractions are quickly computed, this representation has two disadvantages:

Second approach with Stern-Brocot tree

This approach is implemented in the class LighterSternBrocot::Fraction (see also class LighterSternBrocot).

There are two main differences with the class SternBrocot. The first one is that inverses are not stored. With this optimization, there are twice less nodes and each node is lighter. The second one lies in the access to the children of a node. Here, a map type M is provided so that a node \( [u_0; u_1, ..., u_n] \) can access its child node \( [u_0; u_1, ..., u_n, k] \) in the time of the operation M::operator[]. This representation is also different from LightSternBrocot in the sense that nodes have only one set of child nodes and that only fractions greater than 1/1 are stored.

In this representation, the fraction 1/1 has depth 0, like 2/1, 3/1, etc. Furthermore, each fraction \( [u_0,...,u_n] \) has an origin which is the fraction \( [u_0,...,u_{n-1},1] \). It is the top extremity of this branch. The origin has depth n-1 since \( [u_0,...,u_{n-1},1] \) = \( [u_0,...,u_{n-1}+1] \). Inversely a k-child of \( [u_0,...,u_n] \), for k >= 2, is the fraction \( [u_0,...,u_n - 1, k] \). A 1-child of a fraction f is itself, except for the fraction 1/0 where its 1-child is 1/1 by convention.

In practice, also this class has supposedly a better complexity than SternBrocot, it is 1% slower for integers smaller than 10^9 and 5% slower for integers smaller than 10^4. Note however that it takes like 7 times less memory (and asymptotically less when the number of computations tends toward infinity).

This class is not to be instantiated, since it is useless to duplicate it. Use static method LighterSternBrocot::fraction or constructor of LighterSternBrocot::Fraction to obtain your fractions.

#include "DGtal/arithmetic/LighterSternBrocot.h"
...
typedef SternBrocot<DGtal::int64_t, DGtal::int32_t, DGtal::StdMapRebinder> SB;
typedef SB::Fraction Fraction;
typedef Fraction::ConstIterator ConstIterator;
Fraction f1( 117, 45 ); // 117/45
Fraction f2; // null 0/0
f2.push_back( std::make_pair( 2, 0 ) ); // u_0 = 2
f2.push_back( std::make_pair( 1, 1 ) ); // u_1 = 1
f2.push_back( std::make_pair( 2, 2 ) ); // u_2 = 2
// f2 = 8/3
// Displays quotients of f1.
for ( ConstIterator it = f1.begin(), itend = f1.end();
it != itend; ++it )
std::cout << "u_" << (*it).second
<< " = " << (*it).first << std::endl;

You may have a look at testLighterSternBrocot.cpp to see how to use these fractions.

Using irreducible fractions.

Choosing your integers for fractions

Available fractions allows you to choose the ultimate precision you wish to have for your fractions.

You have two different integral types for irreducible fractions: TInteger and TSize. The type TInteger defines the integral type for the numerators and the denominators. The type TSize defines the integral type for the quotients and for the depth of the fraction.

Note
It is worthy to note that quotients are in general much much smaller than numerators and denominators (except for fractions p/1 or 1/q). You may choose thus a smaller type for TSize (like DGtal::int32_t) than for TInteger (for instance DGtal::int64_t).

You may for instance choose your fractions like:

// With SternBrocot
typedef SternBrocot<DGtal::int64_t,DGtal::int32_t>::Fraction Fraction; // not so big fractions
typedef SternBrocot<DGtal::BigInteger,DGtal::BigInteger>::Fraction Fraction; // arbitrary large fractions
// With LighterSternBrocot
typedef LighterSternBrocot<DGtal::int64_t,DGtal::int32_t,DGtal::StdMapRebinder>::Fraction Fraction; // not so big fractions
typedef LighterSternBrocot<DGtal::BigInteger,DGtal::BigInteger,DGtal::StdMapRebinder>::Fraction Fraction; // arbitrary large fractions
Note
If you wish to instantiate fractions with integral types other than DGtal::int32_t, DGtal::int64_t, DGtal::BigInteger, then you must instantiate the singleton of the class. For instance, add somewhere in your code:
// With SternBrocot
typedef SternBrocot<IntegralType1,IntegralType2>::Fraction Fraction; // not so big fractions
template <>
SternBrocot<IntegralType1,IntegralType2>*
SternBrocot<IntegralType1,IntegralType2>::singleton = 0;
Note
In some sense, IntegralType2 should be promotable to IntegralType1. I.e., if t1 is of type IntegralType1 and t2 is of type IntegralType2 then
t1 = (IntegralType1) t2;

should be valid.

Instantiating fractions

You may instantiate a fraction directly by giving the numerator p and denominator q. Note that if \( g = gcd(p,q) \neq 1 \), the instantiated fraction is \( \frac{p'}{q'} \) where \( p'=p/g, q'=q/g \). We follow the example convergents.cpp.

First we give the correct includes:

#include <iostream>
#include "DGtal/arithmetic/LighterSternBrocot.h"

Then the useful types:

typedef LighterSternBrocot<DGtal::int64_t, DGtal::int64_t, StdMapRebinder> SB; // the type of the Stern-Brocot tree
typedef SB::Fraction Fraction; // the type for fractions
typedef Fraction::ConstIterator ConstIterator; // the iterator type for visiting quotients
typedef Fraction::Value Value; // the value of the iterator, a pair (quotient,depth).

The program is given p and q in arguments argv. The fraction is instantiated as follows:

DGtal::int64_t p = atoll( argv[ 1 ] );
DGtal::int64_t q = atoll( argv[ 2 ] );
Fraction f( p, q ); // fraction p/q

Obtaining quotients

You may visit the quotients of a fraction by using the single-pass iterators provided by the class. A fraction is indeed a model of CSinglePassConstRange. The quotients are computed and stored with the call to begin(). We therefore advise you to store the begin() iterator in a variable. Then, all other operations with iterators (copy, displacement, comparison, dereferenciation) are O(1). We display the quotients of the fraction f as follows.

// Visit quotients u_k as pair (u_k,k)
std::cout << "z = ";
ConstIterator itbegin = f.begin(), itend = f.end();
for ( ConstIterator it = itbegin; it != itend; ++it )
{
Value u = *it;
std::cout << ( ( it == itbegin ) ? "[" : "," )
<< u.first;
}
std::cout << "]" << std::endl;
Note
The iterator value is the pair \( (u_i,i) \) where \( 0 \le i \le k \) and k is the depth of the fraction. The type is therefore a std::pair<Size,Size>. This is done for consistency with push_back method.

Obtaining all the convergents of a fraction

One way of obtaining all the convergents of a fraction is simply to create the fraction progressively by adding quotients one at a time. You may use the methods pushBack or push_back to do so. This can be done as follows:

Fraction g; // fraction null, 0/0, invalid
for ( ConstIterator it = itbegin; it != itend; ++it )
{
Value u = *it;
std::cout << "z_" << u.second << " = ";
g.push_back( u ); // add (u_i,i) to existing fractions
std::cout << g.p() << " / " << g.q() << std::endl;
}

You may now execute this program to get the quotients and the convergents of a given fraction. We give below several examples:

# A Fibonacci fraction
$ ./examples/arithmetic/convergents 21 34
z = [0,1,1,1,1,1,1,2]
z_0 = 0 / 1
z_1 = 1 / 1
z_2 = 1 / 2
z_3 = 2 / 3
z_4 = 3 / 5
z_5 = 5 / 8
z_6 = 8 / 13
z_7 = 21 / 34

# Approximations of pi
$ ./examples/arithmetic/convergents 103993 33102
z = [3,7,15,1,292]
z_0 = 3 / 1
z_1 = 22 / 7
z_2 = 333 / 106
z_3 = 355 / 113
z_4 = 103993 / 33102

# A huge fraction
$ ./examples/arithmetic/convergents-biginteger 243224233245235253407096734543059 4324213412343432913758138673203834
z = [0,17,1,3,1,1,12,1,2,33,2,1,1,1,1,49,1,1,1,1,17,34,1,1,304,1,2,1,1,1,2,1,48,1,20,2,3,5,1,1,16,9,1,1,5,1,2,2,7,4,3,1,7,1,1,17,1,1,29,1,12,2,5]
z_0 = 0 / 1
z_1 = 1 / 17
z_2 = 1 / 18
z_3 = 4 / 71
z_4 = 5 / 89
...
z_18 = 23610961 / 419772458
...
z_40 = 832739221613445323225 / 14805030169237188131024
...
z_62 = 243224233245235253407096734543059 / 4324213412343432913758138673203834

Fractions are back insertable

Since fractions have a push_back method, you can use std::back_inserter to create a std::back_insert_iterator on your fraction. The code for computing a fraction from its quotients looks like:

See also
fraction.cpp
Fraction f;
OutputIterator itback = std::back_inserter( f );
for ( Quotient i = 1; i < argc; ++i)
{
Quotient u = atoll( argv[ i ] );
*itback++ = std::make_pair( u, i-1 );
}
std::cout << "z = " << f.p() << " / " << f.q() << std::endl;

which gives

# More precise approximation of pi
$ ./examples/arithmetic/fraction 3 7 15 1 292 1 1 1 2 1 3 1 14
z = 80143857 / 25510582

Computing reduced fractions

If \( z_k \) is a fraction, its 1-reduced is \( z_{k-1} \), its 2-reduced is \( z_{k-2} \) and so on. Reduced are in particular used for computing B├ęzout vectors, or in splitting formulas. You may use the methods reduced (SternBrocot::Fraction::reduced, LighterSternBrocot::Fraction::reduced) or partial (SternBrocot::Fraction::partial, LighterSternBrocot::Fraction::partial) to obtain them.

Assume f is the fraction \( 103993 / 33102 = [3;7,15,1,292] \). Then

Fraction z_3 = f.reduced( 1 ); // [3;7,15,1] but in fact [3;7,16]
Fraction z_2 = f.reduced( 2 ); // [3;7,15]
Note
The depth of a i-reduced fraction is necessarily \( k - i \). It may be \( k - i - 1 \).

You dispose also of the methods getSplit and getSplitBerstel to obtain directly the decomposition of a fraction into two fractions. Their "mediant sum" gives the fraction itself. Fractions are seen as vectors in this case.

See also
SternBrocot::Fraction::getSplit, LighterSternBrocot::Fraction::getSplit, SternBrocot::Fraction::getSplitBerstel, LighterSternBrocot::Fraction::getSplitBerstel.

Getting ancestors in the Stern-Brocot tree

Convergents or partials, and reduced fractions are ancestors of a given fractions. There exist in-between ancestors which are not given by these methods. We list them below with their exact definition, assuming f is \( [u_0,...,u_{k-1},u_k] \):

Inverse of a fraction

You can get the inverse of a fraction in O(1) with the method inverse (SternBrocot::Fraction::inverse, LighterSternBrocot::Fraction::inverse).

Rational approximation of floating-point number

A simple way to get rational approximation of a "real" number is to use the mapping: \( G(x)=\frac{1}{x- \lfloor x \rfloor} \). The successive \( \lfloor x \rfloor \) gives the successive quotients of the rational approximation. You may have a look at example approximation.cpp.

long double epsilon = 1e-14;
long double number0 = strtold( argv[ 1 ], 0 );
long double number = number0;
ASSERT( number >= 0.0 );
Fraction f;
OutputIterator itback = std::back_inserter( f );
Quotient i = 0;
while ( true )
{
long double int_part = floorl( number );
Quotient u = NumberTraits<long double>::castToInt64_t( int_part );
*itback++ = std::make_pair( u, i++ );
long double approx =
( (long double) NumberTraits<Integer>::castToDouble( f.p() ) )
/ ( (long double) NumberTraits<Integer>::castToDouble( f.q() ) );
std::cout << "z = " << f.p() << " / " << f.q()
<< " =~ " << std::setprecision( 16 ) << approx << std::endl;
number -= int_part;
if ( ( (number0 - epsilon ) < approx )
&& ( approx < (number0 + epsilon ) ) ) break;
number = 1.0 / number;
}
std::cout << "z = " << f.p() << " / " << f.q() << std::endl;

Standard arithmetic operations '+', '-', '*', '/'

There are not implemented yet and continued fractions are not the best representation for doing them. Here, you must use numerators and denominators to perform computations with these fractions.

Todo:
Arithmetic operations will soon be implemented with usual arithmetic operations on numerators and denominators. Furthermore, a dedicated system for performing arithmetic operations using only quotients (and guaranteed O(log(k)) complexity) is also under development.