【问题标题】:How to compute the nth root of a very big integer如何计算一个非常大的整数的第 n 个根
【发布时间】:2008-12-10 13:49:19
【问题描述】:

我需要一种方法来计算 Python 中长整数的第 n 个根。

我试过pow(m, 1.0/n),但是不行:

溢出错误:long int 太大而无法转换为浮点数

有什么想法吗?

长整数是指真正的长整数,例如:

11968003966030964356885611480383408833172346450467339251 196093144141045683463085291115677488411620264826942334897996389 485046262847265769280883237649461122479734279424416861834396522 819159219215308460065265520143082728303864638821979329804885526 557893649662037092457130509980883789368448042961108430809620626 059287437887495827369474189818588006905358793385574832590121472 680866521970802708379837148646191567765584039175249171110593159 305029014037881475265618958103073425958633163441030267478942720 703134493880117805010891574606323700178176718412858948243785754 898788359757528163558061136758276299059029113119763557411729353 915848889261125855717014320045292143759177464380434854573300054 940683350937992500211758727939459249163046465047204851616590276 724564411037216844005877918224201569391107769029955591465502737 961776799311859881060956465198859727495735498887960494256488224 613682478900505821893815926193600121890632

【问题讨论】:

  • 正如大卫暗示的那样,pow(n, 1/3) 会给你 n 的三次(即第三)根。
  • 不,它不会,因为 python
  • (但这也不是 OP 想要的)。
  • Py3 没有整数限制......它们可以永远增长,直到内存耗尽。我在我的安装上进行了测试。这是一个解决方案。

标签: python math nth-root


【解决方案1】:

如果它是一个非常大的数字。你可以使用二分搜索。

def find_invpow(x,n):
    """Finds the integer component of the n'th root of x,
    an integer such that y ** n <= x < (y + 1) ** n.
    """
    high = 1
    while high ** n <= x:
        high *= 2
    low = high/2
    while low < high:
        mid = (low + high) // 2
        if low < mid and mid**n < x:
            low = mid
        elif high > mid and mid**n > x:
            high = mid
        else:
            return mid
    return mid + 1

例如:

>>> x = 237734537465873465
>>> n = 5
>>> y = find_invpow(x,n)
>>> y
2986
>>> y**n <= x <= (y+1)**n
True
>>>
>>> x = 119680039660309643568856114803834088331723464504673392511960931441>
>>> n = 45
>>> y = find_invpow(x,n)
>>> y
227661383982863143360L
>>> y**n <= x < (y+1)**n
True
>>> find_invpow(y**n,n) == y
True
>>>

【讨论】:

  • 数字 > ~10**1000 仍然会失败。将 mid = (low + high) // 2 更改为 mid = int((low + high) // 2) + 1 可以解决此问题。
  • 投反对票,因为它有错误。我尝试了find_invpow(64, 3),得到了 3,即使结果应该是 4。
  • 如果你真的想在末尾输出一个整数,那么low = high/2 应该是low = high // 2。
  • 二分法搜索似乎是一种相当缓慢的方法。更喜欢牛顿的迭代,它每次都会使精确数字的数量翻倍。要找到一个非常准确的起始值,请使用浮点数。
【解决方案2】:

Gmpy 是一个 C 编码的 Python 扩展模块,它包装了 GMP 库,为 Python 代码提供快速多精度算术(整数、有理数和浮点数)、随机数生成、高级数论函数等。

包括一个root 函数:

x.root(n):返回一个 2 元素元组 (y,m),使得 y 是 (可能被截断)x 的第 n 个根; m,一个普通的 Python int, 如果根是精确的 (x==y**n),则为 1,否则为 0。n 必须是普通的 Python 整数,>=0。

例如,第 20 个根:

>>> import gmpy
>>> i0=11968003966030964356885611480383408833172346450467339251 
>>> m0=gmpy.mpz(i0)
>>> m0
mpz(11968003966030964356885611480383408833172346450467339251L)
>>> m0.root(20)
(mpz(567), 0)

【讨论】:

  • 使用gmpy2,结果为:'mpz' object has no attribute 'root'。
  • gmpy2 使用基于 MPFR 库的新 mpfr 类型。 gmpy2.root(x, n) -> mpfr 返回 x 的第 n 个根。结果始终是“mpfr”。
  • @Zelphir gmpy2 有 gmpy2.iroot 来计算整数根。
【解决方案3】:

你可以通过避免 while 循环来让它运行得稍微快一点,有利于将 low 设置为 10 ** (len(str(x)) / n) 和 high to low * 10。可能更好的是替换 len( str(x)) 具有按位长度并使用位移。根据我的测试,我估计第一次加速 5%,第二次加速 25%。如果整数足够大,这可能很重要(并且加速可能会有所不同)。不仔细测试就不要相信我的代码。我做了一些基本测试,但可能错过了一个极端情况。此外,这些加速效果因所选数量而异。

如果您使用的实际数据比您在此处发布的数据大得多,那么此更改可能是值得的。

from timeit import Timer

def find_invpow(x,n):
    """Finds the integer component of the n'th root of x,
    an integer such that y ** n <= x < (y + 1) ** n.
    """
    high = 1
    while high ** n < x:
        high *= 2
    low = high/2
    while low < high:
        mid = (low + high) // 2
        if low < mid and mid**n < x:
            low = mid
        elif high > mid and mid**n > x:
            high = mid
        else:
            return mid
    return mid + 1

def find_invpowAlt(x,n):
    """Finds the integer component of the n'th root of x,
    an integer such that y ** n <= x < (y + 1) ** n.
    """
    low = 10 ** (len(str(x)) / n)
    high = low * 10

    while low < high:
        mid = (low + high) // 2
        if low < mid and mid**n < x:
            low = mid
        elif high > mid and mid**n > x:
            high = mid
        else:
            return mid
    return mid + 1

x = 237734537465873465
n = 5
tests = 10000

print "Norm", Timer('find_invpow(x,n)', 'from __main__ import find_invpow, x,n').timeit(number=tests)
print "Alt", Timer('find_invpowAlt(x,n)', 'from __main__ import find_invpowAlt, x,n').timeit(number=tests)

标准 0.626754999161

Alt 0.566340923309

【讨论】:

  • Your second function, find_invpowAlt, gives a wildly wrong answer for x=118997879821732370764604711647724283139870175351576755860556891902958645241483485254092600557474860904935286687480039428945219115513349647465379580432922136155355040992635166676363150438436216219094913514982415747153956476970303302126880391024128871557664284712411567099374094385902892603751471822837746770111, n=3
  • 如果你真的想在末尾输出一个整数,那么low = high/2 应该是low = high // 2。
  • 我故意使用low = high/2 编写它,因为与我提出的改进相比,我直接从Markus's Code 复制了该答案作为对其进行基准测试的一种手段。我会注意到,从那时起,Markus 已经更新了他的代码以更正一个错误,而 tzs 已经回复我的代码,它也有一个错误。但是,我实际上不记得问题的细节,也不记得我的解决方案,所以我不再觉得有能力解决这些缺点。
【解决方案4】:

如果您正在寻找标准的东西,请快速以高精度编写。我会使用小数并将精度 (getcontext().prec) 调整为至少 x 的长度。

代码(Python 3.0)

from decimal import *

x =   '11968003966030964356885611480383408833172346450467339251\
196093144141045683463085291115677488411620264826942334897996389\
485046262847265769280883237649461122479734279424416861834396522\
819159219215308460065265520143082728303864638821979329804885526\
557893649662037092457130509980883789368448042961108430809620626\
059287437887495827369474189818588006905358793385574832590121472\
680866521970802708379837148646191567765584039175249171110593159\
305029014037881475265618958103073425958633163441030267478942720\
703134493880117805010891574606323700178176718412858948243785754\
898788359757528163558061136758276299059029113119763557411729353\
915848889261125855717014320045292143759177464380434854573300054\
940683350937992500211758727939459249163046465047204851616590276\
724564411037216844005877918224201569391107769029955591465502737\
961776799311859881060956465198859727495735498887960494256488224\
613682478900505821893815926193600121890632'

minprec = 27
if len(x) > minprec: getcontext().prec = len(x)
else:                getcontext().prec = minprec

x = Decimal(x)
power = Decimal(1)/Decimal(3)

answer = x**power
ranswer = answer.quantize(Decimal('1.'), rounding=ROUND_UP)

diff = x - ranswer**Decimal(3)
if diff == Decimal(0):
    print("x is the cubic number of", ranswer)
else:
    print("x has a cubic root of ", answer)

回答

x 是 22873918786185635329056863961725521583023133411 的立方数 451452349318109627653540670761962215971994403670045614485973722724603798 107719978813658857014190047742680490088532895666963698551709978502745901 704433723567548799463129652706705873694274209728785041817619032774248488 2965377218610139128882473918261696612098418

【讨论】:

    【解决方案5】:

    哦,对于 大的数字,您可以使用小数模块。

    ns: 你的数字作为字符串

    ns = "11968003966030964356885611480383408833172346450467339251196093144141045683463085291115677488411620264826942334897996389485046262847265769280883237649461122479734279424416861834396522819159219215308460065265520143082728303864638821979329804885526557893649662037092457130509980883789368448042961108430809620626059287437887495827369474189818588006905358793385574832590121472680866521970802708379837148646191567765584039175249171110593159305029014037881475265618958103073425958633163441030267478942720703134493880117805010891574606323700178176718412858948243785754898788359757528163558061136758276299059029113119763557411729353915848889261125855717014320045292143759177464380434854573300054940683350937992500211758727939459249163046465047204851616590276724564411037216844005877918224201569391107769029955591465502737961776799311859881060956465198859727495735498887960494256488224613682478900505821893815926193600121890632"
    from decimal import Decimal
    d = Decimal(ns)
    one_third = Decimal("0.3333333333333333")
    print d ** one_third
    

    答案是:2.287391878618402702753613056E+305

    TZ 指出这不准确……他是对的。这是我的测试。

    from decimal import Decimal
    
    def nth_root(num_decimal, n_integer):
        exponent = Decimal("1.0") / Decimal(n_integer)
        return num_decimal ** exponent
    
    def test():
        ns = "11968003966030964356885611480383408833172346450467339251196093144141045683463085291115677488411620264826942334897996389485046262847265769280883237649461122479734279424416861834396522819159219215308460065265520143082728303864638821979329804885526557893649662037092457130509980883789368448042961108430809620626059287437887495827369474189818588006905358793385574832590121472680866521970802708379837148646191567765584039175249171110593159305029014037881475265618958103073425958633163441030267478942720703134493880117805010891574606323700178176718412858948243785754898788359757528163558061136758276299059029113119763557411729353915848889261125855717014320045292143759177464380434854573300054940683350937992500211758727939459249163046465047204851616590276724564411037216844005877918224201569391107769029955591465502737961776799311859881060956465198859727495735498887960494256488224613682478900505821893815926193600121890632"
        nd = Decimal(ns)
        cube_root = nth_root(nd, 3)
        print (cube_root ** Decimal("3.0")) - nd
    
    if __name__ == "__main__":
        test()
    

    大约差 10**891

    【讨论】:

    • 嗯。这可能有效,但并不准确。使用您的术语,如果答案等于 done_third,那么 (answer3 - d) 应该是多少?
    • 十进制与您需要的一样准确...我的 0.333 字符串只是为了简洁。
    • tz...你是对的。它的方式......哦,好吧。上面的牛顿方法确实很摇滚!
    【解决方案6】:

    可能出于您的好奇心:

    http://en.wikipedia.org/wiki/Hensel_Lifting

    这可能是 Maple 用来实际找到大数的第 n 个根的技术。

    摆出x^n - 11968003.... = 0 mod p的事实,然后从那里开始......

    【讨论】:

      【解决方案7】:

      我想出了自己的答案,它借鉴了@Mahmoud Kassem 的想法,简化了代码,并使其更可重用:

      def cube_root(x):
          return decimal.Decimal(x) ** (decimal.Decimal(1) / decimal.Decimal(3))
      

      我在 Python 3.5.1 和 Python 2.7.8 中测试过,似乎运行良好。

      结果将具有与函数运行的十进制上下文指定的位数相同的位数,默认情况下为 28 位小数。根据decimal 模块中power 函数的文档,“The result is well-defined but only “almost always correctly-rounded”.”。如果需要更准确的结果,可以按如下方式进行:

      with decimal.localcontext() as context:
          context.prec = 50
          print(cube_root(42))
      

      【讨论】:

        【解决方案8】:

        在旧版本的 Python 中,1/3 等于 0。在 Python 3.0 中,1/3 等于 0.33333333333(而1//3 等于 0)。

        因此,要么将代码更改为使用 1/3.0,要么切换到 Python 3.0。

        【讨论】:

        • 据我所知,没有迹象表明它使用了 Python 2。
        • @SolomonUckoI 这个问题和答案都是在 python 3 发布后仅 1 周发布的。
        【解决方案9】:

        我可能会建议四种方法来解决您的任务。首先是基于Binary Search。二是基于Newton's Method。第三是基于Shifting n-th Root Algorithm。第四种是我在picture here中描述的Chord-Tangent方法。

        二分搜索已经在上面的许多答案中实现。在这里我只是介绍一下我自己对它的看法和它的实现。

        作为替代方案,我还实现了优化二进制搜索方法(标记为Opt)。这个方法只是从[hi / 2, hi) 范围开始,如果我们计算第 k 个根,hi 等于 2^(num_bit_length / k)。

        牛顿法在这里是新的,因为我看到它没有在其他答案中实现。它通常被认为比二分搜索更快,尽管我自己在下面的代码中的计时没有显示任何加速。因此,这里的方法仅供参考/兴趣。

        移位方法比优化的二分搜索方法快 30-50%,如果用 C++ 实现应该更快,因为 C++ 具有快速的 64 位算术,该方法部分使用了该方法。

        Chord-Tangent Method 是我on piece of paper 发明的,它受到启发,是对牛顿法的改进。基本上我画了一个Chord 和一个Tangent Line 并找到与水平线y = n 的交点,这两个交点形成了根解(x0, n) 位置的下限和上限近似,其中n = x0 ^ k。这种方法似乎是最快的,而所有其他方法都进行了超过 2000 次迭代,对于 8192 位数字的情况,这种方法只进行了 8 次迭代。所以这种方法比以前的(按速度)换档方法快 200-300x 倍。

        作为示例,我生成了 8192 位大小的非常大的随机整数。并用这两种方法测量求三次根的时间。

        在test() 函数中你可以看到我通过k = 3 作为根的幂(立方根),你可以通过任何幂而不是3。

        Try it online!

        def binary_search(begin, end, f, *, niter = [0]):
            while begin < end:
                niter[0] += 1
                mid = (begin + end) >> 1
                if f(mid):
                    begin = mid + 1
                else:
                    end = mid
            return begin
        
        def binary_search_kth_root(n, k, *, verbose = False):
            # https://en.wikipedia.org/wiki/Binary_search_algorithm
            niter = [0]
            res = binary_search(0, n + 1, lambda root: root ** k < n, niter = niter)
            if verbose:
                print('Binary Search iterations:', niter[0])
            return res
        
        def binary_search_opt_kth_root(n, k, *, verbose = False):
            # https://en.wikipedia.org/wiki/Binary_search_algorithm
            niter = [0]
            hi = 1 << (n.bit_length() // k - 1)
            while hi ** k <= n:
                niter[0] += 1
                hi <<= 1
            res = binary_search(hi >> 1, hi, lambda root: root ** k < n, niter = niter)
            if verbose:
                print('Binary Search Opt iterations:', niter[0])
            return res
        
        def newton_kth_root(n, k, *, verbose = False):
            # https://en.wikipedia.org/wiki/Newton%27s_method
            f = lambda x: x ** k - n
            df = lambda x: k * x ** (k - 1)
            x, px, niter = n, 2 * n, [0]
            while abs(px - x) > 1:
                niter[0] += 1
                px = x
                x -= f(x) // df(x)
            if verbose:
                print('Newton Method iterations:', niter[0])
            mini, minv = None, None
            for i in range(-2, 3):
                v = abs(f(x + i))
                if minv is None or v < minv:
                    mini, minv = i, v
            return x + mini
        
        def shifting_kth_root(n, k, *, verbose = False):
            # https://en.wikipedia.org/wiki/Shifting_nth_root_algorithm
            B_bits = 64
            
            r, y = 0, 0
            B = 1 << B_bits
            Bk_bits = B_bits * k
            Bk_mask = (1 << Bk_bits) - 1
            niter = [0]
            
            for i in range((n.bit_length() + Bk_bits - 1) // Bk_bits - 1, -1, -1):
                alpha = (n >> (i * Bk_bits)) & Bk_mask
                B_y = y << B_bits
                Bk_yk = (y ** k) << Bk_bits
                Bk_r_alpha = (r << Bk_bits) + alpha
                Bk_yk_Bk_r_alpha = Bk_yk + Bk_r_alpha
                beta = binary_search(1, B, lambda beta: (B_y + beta) ** k <= Bk_yk_Bk_r_alpha, niter = niter) - 1
                y, r = B_y + beta, Bk_r_alpha - ((B_y + beta) ** k - Bk_yk)
        
            if verbose:
                print('Shifting Method iterations:', niter[0])
        
            return y
        
        def chord_tangent_kth_root(n, k, *, verbose = False):
            niter = [0]
            hi = 1 << (n.bit_length() // k - 1)
            while hi ** k <= n:
                niter[0] += 1
                hi <<= 1
            f = lambda x: x ** k
            df = lambda x: k * x ** (k - 1)
            # https://i.stack.imgur.com/et9O0.jpg
            x_begin, x_end = hi >> 1, hi
            y_begin, y_end = f(x_begin), f(x_end)
            for icycle in range(1 << 30):
                if x_end - x_begin <= 1:
                    break
                niter[0] += 1
                if 0: # Do Binary Search step if needed
                    x_mid = (x_begin + x_end) >> 1
                    y_mid = f(x_mid)
                    if y_mid > n:
                        x_end, y_end = x_mid, y_mid
                    else:
                        x_begin, y_begin = x_mid, y_mid
                # (y_end - y_begin) / (x_end - x_begin) = (n - y_begin) / (x_n - x_begin) ->
                x_n = x_begin + (n - y_begin) * (x_end - x_begin) // (y_end - y_begin)
                y_n = f(x_n)
                tangent_x = x_n + (n - y_n) // df(x_n) + 1
                
                chord_x = x_n + (n - y_n) * (x_end - x_n) // (y_end - y_n)
                
                assert chord_x <= tangent_x, (chord_x, tangent_x)
                x_begin, x_end = chord_x, tangent_x
                y_begin, y_end = f(x_begin), f(x_end)
                assert y_begin <= n, (chord_x, y_begin, n, n - y_begin)
                assert y_end > n, (icycle, tangent_x - binary_search_kth_root(n, k), y_end, n, y_end - n)
            if verbose:
                print('Chord Tangent Method iterations:', niter[0])
            return x_begin
        
        def test():
            import random, timeit
            
            nruns = 3
            bits = 8192
            n = random.randrange(1 << (bits - 1), 1 << bits)
            
            a = binary_search_kth_root(n, 3, verbose = True)
            b = binary_search_opt_kth_root(n, 3, verbose = True)
            c = newton_kth_root(n, 3, verbose = True)
            d = shifting_kth_root(n, 3, verbose = True)
            e = chord_tangent_kth_root(n, 3, verbose = True)
            assert abs(a - b) <= 0 and abs(a - c) <= 1 and abs(a - d) <= 1 and abs(a - e) <= 1, (a - b, a - c, a - d, a - e)
        
            print()
            print('Binary Search timing:', round(timeit.timeit(lambda: binary_search_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
            print('Binary Search Opt timing:', round(timeit.timeit(lambda: binary_search_opt_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
            print('Newton Method timing:', round(timeit.timeit(lambda: newton_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
            print('Shifting Method timing:', round(timeit.timeit(lambda: shifting_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
            print('Chord Tangent Method timing:', round(timeit.timeit(lambda: chord_tangent_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
        
        if __name__ == '__main__':
            test()
        

        输出:

        Binary Search iterations: 8192
        Binary Search Opt iterations: 2732
        Newton Method iterations: 9348
        Shifting Method iterations: 2752
        Chord Tangent Method iterations: 8
        
        Binary Search timing: 0.506 sec
        Binary Search Opt timing: 0.05 sec
        Newton Method timing: 2.09 sec
        Shifting Method timing: 0.03 sec
        Chord Tangent Method timing: 0.001 sec
        

        【讨论】:

          【解决方案10】:

          尝试将指数转换为浮点数,因为 / 在 Python 中的默认行为是整数除法

          n**(1/float(3))

          【讨论】:

          • n**(1.0/3) 也可以完成这项工作
          • OverflowError: long int too large to convert to float, 它只是不起作用导致无法转换的巨大数字
          • 据我所知,没有迹象表明它使用了 Python 2。
          【解决方案11】:

          好吧,如果你不是特别担心精度,你可以把它转换成一个字符串,去掉一些数字,使用指数函数,然后将结果乘以你切除的根数。

          例如32123约等于32*1000,立方根约等于32*1000的立方根。后者可以用0的个数除以3来计算。

          这避免了使用扩展模块的需要。

          【讨论】:

          • 我担心精​​度,因为我知道这个数字,例如是另一个整数的三次方。 (当然我需要知道这个整数):)
          • 另外,如果数字很大,转换为字符串可能会出现问题。
          • @Attila:不在他在帖子中给出的巨大定义中。此外,不确定他接受输入的确切格式,但截断的基本思想适用于大多数格式。
          • 他有工作代码,只是问题是精度
          猜你喜欢
          • 1970-01-01
          • 2012-02-08
          • 2014-03-24
          • 1970-01-01
          • 2019-02-24
          • 1970-01-01
          • 2014-12-08
          • 1970-01-01
          • 1970-01-01
          相关资源
          最近更新 更多