对于浮点值，tan(x) = 无穷大的不可能性证明（或反例）
Proof of impossibility (or counterexample) for tan(x) = infinity, for floating-point values
C 中 tan 函数族 (tan、tanf、tanl) 的 POSIX page 表示：
If the correct value would cause overflow, a range error shall occur and tan(), tanf(), and tanl() shall return ±HUGE_VAL, ±HUGE_VALF, and ±HUGE_VALL, respectively, with the same sign as the correct value of the function.
然而，实际上在这种情况下非常很难获得infinity/-infinity，因为浮点数需要足够接近 π/2，使其正切大于该类型的最大可表示浮点值。
根据经验，我无法使用标准 glibc 获得这样的结果，即使使用 nextafter 函数来获得最接近 π/2 的可能值（使用 atan(1) * 2 作为“半 π” , 然后从那里开始，在它的左边或右边）。
对所有（32 位）floats 的详尽测试证实这对于 32 位是正确的，至少在我的库中是这样。测试 doubles 和 long doubles 的 π/2 附近表明它们也是这种情况。然而，详尽的测试有点太长了：我需要测试 (2k+1)•π/2 的附近，对于∀k ∈ ℤ.
那么，是否有一些数学或实践论证可以让人们得出这样的结论，至少在“合理正确”的库实现中（例如，像 GNU C library math functions 那样使用一些测量的 ULP 界限），浮点值的精度总是使得 tan(x) 适合这些值的有限表示？换句话说，tan 向无穷大增长的速度不会快于我们接近 π/2 的速度？
请注意，我将 tan(NAN) 和 tan(INFINITY) 排除在讨论之外，因为这些是记录在案的 return NaN 的极端情况。此外，次正规数可能会得到不同的结果，但由于它们只出现在零附近而不是 π/2 附近，我相信我们可以排除它们。
因此，我正在寻找一些表明这不会发生的数学 argument/proof/exhaustive 测试，或者只是一个反例，其中包含任何标准 C 库实现，例如 <math.h> 然后调用 tan 就可以了；但不包括具有非标准 tan 类函数的特定库。
... is there some mathematical or practical argument that allows one to conclude that, ... the precision of floating-point values will always be such that tan(x) will fit in a finite representation of these values?
π 是无理数，some_odd_integer*π/2 - 数学的极点 tangent(some_radian_measure).
所有有限浮点数都是有理数。因此没有 double 是 π/2 的整数倍并且 tan(some_finite_double) 永远不是极点。
some_odd_n*π/2 可能非常接近 double a，例如 6381956970095103 ∗ 2797。如果这是真的，那么 tan(a) 很大。
这接近 K.C. Ng's "ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit" Chapter 2.3。这大致意味着所有 double 的 some_odd_n*π/2 大约 2-62 的最接近大小写差异。这样的值 a 比 some_odd_n*π/2 差 2-62 导致 tan() 为 about 262，完全在 21023.
的 double 范围内
这个关键点是一些FP类型要tan(a)超过FP_MAX，a必须非常接近到some_odd_n*π/2，顺序为 log2(FP_MAX)。对于更宽的 FP 类型，这不太可能 发生，因为精度是线性的，范围与位宽呈指数关系。它可能发生在比 binary16.
窄的人为设计的 FP 类型中
注意：tan() 第一步是将参数减少到主要范围 [-π/2...π/2]。在我看来，这比主要 tan() 评估更难做好。因此，tan() 实施的质量（或缺乏）可能会产生非常大（和不正确）的结果，其值接近 some_odd_n*π/2。
I'm looking for either some mathematical argument/proof/exhaustive test that shows this does not happen,
有关如何确定最坏情况的想法，请参阅上面链接文章中的“使用与 π 相关的连分数”（参考文献 6）。
浮点值的精度由 IEEE 规定，并且只要您使用 IEEE 浮点数就始终相同。特别地，最接近 π/2 的 32 位 IEEE 数字是 1.57079637050628662109375（上图）和 1.5707962512969970703125（下图）。这些值是准确的。这两个数字的正切值在数学上不是特别大，大约为 -2.2E7 和 1.3E7，这与 32 位二进制数的 ~3E38 的 IEEE-754 限制相去甚远。
然而，这并没有说明 tan 对于大参数的行为。首先，绝对不能保证库函数 return 的结果最接近数学上正确的值（ 尤其是 对于大参数）。其次，对于某些奇数 N，N*π/2 的表示可能会非常接近真实的数学值，从而产生超出可表示范围的 tan 值。这两个观察结合起来意味着需要针对所有感兴趣的值详尽地测试 tan 的特定实现。
这里有一个实际演示，即正确舍入的 tan(x) 对于任何有限的 IEEE 754 binary64 double x 永远不会是无限的；事实上，它永远不会超过 2.2e+18。同样的方法也适用于 64 位精度 long double 类型，但搜索时间会稍长一些。它甚至应该适用于 IEEE 754 binary128 和 binary256 浮点格式。除此之外，启发式论据表明 极 不太可能 tan(x) 会溢出任何 IEEE 754 二进制浮点格式，但我不知道有任何证据，并且我怀疑不存在已知证据。
显然，仅考虑正有限浮点数就足够了（此处和下文中，“浮点数”表示 IEEE 754 binary64 浮点值）。对于某些整数 m 和 e 以及 0 < m < 2^53 和 -1074 <= e <= 971.
，每个这样的浮点数都具有 m * 2^e 的形式
所以我们要找到整数n、m和e最小化绝对值|n π/2 - m 2^e|，受约束0 < n、0 < m < 2^53 和 -1074 <= e <= 971。我们实际上只对 n 的 odd 值感兴趣，但最佳 n 已经是奇数：如果它是偶数，那么我们可以将 [=29] 减半=] 并递减 e 以将绝对值减半。此外，很明显检查 e 小于 -53 没有什么意义，因此我们可以进一步限制为 -53 <= e <= 971.
对于 e 的任何特定固定值，我们的问题等同于找到最小化绝对值 |n - m 2^(e+1)/π| 的正整数 m 和 n，服从约束 m < 2^53。如果我们可以有效地解决这个问题，我们可以迭代 e 的 1025 个不同的值来提取整体最佳解决方案。
但是固定指数 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:
的严格下限和上限
# Approximation to π, given by the first 1000 digits after the point. (source:
# https://mathshistory.st-andrews.ac.uk/HistTopics/1000_places/)
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()))
# True value of π is within ERROR of the value given above.
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 and c/d are the closest convergents so far
# sign is 1 if a/b < value < c/d, and -1 otherwise.
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|。如果 returns 两个三元组 (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
最后，我们将所有这些放在一起以找到 π/2.
的奇数倍的整体最佳 binary64 近似值
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，它与 Eric Postpischil 在对原始问题的评论中引用的值 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
因此没有有限 binary64 浮点数的切线超过 2.133485385753704e+18。
C 中 tan 函数族 (tan、tanf、tanl) 的 POSIX page 表示：
If the correct value would cause overflow, a range error shall occur and tan(), tanf(), and tanl() shall return ±HUGE_VAL, ±HUGE_VALF, and ±HUGE_VALL, respectively, with the same sign as the correct value of the function.
然而，实际上在这种情况下非常很难获得infinity/-infinity，因为浮点数需要足够接近 π/2，使其正切大于该类型的最大可表示浮点值。
根据经验，我无法使用标准 glibc 获得这样的结果，即使使用 nextafter 函数来获得最接近 π/2 的可能值（使用 atan(1) * 2 作为“半 π” , 然后从那里开始，在它的左边或右边）。
对所有（32 位）floats 的详尽测试证实这对于 32 位是正确的，至少在我的库中是这样。测试 doubles 和 long doubles 的 π/2 附近表明它们也是这种情况。然而，详尽的测试有点太长了：我需要测试 (2k+1)•π/2 的附近，对于∀k ∈ ℤ.
那么，是否有一些数学或实践论证可以让人们得出这样的结论，至少在“合理正确”的库实现中（例如，像 GNU C library math functions 那样使用一些测量的 ULP 界限），浮点值的精度总是使得 tan(x) 适合这些值的有限表示？换句话说，tan 向无穷大增长的速度不会快于我们接近 π/2 的速度？
请注意，我将 tan(NAN) 和 tan(INFINITY) 排除在讨论之外，因为这些是记录在案的 return NaN 的极端情况。此外，次正规数可能会得到不同的结果，但由于它们只出现在零附近而不是 π/2 附近，我相信我们可以排除它们。
因此，我正在寻找一些表明这不会发生的数学 argument/proof/exhaustive 测试，或者只是一个反例，其中包含任何标准 C 库实现，例如 <math.h> 然后调用 tan 就可以了；但不包括具有非标准 tan 类函数的特定库。
... is there some mathematical or practical argument that allows one to conclude that, ... the precision of floating-point values will always be such that tan(x) will fit in a finite representation of these values?
π 是无理数，some_odd_integer*π/2 - 数学的极点 tangent(some_radian_measure).
所有有限浮点数都是有理数。因此没有 double 是 π/2 的整数倍并且 tan(some_finite_double) 永远不是极点。
some_odd_n*π/2 可能非常接近 double a，例如 6381956970095103 ∗ 2797。如果这是真的，那么 tan(a) 很大。
这接近 K.C. Ng's "ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit" Chapter 2.3。这大致意味着所有 double 的 some_odd_n*π/2 大约 2-62 的最接近大小写差异。这样的值 a 比 some_odd_n*π/2 差 2-62 导致 tan() 为 about 262，完全在 21023.
double 范围内
这个关键点是一些FP类型要tan(a)超过FP_MAX，a必须非常接近到some_odd_n*π/2，顺序为 log2(FP_MAX)。对于更宽的 FP 类型，这不太可能 发生，因为精度是线性的，范围与位宽呈指数关系。它可能发生在比 binary16.
注意：tan() 第一步是将参数减少到主要范围 [-π/2...π/2]。在我看来，这比主要 tan() 评估更难做好。因此，tan() 实施的质量（或缺乏）可能会产生非常大（和不正确）的结果，其值接近 some_odd_n*π/2。
I'm looking for either some mathematical argument/proof/exhaustive test that shows this does not happen,
有关如何确定最坏情况的想法，请参阅上面链接文章中的“使用与 π 相关的连分数”（参考文献 6）。
浮点值的精度由 IEEE 规定，并且只要您使用 IEEE 浮点数就始终相同。特别地，最接近 π/2 的 32 位 IEEE 数字是 1.57079637050628662109375（上图）和 1.5707962512969970703125（下图）。这些值是准确的。这两个数字的正切值在数学上不是特别大，大约为 -2.2E7 和 1.3E7，这与 32 位二进制数的 ~3E38 的 IEEE-754 限制相去甚远。
然而，这并没有说明 tan 对于大参数的行为。首先，绝对不能保证库函数 return 的结果最接近数学上正确的值（ 尤其是 对于大参数）。其次，对于某些奇数 N，N*π/2 的表示可能会非常接近真实的数学值，从而产生超出可表示范围的 tan 值。这两个观察结合起来意味着需要针对所有感兴趣的值详尽地测试 tan 的特定实现。
这里有一个实际演示，即正确舍入的 tan(x) 对于任何有限的 IEEE 754 binary64 double x 永远不会是无限的；事实上，它永远不会超过 2.2e+18。同样的方法也适用于 64 位精度 long double 类型，但搜索时间会稍长一些。它甚至应该适用于 IEEE 754 binary128 和 binary256 浮点格式。除此之外，启发式论据表明 极 不太可能 tan(x) 会溢出任何 IEEE 754 二进制浮点格式，但我不知道有任何证据，并且我怀疑不存在已知证据。
显然，仅考虑正有限浮点数就足够了（此处和下文中，“浮点数”表示 IEEE 754 binary64 浮点值）。对于某些整数 m 和 e 以及 0 < m < 2^53 和 -1074 <= e <= 971.
m * 2^e 的形式
所以我们要找到整数n、m和e最小化绝对值|n π/2 - m 2^e|，受约束0 < n、0 < m < 2^53 和 -1074 <= e <= 971。我们实际上只对 n 的 odd 值感兴趣，但最佳 n 已经是奇数：如果它是偶数，那么我们可以将 [=29] 减半=] 并递减 e 以将绝对值减半。此外，很明显检查 e 小于 -53 没有什么意义，因此我们可以进一步限制为 -53 <= e <= 971.
对于 e 的任何特定固定值，我们的问题等同于找到最小化绝对值 |n - m 2^(e+1)/π| 的正整数 m 和 n，服从约束 m < 2^53。如果我们可以有效地解决这个问题，我们可以迭代 e 的 1025 个不同的值来提取整体最佳解决方案。
但是固定指数 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:
# Approximation to π, given by the first 1000 digits after the point. (source:
# https://mathshistory.st-andrews.ac.uk/HistTopics/1000_places/)
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()))
# True value of π is within ERROR of the value given above.
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 and c/d are the closest convergents so far
# sign is 1 if a/b < value < c/d, and -1 otherwise.
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|。如果 returns 两个三元组 (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
最后，我们将所有这些放在一起以找到 π/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，它与 Eric Postpischil 在对原始问题的评论中引用的值 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
因此没有有限 binary64 浮点数的切线超过 2.133485385753704e+18。