[boost] Re: rational does not seem to throw divide by zero?

From: Reggie Seagraves (reggie_at_[hidden])
Date: 2000-02-09 13:32:09

/*
Here's a couple of additional notes on 'rational'.....

        One place where rationals are often used is in graphics/video data
processing (especially when custom hardware is involved).  I.e., you typically
work in the main cpu with float/double, convert to rationals, then feed them
to a graphics processor.  (Yes, believe it or not, some graphic processors use
integral math, instead of using fixed point math!).  So, with that in mind, I'd
like to propose a couple of additions to your rational class....

Here's a few (sic) additional constructors.....
*/

        explicit rational(float num);		// Construct from a float
        explicit rational(double num);		// Construct from a double
        explicit rational(long double num); 	// Construct from a long double

        rational(float num, long maxD);		// with a maximum denominator
        rational(double num, long maxD);	// with a maximum denominator
        rational(long double num, long maxD);	// with max denominator

/*
The 'with max denominator' constructors specify maximum precision allowed
by the rational (i.e. indirectly, the number of bits used by the numerator and
denominator) which is also useful when dealing with foreign hardware.

and here's a 'utility' routine to convert from floating point values to
rationals that can be used in the new constructors.....
*/

#include <climits>

//-----------------------------------------------------------------
// Given a real number, this finds a fraction to represent it
// whose denominator is not greater than the provided maximum
// value.
//
// If there is a constraint of a maximum numerator, simply pass a
// maximum denominator whose value is maximum numerator/realNum.
//
// This implementation is guaranteed to work for real numbers of
// zero as well as negative numbers, although the range is more
// limited than standard floating point numbers since the maximum
// magnitude must fit into a long integer. However this is more
// accurate than fixed point.
//-----------------------------------------------------------------

void RealToRational(long double realNum, long maxD, long* numerator, long*
denominator) {
        long double n0, n1, n2;					//
numerator and denominator approximations
        long double d0, d1, d2;
        long double intPart;					// Integer
part of real number
        long double fracPart;					// Floating
point part of real number
        long double invertedPart;				// Inverse
of fractional part
        long double fMaxD = maxD;				// Keep as
floating point to avoid conversion later
        int neg = realNum < 0;					// Remember
the sign
        if (neg)
                realNum = -realNum;				// Code
below is for unsigned numbers
        intPart = (long)realNum;				//
Calculate first integer approximation
        fracPart = realNum - intPart;				// and the
fractional part
        n0 = 1;			d0 = 0;				// Initial
values for successive numerator and denominator approximations
        n1 = intPart;	d1 = 1;
        while ((fracPart != 0) && (1 / fracPart < fMaxD)) {	// Stop if
we found exact match, or will overflow
                invertedPart = 1 / fracPart;			// Get the
next inverse
                intPart = long(invertedPart);			// Get its
integer part
                fracPart = invertedPart - intPart;		// and its
fractional part
                n2 = intPart * n1 + n0;				// General
formula for next successive numerator
                d2 = intPart * d1 + d0;				// General
formula for next successive denominator
                if (d2 <= fMaxD && n2 <= LONG_MAX) {		// Are we
still in range?
                        n0 = n1;	d0 = d1;		// Shift
previous values
                        n1 = n2;	d1 = d2;
                } else {
                        intPart = (long)((fMaxD - d0) / d1);	// Try
converging from the other side
                        if (intPart > 0) {
                                n2 = intPart * n1 + n0;
                                if (n2 < LONG_MAX) {
        // Don't overflow numerator
                                        d2 = intPart * d1 + d0;
                                        long double delta1 = realNum - n1 /
d1;	// Last delta (before maxD)
                                        long double delta2 = realNum - n2 /
d2;	// New delta
                                        if (delta1 < 0)
        // get abs(delta1)
                                                delta1 = -delta1;
                                        if (delta2 < 0)
        // get abs(delta2)
                                                delta2 = -delta2;
                                        if (delta2 < delta1) {
        // Do final refinement to get the BEST approximation
                                                n1 = n2;
                                                d1 = d2;
                                        }
                                }
                        }
                        break;
                }
        }
        if (neg)
                n1 = -n1;
        *numerator   = long(n1);				// return
values
        *denominator = long(d1);
}

/*
        I hope you can get this integrated into your current rational....

                Thanks,
                        Reggie
*/