122. 高效指数计算(Efficient exponentiation)

计算\(n^{15}\)最朴素的方式需要\(14\)次乘法,即\(n\times n\times n\cdots\times n=n^{15}\),但使用一种二进制算法,可以只使用\(6\)次乘法:

$$ \begin{array}{ l } n\times n=n^{2}\\ n^{2} \times n^{2} =n^{4}\\ n^{4} \times n^{4} =n^{8}\\ n^{8} \times n^{4} =n^{12}\\ n^{12} \times n^{2} =n^{14}\\ n^{14} \times n=n^{15} \end{array} $$
然而,只使用\(5\)次乘法也是可以的:
$$ \begin{array}{ l } n\times n=n^{2}\\ n^{2} \times n =n^{3}\\ n^{3} \times n^{3} =n^{6}\\ n^{6} \times n^{6} =n^{12}\\ n^{12} \times n^{3} =n^{15} \end{array} $$
\(m(k)\)为计算\(n^k\)所需要的最少乘法次数,比如\(m(15)=5\),对于\(1\le k\le200\),求\(\sum{m(k)}\)

分析:学习过一些算法知识的同学可能都知道可以用二进制的方法快速计算大指数的幂,但正如题目中说明的,二进制算法并不一定是使用乘法次数最少的算法,寻找使用乘法次数最少的问题通常称之为Addition-chain exponentiation问题,这个问题比它初看上去要困难的多,事实上它已经被证明为一个NP完全问题,目前并不存在一个有效的算法,使得我们能够计算出任意幂的最小乘法次数。

在高德纳所著的《计算机程序设计艺术第二卷:半数值算法》第4.6.3节专门讨论了求幂值的加法链算法问题,事实上他给出了一个计算最少乘法次数的递归公式,虽然这个公式只对一千以内的数字有效,但因为我们题目的范围只到两百,所以我们可以使用这个公式,非常快的计算出两百以内每个自然数幂的最少乘法次数。

假设\(m(k)\)表示计算\(n^k\)所需的最少乘法次数,则有:

$$ m(k)=min(m(k-1)+1,m_k)-\delta_k $$

其中\(m_k\)是与\(k\)有关的分段函数:

$$ m_{k} =\begin{cases} +\infty & k\ is\ prime\\ m( p) +m( k/p) & else \end{cases} $$

这里\(p\)表示\(k\)的最小素因子。另\(\delta_k\)也是与\(k\)有关的函数,当\(k\)等于某些特殊的数值时等于一,而在取其它值时则等于零,这些特殊值可见下表:

微信截图_20210716160401

在这个题目中,我们只需要关心两百以内的特殊值。使用以上两个函数,我们就可以用递归公式计算出两百以内的指数幂的最小乘法次数了,代码如下:

# time cost = 15.1 µs ± 12.6 ns

from sympy import isprime
from functools import lru_cache

def p(n):
    for i in range(1,int(n**0.5)+1):
        if n % i == 0 and isprime(i):
            return i

def sigma(n):
    arr = [23,43,59,77,83,107,149,163,165,179]
    if n in arr:
        return 1
    else:
        return 0

@lru_cache(maxsize=None)
def m(k):
    if k == 1:
        return 0
    elif k == 2:
        return 1
    elif isprime(k):
        return m(k-1)+1-sigma(k)
    else:
        return min(m(k-1)+1,m(p(k))+m(k//p(k)))-sigma(k)

def main(N=200):
    res = [m(k) for k in range(1,N+1)]
    return sum(res)