Rational Numbers
January 25, 2011
Scheme provides rational numbers natively, but we will make our own version to show the implementation. We will represent a rational number as a pair with the numerator in the car and denominator in the cdr; our convention will be that the denominator is always positive and the numerator and denominator will have no common factors. We begin with a function frac
that creates a fraction in the appropriate representation:
(define (frac n d)
(if (zero? d) (error 'frac "can't have zero denominator")
(if (negative? d) (frac (- n) (- d))
(let ((g (gcd n d)))
(if (= g 1) (cons n d)
(cons (/ n g) (/ d g)))))))
Then the other functions are from grade-school arithmetic; note that divide
doesn’t check for division by zero because frac
will report any error:
(define (plus x y)
(let ((a (car x)) (b (cdr x)) (c (car y)) (d (cdr y)))
(frac (+ (* a d) (* b c)) (* b d))))
(define (minus x y)
(let ((a (car x)) (b (cdr x)) (c (car y)) (d (cdr y)))
(frac (- (* a d) (* b c)) (* b d))))
(define (times x y)
(let ((a (car x)) (b (cdr x)) (c (car y)) (d (cdr y)))
(frac (* a c) (* b d))))
(define (divide x y)
(let ((a (car x)) (b (cdr x)) (c (car y)) (d (cdr y)))
(frac (* a d) (* b c))))
When my daughter was in fifth grade, she had trouble understanding the method they taught for comparing two fractions, which required the student to multiply each fraction by the denominator of the other fraction, then simplify by removing common factors, and finally to compare numerators. I had trouble understanding the need to simplify the fractions, so I showed her how to compare two fractions by cross-multiplying, and she understood immediately; the next day, when the students had a quiz on comparing fractions, she finished long before any of the other students. At the next parent-teacher conference, the teacher berated me for undermining his teaching, and ordered me not to do so again. Harrumph! Here’s the cross-multiplication method:
(define (less-than? x y)
(let ((a (car x)) (b (cdr x)) (c (car y)) (d (cdr y)))
(< (* a d) (* b c))))
Here are some examples:
> (plus (frac 1 3) (frac -1 7))
(4 . 21)
> (minus (frac 1 3) (frac -1 7))
(10 . 21)
> (times (frac 1 3) (frac -1 7))
(-1 . 21)
> (divide (frac 1 3) (frac -1 7))
(-7 . 3)
> (less-than? (frac 1 3) (frac -1 7))
#f
You can run the program at http://programmingpraxis.codepad.org/EnDFmHgs.
[…] today’s Programming Praxis exercise, our goal is to implement common operations on rational numbers. […]
I believe your multiplication function is wrong. 1/3 * -1/7 should produce -1/21, not -3/7, shouldn’t it?
My Haskell solution (see http://bonsaicode.wordpress.com/2011/01/25/programming-praxis-rational-numbers/ for a version with comments):
Fixed. I must have been asleep when I did that.
Knuth, volume 3, Seminumerical Algorithms, has some simple algorithms for rational arithmetic that computes GCDs on the smallest possible numbers. It complicates matters a bit, but it makes a big difference in some problems; you can find the code in _num.scm in the Gambit source code. (Search for (define-prim (##ratnum.= x y) …)
For example, with your code using Gambit’s bignum operations I get for
(define one (cons 1 1))
(define number-of-loops 1000)
(time
(do ((i 0 (+ i 1))
(r one (plus one (divide one r))))
((= i 1000) (display r)(newline))))
(time
(do ((i 0 (+ i 1))
(r 1 (+ 1 (/ 1 r))))
((= i 1000) (display r)(newline))))
gsi rational
(113796925398360272257523782552224175572745930353730513145086634176691092536145985470146129334641866902783673042322088625863396052888690096969577173696370562180400527049497109023054114771394568040040412172632376 . 70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501)
(time (do ((i 0 (+ i 1)) (r one (plus one (divide one r)))) ((= i 1000) (display r) (newline))))
217 ms real time
210 ms cpu time (210 user, 0 system)
390 collections accounting for 58 ms real time (50 user, 0 system)
108458544 bytes allocated
61 minor faults
no major faults
113796925398360272257523782552224175572745930353730513145086634176691092536145985470146129334641866902783673042322088625863396052888690096969577173696370562180400527049497109023054114771394568040040412172632376/70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501
(time (do ((i 0 (+ i 1)) (r 1 (+ 1 (/ 1 r)))) ((= i 1000) (display r) (newline))))
1 ms real time
10 ms cpu time (0 user, 10 system)
1 collection accounting for 0 ms real time (0 user, 0 system)
411264 bytes allocated
1 minor fault
no major faults
Seems reminiscent of SICP :-)
Since I implemented fractions as pairs, I included two procedures to get the
numerator and denominator. I also included pretty printing of fractions. My
Python solution:
Oops! That was an earlier version. The constructor should be:
A Ruby example
Thanks for the refresher in how to find the GCD, it’s been a long time since I had to do that….
In J, rational numbers are built-in, so this is really a demonstration of J’s extended integer and rational number arithmetic.
In J, a rational number is written as ‘NrD’, where ‘N’ and ‘D’ are integers.
If D is one, the ‘r1’ will be omitted. If D is zero, infinity is represented as “_”.
Negative numbers are indicated with a ‘_’ prefix, eg: ‘_2’ is negative 2.
Now that we now how to read and make rational numbers, let’s do some rational arithmetic in J (a very common thing to do in carpentry).
#include
typedef struct rational
{
int n;
int d;
} rational;
rational normalize(rational r)
{
int a = gcd(r.n, r.d);
r.n/=a;
r.d/=a;
return r;
}
rational add(rational r1, rational r2)
{
rational sum;
sum.d=(r1.d)*(r2.d);
sum.n=(r1.n)*(r2.d)+(r2.n)*(r1.d);
return normalize(sum);
}
rational subtract(rational r1, rational r2)
{
rational sum;
sum.d=(r1.d)*(r2.d);
sum.n=(r1.n)*(r2.d)-(r2.n)*(r1.d);
return normalize(sum);
}
rational multiply(rational r1, rational r2)
{
rational sum;
sum.d=(r1.d)*(r2.d);
sum.n=(r1.n)*(r2.n);
return normalize(sum);
}
rational divide(rational r1, rational r2)
{
rational sum;
sum.d=(r1.d)*(r2.n);
sum.n=(r1.n)*(r2.d);
return normalize(sum);
}
int l_t(rational r1, rational r2)
{
int lcm=(r1.*r2.d)/gcd(r1.d*r2.d);
if(r1.n*(lcm/r1.d)b){
a=a+b;
b=a-b;
a=a-b;
}
if(a==0) return b;
return gcd(b%a, a);
}
void print_rational(rational r)
{
printf(“%d/%d”, r.n, r.d);
}
void read_rational(rational *r)
{
scanf(“%d/%d”, &r->n, &r->d);
}
int main()
{
rational r1, r2;
read_rational(&r1);
read_rational(&r2);
print_rational(divide(r1, r2));
return 0;
}