Package topo :: Package misc :: Module fixedpoint
[hide private]
[frames] | no frames]

Source Code for Module topo.misc.fixedpoint

  1  """ 
  2  FixedPoint objects support decimal arithmetic with a fixed number of 
  3  digits (called the object's precision) after the decimal point.  The 
  4  number of digits before the decimal point is variable & unbounded. 
  5   
  6  The precision is user-settable on a per-object basis when a FixedPoint 
  7  is constructed, and may vary across FixedPoint objects.  The precision 
  8  may also be changed after construction via FixedPoint.set_precision(p). 
  9  Note that if the precision of a FixedPoint is reduced via set_precision, 
 10  information may be lost to rounding. 
 11   
 12  >>> x = FixedPoint("5.55")  # precision defaults to 2 
 13  >>> print x 
 14  5.55 
 15  >>> x.set_precision(1)      # round to one fraction digit 
 16  >>> print x 
 17  5.6 
 18  >>> print FixedPoint("5.55", 1)  # same thing setting to 1 in constructor 
 19  5.6 
 20  >>> repr(x) #  returns constructor string that reproduces object exactly 
 21  "FixedPoint('5.6', 1)" 
 22  >>> 
 23   
 24  When FixedPoint objects of different precision are combined via + - * /, 
 25  the result is computed to the larger of the inputs' precisions, which also 
 26  becomes the precision of the resulting FixedPoint object. 
 27   
 28  >>> print FixedPoint("3.42") + FixedPoint("100.005", 3) 
 29  103.425 
 30  >>> 
 31   
 32  When a FixedPoint is combined with other numeric types (ints, floats, 
 33  strings representing a number) via + - * /, then similarly the computation 
 34  is carried out using-- and the result inherits --the FixedPoint's 
 35  precision. 
 36   
 37  >>> print FixedPoint(1) / 7 
 38  0.14 
 39  >>> print FixedPoint(1, 30) / 7 
 40  0.142857142857142857142857142857 
 41  >>> 
 42   
 43  The string produced by str(x) (implictly invoked by "print") always 
 44  contains at least one digit before the decimal point, followed by a 
 45  decimal point, followed by exactly x.get_precision() digits.  If x is 
 46  negative, str(x)[0] == "-". 
 47   
 48  The FixedPoint constructor can be passed an int, long, string, float, 
 49  FixedPoint, or any object convertible to a float via float() or to a 
 50  long via long().  Passing a precision is optional; if specified, the 
 51  precision must be a non-negative int.  There is no inherent limit on 
 52  the size of the precision, but if very very large you'll probably run 
 53  out of memory. 
 54   
 55  Note that conversion of floats to FixedPoint can be surprising, and 
 56  should be avoided whenever possible.  Conversion from string is exact 
 57  (up to final rounding to the requested precision), so is greatly 
 58  preferred. 
 59   
 60  >>> print FixedPoint(1.1e30) 
 61  1099999999999999993725589651456.00 
 62  >>> print FixedPoint("1.1e30") 
 63  1100000000000000000000000000000.00 
 64  >>> 
 65   
 66  The following Python operators and functions accept FixedPoints in the 
 67  expected ways: 
 68   
 69      binary + - * / % divmod 
 70          with auto-coercion of other types to FixedPoint. 
 71          + - % divmod  of FixedPoints are always exact. 
 72          * / of FixedPoints may lose information to rounding, in 
 73              which case the result is the infinitely precise answer 
 74              rounded to the result's precision. 
 75          divmod(x, y) returns (q, r) where q is a long equal to 
 76              floor(x/y) as if x/y were computed to infinite precision, 
 77              and r is a FixedPoint equal to x - q * y; no information 
 78              is lost.  Note that q has the sign of y, and abs(r) < abs(y). 
 79      unary - 
 80      == != < > <= >=  cmp 
 81      min  max 
 82      float  int  long    (int and long truncate) 
 83      abs 
 84      str  repr 
 85      hash 
 86      use as dict keys 
 87      use as boolean (e.g. "if some_FixedPoint:" -- true iff not zero) 
 88   
 89  Methods unique to FixedPoints: 
 90     .copy()              return new FixedPoint with same value 
 91     .frac()              long(x) + x.frac() == x 
 92     .get_precision()     return the precision(p) of this FixedPoint object 
 93     .set_precision(p)    set the precision of this FixedPoint object 
 94      
 95  Provided as-is; use at your own risk; no warranty; no promises; enjoy! 
 96  """ 
 97   
 98  # Released to the public domain 28-Mar-2001, 
 99  # by Tim Peters (tim.one@home.com). 
100   
101   
102  # 28-Mar-01 ver 0.0,4 
103  #     Use repr() instead of str() inside __str__, because str(long) changed 
104  #     since this was first written (used to produce trailing "L", doesn't 
105  #     now). 
106  # 
107  # 09-May-99 ver 0,0,3 
108  #     Repaired __sub__(FixedPoint, string); was blowing up. 
109  #     Much more careful conversion of float (now best possible). 
110  #     Implemented exact % and divmod. 
111  # 
112  # 14-Oct-98 ver 0,0,2 
113  #     Added int, long, frac.  Beefed up docs.  Removed DECIMAL_POINT 
114  #     and MINUS_SIGN globals to discourage bloating this class instead 
115  #     of writing formatting wrapper classes (or subclasses) 
116  # 
117  # 11-Oct-98 ver 0,0,1 
118  #     posted to c.l.py 
119   
120  __copyright__ = "Copyright (C) Python Software Foundation" 
121  __author__ = "Tim Peters" 
122  __version__ = 0, 1, 0 
123   
124 -def bankersRounding(self, dividend, divisor, quotient, remainder):
125 """ 126 rounding via nearest-even 127 increment the quotient if 128 the remainder is more than half of the divisor 129 or the remainder is exactly half the divisor and the quotient is odd 130 """ 131 c = cmp(remainder << 1, divisor) 132 # c < 0 <-> remainder < divisor/2, etc 133 if c > 0 or (c == 0 and (quotient & 1) == 1): 134 quotient += 1 135 return quotient
136
137 -def addHalfAndChop(self, dividend, divisor, quotient, remainder):
138 """ 139 the equivalent of 'add half and chop' 140 increment the quotient if 141 the remainder is greater than half of the divisor 142 or the remainder is exactly half the divisor and the quotient is >= 0 143 """ 144 c = cmp(remainder << 1, divisor) 145 # c < 0 <-> remainder < divisor/2, etc 146 if c > 0 or (c == 0 and quotient >= 0): 147 quotient += 1 148 return quotient
149 150 # 2002-10-20 dougfort - fake classes for pre 2.2 compatibility 151 try: 152 object 153 except NameError:
154 - class object:
155 pass
156 - def property(x, y):
157 return None
158 159 # The default value for the number of decimal digits carried after the 160 # decimal point. This only has effect at instance initialization. 161 DEFAULT_PRECISION = 2 162
163 -class FixedPoint(object):
164 """Basic FixedPoint object class, 165 The exact value is self.n / 10**self.p; 166 self.n is a long; self.p is an int 167 """ 168 __slots__ = ['n', 'p'] 169
170 - def __init__(self, value=0, precision=None):
171 self.n = self.p = 0 172 if precision == None: 173 precision = DEFAULT_PRECISION 174 175 self.set_precision(precision) 176 p = self.p 177 178 if isinstance(value, type("42.3e5")): 179 n, exp = _string2exact(value) 180 # exact value is n*10**exp = n*10**(exp+p)/10**p 181 effective_exp = exp + p 182 if effective_exp > 0: 183 n = n * _tento(effective_exp) 184 elif effective_exp < 0: 185 n = self._roundquotient(n, _tento(-effective_exp)) 186 self.n = n 187 return 188 189 if isinstance(value, type(42)) or isinstance(value, type(42L)): 190 self.n = long(value) * _tento(p) 191 return 192 193 if isinstance(value, type(self)): 194 temp = value.copy() 195 temp.set_precision(p) 196 self.n, self.p = temp.n, temp.p 197 return 198 199 if isinstance(value, type(42.0)): 200 # XXX ignoring infinities and NaNs and overflows for now 201 import math 202 f, e = math.frexp(abs(value)) 203 assert f == 0 or 0.5 <= f < 1.0 204 # |value| = f * 2**e exactly 205 206 # Suck up CHUNK bits at a time; 28 is enough so that we suck 207 # up all bits in 2 iterations for all known binary double- 208 # precision formats, and small enough to fit in an int. 209 CHUNK = 28 210 top = 0L 211 # invariant: |value| = (top + f) * 2**e exactly 212 while f: 213 f = math.ldexp(f, CHUNK) 214 digit = int(f) 215 assert digit >> CHUNK == 0 216 top = (top << CHUNK) | digit 217 f = f - digit 218 assert 0.0 <= f < 1.0 219 e = e - CHUNK 220 221 # now |value| = top * 2**e exactly 222 # want n such that n / 10**p = top * 2**e, or 223 # n = top * 10**p * 2**e 224 top = top * _tento(p) 225 if e >= 0: 226 n = top << e 227 else: 228 n = self._roundquotient(top, 1L << -e) 229 if value < 0: 230 n = -n 231 self.n = n 232 return 233 234 if isinstance(value, type(42-42j)): 235 raise TypeError("can't convert complex to FixedPoint: " + 236 `value`) 237 238 # can we coerce to a float? 239 yes = 1 240 try: 241 asfloat = float(value) 242 except: 243 yes = 0 244 if yes: 245 self.__init__(asfloat, p) 246 return 247 248 # similarly for long 249 yes = 1 250 try: 251 aslong = long(value) 252 except: 253 yes = 0 254 if yes: 255 self.__init__(aslong, p) 256 return 257 258 raise TypeError("can't convert to FixedPoint: " + `value`)
259
260 - def get_precision(self):
261 """Return the precision of this FixedPoint. 262 263 The precision is the number of decimal digits carried after 264 the decimal point, and is an int >= 0. 265 """ 266 267 return self.p
268
269 - def set_precision(self, precision=DEFAULT_PRECISION):
270 """Change the precision carried by this FixedPoint to p. 271 272 precision must be an int >= 0, and defaults to 273 DEFAULT_PRECISION. 274 275 If precision is less than this FixedPoint's current precision, 276 information may be lost to rounding. 277 """ 278 279 try: 280 p = int(precision) 281 except: 282 raise TypeError("precision not convertable to int: " + 283 `precision`) 284 if p < 0: 285 raise ValueError("precision must be >= 0: " + `precision`) 286 287 if p > self.p: 288 self.n = self.n * _tento(p - self.p) 289 elif p < self.p: 290 self.n = self._roundquotient(self.n, _tento(self.p - p)) 291 self.p = p
292 293 precision = property(get_precision, set_precision) 294
295 - def __str__(self):
296 n, p = self.n, self.p 297 i, f = divmod(abs(n), _tento(p)) 298 if p: 299 frac = repr(f)[:-1] 300 frac = "0" * (p - len(frac)) + frac 301 else: 302 frac = "" 303 return "-"[:n<0] + \ 304 repr(i)[:-1] + \ 305 "." + frac
306
307 - def __repr__(self):
308 return "FixedPoint" + `(str(self), self.p)`
309
310 - def copy(self):
311 return _mkFP(self.n, self.p, type(self))
312 313 __copy__ = copy 314
315 - def __deepcopy__(self, memo):
316 return self.copy()
317 318 # Basic support for pickling
319 - def __getstate__(self):
320 state={} 321 try: 322 for k in self.__slots__: 323 state[k] = getattr(self,k) 324 except AttributeError: 325 pass 326 return state
327
328 - def __setstate__(self,state):
329 for k,v in state.items(): 330 setattr(self,k,v) 331 self.unpickle()
332
333 - def unpickle(self):
334 pass
335 336
337 - def __cmp__(self, other):
338 xn, yn, p = _norm(self, other, FixedPoint=type(self)) 339 return cmp(xn, yn)
340
341 - def __hash__(self):
342 """ Caution! == values must have equal hashes, and a FixedPoint 343 is essentially a rational in unnormalized form. There's 344 really no choice here but to normalize it, so hash is 345 potentially expensive. 346 n, p = self.__reduce() 347 348 Obscurity: if the value is an exact integer, p will be 0 now, 349 so the hash expression reduces to hash(n). So FixedPoints 350 that happen to be exact integers hash to the same things as 351 their int or long equivalents. This is Good. But if a 352 FixedPoint happens to have a value exactly representable as 353 a float, their hashes may differ. This is a teensy bit Bad. 354 """ 355 n, p = self.__reduce() 356 return hash(n) ^ hash(p)
357
358 - def __nonzero__(self):
359 """ Returns true if this FixedPoint is not equal to zero""" 360 return self.n != 0
361
362 - def __neg__(self):
363 return _mkFP(-self.n, self.p, type(self))
364
365 - def __abs__(self):
366 """ Returns new FixedPoint containing the absolute value of this FixedPoint""" 367 if self.n >= 0: 368 return self.copy() 369 else: 370 return -self
371
372 - def __add__(self, other):
373 n1, n2, p = _norm(self, other, FixedPoint=type(self)) 374 # n1/10**p + n2/10**p = (n1+n2)/10**p 375 return _mkFP(n1 + n2, p, type(self))
376 377 __radd__ = __add__ 378
379 - def __sub__(self, other):
380 if not isinstance(other, type(self)): 381 other = type(self)(other, self.p) 382 return self.__add__(-other)
383
384 - def __rsub__(self, other):
385 return (-self) + other
386
387 - def __mul__(self, other):
388 n1, n2, p = _norm(self, other, FixedPoint=type(self)) 389 # n1/10**p * n2/10**p = (n1*n2/10**p)/10**p 390 return _mkFP(self._roundquotient(n1 * n2, _tento(p)), p, type(self))
391 392 __rmul__ = __mul__ 393
394 - def __div__(self, other):
395 n1, n2, p = _norm(self, other, FixedPoint=type(self)) 396 if n2 == 0: 397 raise ZeroDivisionError("FixedPoint division") 398 if n2 < 0: 399 n1, n2 = -n1, -n2 400 # n1/10**p / (n2/10**p) = n1/n2 = (n1*10**p/n2)/10**p 401 return _mkFP(self._roundquotient(n1 * _tento(p), n2), p, type(self))
402
403 - def __rdiv__(self, other):
404 n1, n2, p = _norm(self, other, FixedPoint=type(self)) 405 return _mkFP(n2, p, FixedPoint=type(self)) / self
406
407 - def __divmod__(self, other):
408 n1, n2, p = _norm(self, other, FixedPoint=type(self)) 409 if n2 == 0: 410 raise ZeroDivisionError("FixedPoint modulo") 411 # floor((n1/10**p)/(n2*10**p)) = floor(n1/n2) 412 q = n1 / n2 413 # n1/10**p - q * n2/10**p = (n1 - q * n2)/10**p 414 return q, _mkFP(n1 - q * n2, p, type(self))
415
416 - def __rdivmod__(self, other):
417 n1, n2, p = _norm(self, other, FixedPoint=type(self)) 418 return divmod(_mkFP(n2, p), self)
419
420 - def __mod__(self, other):
421 return self.__divmod__(other)[1]
422
423 - def __rmod__(self, other):
424 n1, n2, p = _norm(self, other, FixedPoint=type(self)) 425 return _mkFP(n2, p, type(self)).__mod__(self)
426
427 - def __float__(self):
428 """Return the floating point representation of this FixedPoint. 429 Caution! float can lose precision. 430 """ 431 n, p = self.__reduce() 432 return float(n) / float(_tento(p))
433
434 - def __long__(self):
435 """EJG/DF - Should this round instead? 436 Note e.g. long(-1.9) == -1L and long(1.9) == 1L in Python 437 Note that __int__ inherits whatever __long__ does, 438 and .frac() is affected too 439 """ 440 answer = abs(self.n) / _tento(self.p) 441 if self.n < 0: 442 answer = -answer 443 return answer
444
445 - def __int__(self):
446 """Return integer value of FixedPoint object.""" 447 return int(self.__long__())
448
449 - def frac(self):
450 """Return fractional portion as a FixedPoint. 451 452 x.frac() + long(x) == x 453 """ 454 return self - long(self)
455
456 - def _roundquotient(self, x, y):
457 """ 458 Divide x by y, 459 return the result of rounding 460 Developers may substitute their own 'round' for custom rounding 461 y must be > 0 462 """ 463 assert y > 0 464 n, leftover = divmod(x, y) 465 return self.round(x, y, n, leftover)
466
467 - def __reduce(self):
468 """ Return n, p s.t. self == n/10**p and n % 10 != 0""" 469 n, p = self.n, self.p 470 if n == 0: 471 p = 0 472 while p and n % 10 == 0: 473 p = p - 1 474 n = n / 10 475 return n, p
476 477 # 2002-10-04 dougfort - Default to Banker's Rounding for backward compatibility 478 FixedPoint.round = bankersRounding 479 480 # return 10L**n 481
482 -def _tento(n, cache={}):
483 """Cached computation of 10**n""" 484 try: 485 return cache[n] 486 except KeyError: 487 answer = cache[n] = 10L ** n 488 return answer
489
490 -def _norm(x, y, isinstance=isinstance, FixedPoint=FixedPoint, 491 _tento=_tento):
492 """Return xn, yn, p s.t. 493 p = max(x.p, y.p) 494 x = xn / 10**p 495 y = yn / 10**p 496 497 x must be FixedPoint to begin with; if y is not FixedPoint, 498 it inherits its precision from x. 499 500 Note that this method is called a lot, so default-arg tricks are helpful. 501 """ 502 assert isinstance(x, FixedPoint) 503 if not isinstance(y, FixedPoint): 504 y = FixedPoint(y, x.p) 505 xn, yn = x.n, y.n 506 xp, yp = x.p, y.p 507 if xp > yp: 508 yn = yn * _tento(xp - yp) 509 p = xp 510 elif xp < yp: 511 xn = xn * _tento(yp - xp) 512 p = yp 513 else: 514 p = xp # same as yp 515 return xn, yn, p
516
517 -def _mkFP(n, p, FixedPoint=FixedPoint):
518 """Make FixedPoint objext - Return a new FixedPoint object with the selected precision.""" 519 f = FixedPoint() 520 #print '_mkFP Debug: %s, value=%s' % (type(f),n) 521 f.n = n 522 f.p = p 523 return f
524 525 # crud for parsing strings 526 import re 527 528 # There's an optional sign at the start, and an optional exponent 529 # at the end. The exponent has an optional sign and at least one 530 # digit. In between, must have either at least one digit followed 531 # by an optional fraction, or a decimal point followed by at least 532 # one digit. Yuck. 533 534 _parser = re.compile(r""" 535 \s* 536 (?P<sign>[-+])? 537 ( 538 (?P<int>\d+) (\. (?P<frac>\d*))? 539 | 540 \. (?P<onlyfrac>\d+) 541 ) 542 ([eE](?P<exp>[-+]? \d+))? 543 \s* $ 544 """, re.VERBOSE).match 545 546 del re 547 548
549 -def _string2exact(s):
550 """Return n, p s.t. float string value == n * 10**p exactly.""" 551 m = _parser(s) 552 if m is None: 553 raise ValueError("can't parse as number: " + `s`) 554 555 exp = m.group('exp') 556 if exp is None: 557 exp = 0 558 else: 559 exp = int(exp) 560 561 intpart = m.group('int') 562 if intpart is None: 563 intpart = "0" 564 fracpart = m.group('onlyfrac') 565 else: 566 fracpart = m.group('frac') 567 if fracpart is None or fracpart == "": 568 fracpart = "0" 569 assert intpart 570 assert fracpart 571 572 i, f = long(intpart), long(fracpart) 573 nfrac = len(fracpart) 574 i = i * _tento(nfrac) + f 575 exp = exp - nfrac 576 577 if m.group('sign') == "-": 578 i = -i 579 580 return i, exp
581
582 -def _test():
583 """Unit testing framework""" 584 fp = FixedPoint 585 o = fp("0.1") 586 assert str(o) == "0.10" 587 t = fp("-20e-2", 5) 588 assert str(t) == "-0.20000" 589 assert t < o 590 assert o > t 591 assert min(o, t) == min(t, o) == t 592 assert max(o, t) == max(t, o) == o 593 assert o != t 594 assert --t == t 595 assert abs(t) > abs(o) 596 assert abs(o) < abs(t) 597 assert o == o and t == t 598 assert t.copy() == t 599 assert o == -t/2 == -.5 * t 600 assert abs(t) == o + o 601 assert abs(o) == o 602 assert o/t == -0.5 603 assert -(t/o) == (-t)/o == t/-o == 2 604 assert 1 + o == o + 1 == fp(" +00.000011e+5 ") 605 assert 1/o == 10 606 assert o + t == t + o == -o 607 assert 2.0 * t == t * 2 == "2" * t == o/o * 2L * t 608 assert 1 - t == -(t - 1) == fp(6L)/5 609 assert t*t == 4*o*o == o*4*o == o*o*4 610 assert fp(2) - "1" == 1 611 assert float(-1/t) == 5.0 612 for p in range(20): 613 assert 42 + fp("1e-20", p) - 42 == 0 614 assert 1/(42 + fp("1e-20", 20) - 42) == fp("100.0E18") 615 o = fp(".9995", 4) 616 assert 1 - o == fp("5e-4", 10) 617 o.set_precision(3) 618 assert o == 1 619 o = fp(".9985", 4) 620 o.set_precision(3) 621 assert o == fp(".998", 10) 622 assert o == o.frac() 623 o.set_precision(100) 624 assert o == fp(".998", 10) 625 o.set_precision(2) 626 assert o == 1 627 x = fp(1.99) 628 assert long(x) == -long(-x) == 1L 629 assert int(x) == -int(-x) == 1 630 assert x == long(x) + x.frac() 631 assert -x == long(-x) + (-x).frac() 632 assert fp(7) % 4 == 7 % fp(4) == 3 633 assert fp(-7) % 4 == -7 % fp(4) == 1 634 assert fp(-7) % -4 == -7 % fp(-4) == -3 635 assert fp(7.0) % "-4.0" == 7 % fp(-4) == -1 636 assert fp("5.5") % fp("1.1") == fp("5.5e100") % fp("1.1e100") == 0 637 assert divmod(fp("1e100"), 3) == (long(fp("1e100")/3), 1)
638 639 if __name__ == '__main__': 640 _test() 641