如何解释计算数的幂的算法?

6

我目前正在阅读斯基纳的《算法设计手册》。

他描述了一种计算数字幂的算法,即计算a^n

他开始时说,最简单的算法就是a*a*a ... *a,因此我们总共需要进行n-1次计算。

然后他继续说,有一种优化方法,要求我们认识到:

n = n/2 + n/2

而且我们也可以说

a^n = ((a^n/2)^2)  (a to the n equals a to the n over 2, squared)

我目前理解的是,从这些方程中,他推导出了一种只需要 O(lg n) 次乘法的算法。

function power(a, n)
  if (n = 0) 
    return(1)

  x = power(a,n/2)

  if (n is even)       
    return(x^2)
  else             
    return(a*x^2)

看起来 x 必须是当前计算出来的值。但是,即使我读了好几遍,我也不明白他是如何从这些方程式中设计出这个算法,甚至是它是如何工作的。有人能解释一下吗?


这是一个递归调用函数power的过程。 - motiur
7个回答

21

这个概念很简单。例如,计算38的值。

你可以使用显然的方法,即 38 = 3 x 3 x 3 x 3 x 3 x 3 x 3 x 3,需要进行7次乘法操作。或者有更好的方法。

假设

  • 如果我们知道34的值,我们只需要进行一次乘法运算就可以计算出38,但我们不知道34的值。
  • 如果我们知道32的值,我们可以在一次乘法运算中计算出34
  • 最后,我们知道32可以轻松地通过3×3来计算。

反推回去,只需要进行3次乘法运算就可以计算出38而不是7次

以下是此过程的清晰视图:

   32 = 3 x 3 = 9
   34 = 32 x 32 = 9 x 9 = 81
   38 = 34 x 34 = 81 x 81 = 6,561

然后,还有另一个问题,如果幂是奇数怎么办。例如:39,该怎么处理?你可以这样做

   39 = 3 x 38   或者
   39 = 3 x 34 x 34

让我们称不断相乘的算法为方法A,称不断将幂除以二的算法为方法B。相比于方法B,方法A有多好?对于像3的8次方这样的小数来说,没有太大的显着改进,即使我们最小化了乘法的数量,但我们也略微增加了除法操作的数量。
所以对于3的8次方,
               乘法      除法
   方法A:        7          0
   方法B:        3          3
然而,对于更大的幂,例如:3的4,294,967,296次方,
               乘法              除法
   方法A:   4,294,967,295          0
   方法B:       32                 32
差异是巨大的。

1
我喜欢你的解释,但我已经接受了mcleod的答案,因为我认为他循序渐进地介绍了基础知识。但现在我已经理解了算法,阅读你的答案非常有帮助。 - patchwork

11

首先,让我们修正您的算法:

function power( a, n )
    if (n = 0) 
        return(1)

    x = power(a,n/2)

    if (n is even) 
        return(x*x)
    else 
        return(a*x*x)

假设你想计算power(2,8),也就是2^8(这里的^当然不是XOR)。

  • 1)(a=2n=8)。 由此得2^8 = (2^4)^2,因此我们要计算x=2^4,然后得到x*x的最终结果。 我们需要再次调用power()来获得xpower(2,4)

  • 2)(a=2n=4)。 由此得2^4 = (2^2)^2,因此我们要计算x=2^2,然后得到x*x的最终结果。 我们需要再次调用power()来获得xpower(2,2)

  • 3)(a=2n=2)。 由此得2^2 = (2^1)^2,因此我们要计算x=2^1,然后得到x*x的最终结果。 我们需要再次调用power()来获得xpower(2,1)

  • 4)(a=2n=1)。 由此得2^1 = (2^0)^2,因此我们要计算x=2^0,然后得到a*x*x的最终结果。 我们需要再次调用power()来获得xpower(2,0)

  • 5)(a=2n=0)。 由于n0,因此2^0 = 1,因此x的值将被返回到步骤#4。

  • 4)(a=2n=1x=1)。 此步骤的最终结果为a*x*x = 2*1*1=2,这是要在步骤#3中分配给x的值。

  • 3)(a=2n=2x=2)。 此步骤的最终结果为x*x = 2*2 = 4。 这是要在步骤#2中分配给x的结果。

  • 2)(a=2n=4x=4)。 此步骤的最终结果为x*x = 4*4 = 16。 这是要在步骤#1中分配给x的结果。

  • 1)(a=2n=8x=16)。 此步骤的最终结果为x*x = 16*16 = 256。 这是由power(2,8)返回的结果。


0
x = power(a, n/2)

将会给你 a^n/2。如果整个语句被平方,得到 (a^n/2)^2。现在如果 n 是奇数,在 n/2 期间,a^1 的幂次将会丢失,所以为了恢复它,需要乘以 a。这是根据所给出的方程。


“n” 怎么可能变成0呢?我对这个算法何时终止感到困惑。 - patchwork
@patchwork:在你的算法中,n / 2 是整数除法,而不是浮点除法。因此,5/2 是2余1,而不是2.5。最终,通过对任何整数进行整数除法,您将达到0。(一些语言使用特殊语法,如双斜杠表示整数除法:n // 2。类型化语言将两个整数类型的除法视为整数除法。)这种整数除以2可能剩下0和1的两个余数是返回幂时分支奇偶数的原因。 - M Oehm

0

函数power(...)的编写方式似乎是为了处理整数除法的影响。请记住,在整数除法规则下,小数部分被舍弃。这只会影响奇数,因为偶数除以二不产生余数。

因此,每当n是一个偶数时,计算出的n/2的值恰好等于n/2,可以采取if的顶部分支:这正是方程所规定的。

每当n是一个奇数时,计算出的n/2的值实际上等于floor(n/2)。换句话说,该语句:

x = power(a,n/2)

已经计算出结果

x = power(a,(n-1)/2)

这意味着您的指数永久损失了1,仅返回x^2将缺少一个a的幂。因此,底部分支会添加回a的幂。

如果您想象一下,计算机不是整数除法,而是能够完全无损实数除法,那么您可以将函数重写为:

function power( a, n )

    if (n = 0) 
        return(1)

    x = power(a,n/2)

    return x^2

你可以轻松地满足自己,这是问题中提供的第二个方程的简单递归实现。


0

这个函数是递归的,也就是说,当被调用时,它会一遍又一遍地调用自己,直到满足某些最终条件为止。在这种情况下,有三个最终条件将停止函数调用自身并返回结果。

如果我是你,我会尝试手动应用算法到几个不同的值上,以便理解。


0

这个公式 a^n = ((a^n/2)^2),你理解的话就是递归算法。

要得到 a^n,你需要先计算 a^(n/2)
要得到 a^(n/2),你需要计算 a^((n/2)/2)
... 以此类推,直到 (n/2/2/...2) 等于 0,这是递归终止条件。

因此,递归算法完全遵循该公式:

要得到 power(a,n),你首先递归计算 power(a,n/2),然后根据 n 是奇数还是偶数来调整返回结果。

关于这种实现方法,也有wikipedia文章。


连续除以2会得到0吗? - patchwork
2
是的,总是这样。当你除以2时,只取除法的整数部分。所以,3 / 2 = 1,而1 / 2 = 0; - kiruwka
如果您考虑到(n/2/2/2.../2)的极限值,那么终止值为1。 - chepner
1
@chepner,不是根据这一行代码:if (n = 0) return(1) - kiruwka
“n/2”永远不会等于零。这是函数中的一个错误/打字错误。当“n”为奇数时,需要减1而不是除以2(我现在意识到,这确实会导致终止条件为0)。 - chepner

0
看到整个操作链作为一行代码(你可以直接用bc运行)确实有助于理解算法在每个阶段要做的事情:
  • 计算星号(*)的总数,即乘法的总数

计算插入符号(^)的个数,以得到总共的平方次数。
  29 ^ 1847  
| expn LSB bitstring := 11101100111
| 29*(29*(29*((29*(29*(((29*(29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 4373  
| expn LSB bitstring := 1010100010001
| 29*((29*((29*((((29*((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 8563  
| expn LSB bitstring := 11001110100001
| 29*(29*(((29*(29*(29*((29*(((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 8689  
| expn LSB bitstring := 10001111100001
| 29*((((29*(29*(29*(29*(29*(((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 9067  
| expn LSB bitstring := 11010110110001
| 29*(29*((29*((29*(29*((29*(29*((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 11197 
| expn LSB bitstring := 10111101110101
| 29*((29*(29*(29*(29*((29*(29*(29*((29*((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 14813 
| expn LSB bitstring := 10111011100111
| 29*((29*(29*(29*((29*(29*(29*(((29*(29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 15749 
| expn LSB bitstring := 10100001101111
| 29*((29*(((((29*(29*((29*(29*(29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 16703 
| expn LSB bitstring := 111111001000001
| 29*(29*(29*(29*(29*(29*(((29*((((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 18919 
| expn LSB bitstring := 111001111001001
| 29*(29*(29*(((29*(29*(29*(29*(((29*(((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 24071 
| expn LSB bitstring := 111000000111101
| 29*(29*(29*(((((((29*(29*(29*(29*((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 24371 
| expn LSB bitstring := 110011001111101
| 29*(29*(((29*(29*(((29*(29*(29*(29*(29*((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 25933 
| expn LSB bitstring := 101100101010011
| 29*((29*(29*(((29*((29*((29*(((29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 26003 
| expn LSB bitstring := 110010011010011
| 29*(29*(((29*(((29*(29*((29*(((29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 27851 
| expn LSB bitstring := 110100110011011
| 29*(29*((29*(((29*(29*(((29*(29*((29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 28097 
| expn LSB bitstring := 100000111011011
| 29*((((((29*(29*(29*((29*(29*((29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 28537 
| expn LSB bitstring := 100111101111011
| 29*(((29*(29*(29*(29*((29*(29*(29*(29*((29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 28603 
| expn LSB bitstring := 110111011111011
| 29*(29*((29*(29*(29*((29*(29*(29*(29*(29*((29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 32869 
| expn LSB bitstring := 1010011000000001
| 29*((29*(((29*(29*(((((((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 33013 
| expn LSB bitstring := 1010111100000001
| 29*((29*((29*(29*(29*(29*((((((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 33107 
| expn LSB bitstring := 1100101010000001
| 29*(29*(((29*((29*((29*(((((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 33529 
| expn LSB bitstring := 1001111101000001
| 29*(((29*(29*(29*(29*(29*((29*((((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 33629 
| expn LSB bitstring := 1011101011000001
| 29*((29*(29*(29*((29*((29*(29*((((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 33809 
| expn LSB bitstring := 1000100000100001
| 29*((((29*((((((29*(((((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 36899 
| expn LSB bitstring := 1100010000001001
| 29*(29*((((29*(((((((29*(((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 38651 
| expn LSB bitstring := 1101111101101001
| 29*(29*((29*(29*(29*(29*(29*((29*(29*((29*(((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 38839 
| expn LSB bitstring := 1110110111101001
| 29*(29*(29*((29*(29*((29*(29*(29*(29*((29*(((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 46093 
| expn LSB bitstring := 1011000000101101
| 29*((29*(29*(((((((29*((29*(29*((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 47279 
| expn LSB bitstring := 1111010100011101
| 29*(29*(29*(29*((29*((29*((((29*(29*(29*((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 47339 
| expn LSB bitstring := 1101011100011101
| 29*(29*((29*((29*(29*(29*((((29*(29*(29*((29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 52183 
| expn LSB bitstring := 1110101111010011
| 29*(29*(29*((29*((29*(29*(29*(29*((29*(((29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 52697 
| expn LSB bitstring := 1001101110110011
| 29*(((29*(29*((29*(29*(29*((29*(29*(((29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 53569 
| expn LSB bitstring := 1000001010001011
| 29*((((((29*((29*((((29*((29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 59471 
| expn LSB bitstring := 1111001000010111
| 29*(29*(29*(29*(((29*(((((29*((29*(29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 63029 
| expn LSB bitstring := 1010110001101111
| 29*((29*((29*(29*((((29*(29*((29*(29*(29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 63667 
| expn LSB bitstring := 1100110100011111
| 29*(29*(((29*(29*((29*((((29*(29*(29*(29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2


  29 ^ 65287 
| expn LSB bitstring := 1110000011111111
| 29*(29*(29*((((((29*(29*(29*(29*(29*(29*(29*(29)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2

这段脚本代码可以自动生成。
gawk 'function ____(__, ___, _) {

    if (_ == "") {

        if ((__ *= (___ = length(__ = (_)__)) \
            * _++ < (__ = int(__))) < ++_)
            return __

        _ ^= (___ += __ + __ < (_+_)^_^_ \
                    ? _ * --___ \
                    : _ * ___ - _^(__ < (_*_^_^_)^_))
        if (__ < _)
            ___ -= __ < (__ += __)
        if (_+_ <= __)
            while ((_+=_+!++___)+_ <= __);

        return (_<__ ? ____(__ -= _, ___, _) \
                   : sprintf("%.*d", ___, !_) !!_
    }
    return !--___ ? !!+__ : (__+=__)<_ ? ____(__, ___, _) !_ \
                                       : ____(__-=_,___,_) !!_
}
BEGIN {
    OFS = ")^2"
} _ = NF {

    printf("\f %11s^%-5s \f\r\t\t| expn LSB bitstring := %*s\f\r\t\t| ",
            __ = $!!_, $!_ = $_, _ = !_, $_ = ____($_)) > ("/dev/stderr")

    print ($!(NF = gsub(/[^01]+/, "")^_ * gsub(_, "(") + gsub(!_,
        __ "*(") * sub(/[^0-9]*$/, "") \
                 * sub(/^$/, !_)^_))) "\n") > ("/dev/stderr") } NF'

========

如果你好奇的话,最小的奇质数的幂与能够由64位双精度类型表示的最大质数之间的关系可以通过以下一系列操作来表示:
       3 ^ 9007199254740881 
    | expn LSB bitstring := 100010011111111111111111
                            11111111111111111111111111111
3*(  (  (  (3*(  (  (3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*
  (3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*
  (3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3*(3
   \
    )^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2 
    )^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2
    )^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2)^2

而且只需要一个1785 TB的驱动器,就可以存储4297万亿位的小数表示。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接