'''
Rational numbers are a representation of the mathematical field of rational
numbers. They have unlimited precision, but many algorithms with rationals
use unbounded space and large amounts of time -- so be careful.

Rationals can be created easily:

>>> a=rational(1)
>>> a
rational(1L)

And are printed in friendly (though misleading ways)
 
>>> print a
1
>>> a/2
rational(1L, 2L)
>>> print a/2
1/2

There are many ways to build rationals:

From integers:
>>> print rational(1)
1

From nominators and denuminators:

>>> print rational(1,2)
1/2

From floats:

>>> print rational(1.3)
5854679515581645/4503599627370496
>>> print float(rational(1.3))
1.3

Not from complex numbers, even if they have 0 imaginary part:

>>> rational(1j)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
  File "Rational.py", line 418, in rational
    raise TypeError('cannot convert arguments')
TypeError: cannot convert arguments

But there is no problem with using floating point literals -- just
protect them by quotes.

>>> print rational("1.3")
13/10
>>> print rational("1.2e-3")
3/2500

From ``rational'' literals
>>> print rational("1/2")
1/2

Or even mix the two
>>> print rational("1.5/2.3")
15/23

Or give them as seperate arguments:
>>> print rational("1.5", "2.3")
15/23


*Note*: The following calculation takes some time

>>> p=0
>>> for i in range(1,1000):
...     p += rational(1)/i
... 

And p is *very* large. Rationals explode quickly in term of space
and (as you have noticed if you tried the above calculation) time.
 
>>> p
rational(53355784417020119952537879239887266136731803921522374204060897246465114565409520646006621457833121819822177013733448523121929191853817388455050878455561717371027592157489651884795570078464505321798967754860322188969172665004633935471096455470633645094270513262722579396248817332458071400971347691033193734596623333937737766140820373673275246317859525956885804716570122271771159715339438239613795876131660183846149167740477557199918997L, 7128865274665093053166384155714272920668358861885893040452001991154324087581111499476444151913871586911717817019575256512980264067621009251465871004305131072686268143200196609974862745937188343705015434452523739745298963145674982128236956232823794011068809262317708861979540791247754558049326475737829923352751796735248042463638051137034331214781746850878453485678021888075373249921995672056932029099390891687487672697950931603520000L)

printing p is not that user friendly, either.

>>> print p
53355784417020119952537879239887266136731803921522374204060897246465114565409520646006621457833121819822177013733448523121929191853817388455050878455561717371027592157489651884795570078464505321798967754860322188969172665004633935471096455470633645094270513262722579396248817332458071400971347691033193734596623333937737766140820373673275246317859525956885804716570122271771159715339438239613795876131660183846149167740477557199918997/7128865274665093053166384155714272920668358861885893040452001991154324087581111499476444151913871586911717817019575256512980264067621009251465871004305131072686268143200196609974862745937188343705015434452523739745298963145674982128236956232823794011068809262317708861979540791247754558049326475737829923352751796735248042463638051137034331214781746850878453485678021888075373249921995672056932029099390891687487672697950931603520000

We can convert p to float. Float conversion always works:

>>> float(p)
7.4844708605503447

even though naive conversion wouldn't:

>>> float(p.n)/float(p.d)
nan

To conserve time, we can trim p: find the closest rational with 
a bounded denominator
>>> p.trim(1L<<32)
rational(12377570989L, 1653767009L)

Subtracting them gives unhelpful results :(

>>> p.trim(1L<<32)-p
rational(-872119393953314540036963267350639988172842400460951889331776938576194338908055712270088035774611947221522363468929809379656051539773356282463641470813867367248573361743120219423901587966011300322682776045585594990374963929108968282002031347349387213913263250396231429588965548616844893074283339877989693620079553025660145795725133735701645191166053131402530951317054053830892681481265653921747553420724047240884052429689973L, 11789482202846854415241649104540583386623406115959677432802223340973331003735668809384664111066145039897075121304472266413879140983139770215358042356369522737429331262552172835890068328534070869038707353333385737580077488092224146480342885552420913445717377586934653792408698993535292433431351520245462060493779590806227120058195691087482315063384946077372429043555366594878942786082108265888537359480000638667859451202061272586716844271680000L)

But we can measure the error as a floating point number to get a rough
estimate:
>>> float(p.trim(1L<<32)-p)
-7.3974359428840768e-20

We can also approximate p: find the closest rational which is closer then
the bound (shifting rationals is the same as multiplying by 2 to the correct
power. Right shifting means negative powers, of course...)

>>> p.approximate(rational(1L)>>32)
rational(860063L, 114913L)
>>> float(p.approximate(rational(1L)>>32))
7.4844708605640786
>>> float(p.approximate(rational(1L)>>32))-p
1.3733902903823036e-11

How many bits is that?

>>> import math
>>> math.log(float(p.approximate(rational(1L)>>32))-p)/math.log(2)
-36.083467374587123

Less then -32, or in other words, the rational number we found is about
2.0**-36 further apart, closer then the bound (2.0**-32).
'''

def _gcd(a, b):
        if a>b:
                b, a = a, b
        if a == 0:
                return b
        while 1:
                c = b % a
                if c == 0:
                        return a
                b = a
                a = c 

def _trim(n, d, max_d):
        if max_d == 1:
                return n/d, 1
        last_n, last_d = 0, 1
        current_n, current_d = 1, 0
        while 1:
                div, mod = divmod(n, d)
                n, d = d, mod
                before_last_n, before_last_d = last_n, last_d
                next_n = last_n + current_n*div
                next_d = last_d + current_d*div
                last_n, last_d = current_n, current_d
                current_n, current_d = next_n, next_d
                if mod == 0 or current_d>=max_d:
                        break
        if current_d == max_d:
                return current_n, current_d
        i = (max_d-before_last_d)/last_d
        alternative_n = before_last_n + i*last_n
        alternative_d = before_last_d + i*last_d
        alternative = _Rational(alternative_n, alternative_d)
        last = _Rational(last_n, last_d)
        num = _Rational(n, d)
        if abs(alternative-num)<abs(last-num):
                return alternative_n, alternative_d 
        else:
                return last_n, last_d 

def _approximate(n, d, err):
        r = _Rational(n, d)
        last_n, last_d = 0, 1
        current_n, current_d = 1, 0
        while 1:
                div, mod = divmod(n, d)
                n, d = d, mod
                next_n = last_n + current_n*div
                next_d = last_d + current_d*div
                last_n, last_d = current_n, current_d
                current_n, current_d = next_n, next_d
                app = _Rational(current_n, current_d)
                if mod == 0 or abs(app-r)<err:
                        break
        return app

import math as _math

def _float_to_ratio(x):
        """\
        x -> (top, bot), a pair of co-prime longs s.t. x = top/bot.

        The conversion is done exactly, without rounding.
        bot > 0 guaranteed.
        Some form of binary fp is assumed.
        Pass NaNs or infinities at your own risk.

        >>> rational(10.0)
        rational(10L, 1L)
        >>> rational(0.0)
        rational(0L)
        >>> rational(-.25)
        rational(-1L, 4L)
        """

        if x == 0:
                return 0L, 1L
        signbit = 0
        if x < 0:
                x = -x
                signbit = 1
        f, e = _math.frexp(x)
        assert 0.5 <= f < 1.0
        # x = f * 2**e exactly

        # Suck up CHUNK bits at a time; 28 is enough so that we suck
        # up all bits in 2 iterations for all known binary double-
        # precision formats, and small enough to fit in an int.
        CHUNK = 28
        top = 0L
        # invariant: x = (top + f) * 2**e exactly
        while f:
                f = _math.ldexp(f, CHUNK)
                digit = int(f)
                assert digit >> CHUNK == 0
                top = (top << CHUNK) | digit
                f = f - digit
                assert 0.0 <= f < 1.0
                e = e - CHUNK
        assert top

        # now x = top * 2**e exactly; fold in 2**e
        r = _Rational(top, 1)
        if e>0:
                r = r << e
        else:
                r = r >> -e
        if signbit:
                return -r
        else:
                return r


class _Rational:

        def __init__(self, n, d):
                if d == 0:
                        return n/d
                n, d = map(long, (n, d))
                if d < 0:
                        n *= -1
                        d *= -1
                f = _gcd(abs(n), d)
                self.n = n/f
                self.d = d/f

        def __repr__(self):
                if self.d == 1:
                        return 'rational(%r)' % self.n
                return 'rational(%(n)r, %(d)r)' % self.__dict__

        def __str__(self):
                if self.d == 1:
                        return str(self.n)
                return '%(n)s/%(d)s' % self.__dict__

        def __coerce__(self, other):
                for int in (type(1), type(1L)):
                        if isinstance(other, int):
                                return self, rational(other)
                if type(other) == type(1.0):
                        return float(self), other
                return NotImplemented

        def __rcoerce__(self, other):
                return coerce(self, other)

        def __add__(self, other):
                return _Rational(self.n*other.d + other.n*self.d, 
                                 self.d*other.d)

        def __radd__(self, other):
                return self+other

        def __mul__(self, other):
                return _Rational(self.n*other.n, self.d*other.d)

        def __rmul__(self, other):
                return self*other

        def inv(self):
                return _Rational(self.d, self.n)

        def __div__(self, other):
                return self*other.inv()

        def __rdiv__(self, other):
                return self.inv()*other

        def __neg__(self):
                return _Rational(-self.n, self.d)

        def __sub__(self, other):
                return self+(-other)

        def __rsub__(self, other):
                return (-self)+other

        def __long__(self):
                if self.d != 1:
                        raise ValueError('cannot convert non-integer')
                return self.n

        def __int__(self):
                return int(long(self))

        def __float__(self):
                # Avoid NaNs like the plague
                if self.d > 1L<<1023:
                        self = self.trim(1L<<1023)
                return float(self.n)/float(self.d)

        def __pow__(self, exp, z=None):
                if z is not None:
                        raise TypeError('pow with 3 args unsupported')
                if isinstance(exp, _Rational):
                        if exp.d == 1:
                                exp = exp.n
                if isinstance(exp, type(1)) or isinstance(exp, type(1L)):
                        if exp < 0:
                                return _Rational(self.d**-exp, self.n**-exp)
                        return _Rational(self.n**exp, self.d**exp)
                return float(self)**exp

        def __cmp__(self, other):
                return cmp(self.n*other.d, self.d*other.n)

        def __hash__(self):
                return hash(self.n)^hash(self.d)

        def __abs__(self):
                return _Rational(abs(self.n), self.d)

        def __complex__(self):
                return complex(float(self))

        def __nonzero__(self):
                return self.n != 0

        def __pos__(self):
                return self

        def __oct__(self):
                return '%s/%s' % (oct(self.n), oct(self.d))

        def __hex__(self):
                return '%s/%s' % (hex(self.n), hex(self.d))

        def __lshift__(self, other):
                if other.d != 1:
                        raise TypeError('cannot shift by non-integer')
                return _Rational(self.n<<other.n, self.d)       

        def __rshift__(self, other):
                if other.d != 1:
                        raise TypeError('cannot shift by non-integer')
                return _Rational(self.n, self.d<<other.n)       

        def trim(self, max_d):
                n, d = self.n, self.d
                if n < 0:
                        n *= -1
                n, d = _trim(n, d, max_d)
                if self.n < 0:
                        n *= -1
                r = _Rational(n, d)
                upwards = self < r
                if upwards:
                        alternate_n = n-1
                else:
                        alternate_n = n+1
                if self == _Rational(alternate_n+n, d*2):
                        new_n = min(alternate_n, n)
                        return _Rational(new_n, d)
                return r

        def approximate(self, err):
                n, d = self.n, self.d
                if n < 0:
                        n *= -1
                app = _approximate(n, d, err)
                if self.n < 0:
                        app *= -1
                return app

def _parse_number(num):
        if '/' in num:
                n, d = num.split('/', 1)
                return _parse_number(n)/_parse_number(d)
        if 'e' in num:
                mant, exp = num.split('e', 1)
                mant = _parse_number(mant)
                exp = long(exp)
                return mant*(rational(10)**rational(exp))
        if '.' in num:
                i, f = num.split('.', 1)
                i = long(i)
                f = rational(long(f), 10L**len(f))
                return i+f
        return rational(long(num))

def rational(n, d=1L):
        if type(n) in (type(''), type(u'')) :
                n = _parse_number(n)
        if type(d) in (type(''), type(u'')) :
                d = _parse_number(d)
        if isinstance(n, type(1.0)):
                n = _float_to_ratio(n)
        if isinstance(d, type(1.0)):
                d = _float_to_ratio(d)
        for arg in (n, d):
                if isinstance(arg, type(1j)):
                        raise TypeError('cannot convert arguments')
        if isinstance(n, _Rational):
                return rational(n.n, n.d*d)
        if isinstance(d, _Rational):
                return rational(n*d.d, d.n)
        return _Rational(n, d)

import __builtin__
__builtin__.rational = rational

def _test():
        import doctest, Rational
        doctest.testmod(Rational)

if __name__ == '__main__':
        _test()