[Python-checkins] CVS: python/nondist/sandbox/rational Rational.py,NONE,1.1
Moshe Zadka
moshez@users.sourceforge.net
Tue, 27 Mar 2001 08:50:12 -0800
Update of /cvsroot/python/python/nondist/sandbox/rational
In directory usw-pr-cvs1:/tmp/cvs-serv2713
Added Files:
Rational.py
Log Message:
First CVS'ed version of Rational.py, prototype for PEP-0239.
--- NEW FILE: Rational.py ---
'''
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(n, 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()