C语言 tan(x)=无穷大的不可能性证明(或反例),对于浮点值

eqqqjvef  于 12个月前  发布在  其他
关注(0)|答案(3)|浏览(97)

在C语言中,tan函数家族(tantanftanl)的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类功能的特定库。

vql8enpb

vql8enpb1#

这里有一个实际的证明,正确舍入的tan(x)对于任何有限的IEEE 754 binary 64 double x都不是无限的;实际上它永远不会超过2.2e+18。同样的方法也适用于64位精度的long double类型,但搜索时间会稍长一些。它甚至可以用于IEEE 754 binary 128和binary 256浮点格式。除此之外,启发式的论点表明tan(x)对于任何IEEE 754二进制浮点格式都不太可能溢出,但我不知道任何证据,我怀疑没有已知的证据存在。
显然,只考虑正有限浮点数就足够了(在这里和下面,“浮点数”意味着IEEE 754二进制64浮点值)。每个这样的浮点数都有m * 2^e的形式,对于一些整数me,以及0 < m < 2^53-1074 <= e <= 971
因此,我们希望找到整数nme,使绝对值|n π/2 - m 2^e|最小化,受到约束0 < n0 < m < 2^53-1074 <= e <= 971。我们实际上只对n的 * 奇 * 值感兴趣,但最佳n已经是奇的了:如果它是偶数,那么我们可以将n减半,并将e递减,以将绝对值减半。此外,很明显,检查e小于-53没有什么意义,所以我们可以进一步限制为-53 <= e <= 971
对于e的任何特定固定值,我们的问题等价于找到使绝对值|n - m 2^(e+1)/π|最小化的正整数mn,受到约束m < 2^53。如果我们能有效地解决这个问题,我们可以通过e1025不同值来提取整体最佳解决方案。
但是固定指数e的子问题可以通过找到2^(e+1) / π的连分式展开来解决:标准数论结果表明,|n - m 2^(e+1) / π|的最小值保证来自于n / m2^(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_LOWERHALF_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

现在是连分数机器。首先,我们给予一个生成函数,当给定一个值的上下界时,该生成函数生成该值的连分数系数序列。本质上,这将分别计算下限和上限的系数,当这些系数匹配时,它们将被提供给调用者。一旦它们不匹配,我们就没有足够的信息来确定下一个系数,函数将在此时引发异常(RuntimeErrorZeroDivisionError)。我们不太关心哪个异常被提出:我们依赖于原始的近似值,以确保足够精确,从而不会遇到异常。

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,在所有nm上最小化|n π/2 - m 2^e|,其中m < 2^53。If返回两个三元组(x, n, error),其中x是浮点值m 2^e,表示为Python floaterror给出近似值中的误差,也转换为浮点数。(所以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的奇数倍的整体最佳binary 64近似。

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

7kjnsjlb

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。这大致意味着对于所有doublesome_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)。

toiithl6

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的所有感兴趣的值进行详尽的测试。

相关问题