【问题标题】:Early termination of fractional exponent calculation?提前终止小数指数计算?
【发布时间】:2021-01-09 21:07:30
【问题描述】:

我需要编写一个函数,该函数取某事物的六次方根(相当于将某事物提高到 1/6 次方),并检查答案是否为整数。我希望这个函数尽可能快并尽可能优化,并且由于这个函数需要运行很多,我认为最好不必计算整个根。

在必须计算整个第六个根之前,我将如何编写一个返回 False(或 0 或等效项)的函数(与语言无关,尽管首选 Python/C/C++)?例如,如果我取 65 的 6 次方根,那么我的函数应该在意识到结果不是 int 后,停止计算并返回 False,而不是首先计算 65 的 6 次方是 2.00517474515 ,然后检查 2.00517474515 是否为int,最后返回False

当然,我问这个问题的印象是提前终止的事情比完整的计算更快,使用类似的东西

print(isinstance(num**(1/6), int))

任何帮助或想法将不胜感激。我也会对可推广到许多分数幂的答案感兴趣,而不仅仅是 x^(1/6)。

【问题讨论】:

  • 由于输入必须是整数,因此一种可能性是对该输入执行素数分解。那么测试和计算就很容易了
  • 您关心的数字的典型大小是多少?你在处理千位整数吗?更大?
  • @john: 六次方,不是六次方!
  • 您是否考虑过为此使用gmpympz.iroot 计算整数根并返回精确度指示;在优化的 C 库中执行此操作可能比在 Python 级别实现的任何早期拒绝算法更快。
  • 单独地,取决于你想优化多少:在做任何数学之前,你可以构造一个六次方的bloom filter——如果最大整数大小是 1e60 然后设置大小是 1e10;具有 ~1% 误报概率的布隆过滤器应该适合 RAM(大约 11GB)。

标签: algorithm math


【解决方案1】:

以下是您可以尝试的一些想法,这些想法可能有助于快速消除非六次方。对于实际的六次方,您最终仍然需要计算六次方。

检查小箱子

如果您给出的数字有合理的小概率(例如,少于 12 位),您可以建立一个小案例表并进行检查。小于 10**12 的六次方只有 100 个。如果你的输入总是更大,那么这个测试的价值就很小,但它仍然是一个非常便宜的测试。

消除小素数

任何小的质因数都必须以 6 的倍数出现。为避免过多的试除,您可以将一些小因数捆绑在一起。 例如,2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 = 223092870,它足够小以适应 Python 中的单个 30 位肢体,因此使用该模数的单个模运算应该很快。

所以给定一个测试号n,计算g = gcd(n, 223092870),如果结果不是1,检查n 是否能被g ** 6 整除。如果不是,n 不是六次方,你就完了。如果n 能被g**6整除,重复n // g**6

检查模124488的值(例如)

如果您执行了上一步,那么此时您的值不能被任何小于25 的素数整除。现在您可以使用精心选择的模数进行模数测试:例如,与124488 = 8 * 9 * 7 * 13 * 19 相对质数的任何六次方都与[1, 15625, 19657, 28729, 48385, 111385] 模数124488 的六个值之一全等。可以使用更大的模数,但代价是必须检查更多可能的残基。

检查是否为正方形

任何六次方都必须是正方形。由于 Python(至少 Python >= 3.8)有一个相当快的内置整数平方根函数,因此在计算完整的六次方根之前检查该值是否为平方是有效的。 (如果它一个平方并且你已经计算了平方根,现在你只需要提取一个立方根而不是六次根。)

使用浮点运算

如果输入不是太大,比如 90 位或更小,并且它是六次方,那么浮点运算就有合理的机会准确地找到六次方。但是,Python 不保证幂运算的准确性,因此值得进行一些额外的检查以确保结果在预期范围内。对于较大的输入,浮点运算获得正确结果的机会较小。 (2**53 + 1)**6 的第六个根不能完全表示为 Python 浮点数(合理假设 Python 的 float 类型与 IEEE 754 binary64 格式匹配),一旦 n 超过 308 位左右,它就太大而无法容纳无论如何都要变成一个浮点数。

使用整数运算

一旦你用尽了所有便宜的技巧,你别无选择,只能计算六次方的底,然后将其与原始数字的六次方进行比较。

这里有一些 Python 代码,将上面列出的所有技巧组合在一起。您应该针对您的特定用例制定自己的时间安排,并选择哪些技巧值得保留,哪些应该调整或丢弃。技巧的顺序也很重要。

from math import gcd, isqrt

# Sixth powers smaller than 10**12.
SMALL_SIXTH_POWERS = {n**6 for n in range(100)}

def is_sixth_power(n):
    """
    Determine whether a positive integer n is a sixth power.

    Returns True if n is a sixth power, and False otherwise.
    """
    # Sanity check (redundant with the small cases check)
    if n <= 0:
        return n == 0

    # Check small cases
    if n < 10**12:
        return n in SMALL_SIXTH_POWERS

    # Try a floating-point check if there's a realistic chance of it working
    if n < 10**90:
        s = round(n ** (1/6.))
        if n == s**6:
            return True
        elif (s - 1) ** 6 < n < (s + 1)**6:
            return False
        # No conclusive result; fall through to the next test.

    # Eliminate small primes
    while True:
        g = gcd(n, 223092870)
        if g == 1:
            break
        n, r = divmod(n, g**6)
        if r:
            return False

    # Check modulo small primes (requires that 
    # n is relatively prime to 124488)
    if n % 124488 not in {1, 15625, 19657, 28729, 48385, 111385}:
        return False

    # Find the square root using math.isqrt, throw out non-squares
    s = isqrt(n)
    if s**2 != n:
        return False

    # Compute the floor of the cube root of s
    # (which is the same as the floor of the sixth root of n).
    # Code stolen from https://stackoverflow.com/a/35276426/270986
    a = 1 << (s.bit_length() - 1) // 3 + 1
    while True:
        d = s//a**2
        if a <= d:
            return a**3 == s
        a = (2*a + d)//3

【讨论】:

  • 这太棒了!愚蠢的我并没有考虑提前计算出我需要的所有能力,然后轻松检查集合中的成员资格。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多