这里有一个实际演示，证明对于任何有限的IEEE 754 binary64 double类型的x，正确舍入的 tan(x) 永远不会是无穷大，事实上它永远不会超过 2.2e+18。相同的方法也适用于64位精度的long double类型，但是搜索需要更长时间。甚至可以对IEEE 754 binary128和binary256浮点格式进行处理。超出这个范围后，启发式论据表明对于任何IEEE 754二进制浮点格式，tan(x)都不可能溢出，但我不知道是否有任何证明，而且我认为不存在已知的证明。
显然，只考虑正有限浮点数就足够了（此处和以下中，“float”指的是IEEE 754 binary64浮点值）。每个这样的 float 都具有形式 m * 2^e ，其中m和e是一些整数，满足 0
因此，我们要找到整数n，m和e，使绝对值 |nπ/2-m2^e| 最小化，同时满足约束条件 0
对于任何特定固定值的 e，我们的问题等价于找到最小化绝对值 |n-m2^(e+1)/π| 的正整数 m 和 n ，其中 m<2^53。如果我们能高效地解决这个问题，我们就可以在1025个不同的 e 值上迭代，提取出总体最佳解决方案。
但是，对于固定的指数
e的子问题可以通过找到
2^(e+1) / π的连分数拓展来解决：标准的数论结果表明，
|n - m 2^(e+1) / π|的最小值保证来自于收敛或半收敛的
n / m至
2^(e+1) / π。我们所要做的就是找到最佳的收敛或半收敛，其分母严格小于
2^53。实际上，连分数拓展允许我们找到既是最佳收敛或半收敛（受制于分母小于
2^53的限制），而且比
|2^(e+1) / π|严格更大，也找到了最佳收敛或半收敛（具有类似限制的分母），其比
|2^(e+1) / π|严格更小，然后需要进行比较以确定哪一个收敛更接近。
下面是一些Python代码，它实现了上述内容，并最终找到了评论中已经提到的相同最佳逼近。首先，我们需要一些稍后会用到的导入：
from fractions import Fraction
from math import floor, ldexp
我们需要一个近似值来表示π/2。实际上，我们需要一对近似值：一个下限和一个上限。我们将使用前1000位π的数字，因为它们有良好的记录，并且在网上易于找到（尽管实际上我们只需要约330个数字）。使用这些数字，我们创建Fraction实例HALF_PI_LOWER和HALF_PI_UPPER，分别代表π/2的紧密下限和上限：
PI_APPROX = Fraction("".join("""
3.14159265358979323846264338327950288419716939937510
58209749445923078164062862089986280348253421170679
82148086513282306647093844609550582231725359408128
48111745028410270193852110555964462294895493038196
44288109756659334461284756482337867831652712019091
45648566923460348610454326648213393607260249141273
72458700660631558817488152092096282925409171536436
78925903600113305305488204665213841469519415116094
33057270365759591953092186117381932611793105118548
07446237996274956735188575272489122793818301194912
98336733624406566430860213949463952247371907021798
60943702770539217176293176752384674818467669405132
00056812714526356082778577134275778960917363717872
14684409012249534301465495853710507922796892589235
42019956112129021960864034418159813629774771309960
51870721134999999837297804995105973173281609631859
50244594553469083026425223082533446850352619311881
71010003137838752886587533208381420617177669147303
59825349042875546873115956286388235378759375195778
18577805321712268066130019278766111959092164201989
""".split()))
ERROR = Fraction("1e-1000")
HALF_PI_LOWER = (PI_APPROX - ERROR) / 2
HALF_PI_UPPER = (PI_APPROX + ERROR) / 2
现在我们来讲一下连分数机制。首先，我们提供了一个生成器函数，当提供该值的下限和上限时，它会生成该值的连分数系数序列。基本上，这个函数会分别计算下限和上限的系数，而只要这些系数匹配，它们就会被返回给调用者。一旦它们不匹配，我们就无法确定下一个系数的值，此时该函数将引发一个异常（可能是RuntimeError或ZeroDivisionError）。我们并不太关心引发哪种异常：我们依赖于我们最初的近似值足够准确，从而避免了遇到该异常的情况。
def continued_fraction_coefficients(lower: Fraction, upper: Fraction):
"""
Generate continued fraction coefficients for a positive number.
Given rational lower and upper bounds for a postive rational or irrational
number, generate the known coefficients of the continued fraction
representation of that number. The tighter the lower and upper bounds, the
more coefficients will be known.
Raises RuntimeError or ZeroDivisionError when there's insufficient
information to determine the next coefficient.
"""
while (q := floor(lower)) == floor(upper):
yield q
lower, upper = 1/(upper - q), 1/(lower - q)
raise RuntimeError("Insufficient resolution")
接下来是一个函数，它使用这个连分数系数序列来产生对应值的最佳下限和最佳上限近似值（对于给定的分母限制）。它从连分数序列中生成连分数逼近，超过分母限制后回溯寻找符合分母限制的最佳边界。
def best_bounds(coefficients, max_denominator):
"""
Find best lower and upper bounds under a given denominator bound.
*coefficients* is the sequence of continued fraction coefficients
for the value we're approximating.
*max_denominator* is the limit on the denominator.
Returns the best lower approximation and the best upper approximation
to the value being approximated, for the given limit on the denominator.
"""
a, b, c, d, sign = 0, 1, 1, 0, 1
for q in coefficients:
a, b, c, d, sign = c, d, a + q * c, b + q * d, -sign
if max_denominator < d:
correction = (max_denominator - d) // b
c, d = c + correction * a, d + correction * b
break
if sign == 1:
return Fraction(a, b), Fraction(c, d)
else:
return Fraction(c, d), Fraction(a, b)
现在我们编写一个专用函数，对于给定的指数
e，在所有
n和
m且
m < 2^53的情况下，最小化
|n π/2 - m 2^e|。该函数返回两个三元组
(x, n, error)，其中
x是浮点值
m 2^e，表示为Python中的
float，
error给出了近似值的误差，也被转换为浮点型。（因此，
error值并不完全准确，但它们足够接近，以便我们找到最佳的整体近似值。）这两个三元组分别表示最佳下限和上限。
def best_for_given_exponent(e):
"""
Find best lower and upper approximations for 2^e / (π/2)
with denominator smaller than 2**53.
Returns two triples of the form (x, n, error), where:
- x is a float close to an integer multiple of π/2
- n is the coefficient of that integer multiple of π/2
- error is the absolute error |x - n π/2|, rounded to the nearest float.
The first triple gives a best lower approximation for a multiple of π/2,
while the second triple gives a best upper approximation.
"""
twopow = Fraction(2)**e
coefficients = continued_fraction_coefficients(twopow/HALF_PI_UPPER, twopow/HALF_PI_LOWER)
best_lower, best_upper = best_bounds(coefficients, max_denominator=2**53 - 1)
n, d = best_lower.as_integer_ratio()
lower_result = ldexp(d, e), n, float(d * twopow - n * HALF_PI_LOWER)
n, d = best_upper.as_integer_ratio()
upper_result = ldexp(d, e), n, float(n * HALF_PI_UPPER - d * twopow)
return lower_result, upper_result
最后，我们将所有内容综合起来，找到一个最佳的binary64逼近值，使其为
奇数个π/2的倍数。
all_results = [
result for e in range(-53, 972)
for result in best_for_given_exponent(e)
]
x, n, error = min(all_results, key=lambda result: result[2])
print(f"Best approximation: {x.hex()} is an approximation to {n}*π/2 with approximate error {error}")
当我在我的电脑上运行这个程序时，大约经过5秒钟的计算，它会输出以下内容：
Best approximation: 0x1.6ac5b262ca1ffp+849 is an approximation to 3386417804515981120643892082331156599120239393299838035242121518428537554064774221620930267583474709602068045686026362989271814411863708499869721322715946622634302011697632972907922558892710830616034038541342154669787134871905353772776431251615694251273653*π/2 with approximate error 4.687165924254628e-19
因此，最糟糕的情况是浮点数 0x1.6ac5b262ca1ffp+849，它匹配了埃里克•波斯特皮斯希尔在原问题的评论中提到的值6381956970095103•2^797：
>>> float.fromhex('0x1.6ac5b262ca1ffp+849') == ldexp(6381956970095103, 797)
True
该值的正切约为±1/4.687165924254628e-19，或2.133485385753704e+18，事实上，这就是我计算机上数学库所产生的值：
>>> from math import tan
>>> tan(float.fromhex('0x1.6ac5b262ca1ffp+849'))
-2.133485385753704e+18
因此，没有任何有限的二进制64位浮点数其正切值超过
2.133485385753704e+18。
tan函数的简化合同，例如：“如果 x 是有限的，则 tan(x) 也是有限的”。 - anol