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:
#!/usr/bin/env python def gcd(a, b): while b: a, b = b, a % b return a def frac(n, d): if d == 0: raise ZeroDivisionError else: sign = -1 if d < 0 else 1 g = gcd(n, d) return (sign * (n // g), abs(d) // g) def num(x): return x[0] def denom(x): return x[1] def plus(x, y): return frac(num(x) * denom(y) + denom(x) * num(y), denom(x) * denom(y)) def minus(x, y): return plus(x, frac(-num(y), denom(y))) def times(x, y): return frac(num(x) * num(y), denom(x) * denom(y)) def divide(x, y): return times(x, frac(denom(y), num(y))) def is_less_than(x, y): return num(x) * denom(y) < denom(x) * num(y) def frac_print(x): print str(num(x)) + "/" + str(denom(x)) return if __name__ == "__main__": frac_print(plus(frac(1, 3), frac(-1, 7))) frac_print(minus(frac(1, 3), frac(-1, 7))) frac_print(times(frac(1, 3), frac(-1, 7))) frac_print(divide(frac(1, 3), frac(-1, 7))) print(is_less_than(frac(1, 3), frac(-1, 7)))Oops! That was an earlier version. The constructor should be:
def frac(n, d): if d == 0: raise ZeroDivisionError else: sign = -1 if d < 0 else 1 ad = abs(d) g = gcd(n, ad) return (sign * (n // g), ad // g)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….
class RationalNumber include Comparable public attr_reader :denominator, :numerator def initialize(numerator, denominator = 1) raise "Division by zero" unless denominator != 0 if (denominator < 0) numerator = -numerator denominator = -denominator end @numerator = numerator.to_i @denominator = denominator.to_i reduce! end def -@ RationalNumber.new(-@numerator, @denominator) end def * other RationalNumber.new(@numerator * other.numerator, @denominator * other.denominator) end def / other self * other.invert end def + other RationalNumber.new(@numerator * other.denominator + other.numerator * @denominator, @denominator * other.denominator) end def - other self + (-other) end def <=> other @numerator * @denominator <=> other.numerator * other.denominator end protected def invert RationalNumber.new(@denominator, @numerator) end def reduce! gdc = gdc(@numerator, @denominator) @denominator /= gdc @numerator /= gdc end def gdc(a, b) a, b = b, a % b while b > 0 a end endIn 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).
do forever say 'Bitte Rechnung eingeben, z.B. 1/2 * 1/3' parse pull f1 op f2 if strip(f1) == '' & strip(op) == '' & strip(f2) == '' then leave say f1||' '||op||' '||f2||' = 'rat_rech(op,f1,f2) end exit rat_rech: parse arg op,f1,f2 parse value f1 with z1'/'n1 parse value f2 with z2'/'n2 if n1 == '' then n1 = 1 if n2 == '' then n2 = 1 select when op == '+' | op == '-' then do n = n1 * n2 interpret 'z = (z1 * n2) 'op' (z2 * n1)' end when op == '*' then do n = n1 * n2 z = z1 * z2 end when op == '/' then do n = n1 * z2 z = z1 * n2 end end return kuerzen(n,z) kuerzen: /* Zaehler und Teiler durch GGT teilen */ parse arg y,x z = 0 ggt = euklid(x,y) x = x / ggt y = y / ggt if x > y then do z = x % y x = x // y end /* Ganzzahlen separat */ if z > 0 then return strip(z||' '||x||'/'||y) if y == 1 then return strip(x) return strip(x||'/'||y) euklid: /* gcd */ parse arg a,b if a == 0 then return b do while b \= 0 if a > b then a = a - b else b = b - a end return a#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;
}