619 lines
19 KiB
Python
619 lines
19 KiB
Python
#!/usr/bin/env python
|
|
"""
|
|
FixedPoint objects support decimal arithmetic with a fixed number of
|
|
digits (called the object's precision) after the decimal point. The
|
|
number of digits before the decimal point is variable & unbounded.
|
|
|
|
The precision is user-settable on a per-object basis when a FixedPoint
|
|
is constructed, and may vary across FixedPoint objects. The precision
|
|
may also be changed after construction via FixedPoint.set_precision(p).
|
|
Note that if the precision of a FixedPoint is reduced via set_precision,
|
|
information may be lost to rounding.
|
|
|
|
>>> x = FixedPoint("5.55") # precision defaults to 2
|
|
>>> print x
|
|
5.55
|
|
>>> x.set_precision(1) # round to one fraction digit
|
|
>>> print x
|
|
5.6
|
|
>>> print FixedPoint("5.55", 1) # same thing setting to 1 in constructor
|
|
5.6
|
|
>>> repr(x) # returns constructor string that reproduces object exactly
|
|
"FixedPoint('5.6', 1)"
|
|
>>>
|
|
|
|
When FixedPoint objects of different precision are combined via + - * /,
|
|
the result is computed to the larger of the inputs' precisions, which also
|
|
becomes the precision of the resulting FixedPoint object.
|
|
|
|
>>> print FixedPoint("3.42") + FixedPoint("100.005", 3)
|
|
103.425
|
|
>>>
|
|
|
|
When a FixedPoint is combined with other numeric types (ints, floats,
|
|
strings representing a number) via + - * /, then similarly the computation
|
|
is carried out using-- and the result inherits --the FixedPoint's
|
|
precision.
|
|
|
|
>>> print FixedPoint(1) / 7
|
|
0.14
|
|
>>> print FixedPoint(1, 30) / 7
|
|
0.142857142857142857142857142857
|
|
>>>
|
|
|
|
The string produced by str(x) (implictly invoked by "print") always
|
|
contains at least one digit before the decimal point, followed by a
|
|
decimal point, followed by exactly x.get_precision() digits. If x is
|
|
negative, str(x)[0] == "-".
|
|
|
|
The FixedPoint constructor can be passed an int, long, string, float,
|
|
FixedPoint, or any object convertible to a float via float() or to a
|
|
long via long(). Passing a precision is optional; if specified, the
|
|
precision must be a non-negative int. There is no inherent limit on
|
|
the size of the precision, but if very very large you'll probably run
|
|
out of memory.
|
|
|
|
Note that conversion of floats to FixedPoint can be surprising, and
|
|
should be avoided whenever possible. Conversion from string is exact
|
|
(up to final rounding to the requested precision), so is greatly
|
|
preferred.
|
|
|
|
>>> print FixedPoint(1.1e30)
|
|
1099999999999999993725589651456.00
|
|
>>> print FixedPoint("1.1e30")
|
|
1100000000000000000000000000000.00
|
|
>>>
|
|
|
|
The following Python operators and functions accept FixedPoints in the
|
|
expected ways:
|
|
|
|
binary + - * / % divmod
|
|
with auto-coercion of other types to FixedPoint.
|
|
+ - % divmod of FixedPoints are always exact.
|
|
* / of FixedPoints may lose information to rounding, in
|
|
which case the result is the infinitely precise answer
|
|
rounded to the result's precision.
|
|
divmod(x, y) returns (q, r) where q is a long equal to
|
|
floor(x/y) as if x/y were computed to infinite precision,
|
|
and r is a FixedPoint equal to x - q * y; no information
|
|
is lost. Note that q has the sign of y, and abs(r) < abs(y).
|
|
unary -
|
|
== != < > <= >= cmp
|
|
min max
|
|
float int long (int and long truncate)
|
|
abs
|
|
str repr
|
|
hash
|
|
use as dict keys
|
|
use as boolean (e.g. "if some_FixedPoint:" -- true iff not zero)
|
|
|
|
Methods unique to FixedPoints:
|
|
.copy() return new FixedPoint with same value
|
|
.frac() long(x) + x.frac() == x
|
|
.get_precision() return the precision(p) of this FixedPoint object
|
|
.set_precision(p) set the precision of this FixedPoint object
|
|
|
|
Provided as-is; use at your own risk; no warranty; no promises; enjoy!
|
|
"""
|
|
|
|
# Released to the public domain 28-Mar-2001,
|
|
# by Tim Peters (tim.one@home.com).
|
|
|
|
|
|
# 28-Mar-01 ver 0.0,4
|
|
# Use repr() instead of str() inside __str__, because str(long) changed
|
|
# since this was first written (used to produce trailing "L", doesn't
|
|
# now).
|
|
#
|
|
# 09-May-99 ver 0,0,3
|
|
# Repaired __sub__(FixedPoint, string); was blowing up.
|
|
# Much more careful conversion of float (now best possible).
|
|
# Implemented exact % and divmod.
|
|
#
|
|
# 14-Oct-98 ver 0,0,2
|
|
# Added int, long, frac. Beefed up docs. Removed DECIMAL_POINT
|
|
# and MINUS_SIGN globals to discourage bloating this class instead
|
|
# of writing formatting wrapper classes (or subclasses)
|
|
#
|
|
# 11-Oct-98 ver 0,0,1
|
|
# posted to c.l.py
|
|
|
|
__copyright__ = "Copyright (C) Python Software Foundation"
|
|
__author__ = "Tim Peters"
|
|
__version__ = 0, 1, 0
|
|
|
|
def bankersRounding(self, dividend, divisor, quotient, remainder):
|
|
"""
|
|
rounding via nearest-even
|
|
increment the quotient if
|
|
the remainder is more than half of the divisor
|
|
or the remainder is exactly half the divisor and the quotient is odd
|
|
"""
|
|
c = cmp(remainder << 1, divisor)
|
|
# c < 0 <-> remainder < divisor/2, etc
|
|
if c > 0 or (c == 0 and (quotient & 1) == 1):
|
|
quotient += 1
|
|
return quotient
|
|
|
|
def addHalfAndChop(self, dividend, divisor, quotient, remainder):
|
|
"""
|
|
the equivalent of 'add half and chop'
|
|
increment the quotient if
|
|
the remainder is greater than half of the divisor
|
|
or the remainder is exactly half the divisor and the quotient is >= 0
|
|
"""
|
|
c = cmp(remainder << 1, divisor)
|
|
# c < 0 <-> remainder < divisor/2, etc
|
|
if c > 0 or (c == 0 and quotient >= 0):
|
|
quotient += 1
|
|
return quotient
|
|
|
|
# 2002-10-20 dougfort - fake classes for pre 2.2 compatibility
|
|
try:
|
|
object
|
|
except NameError:
|
|
class object:
|
|
pass
|
|
def property(x, y):
|
|
return None
|
|
|
|
# The default value for the number of decimal digits carried after the
|
|
# decimal point. This only has effect at compile-time.
|
|
DEFAULT_PRECISION = 2
|
|
|
|
class FixedPoint(object):
|
|
"""Basic FixedPoint object class,
|
|
The exact value is self.n / 10**self.p;
|
|
self.n is a long; self.p is an int
|
|
"""
|
|
__slots__ = ['n', 'p']
|
|
def __init__(self, value=0, precision=DEFAULT_PRECISION):
|
|
self.n = self.p = 0
|
|
self.set_precision(precision)
|
|
p = self.p
|
|
|
|
if isinstance(value, type("42.3e5")):
|
|
n, exp = _string2exact(value)
|
|
# exact value is n*10**exp = n*10**(exp+p)/10**p
|
|
effective_exp = exp + p
|
|
if effective_exp > 0:
|
|
n = n * _tento(effective_exp)
|
|
elif effective_exp < 0:
|
|
n = self._roundquotient(n, _tento(-effective_exp))
|
|
self.n = n
|
|
return
|
|
|
|
if isinstance(value, type(42)) or isinstance(value, type(42L)):
|
|
self.n = long(value) * _tento(p)
|
|
return
|
|
|
|
if isinstance(value, type(self)):
|
|
temp = value.copy()
|
|
temp.set_precision(p)
|
|
self.n, self.p = temp.n, temp.p
|
|
return
|
|
|
|
if isinstance(value, type(42.0)):
|
|
# XXX ignoring infinities and NaNs and overflows for now
|
|
import math
|
|
f, e = math.frexp(abs(value))
|
|
assert f == 0 or 0.5 <= f < 1.0
|
|
# |value| = 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: |value| = (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
|
|
|
|
# now |value| = top * 2**e exactly
|
|
# want n such that n / 10**p = top * 2**e, or
|
|
# n = top * 10**p * 2**e
|
|
top = top * _tento(p)
|
|
if e >= 0:
|
|
n = top << e
|
|
else:
|
|
n = self._roundquotient(top, 1L << -e)
|
|
if value < 0:
|
|
n = -n
|
|
self.n = n
|
|
return
|
|
|
|
if isinstance(value, type(42-42j)):
|
|
raise TypeError("can't convert complex to FixedPoint: " +
|
|
`value`)
|
|
|
|
# can we coerce to a float?
|
|
yes = 1
|
|
try:
|
|
asfloat = float(value)
|
|
except:
|
|
yes = 0
|
|
if yes:
|
|
self.__init__(asfloat, p)
|
|
return
|
|
|
|
# similarly for long
|
|
yes = 1
|
|
try:
|
|
aslong = long(value)
|
|
except:
|
|
yes = 0
|
|
if yes:
|
|
self.__init__(aslong, p)
|
|
return
|
|
|
|
raise TypeError("can't convert to FixedPoint: " + `value`)
|
|
|
|
def get_precision(self):
|
|
"""Return the precision of this FixedPoint.
|
|
|
|
The precision is the number of decimal digits carried after
|
|
the decimal point, and is an int >= 0.
|
|
"""
|
|
|
|
return self.p
|
|
|
|
def set_precision(self, precision=DEFAULT_PRECISION):
|
|
"""Change the precision carried by this FixedPoint to p.
|
|
|
|
precision must be an int >= 0, and defaults to
|
|
DEFAULT_PRECISION.
|
|
|
|
If precision is less than this FixedPoint's current precision,
|
|
information may be lost to rounding.
|
|
"""
|
|
|
|
try:
|
|
p = int(precision)
|
|
except:
|
|
raise TypeError("precision not convertable to int: " +
|
|
`precision`)
|
|
if p < 0:
|
|
raise ValueError("precision must be >= 0: " + `precision`)
|
|
|
|
if p > self.p:
|
|
self.n = self.n * _tento(p - self.p)
|
|
elif p < self.p:
|
|
self.n = self._roundquotient(self.n, _tento(self.p - p))
|
|
self.p = p
|
|
|
|
precision = property(get_precision, set_precision)
|
|
|
|
def __str__(self):
|
|
n, p = self.n, self.p
|
|
i, f = divmod(abs(n), _tento(p))
|
|
if p:
|
|
frac = repr(f)[:-1]
|
|
frac = "0" * (p - len(frac)) + frac
|
|
else:
|
|
frac = ""
|
|
return "-"[:n<0] + \
|
|
repr(i)[:-1] + \
|
|
"." + frac
|
|
|
|
def __repr__(self):
|
|
return "FixedPoint" + `(str(self), self.p)`
|
|
|
|
def copy(self):
|
|
return _mkFP(self.n, self.p, type(self))
|
|
|
|
__copy__ = copy
|
|
|
|
def __deepcopy__(self, memo):
|
|
return self.copy()
|
|
|
|
def __cmp__(self, other):
|
|
xn, yn, p = _norm(self, other, FixedPoint=type(self))
|
|
return cmp(xn, yn)
|
|
|
|
def __hash__(self):
|
|
""" Caution! == values must have equal hashes, and a FixedPoint
|
|
is essentially a rational in unnormalized form. There's
|
|
really no choice here but to normalize it, so hash is
|
|
potentially expensive.
|
|
n, p = self.__reduce()
|
|
|
|
Obscurity: if the value is an exact integer, p will be 0 now,
|
|
so the hash expression reduces to hash(n). So FixedPoints
|
|
that happen to be exact integers hash to the same things as
|
|
their int or long equivalents. This is Good. But if a
|
|
FixedPoint happens to have a value exactly representable as
|
|
a float, their hashes may differ. This is a teensy bit Bad.
|
|
"""
|
|
n, p = self.__reduce()
|
|
return hash(n) ^ hash(p)
|
|
|
|
def __nonzero__(self):
|
|
""" Returns true if this FixedPoint is not equal to zero"""
|
|
return self.n != 0
|
|
|
|
def __neg__(self):
|
|
return _mkFP(-self.n, self.p, type(self))
|
|
|
|
def __abs__(self):
|
|
""" Returns new FixedPoint containing the absolute value of this FixedPoint"""
|
|
if self.n >= 0:
|
|
return self.copy()
|
|
else:
|
|
return -self
|
|
|
|
def __add__(self, other):
|
|
n1, n2, p = _norm(self, other, FixedPoint=type(self))
|
|
# n1/10**p + n2/10**p = (n1+n2)/10**p
|
|
return _mkFP(n1 + n2, p, type(self))
|
|
|
|
__radd__ = __add__
|
|
|
|
def __sub__(self, other):
|
|
if not isinstance(other, type(self)):
|
|
other = type(self)(other, self.p)
|
|
return self.__add__(-other)
|
|
|
|
def __rsub__(self, other):
|
|
return (-self) + other
|
|
|
|
def __mul__(self, other):
|
|
n1, n2, p = _norm(self, other, FixedPoint=type(self))
|
|
# n1/10**p * n2/10**p = (n1*n2/10**p)/10**p
|
|
return _mkFP(self._roundquotient(n1 * n2, _tento(p)), p, type(self))
|
|
|
|
__rmul__ = __mul__
|
|
|
|
def __div__(self, other):
|
|
n1, n2, p = _norm(self, other, FixedPoint=type(self))
|
|
if n2 == 0:
|
|
raise ZeroDivisionError("FixedPoint division")
|
|
if n2 < 0:
|
|
n1, n2 = -n1, -n2
|
|
# n1/10**p / (n2/10**p) = n1/n2 = (n1*10**p/n2)/10**p
|
|
return _mkFP(self._roundquotient(n1 * _tento(p), n2), p, type(self))
|
|
|
|
def __rdiv__(self, other):
|
|
n1, n2, p = _norm(self, other, FixedPoint=type(self))
|
|
return _mkFP(n2, p, FixedPoint=type(self)) / self
|
|
|
|
def __divmod__(self, other):
|
|
n1, n2, p = _norm(self, other, FixedPoint=type(self))
|
|
if n2 == 0:
|
|
raise ZeroDivisionError("FixedPoint modulo")
|
|
# floor((n1/10**p)/(n2*10**p)) = floor(n1/n2)
|
|
q = n1 / n2
|
|
# n1/10**p - q * n2/10**p = (n1 - q * n2)/10**p
|
|
return q, _mkFP(n1 - q * n2, p, type(self))
|
|
|
|
def __rdivmod__(self, other):
|
|
n1, n2, p = _norm(self, other, FixedPoint=type(self))
|
|
return divmod(_mkFP(n2, p), self)
|
|
|
|
def __mod__(self, other):
|
|
return self.__divmod__(other)[1]
|
|
|
|
def __rmod__(self, other):
|
|
n1, n2, p = _norm(self, other, FixedPoint=type(self))
|
|
return _mkFP(n2, p, type(self)).__mod__(self)
|
|
|
|
def __float__(self):
|
|
"""Return the floating point representation of this FixedPoint.
|
|
Caution! float can lose precision.
|
|
"""
|
|
n, p = self.__reduce()
|
|
return float(n) / float(_tento(p))
|
|
|
|
def __long__(self):
|
|
"""EJG/DF - Should this round instead?
|
|
Note e.g. long(-1.9) == -1L and long(1.9) == 1L in Python
|
|
Note that __int__ inherits whatever __long__ does,
|
|
and .frac() is affected too
|
|
"""
|
|
answer = abs(self.n) / _tento(self.p)
|
|
if self.n < 0:
|
|
answer = -answer
|
|
return answer
|
|
|
|
def __int__(self):
|
|
"""Return integer value of FixedPoint object."""
|
|
return int(self.__long__())
|
|
|
|
def frac(self):
|
|
"""Return fractional portion as a FixedPoint.
|
|
|
|
x.frac() + long(x) == x
|
|
"""
|
|
return self - long(self)
|
|
|
|
def _roundquotient(self, x, y):
|
|
"""
|
|
Divide x by y,
|
|
return the result of rounding
|
|
Developers may substitute their own 'round' for custom rounding
|
|
y must be > 0
|
|
"""
|
|
assert y > 0
|
|
n, leftover = divmod(x, y)
|
|
return self.round(x, y, n, leftover)
|
|
|
|
def __reduce(self):
|
|
""" Return n, p s.t. self == n/10**p and n % 10 != 0"""
|
|
n, p = self.n, self.p
|
|
if n == 0:
|
|
p = 0
|
|
while p and n % 10 == 0:
|
|
p = p - 1
|
|
n = n / 10
|
|
return n, p
|
|
|
|
# 2002-10-04 dougfort - Default to Banker's Rounding for backward compatibility
|
|
FixedPoint.round = bankersRounding
|
|
|
|
# return 10L**n
|
|
|
|
def _tento(n, cache={}):
|
|
"""Cached computation of 10**n"""
|
|
try:
|
|
return cache[n]
|
|
except KeyError:
|
|
answer = cache[n] = 10L ** n
|
|
return answer
|
|
|
|
def _norm(x, y, isinstance=isinstance, FixedPoint=FixedPoint,
|
|
_tento=_tento):
|
|
"""Return xn, yn, p s.t.
|
|
p = max(x.p, y.p)
|
|
x = xn / 10**p
|
|
y = yn / 10**p
|
|
|
|
x must be FixedPoint to begin with; if y is not FixedPoint,
|
|
it inherits its precision from x.
|
|
|
|
Note that this method is called a lot, so default-arg tricks are helpful.
|
|
"""
|
|
assert isinstance(x, FixedPoint)
|
|
if not isinstance(y, FixedPoint):
|
|
y = FixedPoint(y, x.p)
|
|
xn, yn = x.n, y.n
|
|
xp, yp = x.p, y.p
|
|
if xp > yp:
|
|
yn = yn * _tento(xp - yp)
|
|
p = xp
|
|
elif xp < yp:
|
|
xn = xn * _tento(yp - xp)
|
|
p = yp
|
|
else:
|
|
p = xp # same as yp
|
|
return xn, yn, p
|
|
|
|
def _mkFP(n, p, FixedPoint=FixedPoint):
|
|
"""Make FixedPoint objext - Return a new FixedPoint object with the selected precision."""
|
|
f = FixedPoint()
|
|
#print '_mkFP Debug: %s, value=%s' % (type(f),n)
|
|
f.n = n
|
|
f.p = p
|
|
return f
|
|
|
|
# crud for parsing strings
|
|
import re
|
|
|
|
# There's an optional sign at the start, and an optional exponent
|
|
# at the end. The exponent has an optional sign and at least one
|
|
# digit. In between, must have either at least one digit followed
|
|
# by an optional fraction, or a decimal point followed by at least
|
|
# one digit. Yuck.
|
|
|
|
_parser = re.compile(r"""
|
|
\s*
|
|
(?P<sign>[-+])?
|
|
(
|
|
(?P<int>\d+) (\. (?P<frac>\d*))?
|
|
|
|
|
\. (?P<onlyfrac>\d+)
|
|
)
|
|
([eE](?P<exp>[-+]? \d+))?
|
|
\s* $
|
|
""", re.VERBOSE).match
|
|
|
|
del re
|
|
|
|
|
|
def _string2exact(s):
|
|
"""Return n, p s.t. float string value == n * 10**p exactly."""
|
|
m = _parser(s)
|
|
if m is None:
|
|
raise ValueError("can't parse as number: " + `s`)
|
|
|
|
exp = m.group('exp')
|
|
if exp is None:
|
|
exp = 0
|
|
else:
|
|
exp = int(exp)
|
|
|
|
intpart = m.group('int')
|
|
if intpart is None:
|
|
intpart = "0"
|
|
fracpart = m.group('onlyfrac')
|
|
else:
|
|
fracpart = m.group('frac')
|
|
if fracpart is None or fracpart == "":
|
|
fracpart = "0"
|
|
assert intpart
|
|
assert fracpart
|
|
|
|
i, f = long(intpart), long(fracpart)
|
|
nfrac = len(fracpart)
|
|
i = i * _tento(nfrac) + f
|
|
exp = exp - nfrac
|
|
|
|
if m.group('sign') == "-":
|
|
i = -i
|
|
|
|
return i, exp
|
|
|
|
def _test():
|
|
"""Unit testing framework"""
|
|
fp = FixedPoint
|
|
o = fp("0.1")
|
|
assert str(o) == "0.10"
|
|
t = fp("-20e-2", 5)
|
|
assert str(t) == "-0.20000"
|
|
assert t < o
|
|
assert o > t
|
|
assert min(o, t) == min(t, o) == t
|
|
assert max(o, t) == max(t, o) == o
|
|
assert o != t
|
|
assert --t == t
|
|
assert abs(t) > abs(o)
|
|
assert abs(o) < abs(t)
|
|
assert o == o and t == t
|
|
assert t.copy() == t
|
|
assert o == -t/2 == -.5 * t
|
|
assert abs(t) == o + o
|
|
assert abs(o) == o
|
|
assert o/t == -0.5
|
|
assert -(t/o) == (-t)/o == t/-o == 2
|
|
assert 1 + o == o + 1 == fp(" +00.000011e+5 ")
|
|
assert 1/o == 10
|
|
assert o + t == t + o == -o
|
|
assert 2.0 * t == t * 2 == "2" * t == o/o * 2L * t
|
|
assert 1 - t == -(t - 1) == fp(6L)/5
|
|
assert t*t == 4*o*o == o*4*o == o*o*4
|
|
assert fp(2) - "1" == 1
|
|
assert float(-1/t) == 5.0
|
|
for p in range(20):
|
|
assert 42 + fp("1e-20", p) - 42 == 0
|
|
assert 1/(42 + fp("1e-20", 20) - 42) == fp("100.0E18")
|
|
o = fp(".9995", 4)
|
|
assert 1 - o == fp("5e-4", 10)
|
|
o.set_precision(3)
|
|
assert o == 1
|
|
o = fp(".9985", 4)
|
|
o.set_precision(3)
|
|
assert o == fp(".998", 10)
|
|
assert o == o.frac()
|
|
o.set_precision(100)
|
|
assert o == fp(".998", 10)
|
|
o.set_precision(2)
|
|
assert o == 1
|
|
x = fp(1.99)
|
|
assert long(x) == -long(-x) == 1L
|
|
assert int(x) == -int(-x) == 1
|
|
assert x == long(x) + x.frac()
|
|
assert -x == long(-x) + (-x).frac()
|
|
assert fp(7) % 4 == 7 % fp(4) == 3
|
|
assert fp(-7) % 4 == -7 % fp(4) == 1
|
|
assert fp(-7) % -4 == -7 % fp(-4) == -3
|
|
assert fp(7.0) % "-4.0" == 7 % fp(-4) == -1
|
|
assert fp("5.5") % fp("1.1") == fp("5.5e100") % fp("1.1e100") == 0
|
|
assert divmod(fp("1e100"), 3) == (long(fp("1e100")/3), 1)
|
|
|
|
if __name__ == '__main__':
|
|
_test()
|
|
|