'''
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
>>> print rational("-1.3")
-13/10
>>> print rational("-.25")
-1/4

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)
Traceback (most recent call last):
  File "<string>", line 1, in ?
OverflowError: long int too large to convert to float

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-020

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.3733999215941832e-011

How many bits is that?

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

Less than -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)
    >>> rational(0.0)
    rational(0L)
    >>> rational(-.25)
    rational(-1L, 4L)
    """

    if x == 0:
        return _Rational(0, 1)
    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(object):

    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, o):
        other = rational(o)
        return _Rational(self.n*other.d + other.n*self.d, 
                         self.d*other.d)

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

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

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

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

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

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

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

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

    def __rsub__(self, other):
        return (-self)+rational(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), type(1L))):
            if exp < 0:
                return _Rational(self.d**-exp, self.n**-exp)
            return _Rational(self.n**exp, self.d**exp)
        return float(self)**float(exp)

    def __rpow__(self, base, z=None):
        return rational(base) ** self

    def __cmp__(self, o):
        other = rational(o)
        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, o):
        other = rational(o)
        if other.d != 1:
            raise TypeError('cannot shift by non-integer')
        return _Rational(self.n<<other.n, self.d)   

    def __rshift__(self, o):
        other = rational(o)
        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 num.strip().startswith("-"):
        num = num.strip()[1:]
        sign = -1
    else:
        sign = 1
    if '/' in num:
        n, d = num.split('/', 1)
        return sign * _parse_number(n) / _parse_number(d)
    if 'e' in num:
        mant, exp = num.split('e', 1)
        mant = _parse_number(mant)
        exp = long(exp)
        return sign * mant * (rational(10)**rational(exp))
    if '.' in num:
        i, f = num.split('.', 1)
        if i == "":
            i = 0
        else:
            i = long(i)
        f = rational(long(f), 10L**len(f))
        return sign * (i+f)
    return rational(sign * 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):
        if d == 1L:
            return n
        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()