# Source code for topo.misc.fixedpoint

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

__author__ = "Tim Peters"
__version__ = 0, 1, 0

[docs]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

[docs]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 instance initialization.
DEFAULT_PRECISION = 2

[docs]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=None):
self.n = self.p = 0
if precision == None:
precision = DEFAULT_PRECISION

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)

[docs]    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

[docs]    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()

# Basic support for pickling
def __getstate__(self):
state={}
try:
for k in self.__slots__:
state[k] = getattr(self,k)
except AttributeError:
pass
return state

def __setstate__(self,state):
for k,v in state.items():
setattr(self,k,v)
self.unpickle()

def unpickle(self):
pass

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

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

def __sub__(self, other):
if not isinstance(other, type(self)):
other = type(self)(other, self.p)

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
"""
if self.n < 0:

def __int__(self):
"""Return integer value of FixedPoint object."""
return int(self.__long__())

[docs]    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

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()