prePEP: Decimal data type
Bengt Richter
bokr at oz.net
Sat Nov 8 17:42:13 EST 2003
On Thu, 6 Nov 2003 10:43:52 -0500, "Tim Peters" <tim.one at comcast.net> wrote:
[...]
>
>It has to support some half dozen specific rounding modes already (the spec
>requires them). The mandated modes have been reviewed by many people now
>over several years, so it's tempting to call YAGNI on adding more rounding
>complexity. My FixedPoint does support user-defined rounding modes; alas,
>all evidence to date suggests I'm the only person to date able to define one
>correctly <wink -- but doing rounding exactly right requires exact
>understanding of every stinking detail of every stinking endcase -- it's not
>easy the first time someone tries it>.
>
Following is an exact decimal/<somewhat>rational module hack that implements rounding
when the .round() method is called. This is not too elegant, but it might make
an interesting toy to cross-check results with. 30 years of monthly compound interest
does generate a fair amount of digits, but it can handle it. And you can check the _exact_
difference from doing it floating point, and perhaps demonstrate when/how it's safe
to use hardware floating point (the below has a software decimal floating point).
About the 30 years interest,
>>> from exactdec import ED
>>> monthly5 = 1 + ED( '.05/12')
>>> monthly5
ED('3.0125 / 3')
>>> start = ED(1)
>>> start
ED('1')
>>> fm5 = 1.0+ 0.05/12
>>> fm5
1.0041666666666667
>>> monthly5.round(20)
ED('1.00416666666666666667')
>>> fstart = 1.0
>>> for i in xrange(30*12):
... fstart *= fm5
... start *= monthly5
...
>>> fstart
4.4677443140061053
>>> start.round(20)
ED('4.46774431400613221243')
>>> err = ED(fstart,'all')-start
>>> err
ED('-1.56304148231013236986795175352733832996669654337153733730960097118272041942305024046142360
963012057392399385808513132735004380045887019491590594848567543746099279247782492404211092305697
981942603081504204074464788922982296549717752974330439000221572233435214122399335151325954144618
900462032107403557206447160788112551619335537177864834308165003815651982743545600184830416792430
006198225629411824849523906416912497782230577924310271500267297956136382334806613723525165395284
535900838608901737436464025004714288622793588896506577479134017475377641227803205011713549091096
133145948633022116071771866555211557151726225215476514565787630190565855468620450845083504202642
689190693005234276689807253839309304532499468247786496752880094937190474830588626230621777889479
703217124659874396057456877927948952659639736522210676920364370104538827117231761130501253771264
878907949516858455102771211294274323473541556175745838790698849963815157512517119109683330894122
419465602517028421797269789171429961866520967620825108423637749408425943209381542284540329915368
857997235383065705584065571271008721284350021123549868449933295248144722514707438777251932770746
090987731235646452226688847245082616867972045549031014767048143096992525819923039646122088361644
579918406562388873940433411106385019399448321943189013300306911684700844253559215810364410055218
335449640228259572789992966032626222008918446886031635163380806543318900671084821575936160455692
493336358710092579196071811371493642452908504603585136777462572632798060122885204925037488881947
910781074947690321238995118414971674525304479175247251987457275390625e158 / 58029883553010642635
729412690188122648867574983421180398276039568933683892061549947054264296629049795373656629387735
93537918329060451697709838267789886857631025156060199201')
>>> err.round(50)
ED('-2.693511319701830367844690393524096646e-14')
>>> err.round(49)
ED('-2.69351131970183036784469039352409665e-14')
>>> err.round(48)
ED('-2.6935113197018303678446903935240966e-14')
But there is a lot of digits in the exact result, and the percentage error of the floating
point isn't so bad:
>>> (err/start).round(30)
ED('-6.028794690103059e-15')
====< exactdec.py >====================================================================
# exactdec.py -- exact decimal/rational class
# 20031107 14:54:21 alpha -- bokr
# 20031108 09:42:18 added contruction from float with all info or specified rounding -- bokr
# 20031108 13:58:10 added more docs and interactive calculator thing -- bokr
#
# WARNING: Not signigicantly tested. Not optimized. No warranty. Use at your own risk. And,
# if you run this program interactively, it passes your raw_input to eval unvalidated!
#
# The ED class preserves accuracy exactly by including a denominator for a kind of
# hybrid rational. The info defining the value can be retrieved via ED.astuple
# as a 3-tuple of (numerator, denominator, power_of_10) for the exact value
# (numerator*10**power_of_10)/denominator.
#
# Currently, rounding is a separate operation that returns an exact result. Only rounding
# away from zero on halves is implemented (and not very well tested ;-)
#
# Examples of valid string literals being passed to the constructor are:
# ED('-12.34e-56/78') -- which by the way currently reprs as the exactly equivalent
# normalized ED('-6.17e-56 / 39')
# ED('.1') which reprs as ED('.1')
# Converting the same floating point value preserving total or otherwise specified
# accuracy get you:
# ED(.1, 'all') => ED('1.000000000000000055511151231257827021181583404541015625e-1')
# ED(.1, 24) => ED('1.00000000000000005551115e-1')
# ED(.1, 23) => ED('1.0000000000000000555112e-1')
# ED(.1, 22) => ED('1.000000000000000055511e-1')
# ED(.1, 20) => ED('0.10000000000000000555')
# ED(.1, 19) => ED('0.1000000000000000056')
# ED(.1, 18) => ED('0.100000000000000006')
# ED(.1, 17) => ED('0.10000000000000001')
# ED(.1, 16) => ED('0.1')
#
# Integers and longs are accepted as accurate, and value will be preserved, though
# representation may be normalized.
#
# ED(2L**55) => ED('36028797018963968')
# ED(20L**55) => ED('3.6028797018963968e71')
# ED(10L**55) => ED('1.0e55')
# The inverse can be calculated accurately, and mutiplied out again:
# 1 / ED(10L**55) => ED('1.0e-55')
# 1 / ED(20L**55) => ED('2.77555756156289135105907917022705078125e-72')
# (1 / ED(20L**55)) * ED(2**55) => ED('1.0e-55')
# (1 / ED(20L**55)) * ED(2**55) * ED('1e55') => ED('1')
#
#############################################################################
#
# The literal for an exact value is a string
import re
rxo = re.compile(r'\s*([+-]?\d*)(\.(\d*))?(e([+-]?\d+))?(\s*/\s*([+-]?\d+))?')
# example
# rxo.search('12.34e-45 / 678' ).groups()[::2]
# ('12', '34', '-45', '678')
def fracnorm(num,den):
"""make numerator carry sign and return equivalent num,den with no common divisor"""
if num==0: return 0,1
elif den==0: raise ValueError,'zero fracnorm den!'
if den<0: num=-num; den=-den
n=abs(num); d=den
if n<d: n,d = d,n
while d != 0: n, d = d, n % d
return num//n, den//n
class ED(object):
"""exact decimal"""
def __init__(self, v=0, decimals=None):
"""
decimals is only mandatory for exactifying float, done by rounding
to a least significant digit of 10**-decimals (decimals may be <=0),
but may be specified as 'all' to capture all bits of float as an accurate value.
"""
if isinstance(v, (int,long)):
self.num = v
self.den = 1
self.p10=0
elif isinstance(v, str):
"""r'\s*([+-]?\d*)(\.(\d*))?(e([+-]?\d+))?(\s*/\s*([+-]?\d+))?'"""
# XXX needs validation checks
# <i1>.<i2>e<p10> / <den>
i1,i2,p10,den = rxo.search(v).groups()[::2]
if den is None: den = 1
else: den = int(den)
if p10 is None: p10=0
else: p10 = int(p10)
if i2 is None: i2=''
p10 -= len(i2)
self.num = int(i1+i2)
self.den = den
self.p10 = p10
elif isinstance(v, float):
if decimals is None or (not isinstance(decimals, int) and decimals!='all'):
raise ValueError(
"Specify decimals for least significant digit of 10**(-decimals)\n"
"(decimals may also be specified as 'all' to capture all bits of float)")
sign = v<0
if sign: v = -v
acc = type(self)(0)
p2=1; shift2=2**27 # usual double has 53 bits
while v:
i = int(v) # depend on accurate truncation and promotion to long
acc = acc + type(self)((i,p2,0))
v = (v-i)*shift2
p2 *= shift2
if sign: acc.num = -acc.num
self.num = acc.num
self.den = acc.den
self.p10 = acc.p10
elif isinstance(v, type(self)):
self.num = v.num
self.den = v.den
self.p10 = v.p10
elif isinstance(v, tuple):
self.num = v[0]
self.den = v[1]
self.p10 = v[2]
else:
raise ValueError, 'Exactdec constructor does not accept %r' %type(v)
self.normalize()
if decimals !='all' and decimals!=None:
self.num, self.den, self.p10 = self.round(decimals).astuple()
def __repr__(self):
self.normalize()
num = self.num; den=self.den; p10=self.p10
s = str(abs(num))
sign = '-'[:num<0] # '' or '-'
dens = den != 1 and ' / %s'%den or ''
# normalize to x.xxxxe<p10> if total string > 20 long
zs = len(s)+p10 # leading digits or leading-zero deficit if negative
if (zs<0 and -zs or p10)+len(s)>20:
return "ED('%s%s.%se%s%s')" %(sign, s[0], s[1:].rstrip('0') or '0', zs-1, dens)
zs = zs<0 and '0'*-zs or ''
s = zs+s
if p10>=0:
return "ED('%s%s%s%s')" %(sign, s, '0'*p10, dens)
elif zs>0:
return "ED('%s%s.%s%s')" %(sign, s[:p10] or '0', s[p10:], dens)
else:
return "ED('%s%se%s%s')"%(sign,num,p10,dens)
def __mul__(self,other):
"""multiply exactly"""
other = type(self)(other)
num = self.num * other.num
den = self.den * other.den
p10 = self.p10 + other.p10
return type(self)((num,den,p10)).normalize()
__rmul__ = __mul__
def __add__(self, other):
"""add exactly"""
other = type(self)(other)
num = self.num; p10 = self.p10; den = self.den
numo = other.num; p10o = other.p10; deno = other.den
dshift = p10-p10o
if dshift:
if dshift<0:
p10o += dshift
numo *= 10**-dshift
else:
p10 -= dshift
num *= 10**dshift
numsum = num*deno+numo*den
densum = den*deno
return type(self)((numsum, densum, p10)).normalize()
__radd__ = __add__
def __sub__(self, other):
"""substract exactly (negates other and adds)"""
other = type(self)(other)
other.num = -other.num
return self + other
def __rsub__(self, other):
"""add as right operand"""
return type(self)(other) - self
def __div__(self,other):
"""divide exactly -- multiplies inverse"""
other = type(self)(other)
num = self.num * other.den
den = self.den * other.num
p10 = self.p10 - other.p10
return type(self)((num,den,p10)).normalize()
def __rdiv__(self, other):
"""divide in role of right arg"""
return type(self)(other)/self
def normalize(self):
"""
make num part carry sign
remove common factors of num, den parts
convert factors of 10 to exponent count
convert den factors of 2 to 10/5 and remove 10's to exponent
"""
self.num, self.den = fracnorm(self.num, self.den)
while self.num and self.num%10==0:
self.num //= 10
self.p10+=1
while not self.den%2:
self.num *=5
self.den //=2
self.p10 -=1
while not self.den%5:
self.num *=2
self.den //=5
self.p10 -=1
return self
def round(self, ndec=2):
"""
Return rounded value as exact decimal.
Round away from zero on fractions of .5 and more
XXX needs alternative rounding choices.
"""
ret = type(self)(self)
sign = ret.num<0
ret.num = abs(ret.num)
ret = ret + type(self)((5,1,-ndec-1))
psh = -(ret.p10+ndec)
ret.num //= (ret.den*10**psh)
ret.den = 1
ret.p10 += psh
if sign: ret.num = -ret.num
ret.normalize()
return ret
def astuple(self): self.normalize(); return self.num, self.den, self.p10
def __eq__(self, other): return self.astuple() == type(self)(other).astuple()
def __lt__(self, other): return (self-other).num<0
def __nonzero__(self): return self.num!=0
def __int__(self): return (self.num*10**self.p10)//self.den
def test():
import sys
for x in xrange(-9,10):
sys.stdout.write('+')
for y in xrange(-9,10):
if not y: continue # skip zero denom
for p10 in xrange(-9,10):
s='%se%s/%s'%(x,p10,y)
assert ED(s) == ED((x,y,p10))
print
if __name__ == '__main__':
print """
Enter 'test' to run test (such as it is), otherwise
Enter arguments to be passed to the ED constructor, or
prefix with *, /, +, or - to operate with previous result
Enter r <number of decimals> to round previous result.
Note that your input is passed to eval unvalidated, so
do NOT let this run where that would be dangerous!!
Enter quotes to make string literals for ED constructor.
To pass a float per se, specify decimals, e.g.,
> 0.1, 'all'
ED('0.1000000000000000055511151231257827021181583404541015625')
> 0.1, 22
ED('0.1000000000000000055511')
> 0.1, 18
ED('0.100000000000000006')
> 0.1, 16
ED('0.1')
To constrast, a string literal (second below) is accurate:
> 0.1, 'all'
ED('0.1000000000000000055511151231257827021181583404541015625')
> - '0.1'
ED('5.5511151231257827021181583404541015625e-18')
Press Enter only, to Exit."""
value = ED(0)
while 1:
try:
cmd = raw_input('%r\n> '%value).strip()
if not cmd: break
if cmd == 'test': test(); continue
if not cmd[0] in '*/+-r':
value = eval('ED(%s)'%cmd)
else:
if cmd[0]=='r':
value = eval('%s.round(%s)'%(value, cmd[1:].strip()))
else:
value = eval('%s %s ED(%s)'%(value, cmd[0], cmd[1:].strip()))
except (EOFError, KeyboardInterrupt), e:
raise SystemExit, 'Exiting due to %s'%e.__class__.__name__
except Exception, e:
print '%s: %s'%(e.__class__.__name__, e)
=======================================================================================
Regards,
Bengt Richter
More information about the Python-list
mailing list