/*
Fraction class by M Phillips - 2005.

This code is provided as is with no warranties or guarantees of
any kind.

Please send an email to M Phillips (mbp2@i4free.co.nz)
  - if you use this file in a released product, or
  - if you find any bugs, or
  - if you have any suggestions
*/

#ifndef FRACTION_H
#define FRACTION_H

#include <sstream>
#include <iostream>
#include <cmath>
#include <limits>
#include <utility>
#include <iomanip>

// compile-time asserts (failure results in error C2118: negative subscript)
#ifndef C_ASSERT
#define ASSERT_CONCAT_(a, b) a##b
#define ASSERT_CONCAT(a, b) ASSERT_CONCAT_(a, b)
#define C_ASSERT(e) typedef char ASSERT_CONCAT(assert_line_, __LINE__)[(e)?1:-1]
#endif

// This class offers weak exception guarantees

template <class SINT_N, class SINT_2N>
class fraction {
	C_ASSERT(sizeof(SINT_2N) == 2 * sizeof(SINT_N));

protected:
	enum {
		MIN_DENOM = -SINT_N(~0UL>>1)-1,
		MAX_DENOM = SINT_N(~0UL>>1)
	};

	SINT_N numer, denom;

	template <typename T> void internalInitialiseFromFloating(T f) {
		if (f > MAX_DENOM) {
			// + infinity
			numer = 1;
			denom = 0;
			return;
		}
		if (f < MIN_DENOM) {
			// - infinity
			numer = -1;
			denom = 0;
			return;
		}
		if (f == ::floor(f)) {
			// whole number
			numer = SINT_N(::floor(f));
			denom = 1;
			return;
		}
		int sign = 1;
		if (f < T(0.0)) {
			sign = -1;
			f = -f;
		}
		if (f < T(1.0)/MAX_DENOM) {
			// too small to represent -> 0
			numer = 0;
			denom = 1;
			return;
		}
		T z = f;
		T prevNumer = 0;
		T prevDenom = 0;
		T fractDenom = 1;
		T fractNumer = 0;
		do {
			z = T(1.0)/(z - ::floor(z));
			T scratch = fractDenom;
			prevNumer = fractNumer;
			fractDenom = fractDenom * ::floor(z) + prevDenom;
			prevDenom = scratch;
			fractNumer = ::floor(f * fractDenom + T(0.5));
		} while (fractNumer < MAX_DENOM && fractDenom < MAX_DENOM &&
			::fabs(f - fractNumer/fractDenom) >= std::numeric_limits<T>::epsilon()*32 && z != ::floor(z));
		if (fractNumer < MAX_DENOM && fractDenom < MAX_DENOM) {
			numer = SINT_N(::floor(fractNumer)) * sign;
			denom = SINT_N(::floor(fractDenom));
		} else {
			numer = SINT_N(::floor(prevNumer)) * sign;
			denom = SINT_N(::floor(prevDenom));
		}
	}

	template <typename ostream_t>
	static ostream_t& output(ostream_t &os, const fraction &f) {
		typedef typename ostream_t::char_type char_t;

		std::ios::fmtflags fmtfl(os.flags());

		SINT_N base = 10;
		if (os.flags() & std::ios::hex) {
			base = 16;
		} else if (os.flags() & std::ios::oct)
			base = 8;

		std::streamsize width = os.width()-1;

		if (os.flags() & std::ios::left) {
			if (f.numer < 0 || os.flags() & std::ios::showpos)
				--width;
			SINT_N val = f.numer;
			do {
				val /= base;
				--width;
			} while (val != 0);

			os << std::setw(0) << f.numer << (char_t)'/' << std::setw(width);
		} else if (os.flags() & std::ios::internal) {

			os << std::setw((width+1)/2) << f.numer << (char_t)'/' << std::right << std::setw(width/2);
		} else {
			SINT_N val = f.denom;
			do {
				val /= base;
				--width;
			} while (val != 0);

			os << std::setw(width) << f.numer << std::setw(0) << (char_t)'/';
		}
		os.unsetf(std::ios::showpos);
		os << f.denom;
		os.flags(fmtfl);
		return os;
	}
	
	template <typename istream_t>
	static istream_t& input(istream_t &is, fraction &val) {
		typedef typename istream_t::char_type char_t;
		SINT_N numer, denom;
		std::ios::iostate except(is.exceptions());
		is.exceptions(std::ios::goodbit);

		if (is >> numer) {
			char_t ch;
			is >> ch;
			if (ch != (char_t)'/') {
				is.putback(ch);
				is.setstate(std::ios::failbit);
			} else if (!is.fail()) {
				if (is >> denom)
					val = fraction(numer, denom);
			}
		}

		is.exceptions(except);
		return is;
	}

public:
	fraction() : numer(0), denom(0) {}	// initialise to nan
	fraction(SINT_N x) : numer(x), denom(1) {}
/*
	fraction(char				x) : numer(x), denom(1) {}
	fraction(signed char		x) : numer(x), denom(1) {}
	fraction(short				x) : numer(x), denom(1) {}
	fraction(int				x) : numer(x), denom(1) {}
	fraction(long				x) : numer(x), denom(1) {}
	fraction(long long			x) : numer(x), denom(1) {}
	fraction(unsigned char		x) : numer(x), denom(1) {}
	fraction(unsigned short		x) : numer(x), denom(1) {}
	fraction(unsigned int		x) : numer(x), denom(1) {}
	fraction(unsigned long		x) : numer(x), denom(1) {}
	fraction(unsigned long long x) : numer(x), denom(1) {}
*/
	template <class T> fraction(const T &x, const T &y) {
		simplify(x, y);
	}

	fraction(float f) {
		internalInitialiseFromFloating(f);
	}
	fraction(double d) {
		internalInitialiseFromFloating(d);
	}
	fraction(long double ld) {
		internalInitialiseFromFloating(ld);
	}

	friend void swap(fraction &a, fraction &b) {
		using std::swap;
		swap(a.numer, b.numer);
		swap(a.denom, b.denom);
	}

	// These offer only weak exception guarantees, as do member functions below which use them
	template <class T> fraction &simplify(const T &x, const T &y) {
		if (y == 0) {
			// + or - infinity, ensure numer is +1 or -1 for equality tests
			if (x < 0)
				numer = -1;
			else if (x > 0)
				numer = 1;
			else
				numer = 0;	// 0/0 = NAN!
			denom = 0;
		} else if (x == 0) {
			// ensure zero is always represented as 0/1, not 0/n, etc for equality tests
			numer = 0;
			denom = 1;
		} else if (y == 1) {
			// special (faster) case for whole numbers
			numer = static_cast<SINT_N>(x);
			denom = 1;
		} else {
			// gcd algorithm
			T a = x, b = y;
			do {
				T z = b;
				b = a % b;
				a = z;
			} while (b != 0);
			// divide by simplification factor
			numer = static_cast<SINT_N>(x / a);
			denom = static_cast<SINT_N>(y / a);
			// check for negatives, and ensure only numer can be negative (for equality tests and output formatting)
			if (denom < 0) {
				denom = -denom;
				numer = -numer;
			}
		}
		return *this;
	}

	friend std::ostream& operator << (std::ostream &os, const fraction &f) {
		return output(os, f);
	}
	friend std::wostream& operator << (std::wostream &wos, const fraction &f) {
		return output(wos, f);
	}

	friend std::istream& operator >> (std::istream &is, fraction &val) {
		return input(is, val);
	}
	friend std::wistream& operator >> (std::wistream &wis, fraction &val) {
		return input(wis, val);
	}

	// End of only weak exception guarantee

	// misc functions
	friend const fraction abs(const fraction &a) {
		fraction b;
		b.numer = a.numer < 0 ? -a.numer : a.numer;
		b.denom = a.denom;
		return b;
	}
	friend const fraction sqr(const fraction &a) {
		fraction b;
		b.numer = a.numer*a.numer;
		b.denom = a.denom*a.denom;
		return b;
	}
	friend const fraction pow(fraction a, int b) {
		if (b < 0)
			return (fraction)0;
		if (a.numer == a.denom)
			return a;
		fraction result(1);
		while (b) {
			if (b & 1) {
				result.numer *= a.numer;
				result.denom *= a.denom;
			}
			b >>= 1;
			a.numer *= a.numer;
			a.denom *= a.denom;
		}
		return result;
	}
	static const fraction infinity() {
		return fraction(1, 0);
	}

	std::string toString() const {
		return *this;
	}

	std::wstring toWString() const {
		return *this;
	}

	operator std::string() const {
		std::stringstream ss;
		ss << *this;
		return ss.str();
	}

	operator std::wstring() const {
		std::wstringstream wss;
		wss << *this;
		return wss.str();
	}

	bool isValid() const {
		fraction f;
		return ((numer != 0 || denom == 1)
			&& (denom != 0 || numer == 1 || numer == -1)
			&& simplify(numer, denom).numer == numer);
	}

	// quantisation functions
	friend SINT_N trunc(const fraction &f) {
		return f.numer / f.denom;
	}
	friend SINT_N floor(const fraction &f) {
		return ((f.numer >= 0) ? f.numer : f.numer-f.denom+1) / f.denom;
	}
	friend SINT_N ceil(const fraction &f) {
		return ((f.numer >= 0) ? f.numer+f.denom-1 : f.numer) / f.denom;
	}
	friend SINT_N round(const fraction &f) {
		return ((f.numer >= 0) ? f.numer+(f.denom>>1) : f.numer-(f.denom>>1)) / f.denom;
	}

	operator SINT_N() const { return numer / denom; }
/*
	operator char() const				{ return numer / denom; }
	operator signed char() const		{ return numer / denom; }
	operator short() const				{ return numer / denom; }
	operator int() const				{ return numer / denom; }
	operator long() const				{ return numer / denom; }
	operator long long() const			{ return numer / denom; }
	operator unsigned char() const		{ return numer / denom; }
	operator unsigned short() const		{ return numer / denom; }
	operator unsigned int() const		{ return numer / denom; }
	operator unsigned long() const		{ return numer / denom; }
	operator unsigned long long() const	{ return numer / denom; }
*/
	operator float() const			{ return (float)numer / (float)denom; }
	operator double() const			{ return (double)numer / (double)denom; }
	operator long double() const	{ return (long double)numer / (long double)denom; }

	SINT_N getNumerator() const {
		return numer;
	}

	SINT_N getDenominator() const {
		return denom;
	}

	// 'invert' is the same as performing x = 1/x
	fraction& invert() {
		using std::swap;
		swap(numer, denom);
		if (denom == 0) {
			// + or - infinity, ensure numer is +1 or -1 for equality tests
			if (numer < 0)
				numer = -1;
			else if (numer > 0)
				numer = 1;
		} else if (numer == 0) {
			// ensure zero is always represented as 0/1, not 0/n, etc for equality tests
			denom = 1;
		} else if (denom < 0) {
			denom = -denom;
			numer = -numer;
		}
		return *this;
	}

	// modifying binary operators
	fraction& operator += (const fraction &a) {
		if (denom == a.denom)
			return simplify(numer + a.numer, denom);
		else
			return simplify(SINT_2N(numer)*SINT_2N(a.denom) + SINT_2N(a.numer)*SINT_2N(denom), SINT_2N(denom) * SINT_2N(a.denom));
	}
	fraction& operator -= (const fraction &a) {
		if (denom == a.denom)
			return simplify(numer - a.numer, denom);
		else
			return simplify(SINT_2N(numer)*SINT_2N(a.denom) - SINT_2N(a.numer)*SINT_2N(denom), SINT_2N(denom) * SINT_2N(a.denom));
	}
	fraction& operator *= (const fraction &a) {
		return simplify(SINT_2N(numer)*SINT_2N(a.numer), SINT_2N(denom) * SINT_2N(a.denom));
	}
	fraction& operator /= (const fraction &a) {
		return simplify(SINT_2N(numer)*SINT_2N(a.denom), SINT_2N(denom) * SINT_2N(a.numer));
	}

	const fraction& operator>>= (int shift) {
		return simplify(SINT_2N(numer), SINT_2N(denom)<<shift);
	}
	const fraction& operator<<= (int shift) {
		return simplify(SINT_2N(numer)<<shift, SINT_2N(denom));
	}

	// modifying unary operators
	const fraction operator++(int) {	// Post increment operator
		fraction a = (const fraction&)*this;
		numer+=denom;
		return a;
	}
	const fraction operator--(int) {	// Post decrement operator
		fraction a = (const fraction&)*this;
		numer-=denom;
		return a;
	}
	const fraction& operator++() {	// Pre increment operator
		numer+=denom;
		return *this;
	}
	const fraction& operator--() {	// Pre decrement operator
		numer-=denom;
		return *this;
	}

	// binary operators
	friend const fraction operator + (fraction a, const fraction &b) {
		return a += b;
	}
	friend const fraction operator - (fraction a, const fraction &b) {
		return a -= b;
	}
	friend const fraction operator * (fraction a, const fraction &b) {
		return a *= b;
	}
	friend const fraction operator / (fraction a, const fraction &b) {
		return a /= b;
	}
	friend bool operator < (const fraction &a, const fraction &b) {
		return SINT_2N(a.numer) * SINT_2N(b.denom) < SINT_2N(a.denom) * SINT_2N(b.numer);
	}
	friend bool operator > (const fraction &a, const fraction &b) {
		return SINT_2N(a.numer) * SINT_2N(b.denom) > SINT_2N(a.denom) * SINT_2N(b.numer);
	}
	friend bool operator <=(const fraction &a, const fraction &b) {
		return SINT_2N(a.numer) * SINT_2N(b.denom) <=SINT_2N(a.denom) * SINT_2N(b.numer);
	}
	friend bool operator >=(const fraction &a, const fraction &b) {
		return SINT_2N(a.numer) * SINT_2N(b.denom) >=SINT_2N(a.denom) * SINT_2N(b.numer);
	}
	friend bool operator ==(const fraction &a, const fraction &b) {
		return a.numer == b.numer && b.denom == a.denom;
	}
	friend bool operator !=(const fraction &a, const fraction &b) {
		return a.numer != b.numer || b.denom != a.denom;
	}

	const fraction operator >> (int shift) const {
		return fraction(SINT_2N(numer), SINT_2N(denom)<<shift);
	}
	const fraction operator << (int shift) const {
		return fraction(SINT_2N(numer)<<shift, SINT_2N(denom));
	}

	// unary operators
	bool operator ! () const {
		return !numer;
	}
	const fraction& operator + () const {
		return *this;
	}
	const fraction operator - () const {
		fraction a;
		a.numer = -numer;
		a.denom = denom;
		return a;
	}

	// base-type binary operators
	friend const fraction operator + (const SINT_N &a, const fraction &b) {
		if (b.denom == 1)
			return fraction(a + b.numer, b.denom);
		else
			return fraction(SINT_2N(a)*SINT_2N(b.denom) + SINT_2N(b.numer), SINT_2N(b.denom));
	}
	friend const fraction operator - (const SINT_N &a, const fraction &b) {
		if (b.denom == 1)
			return fraction(a - b.numer, b.denom);
		else
			return fraction(SINT_2N(a)*SINT_2N(b.denom) - SINT_2N(b.numer), SINT_2N(b.denom));
	}
	friend const fraction operator * (const SINT_N &a, const fraction &b) {
		return fraction(SINT_2N(a) * SINT_2N(b.numer), SINT_2N(b.denom));
	}
	friend const fraction operator / (const SINT_N &a, const fraction &b) {
		return fraction(SINT_2N(b.numer), SINT_2N(a) * SINT_2N(b.denom));
	}
	friend bool operator < (const SINT_N &a, const fraction &b) {
		return SINT_2N(a)*SINT_2N(b.denom) < SINT_2N(b.numer);
	}
	friend bool operator > (const SINT_N &a, const fraction &b) {
		return SINT_2N(a)*SINT_2N(b.denom) > SINT_2N(b.numer);
	}
	friend bool operator <=(const SINT_N &a, const fraction &b) {
		return SINT_2N(a)*SINT_2N(b.denom) <=SINT_2N(b.numer);
	}
	friend bool operator >=(const SINT_N &a, const fraction &b) {
		return SINT_2N(a)*SINT_2N(b.denom) >=SINT_2N(b.numer);
	}
	friend bool operator ==(const SINT_N &a, const fraction &b) {
		return SINT_2N(a)*SINT_2N(b.denom) ==SINT_2N(b.numer);
	}
	friend bool operator !=(const SINT_N &a, const fraction &b) {
		return SINT_2N(a)*SINT_2N(b.denom) !=SINT_2N(b.numer);
	}
	friend const fraction operator + (const fraction &a, const SINT_N &b) {
		if (a.denom == 1)
			return fraction(a.numer + b, a.denom);
		else
			return fraction(SINT_2N(a.numer) + SINT_2N(b)*SINT_2N(a.denom), SINT_2N(a.denom));
	}
	friend const fraction operator - (const fraction &a, const SINT_N &b) {
		if (a.denom == 1)
			return fraction(a.numer - b, a.denom);
		else
			return fraction(SINT_2N(a.numer) + SINT_2N(b)*SINT_2N(a.denom), SINT_2N(a.denom));
	}
	friend const fraction operator * (const fraction &a, const SINT_N &b) {
		return fraction(SINT_2N(a.numer) * SINT_2N(b), SINT_2N(a.denom));
	}
	friend const fraction operator / (const fraction &a, const SINT_N &b) {
		return fraction(SINT_2N(a.numer), SINT_2N(a.denom) * SINT_2N(b));
	}
	friend bool operator < (const fraction &a, const SINT_N &b) {
		return SINT_2N(a.numer) < SINT_2N(b)*SINT_2N(a.denom);
	}
	friend bool operator > (const fraction &a, const SINT_N &b) {
		return SINT_2N(a.numer) > SINT_2N(b)*SINT_2N(a.denom);
	}
	friend bool operator <=(const fraction &a, const SINT_N &b) {
		return SINT_2N(a.numer) <=SINT_2N(b)*SINT_2N(a.denom);
	}
	friend bool operator >=(const fraction &a, const SINT_N &b) {
		return SINT_2N(a.numer) >=SINT_2N(b)*SINT_2N(a.denom);
	}
	friend bool operator ==(const fraction &a, const SINT_N &b) {
		return SINT_2N(a.numer) ==SINT_2N(b)*SINT_2N(a.denom);
	}
	friend bool operator !=(const fraction &a, const SINT_N &b) {
		return SINT_2N(a.numer) !=SINT_2N(b)*SINT_2N(a.denom);
	}

	// float binary operators
	friend const fraction operator + (float a, const fraction &b) { return fraction(a) + b; }
	friend const fraction operator - (float a, const fraction &b) { return fraction(a) - b; }
	friend const fraction operator * (float a, const fraction &b) { return fraction(a) * b; }
	friend const fraction operator / (float a, const fraction &b) { return fraction(a) / b; }
	friend const fraction operator + (const fraction &a, float b) { return a + fraction(b); }
	friend const fraction operator - (const fraction &a, float b) { return a - fraction(b); }
	friend const fraction operator * (const fraction &a, float b) { return a * fraction(b); }
	friend const fraction operator / (const fraction &a, float b) { return a / fraction(b); }
	friend bool operator < (float a, const fraction &b) { return a * (float)b.denom < (float)b.numer; }
	friend bool operator > (float a, const fraction &b) { return a * (float)b.denom > (float)b.numer; }
	friend bool operator <=(float a, const fraction &b) { return a * (float)b.denom <=(float)b.numer; }
	friend bool operator >=(float a, const fraction &b) { return a * (float)b.denom >=(float)b.numer; }
	friend bool operator ==(float a, const fraction &b) { return a * (float)b.denom ==(float)b.numer; }
	friend bool operator !=(float a, const fraction &b) { return a * (float)b.denom !=(float)b.numer; }
	friend bool operator < (const fraction &a, float b) { return (float)a.numer < (float)a.denom * b; }
	friend bool operator > (const fraction &a, float b) { return (float)a.numer > (float)a.denom * b; }
	friend bool operator <=(const fraction &a, float b) { return (float)a.numer <=(float)a.denom * b; }
	friend bool operator >=(const fraction &a, float b) { return (float)a.numer >=(float)a.denom * b; }
	friend bool operator ==(const fraction &a, float b) { return (float)a.numer ==(float)a.denom * b; }
	friend bool operator !=(const fraction &a, float b) { return (float)a.numer !=(float)a.denom * b; }

	// double binary operators
	friend const fraction operator + (double a, const fraction &b) { return fraction(a) + b; }
	friend const fraction operator - (double a, const fraction &b) { return fraction(a) - b; }
	friend const fraction operator * (double a, const fraction &b) { return fraction(a) * b; }
	friend const fraction operator / (double a, const fraction &b) { return fraction(a) / b; }
	friend const fraction operator + (const fraction &a, double b) { return a + fraction(b); }
	friend const fraction operator - (const fraction &a, double b) { return a - fraction(b); }
	friend const fraction operator * (const fraction &a, double b) { return a * fraction(b); }
	friend const fraction operator / (const fraction &a, double b) { return a / fraction(b); }
	friend bool operator < (double a, const fraction &b) { return a * (double)b.denom < (double)b.numer; }
	friend bool operator > (double a, const fraction &b) { return a * (double)b.denom > (double)b.numer; }
	friend bool operator <=(double a, const fraction &b) { return a * (double)b.denom <=(double)b.numer; }
	friend bool operator >=(double a, const fraction &b) { return a * (double)b.denom >=(double)b.numer; }
	friend bool operator ==(double a, const fraction &b) { return a * (double)b.denom ==(double)b.numer; }
	friend bool operator !=(double a, const fraction &b) { return a * (double)b.denom !=(double)b.numer; }
	friend bool operator < (const fraction &a, double b) { return (double)a.numer < (double)a.denom * b; }
	friend bool operator > (const fraction &a, double b) { return (double)a.numer > (double)a.denom * b; }
	friend bool operator <=(const fraction &a, double b) { return (double)a.numer <=(double)a.denom * b; }
	friend bool operator >=(const fraction &a, double b) { return (double)a.numer >=(double)a.denom * b; }
	friend bool operator ==(const fraction &a, double b) { return (double)a.numer ==(double)a.denom * b; }
	friend bool operator !=(const fraction &a, double b) { return (double)a.numer !=(double)a.denom * b; }
};

#endif