求解计算a^b的快速幂函数的循环不变式

函数背景 #

给定的R语言函数是快速幂算法的实现，用于计算实数a的自然数b次幂a^b：

power <- function(a, b){
     c <- 1
     while(b > 0){
          if(b %% 2 == 1){
               c <- c * a
          }
          b <- floor(b / 2)
          a <- a * a
     }
     return c
}

核心循环不变式 #

该循环的关键不变式为：在每次循环迭代的起始时刻，c * a^b = a₀^b₀，其中a₀和b₀是函数调用时传入的初始a和b值。

验证不变式的有效性 #

1. 初始状态
   循环开始前，c=1，a=a₀，b=b₀，代入得：1 * a₀^b₀ = a₀^b₀，等式成立。

2. 迭代过程验证
   假设迭代前的状态满足c_old * a_old^b_old = a₀^b₀，分两种情况验证迭代后的状态：

   • 当b_old为奇数时
     执行c = c_old * a_old，随后b = floor(b_old/2) = (b_old-1)/2，a = a_old^2。
     代入新状态计算：

     c_new * a_new^b_new = (c_old*a_old) * (a_old²)^((b_old-1)/2)
     = c_old*a_old * a_old^(b_old-1)
     = c_old*a_old^b_old
     = a₀^b₀
     

     等式依然成立。

   • 当b_old为偶数时
     不修改c，直接执行b = floor(b_old/2) = b_old/2，a = a_old^2。
     代入新状态计算：

     c_old * a_new^b_new = c_old * (a_old²)^(b_old/2)
     = c_old * a_old^b_old
     = a₀^b₀
     

     等式依然成立。

3. 循环终止时
   当b > 0不成立时，b=0，此时c * a^0 = c * 1 = c，根据不变式可得c = a₀^b₀，与函数返回值一致，符合算法预期。

结合你观察到的规律 #

你提到第k次迭代时a = a₀^(2^k)，这和不变式完全契合：每次迭代a平方一次，k次迭代后确实是初始a的2^k次方；而c则负责累积b二进制表示中为1的位对应的a的幂次，两者共同维持c * a^b = 初始结果的不变关系。

内容的提问来源于stack exchange，提问作者Jeremy