在C语言中,tan
函数家族(tan
,tanf
,tanl
)的POSIX page表示:
如果正确的值会导致溢出,则应发生范围错误,并且tan()、tanf()和tanl()应分别返回±HUGE_VAL、±HUGE_VALF和±HUGE_VALL,符号与函数的正确值相同。
然而,在实际中,在这种情况下,实际上很难获得+∞或−∞作为计算结果,因为浮点数需要非常接近π/2的奇数倍,以至于它的正切将大于该类型的最大可表示浮点值。由于可表示值的间隔,我们不能总是像我们希望的那样接近特定的真实的数。
从经验上讲,我无法用标准的glibc得到这样的结果,即使使用nextafter
函数来获得最接近π/2的值(使用atan(1) * 2
作为“半π",并从那里开始,要么向左,要么向右)。
对所有(32位)float
s的详尽测试证实了这对于32位是正确的,至少在我的库中是这样。测试double
s和long double
s的π/2附近表明,这也是它们的情况。然而,详尽的测试有点太长了:我需要测试(2k+1)·π/2的邻域,对于所有k ∈ π,使得(2k+1)·π/2在浮点有限范围内。
那么,是否有一些数学或实际的论据可以让人们得出这样的结论:至少在“合理准确”的库实现中(* 例如 *,像对GNU C library math functions所做的那样,使用一些测量的ULP边界),浮点值的精度将总是使得tan(x)
将适合这些值的有限表示?换句话说,tan
向无穷大的增长速度不会比我们接近π/2的速度快?
请注意,我在讨论中排除了tan(NAN)
和tan(INFINITY)
,因为这些都是返回NaN的有文档记录的极端情况。同样,对于次正规数也可能得到不同的结果,但由于它们只出现在零附近而不是π/2附近,我相信我们可以排除它们。
因此,我正在寻找一些数学论证/证明/详尽的测试来证明这不会发生,或者只是一个反例,在任何标准C库实现中,包括<math.h>
并调用tan
就可以做到这一点;但不包括具有非标准X1 M15 N1 X类功能的特定库。
3条答案
按热度按时间vql8enpb1#
这里有一个实际的证明,正确舍入的
tan(x)
对于任何有限的IEEE 754 binary 64 doublex
都不是无限的;实际上它永远不会超过2.2e+18
。同样的方法也适用于64位精度的long double
类型,但搜索时间会稍长一些。它甚至可以用于IEEE 754 binary 128和binary 256浮点格式。除此之外,启发式的论点表明tan(x)
对于任何IEEE 754二进制浮点格式都不太可能溢出,但我不知道任何证据,我怀疑没有已知的证据存在。显然,只考虑正有限浮点数就足够了(在这里和下面,“浮点数”意味着IEEE 754二进制64浮点值)。每个这样的浮点数都有
m * 2^e
的形式,对于一些整数m
和e
,以及0 < m < 2^53
和-1074 <= e <= 971
。因此,我们希望找到整数
n
,m
和e
,使绝对值|n π/2 - m 2^e|
最小化,受到约束0 < n
,0 < m < 2^53
和-1074 <= e <= 971
。我们实际上只对n
的 * 奇 * 值感兴趣,但最佳n
已经是奇的了:如果它是偶数,那么我们可以将n
减半,并将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代码,并最终找到了与注解中提到的相同的最佳近似。首先,我们稍后需要一些导入:
我们需要
π/2
的近似值。实际上,我们需要一对近似:下界和上界。我们将使用π
的前1000位,因为它们有很好的文档记录,并且很容易在网上找到(尽管实际上我们实际上只需要大约330位)。使用这些数字,我们创建Fraction
示例HALF_PI_LOWER
和HALF_PI_UPPER
,表示π/2
的严格下限和上限:现在是连分数机器。首先,我们给予一个生成函数,当给定一个值的上下界时,该生成函数生成该值的连分数系数序列。本质上,这将分别计算下限和上限的系数,当这些系数匹配时,它们将被提供给调用者。一旦它们不匹配,我们就没有足够的信息来确定下一个系数,函数将在此时引发异常(
RuntimeError
或ZeroDivisionError
)。我们不太关心哪个异常被提出:我们依赖于原始的近似值,以确保足够精确,从而不会遇到异常。接下来,这里有一个函数,它使用这个连分数系数序列来产生对应值的最佳下近似值和最佳上近似值(对于给定的分母限制)。它从连分数序列生成连续的收敛,并且在超过分母限制之后,它回溯以找到适合分母限制的最佳界限。
现在我们编写一个专用函数,对于固定指数
e
,在所有n
和m
上最小化|n π/2 - m 2^e|
,其中m < 2^53
。If返回两个三元组(x, n, error)
,其中x
是浮点值m 2^e
,表示为Pythonfloat
,error
给出近似值中的误差,也转换为浮点数。(所以error
的值并不完全精确,但它们足够接近,可以让我们找到最佳的整体近似值。)两个三元组代表最佳的下限和上限。最后,我们把这些放在一起,找到
π/2
的奇数倍的整体最佳binary 64近似。当我在我的机器上运行这个程序时,它给出了以下输出(经过大约5秒的计算):
因此,总体上最差的情况是浮点数
0x1.6ac5b262ca1ffp+849
,它与Eric Postpischil在对原始问题的评论中引用的值6381956970095103•2^797
相匹配:这个值的正切值大约是
±1/4.687165924254628e-19
,或者2.133485385753704e+18
,实际上这就是我机器上的数学库产生的值:因此,没有一个有限binary64浮点数的正切值超过
2.133485385753704e+18
。7kjnsjlb2#
.是否有一些数学或实际的论点,让一个结论,.浮点值的精度将总是这样,tan(x)将适合这些值的有限表示?
π是一个无理数,就像
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()
为 * 约 * 262,正好在double
的21023范围内。这个关键点是,对于某些FP类型,
tan(a)
超过FP_MAX,a
必须 * 非常接近 *some_odd_n*π/2
,大约为log 2(FP_MAX)。对于更宽的FP类型,这不太可能发生,因为精度是线性的,范围与位宽呈指数关系。这可能发生在比binary16窄的人工FP类型上。注意:
tan()
第一步是将参数减少到主范围[-π/2... π/2]。这一点,IMO,比主要的tan()
评估更难做好。因此,tan()
实现的质量(或缺乏)可能会在some_odd_n*π/2
附近产生更大(和不正确)的结果。我在寻找一些数学论证/证明/详尽的测试来证明这不会发生,
有关如何确定最差情况的想法,请参见上述链接文章中的“使用与π相关的连分数”(参考文献6)。
toiithl63#
浮点值的精度由IEEE规定,并且只要您使用IEEE浮点数,精度就始终相同。特别是,最接近π/2的32位IEEE数是1.57079637050628662109375(上图)和1.5707962512969970703125(下图)。这些值是精确的。这两个数字的正切的数学正确值不是特别大,大约是-2.2E7和1.3E7,这与IEEE-754对32位二进制数的限制~ 3E 38相差甚远。
然而,这并没有说明
tan
对于大参数的行为。首先,绝对不能保证库函数将返回最接近数学正确值的结果(* 特别是 * 对于大参数)。第二,对于某个奇数N,N*π/2的表示可能足够接近真实的数学值,从而产生一个在可表示范围之外的值tan
。这两个观察结果结合起来意味着需要对tan
的所有感兴趣的值进行详尽的测试。