“就本质来说,浮点算术是不精确的,而且程序员们很容易滥用它,从而使计算的结果几乎全部由噪声组成”
–Donald Knuth(《程序设计艺术》( 第二卷) 半数值算法)
一. 舍入困惑
Python2中,round函数使用靠 近 最 近 和 等 距 远 离 0 ‾ \underline{靠近最近和等距远离0}
靠近最近和等距远离0
(ROUND_HALF_UP)策略,是通常习惯使用的四舍五入模式。
(Values are rounded to the closest multiple of 10 to the power minus ndigits; if two multiples are equally close, rounding is done away from 0 )
而在Python3中,round函数的取舍方式使用靠 近 最 近 和 等 距 靠 近 偶 数 ‾ \underline{靠近最近和等距靠近偶数}
靠近最近和等距靠近偶数
(ROUND_HALF_EVEN)策略,
(values are rounded to the closest multiple of 10 to the power minus ndigits; if two multiples are equally close, rounding is done toward the even choice)
同时,float类型采用双精度二进制存储(参照IEEE754标准),round函数使用的是二进制存储值,在舍入时会对结果产生影响,而round本身没有使用四舍五入规范,就造成了一些困惑。
不好理解的舍入:
>>> round(1.205, 2) | |
1.21 | |
...... | |
>>> round(1.245, 2) | |
1.25 | |
>>> round(1.255, 2) | |
1.25 | |
>>> round(1.265, 2) | |
1.26 | |
...... | |
>>> round(1.295, 2) | |
1.29 | |
Numpy的round出现同样问题:
>>> numpy.round(1.255, 2) | |
1.25 | |
>>> numpy.round(1.265, 2) | |
1.26 |
大多数浮点数出现“尾部乱码”与“999变异“:
>>> format(1.215, '.52f') | |
'1.2150000000000000799360577730112709105014801025390625' | |
>>> round(1.215, 2) | |
1.22 | |
>>>formatl(1.275, '.52f') | |
'1.2749999999999999111821580299874767661094665527343750' | |
>>> round(1.275, 2) | |
1.27 |
似乎并不是全部:
>>> format(1.125, '.52f') | |
'1.1250000000000000000000000000000000000000000000000000' | |
>>> round(1.125, 2) | |
1.12 | |
>>> format(0.5, '.52f') | |
'0.5000000000000000000000000000000000000000000000000000' | |
>>> round(0.5, 0) | |
0 |
当需要掌握精确的计算结果时,round()似乎成为头疼的问题。
本来一个简单的问题,为什么变得难以理解,难道是Python的round策略出了问题?
作为一个专业级语言平台,Python显然不会没有想到这样的问题。
其实,浮点数问题并非我们想象的那样简单。
二. 问题分析
就浮点数的精度问题,正如计算大师Knuth所说,这是一个超出我们想象的复杂问题,其结果往往使我们很颓丧,甚至不敢再相信计算机。要掌握其计算规律,确实需要一些较为耐心和深入的探讨。在浮点数计算方面,我们有很多错觉,往往会将计算机的有限计算与代数中的计算相混淆,以至于对很多结果产生迷惑。对于一个认真的程序员来说,如果想知道浮点数运算结果的确切性,需要了解计算机对浮点数的表示和处理方式。想得到比较彻底的答案,需要具体研究一下IEEE 754标准,或研读一下Kuth的《计算机程序设计艺术》。
针对Python中浮点数float的四舍五入问题, 我们需要理解两个方面:二进制表示方式、十进制舍入策略。
$\underline{浮点数类型float与round的设计}$
- float采用IEEE754格式,使用二进制编码存储浮点数,与其他计算机语言并无太多差异。float类型是双精度浮点数。
- round取舍方式是靠近最近和偶数,这个策略符合大规模计算的总体逼近优化原则,未采用通常的四舍五入策略。
- 对于高精度运算和十进制小数精确表示,Python提供了专门的模块decimal,并提供了可选择的舍入策略,包括四舍五入。
$\underline{float类型的二进制表示}$
- float采用二进制编码描述浮点数。在二进制表示中,大多数有限位十进制小数无法使用二进制进行有限位精确表示。也就是说,有限位数的十进制小数,往往会变为无限位数的二进制小数。 事实上,分母中含有非2质数因子的分数,都不能使用有限位二进制小数表示。而十进制小数分母中含有质数因子5,如果约分后分母中仍然含有因子5,就会变成无限位二进制小数。
- 对不能使用有限位二进制小数表示的十进制有限位小数,在系统中存储的是这些十进制浮点数的近似值。在近似值中,分为进位和截断两种类型,近似误差一般在${10}^{-17}$左右。进位近似值大于原值,截断近似值小于原值,所以对小数位较小的数值(如1.215被进位,1.275被截断),进位近似值会出现尾部的增加值(上面说的尾部乱码),截断近似值小于原值,会出现”999...“的近似值现象。
- 表示为二进制近似值后,Python系统在进行round计算时,使用近似值,不是使用原值。
【例1】十进制小数0.1无法使用有限位二进制精确表示,同一原值的近似值都是相等的。
从十进制到二进制 :$(0.1){10} = (0.0001100110011001100110011001100110011001100110011001...){2}$
在系统中,原值转为近似值:
>> a = 1.2355000000000000426325641456060111522674560546875 | |
>> b = 1.2355 | |
>> a == b | |
True |
【例2】经过二进制表示转换后,原值已经不存在,但有时系统显示会产生误会。
我们赋值给Python的一个数值,他会依据浮点数标准,将其转换为二进制存储的近似值。于是,系统内已经不是我们原来的实际数值。而对于系统来说,他把这两个数值认为是一样的。在实际运算中,他使用的是那个近似值。而容易迷惑的是,为了显示简捷,Python显示给用户的有时还是原值。
>>> x = 0.1 | |
>>> print(x) | |
0.1 | |
>>> format(x, '.20f') | |
'0.10000000000000000555' | |
>>> format(x, '.50f') | |
'0.10000000000000000555111512312578270211815834045410' |
$\underline{内置函数round的舍入分析}$
<font face="黑体" color="Crimson">round(number, ndigits)
【values are rounded to the closest multiple of 10 to the power minus ndigits; if two multiples are equally close, rounding is done toward the even choice (so, for example, both round(0.5) and round(-0.5) are 0, and round(1.5) is 2). 】
对一个给定的浮点数$v$,这里不妨只考虑小于1的小数,即浮点数的有效尾数部分,忽略符号和指数部分。
其十进制表示为:
$v{10} = d{1} \cdots d{m} =d{1}10^{-1}+\cdots+d_{m}10^{-m}$
在Python中会使用float类型处理,使用与其最接近的二进制值存储表示:
$v{2} = b{1} \cdots b{n} = b{1}2^{-1}+\cdots+b_{n}2^{-n}$
这个二进制小数近似值$v_{2}$还原为10进制时,一般是不等于原值的,除非原值可以表示为2的有限负幂和。
如:0.125 = $0\frac{1}{2}+0\frac{1}{4}+1*\frac{1}{8} = 0.001$
$\underline {round进行小数位的舍入时, 实际上是对 v_{2}的十进制值进行处理}$
使用round对十进制数$v{10}$的第k位四舍五入时,实际上转为对$v{2}$的十进制值的处理,Python3目前的处理原则是:
(1)$v_{2}$第k+1位是0-4或6-9时,舍入是非常明确的,可以称之位4舍6入。
(2)$v{2}$第k+1位是5时,要看$v{2}$与$d{1}...d{k}$和$d{1}...(d{k}+1)$之间的距离:
- -【1】 如果两个距离不相等,结果取最近的取值。
- -【2】 如果两个距离相等,采用靠近偶数原则,即,如果$d{k}$是偶数就取$d{1}...d{k}$,否则取$d{1}...(d_{k}+1)$。
<font face="黑体" color=#0099ff>看下面两个浮点数保留两位小数的四舍五入:
- 1.275 的二进制近似表示值: $v_{2}(1.275)$ = 1.274999999999999911182158029987476766109466552734375 计算1.27和1.28哪一个离$v_2(1.275)$更近, 显然1.27更近,于是:round(1.275, 2) = 1.27round(1.215, 2) = 1.22<font face="黑体" color=#0099ff>两边一样近的情况(原数值是2的幂次的有限表示)向偶数靠近, 如1.125 = 1 + 0$\frac{1}{2}$ + 0$\frac{1}{4}$ + 1$\frac{1}{8}$ = $1 + 12^{-3}$ ,结果为:
- 1.215 的二进制近似表示值: $v_{2}(1.215)$ = 1.2150000000000000799360577730112709105014801025390625 由于$v_2(1.215)$与1.22比1.21更近,结果便是:
>>> round(1.125, 2) | |
1.2 | |
>>> round(1.5, 0) | |
2.0 | |
>>> round(0.5, 0) | |
0.0 |
看来靠round解决四舍五入是不行的,因为其设计目标不在于此,简单调整也满足不了四舍五入的要求。
逼近舍入有利于数据分析的精确性,是一个误差最小策略。
同时,从用户的角度来看,round也受二进制表示的影响。只考虑四舍五入问题的话,一定精度范围内
仅与舍入规则有关。但用户给出的值首先要转换到双精度近似值,round的规则用于这个近似值四舍五入,
就要考虑精度范围。这个范围对float 来说,就是52位二进制存储精度,即十进制的17位小数有效位之内。
这个有效位数包括整数部分位数。
三. 解决方法
字符串转换方法(不支持位数为负数,其他问题有待检验...):
def round45s(v, d=0): | |
vl = str(v).split('.') | |
sign = -1 if v < 0 else 1 | |
if len(vl) == 2: | |
if len(vl[1]) > d: | |
if vl[1][d] >= '5': | |
if d > 0: | |
vp = (eval(vl[1][:d]) + 1)/10**d * sign | |
return int(v)+vp | |
else: | |
return int(v)+sign | |
else: | |
if d > 0: | |
return int(v) + eval(vl[1][:d])/10**d * sign | |
else: | |
return int(v) | |
else: | |
return v | |
return int(v) |
改进的方法:
这个方法将原值转换为略大的十进制值,从而使输入值的有限小数位数字不会发生变化(后面是0值,直到15位),避免出现"999变异"。但受到双精度二进制近似存储影响,只能在十进制有效位数15位以内使用(digits<15)。需要注意,整数部分的位数也考虑在有效位数之内。
def round45r(number, digits=0): | |
int_len = len(str(int(abs(number)))) | |
signal_ = 1 if number >= 0 else -1 | |
err_place = 16 - int_len - 1 | |
if err_place > 0: | |
err_ = 10**-err_place | |
return round(number + err_ * signal_, digits) | |
else: | |
raise NotImplemented # 受到float表示精度的限制! |
round45r() 对负数和整数也有效,即支持v, d为负数的情况:
1.205, 2) | round45r(-|
-1.210000000000002 # 在16位补误差,保障前面的数字不会变化 | |
123.205, -2) | round45r(|
100.00000000000001 | |
153.205, -2) | round45r(|
200.0 | |
153.205, -2) | round45r(-|
-200.0 |
如果运行时间可以承受,尽量考虑使用高精度十进制表示模块decimal(精度和舍入策略都可以控制):
from decimal import Decimal, Context, ROUND_HALF_UP | |
def roundl45d(v, d): | |
if int(v) > 0: | |
d = d + str(v).find('.') # 有效位包括整数部分 | |
return float(Context(prec=d, rounding=ROUND_HALF_UP).create_decimal(str(v))) | |
1.205) | decimal45(|
1.21 | |
1.255) | decimal45(|
1.26 |
效率测试结果:
1000000) | test_round45s(number=|
6.26826286315918e-06 | |
1000000) | test_round45r(number=|
1.287583589553833e-06 | |
1000000) | test_round45d(number=|
1.7323946952819824e-06 |
round45s比round45r有5倍的运算速度差异。
round45d与round45r基本在一个水平。
四. 进一步思考
- 实现更高效的方法,应该考虑使用c模块去写round45。 >>> x = round45(1.275, 2); print(x) 1.28 >>> format(x, '.30f') '1.280000000000000026645352591004'
- 如果计算结果仍然使用浮点类型float表示,其值还是Python中存储的二进制双精度近似数。
- 要精确进行浮点数运算,建议使用decimal模块,并通过字符串进行赋值,并根据计算需要设置精度和舍入策略。
>>> from decimal import Decimal, gecontext, ROUND_HALF_UP | |
>>> Decimal(‘1.675’) | |
Decimal(‘1.675’) # 使用字符串得到精确值 | |
>>> getcontext().prec = 52 # 设置精度为52位(这里指十进制) | |
>>> Decimal(‘1’) / Decimal(str(2**52)) | |
Decimal('2.220446049250313080847263336181640625E-16') # 按照要求精确表示有效位数字(实际有效位数37位) | |
>>> Decimal('0.1') + Decimal('0.2') | |
Decimal('0.3') # 精度范围内,准确运算 | |
>>> getcontext().prec = 3 | |
>>> getcontext().rounding = ROUND_HALF_UP | |
>>> Decimal('1.275') | |
Decimal('1.28') |
了解一下decimal的八种舍入策略
在decimal中,实现了八种rounding策略,计算时可以根据需要选取。
这些应该是科学计算和数据处理中经常使用的,也基本已经成为标准。
有些文章也谈到了Java BigDecimal的舍入策略,基本与Python decimal 类同,只是decimal多了一个ROUND_05UP。
通过decimal运算环境Context管理,设置取舍精度和策略
# 获取运算环境 | tc = getcontext()|
tc.prec | |
28 # 缺省精度 | |
tc.rounding | |
decimal.ROUND_HALF_EVEN # 缺省策略为ROUND_HALF_EVEN | |
5 # 设置5位有效数字精度 | tc.prec =|
# 设置为新的舍入策略 | tc.rounding = decimal.ROUND_HALF_UP
decimal运算中的八种舍入策略:
1) ROUND_CEILING 向正无穷(Infinity)靠近
tc.rounding = decimal.ROUND_CEILING | |
1.12345’) | tc.create_decimal(‘|
Decimal('1.1235') # 正数时靠近较大方向 | |
1.12345’) | tc.create_decimal(‘-|
Decimal('-1.1234') # 负数时靠近绝对值较小方向 |
2) ROUND_DOWN 向0靠近
tc.rounding = decimal.ROUND_DOWN | |
1.12345’) | tc.create_decimal(‘|
Decimal('1.1234') # 正数时是向下取整 | |
1.12345’) | tc.create_decimal(‘-|
Decimal('-1.1234') # 负数时靠近绝对值较小方向 |
3) ROUND_FLOOR 向负无穷(-Infinity)靠近
tc.rounding = decimal.ROUND_FLOOR | |
1.12345’) | tc.create_decimal(‘|
Decimal('1.1234') # 正数时是向下取整 | |
1.12345’) | tc.create_decimal(‘-|
Decimal('-1.1235') # 负数时靠近绝对值较大方向 |
4) ROUND_HALF_DOWN 向最接近的近似值靠近,两边相等时靠近0方向
tc.rounding = decimal.ROUND_HALF_DOWN | |
1.12346’) | tc.create_decimal(‘|
Decimal('1.1235') # 两端不相等,为4舍6入 | |
1.12346’) | tc.create_decimal(‘-|
Decimal('-1.1235') # 两端不相等,为4舍6入 | |
1.12345’) | tc.create_decimal(‘|
Decimal('1.1234') # 两端相等,正数时是向下取整 | |
1.12345’) | tc.create_decimal(‘-|
Decimal('-1.1234') # 两端相等,负数时靠近绝对值较小方向 |
5) ROUND_HALF_EVEN 向最接近的近似值靠近;两边相等时,前面是奇数进位,偶数不进位
tc.rounding = decimal.ROUND_EVEN | |
1.12346’) | tc.create_decimal(‘|
Decimal('1.1235') # 两端不相等,为4舍6入 | |
1.12346’) | tc.create_decimal(‘-|
Decimal('-1.1235') # 两端不相等,为4舍6入 | |
1.12345’) | tc.create_decimal(‘|
Decimal('1.1234') # 两端相等,向偶数靠近 | |
1.12335’) | tc.create_decimal(‘-|
Decimal('-1.1234') # 两端相等,向偶数靠近 |
6) ROUND_HALF_UP 向最接近的近似值靠近,两边相等时靠近远离0的方向
tc.rounding = decimal.ROUND_HALF_UP | |
1.12346’) | tc.create_decimal(‘|
Decimal('1.1235') # 两端不相等,为4舍6入 | |
1.12346’) | tc.create_decimal(‘-|
Decimal('-1.1235') # 两端不相等,为4舍6入 | |
1.12345’) | tc.create_decimal(‘|
Decimal('1.1235') # 两端相等,正数向靠Infinity近 | |
1.12345’) | tc.create_decimal(‘-|
Decimal('-1.1235') # 两端相等,负数向-Infinity靠近 |
7) ROUND_UP 靠近远离0的方向
tc.rounding = decimal.ROUND_UP | |
1.12345’) | tc.create_decimal(‘|
Decimal('1.1235') # 正数时向上取整 | |
1.12345’) | tc.create_decimal(‘-|
Decimal('-1.1235') # 负数时向下取整 |
8) ROUND_05UP 如果向0靠近取舍后保留小数的最后一位是0或5,就向远离0靠近,否则就向0靠近。
tc.rounding = decimal.ROUND_05UP | |
1.12343’) | tc.create_decimal(‘|
Decimal('1.1234') # 向0靠近取舍后最后一位不是0或5 | |
1.12351’) | tc.create_decimal(‘|
Decimal('1.1236') # 向0靠近取舍后最后一位是5,选择远离0 | |
1.12355’) | tc.create_decimal(‘-|
Decimal('-1.1236') # 向0靠近取舍后最后一位是5,选择远离0 | |
1.12305’) | tc.create_decimal(‘-|
Decimal('-1.1231') # 向0靠近取舍后最后一位是0,选择远离0 |