如何在双双算术中正确四舍五入到最接近的整数

问：

我必须使用 Python3（PyPy 实现）分析大量数据，在其中我对相当大的浮点数进行一些操作，并且必须检查结果是否足够接近整数。

举例来说，假设我正在生成随机的数字对，并检查它们是否形成毕达哥拉斯三元组（是具有整数边的直角三角形的边）：

from math import hypot
from pprint import pprint
from random import randrange
from time import time

def gen_rand_tuples(start, stop, amount):
    '''
    Generates random integer pairs and converts them to tuples of floats.
    '''
    for _ in range(amount):
        yield (float(randrange(start, stop)), float(randrange(start, stop)))

t0 = time()
## Results are those pairs that results in integer hypothenuses, or
## at least very close, to within 1e-12.
results = [t for t in gen_rand_tuples(1, 2**32, 10_000_000) if abs((h := hypot(*t)) - int(h)) < 1e-12]
print('Results found:')
pprint(results)
print('finished in:', round(time() - t0, 2), 'seconds.')

运行它我得到了：

Python 3.9.17 (a61d7152b989, Aug 13 2023, 10:27:46)
[PyPy 7.3.12 with GCC 13.2.1 20230728 (Red Hat 13.2.1-1)] on linux
Type "help", "copyright", "credits" or "license()" for more information.
>>> 
===== RESTART: /home/user/Downloads/pythagorean_test_floats.py ====
Results found:
[(2176124225.0, 2742331476.0),
 (342847595.0, 3794647043.0),
 (36.0, 2983807908.0),
 (791324089.0, 2122279232.0)]
finished in: 2.64 seconds.

有趣的是，它运行得很快，在 2 秒多一点的时间内处理了 1000 万个数据点，我甚至找到了一些匹配的数据。hypothenuse 显然是整数：

>>> pprint([hypot(*x) for x in results])
[3500842551.0, 3810103759.0, 2983807908.0, 2265008378.0]

但事实并非如此，如果我们使用十进制任意精度模块检查结果，我们会发现结果实际上并不够接近整数：

>>> from decimal import Decimal
>>> pprint([(x[0]*x[0] + x[1]*x[1]).sqrt() for x in (tuple(map(Decimal, x)) for x in results)])
[Decimal('3500842551.000000228516418075'),
 Decimal('3810103758.999999710375341513'),
 Decimal('2983807908.000000217172157183'),
 Decimal('2265008377.999999748566051441')]

所以，我认为问题在于数字足够大，以至于落在 python 浮点数缺乏精度的范围内，因此返回了误报。

现在，我们可以将程序更改为在任何地方使用任意精度小数：

from decimal import Decimal
from pprint import pprint
from random import randrange
from time import time

def dec_hypot(x, y):
    return (x*x + y*y).sqrt()

def gen_rand_tuples(start, stop, amount):
    '''
    Generates random integer pairs and converts them to tuples of decimals.
    '''
    for _ in range(amount):
        yield (Decimal(randrange(start, stop)), Decimal(randrange(start, stop)))

t0 = time()
## Results are those pairs that results in integer hypothenuses, or
## at least very close, to within 1e-12.
results = [t for t in gen_rand_tuples(1, 2**32, 10_000_000) if abs((h := dec_hypot(*t)) - h.to_integral_value()) < Decimal(1e-12)]
print('Results found:')
pprint(results)
print('finished in:', round(time() - t0, 2), 'seconds.')

现在，我们没有得到任何误报，但我们对性能造成了很大的影响。以前需要 2 秒多一点，现在需要 100 多秒。小数点似乎对 JIT 不友好：

====== RESTART: /home/user/Downloads/pythagorean_test_dec.py ======
Results found:
[]
finished in: 113.82 seconds.

我找到了这个问题的答案，CPython 和 PyPy Decimal 操作性能，建议使用双倍精度数字作为小数的更快、JIT 友好的替代方案，以获得比内置浮点数更好的精度。所以我pip安装了doubledouble第三方模块，并相应地更改了程序：

from doubledouble import DoubleDouble
from decimal import Decimal
from pprint import pprint
from random import randrange
from time import time

def dd_hypot(x, y):
    return (x*x + y*y).sqrt()

def gen_rand_tuples(start, stop, amount):
    for _ in range(amount):
        yield (DoubleDouble(randrange(start, stop)), DoubleDouble(randrange(start, stop)))

t0 = time()
print('Results found:')
results = [t for t in gen_rand_tuples(1, 2**32, 10_000_000) if abs((h := dd_hypot(*t)) - int(h)) < DoubleDouble(1e-12)]
pprint(results)
print('finished in:', round(time() - t0, 2), 'seconds.')

但是我收到这个错误：

======= RESTART: /home/user/Downloads/pythagorean_test_dd.py ======
Results found:
Traceback (most recent call last):
  File "/home/user/Downloads/pythagorean_test_dd.py", line 24, in <module>
    results = [t for t in gen_rand_tuples(1, 2**32, 10_000_000) if abs((h := dd_hypot(*t)) - int(h)) < DoubleDouble(1e-12)]
  File "/home/user/Downloads/pythagorean_test_dd.py", line 24, in <listcomp>
    results = [t for t in gen_rand_tuples(1, 2**32, 10_000_000) if abs((h := dd_hypot(*t)) - int(h)) < DoubleDouble(1e-12)]
TypeError: int() argument must be a string, a bytes-like object or a number, not 'DoubleDouble'

我认为问题是该模块没有指定转换或舍入到最接近的整数方法。我能写的最好的是一个非常人为的“int”函数，它通过对字符串和小数进行往返并返回 DoubleDouble 来将 double-double 四舍五入到最接近的整数：

def contrived_int(dd):
    rounded_str = (Decimal(dd.x) + Decimal(dd.y)).to_integral_value()
    hi = float(rounded_str)
    lo = float(Decimal(rounded_str) - Decimal(hi))
    return DoubleDouble(hi, lo)

但它非常迂回，违背了回避小数的目的，并使程序甚至比全十进制版本更慢。

然后我问，有没有一种快速的方法可以将双精度数字直接舍入到最接近的整数，而无需中间步骤通过小数或字符串？