Rational.ck

From CSWiki
Revision as of 14:33, 15 December 2009 by Rdpoor (talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Implementation of a Rational_Number_Package

See additional documentation in Rational_Number_Package.

/*
File: rational.ck -- represent numbers as quotient of two integers.

=== Overview
The Rational class lets you represent numbers as a quotient of two
integers.  One advantage of rational representation is that precision
is maintained across arithmetic operations: 2/3 is really 2/3 and not
.6666667.

=== Rational approximation
Of special note is the internal method that converts a floating point
value to a rational.  We use a modified version of Farey's technique,
which is efficient and numerically stable.  The method allows you to
specify the largest permissible denominator used to approximate the
floating point value.  Even limiting the denominator to less than 1000
produces rational approximations that are very close to their
irrational counterparts

=== Suggestions, feedback, bugfixes or bug reports are welcome:
rdpoor (at) gmail (dot) com

=== Revision History:
rdpoor	15 Nov 2009: initial version.

Copyright (c) 2009 Robert Poor

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.
*/

public class Rational {

    // ================================================================
    // class constants

    -1 => static int CMP_LT;
    0 => static int CMP_EQ;
    1 => static int CMP_GT;

    10000 => static int DEFAULT_DLIMIT;

    // ================================================================
    // constructors
    //
    // Implementation note: We assume (and enforce) that Rational values
    // are always normalized.

    fun static Rational create(int i) {
	return create(i, 1);
    }
    fun static Rational create(float v) {
	return create(v, DEFAULT_DLIMIT);
    }
    fun static Rational create(float v, int dlimit) {
	return farey(v, dlimit);
    }
    fun static Rational create(Rational other) {
	return create(other.num(), other.den());
    }
    fun static Rational create(int n, int d) {
	return (new Rational).set_(n, d).norm_();
    }

    // ================================================================
    // (private) instance variables

    int _num;
    int _den;

    // ================================================================
    // (public) instance methods
    //
    // NOTE: methods names that end in '_' signify destructive operations:
    // they modify 'this'.  All others leave 'this' untouched.

    fun int num() { return _num; }
    fun int den() { return _den; }

    fun Rational set_(Rational other) { 
	return set_(other.num(), other.den());
    }

    fun Rational set_(int n, int d) {
	n => _num;
	d => _den;
	return this;
    }

    fun float toFloat() { return (num() $ float) / (den() $ float); }

    fun string toString() { return num() + "/" + den(); }

    // ================
    // operations with integer results

    // compare receiver against 0
    fun int signum() {	return ((num()==0)?CMP_EQ:(num()<0)?CMP_LT:CMP_GT); }

    // compare receiver against integer n
    fun int cmp(int i) {
	num() - (i * den()) => int n;
	return (n==0)?CMP_EQ:((n<0)?CMP_LT:CMP_GT);
    }

    // compare receiver against other Rational
    fun int cmp(Rational other) { return sub(other).signum(); }

    // trunc() returns the largest integer value no greater in
    // magnitude than the receiver ("round towards zero")
    //
    fun int trunc() { return num() / den(); }
    
    // floor() returns the largest integer not greater than the
    // receiver ("round towards minus infinity").
    //
    // Implementation note: floor(x) is the same as trunc(x), except
    // when the residual (x - trunc(x)) is less than zero, in which
    // case the result is trunc(x) - 1.
    //
    // Note that "((x - trunc(x)) < 0)" and "(x < trunc(x))" are
    // equivalent.  We use the latter form to take advantage of the
    // more efficient integer cmp() function.
    fun int floor() {
	trunc() => int trunc;
	return (cmp(trunc) == CMP_LT)?trunc-1:trunc;
    }

    // ceil() returns the smallest integer not less than the receiver
    // ("round towards plus infinity").
    //
    // Implementation note: floor(x) is the same as trunc(x), except
    // when the residual (x - trunc(x)) is greater than zero, in which
    // case the result is trunc(x) + 1.
    fun int ceil() {
	trunc() => int trunc;
	return (cmp(trunc) == CMP_GT)?trunc+1:trunc;
    }

    // round() returns the closest integer to the receiver, rounding
    // to even when the receiver is halfway between two integers.
    //
    // Implementation note: round(x) is the same as trunc(x), except
    // when the residual (x - trunc(x)) is less or equal to than -.5,
    // in which case the result is trunc(x)-1, or if the residual is
    // greater than or equal to +.5, in which case the result is
    // trunc(x)+1.
    //
    // We scale the receiver as well as the residual by a factor of 2
    // so we can use the more efficient integer cmp() operations.
    fun int round() {
	mul(2) @=> Rational @ twothis;
	trunc() => int trunc;
	if (twothis.cmp(2*(trunc-1)) != CMP_GT) {
	    return trunc-1;	// residual <= -0.5
	} else if (twothis.cmp(2*(trunc+1)) != CMP_LT) {
	    return trunc+1;	// residual >= 0.5
	} else {
	    return trunc;	// -0.5 < residual < 0.5
	}
    }

    // Truncating division.
    fun int div(Rational other) {
	return (num() * other.den()) / (den() * other.num());
    }

    fun int div(int i) {
	return num() / (den() * i);
    }

    fun Rational mod(Rational other) {
	div(other) => int div;
	return this.sub(other.mul(div));
    }

    fun Rational mod(int i) {
	div(i) => int div;
	return this.sub(i * div);
    }

    // ================
    // operations with Rational results

    // Normalize a Rational.  Guarantess that ratio is in 
    // lowest form (## what's the term for that? ##) and
    // that the denominator is always positive.
    fun Rational norm() { return norm(new Rational); }
    fun Rational norm_() { return norm(this); }
    fun Rational norm(Rational result) {
	num() => int num;
	den() => int den;

	if (num == 0) {
	    return result.set_(0, 1);
	} else if (den == 0) {
	    return result.set_(1, 0);
	}
	gcd(num, den) => int gcd;
	if (((gcd > 0) && (den < 0)) || ((gcd < 0) && (den > 0))){
	    -num => num;
	    -den => den;
	}
	return result.set_(num/gcd, den/gcd);
    }

    // result = abs(this)
    fun Rational abs() { return abs(new Rational); }
    fun Rational abs_() { return abs(this); }
    fun Rational abs(Rational result) {
	result.set_(this);
	if (signum() == -1) result.negate_();
	return result;
    }

    // result = - this
    fun Rational negate() { return negate(new Rational); }
    fun Rational negate_() { return negate(this); }
    fun Rational negate(Rational result) {
	return result.set_(-num(), den());
    }

    // result = 1 / this
    fun Rational reciprocal() { return reciprocal(new Rational); }
    fun Rational reciprocal_() { return reciprocal(this); }
    fun Rational reciprocal(Rational result) { 
	return result.set_(den(), num());
    }

    // result = this + other
    fun Rational add(Rational other) { return add(other, new Rational); }
    fun Rational add_(Rational other) { return add(other, this); }
    fun Rational add(Rational other, Rational result) {
	return result.set_(num() * other.den() + other.num() * den(), 
			   den() * other.den()).norm_();
    }
    fun Rational add(int i) { return add(i, new Rational); }
    fun Rational add_(int i) { return add(i, this); }
    fun Rational add(int i, Rational result) {
	return result.set_(num() + (i * den()), den()).norm_();
    }

    // result = this - other
    fun Rational sub(Rational other) { return sub(other, new Rational); }
    fun Rational sub_(Rational other) { return sub(other, this); }
    fun Rational sub(Rational other, Rational result) {
	return result.set_(num() * other.den() - other.num() * den(), 
			   den() * other.den()).norm_();
    }
    fun Rational sub(int i) { return sub(i, new Rational); }
    fun Rational sub_(int i) { return sub(i, this); }
    fun Rational sub(int i, Rational result) {
	return result.set_(num() - (i * den()), den()).norm_();
    }

    // result = this * other
    fun Rational mul(Rational other) { return mul(other, new Rational); }
    fun Rational mul_(Rational other) { return mul(other, this); }
    fun Rational mul(Rational other, Rational result) {
	return result.set_(num() * other.num(), den() * other.den()).norm_();
    }
    fun Rational mul(int i) { return mul(i, new Rational); }
    fun Rational mul_(int i) { return mul(i, this); }
    fun Rational mul(int i, Rational result) {
	return result.set_(num() * i, den()).norm_();
    }
	
    // result = this / other
    fun Rational quo(Rational other) { return quo(other, new Rational); }
    fun Rational quo_(Rational other) { return quo(other, this); }
    fun Rational quo(Rational other, Rational result) {
	return result.set_(num() * other.den(), den() * other.num()).norm_();
    }
    fun Rational quo(int i) { return quo(i, new Rational); }
    fun Rational quo_(int i) { return quo(i, this); }
    fun Rational quo(int i, Rational result) {
	return result.set_(num(), den() * i).norm_();
    }
	

    // ================================================================
    // private methods

    // Return the closest Rational that approximates float value v
    // with a denominator of less than or equal to dlim using Farey's
    // technique.
    //
    // Gratefully converted from vegaseat's code at:
    //   http://www.python-forum.org/pythonforum/viewtopic.php?f=2&t=5692
    // which, in turn was derived from:
    //   http://en.wikipedia.org/wiki/Farey_series
    // 
    fun static Rational farey(float v, int dlimit) {
	if (v < 0) return farey(-v, dlimit).negate_();

	Rational lower, upper, mediant;
	lower.set_(0, 1);
	upper.set_(1, 0);
	while (true) {
	    mediant.set_(lower.num() + upper.num(), lower.den() + upper.den());
	    if (v * mediant.den() > mediant.num()) {
		if (dlimit < mediant.den()) return upper;
		lower.set_(mediant);
	    } else if (v * mediant.den() == mediant.num()) {
		if (dlimit >= mediant.den()) return mediant;
		if (lower.den() < upper.den()) return lower;
		return upper;
	    } else {
		if (dlimit < mediant.den()) return lower;
		upper.set_(mediant);
	    }
	}
    }

    // Return GCD of integers a and b using Euclid's method
    //
    fun static int gcd(int a, int b) {
	if (b == 0) 
	    return a;
	else 
	    return gcd(b, a % b);
    }


    // Return LCM of integers a and b.
    // 
    // Implementation note: To avoid overflow, we distribute the
    // division of gcd rather than returning (a*b)/gcd
    // 
    fun static int lcm(int a, int b) {
	gcd(a, b) => int gcd;
	return (a/gcd) * (b/gcd);
    }

} // class Rational