GCD and fraction class
published: Wed, 21-Apr-2004 | updated: Wed, 21-Apr-2004
Ah, there's a fun little algorithm coming up; one of the oldest known to man.
One of the things I like doing when I have some spare time is to peruse the list of search phrases that are used to access my site. Some of them are peculiar and some relate to me area of expertise. From April's collection of search phrases, I've selected a few that could be fleshed out in an article and so I'll be writing about them over the next few days.
Today we're going to look at the GCD algorithm to calculate the Greatest Common Divisor of two integers. This was first described by Euclid way back when in 300 B.C.E. or thereabouts. Essentially it says, for positive integers, that the GCD of two integers a and b, where a is greater than b, is equal to the GCD of b and (a mod b). This is a recursive definition, so we have to pay attention to the terminating cases. First, if a mod b is 1, the GCD of a and b is 1; a and b are called co-prime. Second, if a mod b is 0, then b divides a without remainder and so the GCD of a and b is b. I won't go into a proof here (if you want to see it, there are plenty of places on the Internet that explain it).
Let's see an example of one of these calculations by hand:
GCD(225, 141) = GCD(141, 84) since 225 mod 141 is 84
= GCD(84, 57) since 141 mod 84 is 57
= GCD(57, 27) since 84 mod 57 is 27
= GCD(27, 3) since 57 mod 27 is 3
= 3 since 27 mod 3 is 0.
A straight conversion of the definition into C# becomes this:
static public int GcdPrim(int a, int b) {
int rem = a % b;
if (rem == 0)
return b;
if (rem == 1)
return 1;
return GcdPrim(b, rem);
}
static public int wGcd(int a, int b) {
if ((a <= 0) || (b <= 0))
return 0;
if (a < b)
return GcdPrim(b, a);
else
return GcdPrim(a, b);
}
Notice how the function you call, Gcd(), prepares the data for the recursive routine. This is a common pattern when using recursive routines like this.
However, we note that the recursive call is at the end of the GcdPrim() recursive function, and so we can replace it with a loop and rearrange a little.
static private int GcdPrim(int a, int b) {
int rem = a % b;
while (rem > 1) {
a = b;
b = rem;
rem = a % b;
}
if (rem == 0)
return b;
return 1;
}
It's a pretty easy algorithm, so instead of just stopping there, let's write a fraction class and use it for real. Our fraction class will perform all the simple arithmetic operations: +, -, *, /. From your school days, you'll recollect the formulae for the basic operations:
a/b + c/d = (ad + bc)/bd a/b - c/d = (ad - bc)/bd a/b * c/d = ac/bd a/b / c/d = ad/bc
But of course we have to reduce the results to their simplest form. To do that, we have to find the GCD of the numerator and denominator and divide both by this. Here's the code in C#; note how I've taken into account the sign of the fractions in the private reduce() method. Notice also that the Fraction type is an immutable type. Once you create a Fraction instance with a particular value, you can't change that value.
using System;
namespace JmBucknall.Fractions {
public struct Fraction {
private int d;
private int n;
public Fraction(int numerator, int denominator) {
if (denominator == 0)
throw new ArgumentException("Denominator of fraction cannot be zero", "denominator");
this.n = numerator;
this.d = denominator;
reduce();
}
public int Denominator {
get { return d; }
}
public int Numerator {
get { return n; }
}
static private int GcdPrim(int a, int b) {
int rem = a % b;
while (rem > 1) {
a = b;
b = rem;
rem = a % b;
}
if (rem == 0)
return b;
return 1;
}
static private int Gcd(int a, int b) {
if (a < b)
return GcdPrim(b, a);
return GcdPrim(a, b);
}
private void reduce() {
bool isNegative = (n < 0) ^ (d < 0);
n = System.Math.Abs(n);
d = System.Math.Abs(d);
int gcd = Gcd(n, d);
n /= gcd;
d /= gcd;
if (isNegative)
n = -n;
}
// unary minus operator
static public Fraction operator- (Fraction f) {
return new Fraction(-f.n, f.d);
}
// addition operator
static public Fraction operator+ (Fraction f, Fraction g) {
int num = (f.n * g.d) + (f.d * g.n);
int denom = f.d * g.d;
return new Fraction(num, denom);
}
// subtraction operator
static public Fraction operator- (Fraction f, Fraction g) {
int num = (f.n * g.d) - (f.d * g.n);
int denom = f.d * g.d;
return new Fraction(num, denom);
}
// multiplication operator
static public Fraction operator* (Fraction f, Fraction g) {
int num = f.n * g.n;
int denom = f.d * g.d;
return new Fraction(num, denom);
}
// division operator
static public Fraction operator/ (Fraction f, Fraction g) {
int num = f.n * g.d;
int denom = f.d * g.n;
return new Fraction(num, denom);
}
}
}
Of course in writing such a type, we have to provide a little more plumbing. For a start we should override the Equals() and GetHashCode() methods as well as the equality and inequality operators. Similarly we override the ToString() method.
// overridden Equals method
public override bool Equals(object obj) {
if (obj == null)
return false;
if (obj.GetType() == typeof(Fraction))
return this == (Fraction) obj;
return false;
}
// overridden GetHashCode method
public override int GetHashCode() {
return (n << 3) + d;
}
// overidden ToString method
public override string ToString() {
return n.ToString() + '/' + d.ToString();
}
// equality operator
static public bool operator== (Fraction f, Fraction g) {
return (f.n == g.n) && (f.d == g.d);
}
// inequality operator
static public bool operator!= (Fraction f, Fraction g) {
return (f.n != g.n) || (f.d != g.d);
}
Finally, we should add a few conversion operators so that we can convert to and from integers and double values. The only one that is accurate and doesn't lose any information is the conversion from an integer to a fraction (and so it's marked as implicit), all the others will lose some information to a certain extent (and so are marked explicit). FOr example, in converting to an integer, we'll lose the fractional part; in converting to a double, we'll get a good approximation most of the time (the only times we will get an exact conversion is if the fraction's denominator is a power of 2); and in converting from a double we'll be calculating the most "reasonable" fraction from the double value.
// implicit conversion from int to Fraction
static public implicit operator Fraction(int n) {
return new Fraction(n, 1);
}
// explicit conversion from double to Fraction
static public explicit operator Fraction(double value) {
// uses the Continued Fraction algorithm
bool isNegative = (value < 0.0);
value = Math.Abs(value);
int prevPrevN = 0;
int prevPrevD = 1;
int prevN = 1;
int prevD = 0;
int a = (int) Math.Floor(value);
double x = value - Math.Floor(value);
int n = (a * prevN) + prevPrevN;
int d = (a * prevD) + prevPrevD;
double error = Math.Abs(value - (n / d));
while (error > 1e-6) {
prevPrevN = prevN;
prevPrevD = prevD;
prevD = d;
prevN = n;
x = 1.0 / x;
a = (int) Math.Floor(x);
x = x - Math.Floor(x);;
n = (a * prevN) + prevPrevN;
d = (a * prevD) + prevPrevD;
error = Math.Abs(value - ((double) n / d));
}
if (isNegative)
n = -n;
return new Fraction(n, d);
}
// explicit conversion from Fraction to double
static public explicit operator double(Fraction f) {
return f.n / f.d;
}
// explicit conversion from Fraction to int
static public explicit operator int(Fraction f) {
return (int) Math.Round((double) f);
}
The funky code in the conversion from a double to a fraction is an implementation of the Continued Fraction algorithm. This essentially uses the fact that a rational value can always be expressed as a finite continued fraction, and gives the method for calculating that continued fraction (and from that we can work out the fraction itself).
So let's talk a little about continued fractions. In the 19th century these were the subject of research for many mathematicians, and many important theorems were proved using them. but these days they've fallen a bit out of fashion. When I did my Maths degree, there was no course on continued fractions, nor were they used anywhere else. I did write a paper on them, however, for part of my degree.
A continued fraction looks like this:
a0 + b0
----------------------------
a1 + b1
-----------------------
a2 + b2
------------------
a3 + b3
-------------
a4 + b4
--------
a5 + ...
It looks wacky and it sometimes hurts just thinking about it (the ellipsis at the end is yet more continued fraction). Continued fractions can be finite (in which case the fraction ends somewhere), or they can be infinite. Lots of times, the b's are all assumed or forced to be one, and in this case you can write down the quotients of the continued fraction as if it were an array:
[a0, a1, a2, a3, a4, a5, ...]
When a continued fraction is finite, it represents a rational number, or a fraction. Infinite continued fractions with repeating quotients represent solutions to quadratic equations. Here's a couple of samples of continued fractions and what they equal.
3/8 = 0.375 = [0, 2, 1, 2] SQRT(2) = 1.4142135 = [1, 2, 2, 2, 2, 2, ...] Golden ratio = (1 + SQRT(5))/2 = [1, 1, 1, 1, 1, 1, 1, ....] e = 2.7182818 = [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, ...]
How do you evaluate a continued fraction? For a finite one we could evaluate it form the bottom up, but for an infinite one, we'd be out of luck. Consequently, the accepted way is to calculate the convergents of a continued fraction.
Using the continued fraction for e we proceed as follows. The first convergent is just the first quotient of the fraction, 2. The next convergent is the value of [2, 1], which, if you write it out, is 3. The next convergent is 2.6667. To see how to work this out longhand, follow along:
[2, 1, 2] = 2 + 1
-----
1 + 1
-
2
= 2 + 1 / (3/2)
= 2 + 2/3
= 2.6667
The next convergent is 2.75, the next 2.7143, the next 2.7188, and as you can see the convergents are slowly converging to e, 2.7183.
This calculation is still awkward, so it's a good job that convergents can be calculated using a simple iterative algorithm:
Assume the continued fraction is [a0, a1, a2, a3, ...] The convergents are the fractions Ni / Di, as i takes the values 0, 1, 2, 3, etc, where N0 = a0, D0 = 1 N1 = (a1 * N0) + 1, D1 = (a1 * D0) N2 = (a2 * N1) + N0, D2 = (a2 * D1) + D0 N3 = (a3 * N2) + N1, D3 = (a3 * D2) + D1 ...
So the next value for N is the next quotient multiplied by the previous value for N, with the previous previous value of N added in. Similarly for D. If you look back you'll see this algorithm embedded in the code, using variables like prevPrevN, prevN, and so forth, to hold previous values of the numerator N and denominator D.
Next we should show how to calculate the continued fraction for a floating-point value. Essentially it's very simple. Take the number, remove the integer part. This is the first quotient (a0) of the continued fraction. Find the inverse or reciprocal of the fractional part. Remove the integer part for the next quotient of the continued fraction, calculate the inverse of the factional part. Continue like this for as long as you need to.
The code above in the conversion method essentially merges these two algorithms in one loop, terminating the loop once the latest convergent is within a small tolerance to the original number.
So, in this article we've gone from Euclid's GCD algorithm, to a full- blown Fraction class, and then discussed an algorithm for calculating the most reasonable fraction given a double value. Have fun with it!
Update
I forgot the other comparison operators. Here's the implementation:
// less than operator
static public bool operator< (Fraction f, Fraction g) {
return (f.n * g.d) < (f.d * g.n);
}
// less than or equal operator
static public bool operator<= (Fraction f, Fraction g) {
return (f.n * g.d) <= (f.d * g.n);
}
// greater than operator
static public bool operator> (Fraction f, Fraction g) {
return (f.n * g.d) > (f.d * g.n);
}
// greater than operator
static public bool operator>= (Fraction f, Fraction g) {
return (f.n * g.d) >= (f.d * g.n);
}
