
[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