JOIN
Get Time
forums   
Search | Watch Thread  |  My Post History  |  My Watches  |  User Settings
View: Flat  | Threaded  | Tree
Previous Thread  |  Next Thread
1.3. Using Rational Arithmetic to Avoid Floating Point Imprecision | Reply
Using Rational Arithmetic to Avoid Floating Point Imprecision

This is a rewrite of this recipe by Nemui.

Problem

As discussed in previous recipe, floating point numbers can be dangerous from imprecision point of view. Can they be avoided completely?

Solution

If the problem requires calculating square roots, exponents, logarithms etc., floating-point numbers can't be avoided. However, if all numbers in your calculations originate from arithmetics, you can eliminate floating-point numbers by using rational arithmetics and processing all numbers as fractions.

A common fraction consists of numerator N and denominator D (both integer), and represents a number N/D. The operations on common fractions are done as follows:

- addition/subtraction: (N1/D1) +/- (N2/D2) = (N1*D2 +/- N2*D1)/(D1*D2)
- multiplication: N1/D1 * N2/D2 = (N1*N2)/(D1*D2)
- division: (N1/D1) / (N2/D2) = (N1*D2)/(D1*N2)
- comparison: (N1/D1) ? (N2/D2) ~ N1*D2 ? N2*D1 (as long as D1*D2 is positive)

This simple technique can be surprisingly useful, especially if the problem involves a lot of math, and you have a pre-written class to handle fractions.

Discussion

Note that the more operations you do, the larger numerator and denominator of resulting fractions become. This might not be a problem if you're using just a few operations, but pre-written fractions class should take care of this by dividing both of them by their greatest common divisor - just in case. It's also a good idea to use long integers for storing numerator and denominator to be even safer; C++ users can implement Fraction class as a template and use whichever type fits the specific problem better.

Another thing to take care of is the sign of denominator: keeping it positive simplifies the things a lot, and requires only the line
if (d<0) { d=-d; n=-n; }

in implementation of division, since this is the only operation that doesn't necessarily keep the denominator positive.

Here is a sample implementation of fractions-handling class in C++. Note that it uses operators overloading; in Java you'll have to implement class methods. Another thing to note is shortening the definition of similar operands by using macros - this is handy when you don't want your pre-written code to violate Extra Code Rule.

#include <sstream>
#include <string>
 
using namespace std;
 
class Fraction {
    public:
    int N,D;
    int gcd(int a, int b) { return (a==0?b:gcd(b%a,a)); }
    void reduce() {
        int g = gcd(N,D);
        N/=g;
        D/=g;
    }
    Fraction(int N1 = 0, int D1 = 1) {
        N=N1;
        D=D1;
        reduce();
    }
    string toString() {
        stringstream A;
        A << N << "/" << D;
        return A.str();
    }
    double toDouble() {
        return ((double)(N))/D;
    }
};
 
//redefine operators
#define MATH(OP, NUM, DEN) Fraction operator OP (Fraction f1, Fraction f2) { return Fraction(NUM, DEN); }
MATH(+, f1.N*f2.D+f2.N*f1.D, f1.D*f2.D)
MATH(-, f1.N*f2.D-f2.N*f1.D, f1.D*f2.D)
MATH(*, f1.N*f2.N, f1.D*f2.D)
MATH(/, f1.N*f2.D, f1.D*f2.N)
 
#define ARITH(OP) Fraction& operator OP##=(Fraction&f1, Fraction f2) { return f1 = f1 OP f2; }
ARITH(+) ARITH(-) ARITH(*) ARITH(/)
 
#define COMP(OP) bool operator OP (Fraction f1, Fraction f2) { return (f1-f2).N OP 0; }
COMP(<) COMP(<=) COMP(>) COMP(>=) COMP(==) COMP(!=)


Now let's have a look at how rational arithmetics works in SRM problems.

ContiguousSubsequences
This problem is pretty straightforward, and possible floating-point imprecision is its only trick. We need quite a few operations, and the return value is based not on the value of the fractions but on the way there were obtained, so using the pre-written fractions class is not necessary:

public class ContiguousSubsequences {
	public int[] findMaxAverage(int[] a, int K) {
		int[] best = {0, K};
		int[] sum = new int[a.length+1];
		//sum[i] = a[0]+...+a[i-1]
		//sum[j]-sum[i] = a[i]+...+a[j-1]
		int i,j;
		for (i=1; i<=a.length; i++)
			sum[i] = sum[i-1]+a[i-1];
		for (i=0; i<a.length; i++)
		for (j=i+K; j<a.length+1; j++)
			if ((sum[j]-sum[i]) * (best[1]-best[0]) > (sum[best[1]]-sum[best[0]]) * (j-i) || 
			    (sum[j]-sum[i]) * (best[1]-best[0]) == (sum[best[1]]-sum[best[0]]) * (j-i) && 
			    (j-i)>(best[1]-best[0]) )
            {   best[0] = i;
                best[1] = j;
            }
		best[1]--;
		return best;
	}
}


On the contrary, in problem like BComputation there are a lot of computations, and all of them should be done in fractions, so a pre-written class helps a lot. Here is a possible solution:

#include <sstream>
#include <string>
 
using namespace std;
 
class Fraction {
    public:
    int N,D;
    int gcd(int a, int b) { return (a==0?b:gcd(b%a,a)); }
    void reduce() {
        int g = gcd(N,D);
        N/=g;
        D/=g;
		if (D<0) { D*=-1; N*=-1; }
    }
    Fraction(int N1 = 0, int D1 = 1) {
        N=N1;
        D=D1;
        reduce();
    }
    string toString() {
        stringstream A;
        A << N << "/" << D;
        return A.str();
    }
};
 
//redefine operators
#define MATH(OP, NUM, DEN) Fraction operator OP (Fraction f1, Fraction f2) { return Fraction(NUM, DEN); }
MATH(+, f1.N*f2.D+f2.N*f1.D, f1.D*f2.D)
MATH(-, f1.N*f2.D-f2.N*f1.D, f1.D*f2.D)
MATH(*, f1.N*f2.N, f1.D*f2.D)
MATH(/, f1.N*f2.D, f1.D*f2.N)
 
#define ARITH(OP) Fraction& operator OP##=(Fraction&f1, Fraction f2) { return f1 = f1 OP f2; }
ARITH(+) ARITH(-) ARITH(*) ARITH(/)
 
class BComputation {
	public:
	string getValue(int B0, int pos)  {
		int i,j;
		int C[20][20];
		memset(C, 0, sizeof(C));
		for(i=0; i<20; i++)
		{	C[i][0]=C[i][i]=1;
			for(j=1; j<i; j++)
				C[i][j]=C[i-1][j-1]+C[i-1][j];
		}
		Fraction sum, B[20];
		B[0] = B0;
		for (i=1; i<=pos; i++)
		{	sum = 0;
			for (j=0; j<i; j++)
				sum += B[j] * C[i+1][j];
			B[i] = sum / (-C[i+1][i]);
		}
		return B[pos].toString();
	}
};


It looks rather lengthy, but once you have pre-written code for both calculating binomial coefficients and processing fractions, the actual solution takes about 10 lines. Note that integers are converted into Fractions automatically, without explicit conversion.
Re: 1.3. Using Rational Arithmetic to Avoid Floating Point Imprecision (response to post by Nickolas) | Reply
Since the nominators and denominators can get very large, I think it would be safer to use long longs for them (or probably even better, make it possible to write it in a way which allows changing back to ints when time/space const turns out to be important).

For completeness, here is my implementation of fractions: (ll is long long, ld is long double)

//@? reduceGCD fll
template<class T> void reduceGCD(T& l, T& m) {
  T a = gcd(l,m); if(a<0) a=-a; if(m<0) a=-a; l/=a; m/=a; if(m==0&&l<0) l=-l;
  }
 
//@? fll
#define O(P,v) operator P (fll x, fll y) {return v;}
#define F(P,d) fll& operator P##= (fll&x, fll y) {d; return x;} fll O(P, x P##= y)
#define D(P) bool O(P, (x-y).l P 0)
 
struct fll {
  ll l; ll m;
  void red() {reduceGCD(l,m);}
  void unr(ll a) {l*=a; m*=a;}
  fll(ll x = 0, ll y = 1) {l=x; m=y; red();}
  };
 
F(+, ll a = gcd(x.m, y.m); ll b=x.m; x.unr(y.m/a); y.unr(b/a); x.l += y.l; x.red())
F(-, y.l = -y.l; x+=y)
F(*, swap(x.m, y.m); x.red(); y.red(); x.l *= y.l; x.m *= y.m; x.red())
F(/, swap(x.m, y.l); x.red(); y.red(); x.l *= y.m; x.m *= y.l; x.red())
 
D(<) D(<=) D(>) D(>=) D(==) D(!=)
 
//@? ftd
ld ftd(fll x) {return ld(x.l)/x.m;}
 
//@? fllts
string llts(ll x) {char buf[30]; sprintf(buf, "%lld", x); return buf;}
string fllts(fll x) {return llts(x.l) + "/" + llts(x.m);}
Re: 1.3. Using Rational Arithmetic to Avoid Floating Point Imprecision (response to post by Eryx) | Reply
Agreed, added a note about using long integers. Thanks!
RSS