[Python-Dev] Fixed-decimal types

[Tim Peters]

[Guido]
> ...
> Not arguing for this interpretation, just indicating that doing
> fixed precision arithmetic right is hard.

It's not so much hard as it is arbitrary.  The floating-point world is
standardized now, but the fixed-point world remains a mish-mash of
incompatible legacy schemes carried across generations of products for no
reason other than product-specific compatibility.  So despite that
fixed-point has a specialty audience, whatever rules Python chooses will
leave it incompatible with much of that audience's (mixed!) expectations.

If fixed-point is needed, and my FixedPoint.py isn't good enough (all other
fixed point pkgs I've seen for Python were braindead), then it should be
implemented such that developers can control both rounding and precision
propagation.  I'll attach suitable kernels; they haven't been tested but any
bugs discovered will be trivial to fix (there are no difficulties here, but
typos are likely); the kernels supply the bulk of what's required, whether
implemented in Python or C; various packages can wrap them to supply
whatever policies they like; see FixedPoint.py for exact string<->FixedPoint
and exact float->FixedPoint conversions; and that's the end of my
involvement in fixed-point <wink>.

Python should certainly *not* add a "scale factor" to its current long
implementation; fixed-point should be a distinct type, as scale-factor
fiddling is clumsy and pervasive (long arithmetic is challenging enough to
get correct and quick without this obfuscating distraction; and by leaving
scale factors out of it, it's much easier to plug in alternative bigint
implementations (like GMP)).

One other point:  some people are going to want BCD (binary-coded decimal),
which suffers the same mish-mash of legacy policies, but with a different
data representation.  The point is that many commercial applications spend
much more time doing I/O conversions than arithmetic, and BCD accepts slow
arithmetic (in the absence of special HW support) in return for fast scaling
& I/O conversion.

Forgetting the database-heads for a moment, decimal *floating*-point is what
calculators do, so that's what "real people" are most comfortable with.  The
IEEE-854 std (IEEE-754's younger and friendlier brother) specifies that
completely.  Add a means to boost "global" precision (a la REXX), and it's a
powerful tool even for experts (benefits approximating those of unbounded
rational arithmetic but with bounded & user-controllable expense).

can-never-have-too-many-numeric-types-but-always-have-
    too-few-literal-notations-ly y'rs  - tim


# Kernels for fixed-point decimal arithmetic.

# _add, _sub, _mul, _div all have arglist
#     n1, p1, n2, p2, p, round=DEFAULT_ROUND
# n1 and n2 are longs; p1, p2 and p ints >= 0.
# The inputs are exactly n1/10**p1 and n2/10**p2.
#
# The return value is the integer n such that n/10**p is the best
# approximation to the infinite-precision result.  In other words, p1
# and p2 are the input precisions and p is the desired output
# precision, where precision is the # of digits *after* the decimal
# point.
#
# What "best approximation" means is determined by the round function.
# In many cases rounding isn't required, but when it is
#     round(top, bot)
# is returned.  top and bot are longs, with bot > 0 guaranteed.  The
# infinite-precision result is top/bot.  round must return an integer
# (long) approximation to top/bot, using whichever rounding discipline
# you want.  By default, IEEE round-to-nearest/even is used; see the
# _roundXXX functions for examples of suitable rounding functions.
#
# Note:  The only code here that knows we're working in decimal is
# function _tento; simply change the "10L" in that to do fixed-point
# arithmetic in some other base.
#
# Example:
#
# >>> r7 = _div(1L, 0, 7L, 0, 20)  # 1/7
# >>> r7
# 14285714285714285714L
# >>> r5 = _div(1L, 0, 5L, 0, 20)  # 1/5
# >>> r5
# 20000000000000000000L
# >>> sum = _add(r7, 20, r5, 20, 20)  # 1/7 + 1/5 = 12/35
# >>> sum
# 34285714285714285714L
# >>> _mul(sum, 20, 35L, 0, 20)
# 1199999999999999999990L
# >>> _mul(sum, 20, 35L, 0, 18)
# 12000000000000000000L
# >>> _mul(sum, 20, 35L, 0, 0)
# 12L
# >>>

###################################################################
# Sample rounding functions.
###################################################################

# Round to minus infinity.

def _roundminf(top, bot):
    assert bot > 0
    return top / bot

# Round to plus infinity.

def _roundpinf(top, bot):
    assert bot > 0
    q, r = divmod(top, bot)
    # answer is exactly q + r/bot; and 0 <= r < bot
    if r:
        q = q + 1
    return q

# IEEE nearest/even rounding (closest integer; in case of tie closest
# even integer).

def _roundne(top, bot):
    assert bot > 0
    q, r = divmod(top, bot)
    # answer is exactly q + r/bot; and 0 <= r < bot
    c = cmp(r << 1, bot)
    # c < 0 <-> r < bot/2, etc
    if c > 0 or (c == 0 and (q & 1) == 1):
        q = q + 1
    return q

# "Add a half and chop" rounding (remainder < 1/2 toward 0; remainder
# >= half away from 0).

def _roundhalf(top, bot):
    assert bot > 0
    q, r = divmod(top, bot)
    # answer is exactly q + r/bot; and 0 <= r < bot
    c = cmp(r << 1, bot)
    # c < 0 <-> r < bot/2, etc
    if c > 0 or (c == 0 and q >= 0):
        q = q + 1
    return q

# Round toward 0 (throw away remainder).

def _roundchop(top, bot):
    assert bot > 0
    q, r = divmod(top, bot)
    # answer is exactly q + r/bot; and 0 <= r < bot
    if r and q < 0:
        q = q + 1
    return q

###################################################################
# Kernels for + - * /.
###################################################################

DEFAULT_ROUND = _roundne

def _add(n1, p1, n2, p2, p, round=DEFAULT_ROUND):
    assert p1 >= 0
    assert p2 >= 0
    assert p >= 0
    # (n1/10**p1 + n2/10**p2) * 10**p ==
    # (n1*10**(max-p1) + n2*10**(max-p2))/10**max * 10**p
    max = p1    # until proven otherwise
    if p1 < p2:
        n1 = n1 * _tento(p2 - p1)
        max = p2
    elif p2 < p1:
        n2 = n2 * _tento(p1 - p2)
    n3 = n1 + n2
    p3 = p - max
    if p3 > 0:
        n3 = n3 * _tento(p3)
    elif p3 < 0:
        n3 = round(n3, _tento(-p3))
    return n3

def _sub(n1, p1, n2, p2, p, round=DEFAULT_ROUND):
    assert p1 >= 0
    assert p2 >= 0
    assert p >= 0
    return _add(n1, p1, -n2, p2, p, round)

def _mul(n1, p1, n2, p2, p, round=DEFAULT_ROUND):
    assert p1 >= 0
    assert p2 >= 0
    assert p >= 0
    # (n1/10**p1 * n2/10**p2) * 10**p ==
    # (n1*n2)/10**(p1+p2) * 10**p
    n3 = n1 * n2
    p3 = p - p1 - p2
    if p3 > 0:
        n3 = n3 * _tento(p3)
    elif p3 < 0:
        n3 = round(n3, _tento(-p3))
    return n3

def _div(n1, p1, n2, p2, p, round=DEFAULT_ROUND):
    assert p1 >= 0
    assert p2 >= 0
    assert p >= 0
    if n2 == 0:
        raise ZeroDivisionError("scaled integer")
    # (n1/10**p1 / n2/10**p2) * 10**p ==
    # (n1/n2) * 10**(p2-p1+p)
    p3 = p2 - p1 + p
    if p3 > 0:
        n1 = n1 * _tento(p3)
    elif p3 < 0:
        n2 = n2 * _tento(-p3)
    if n2 < 0:
        n1 = -n1
        n2 = -n2
    return round(n1, n2)

def _tento(i, _cache={}):
    assert i >= 0
    try:
        return _cache[i]
    except KeyError:
        answer = _cache[i] = 10L ** i
        return answer