快速幂及其应用

　　计算指数函数a^p（p>0）在p点的值，朴素算法将a累乘p次得出结果，时间复杂度为O(p)，当指数p很大时（>1e7）算法在机器上需要花大量运行时间。

　　改进朴素算法，采用分治策略：当p为偶数时，将a^p分解为(a^p/2)*(a^p/2)两因子，而a^p/2又可继续做同样的划分，直到降为容易计算的低次幂；当p为奇数时，从式中提出因子a，剩下的a^p-1就化为了偶指数形式。

　　形式化描述为：

egin{equation}
a(p) = left{
egin{array}{**l**}
a(frac{p}{2}) cdot a(frac{p}{2}) & p=2k, k in N^+ \
a(frac{p-1}{2}) cdot a(frac{p-1}{2}) cdot a & p=2k+1, k in N \
1 & p=0
end{array}
ight.
end{equation}

　　设T(p)为算法在规模p下花费的时间，则根据递推方程有：

egin{equation}
T(p) = left{
egin{array}{**l**}
T(lfloor frac{p}{2} floor) + O(1) & p>0 \
O(1) & p=0
end{array}
ight.
end{equation}

　　根据Master Theory可得$ T(p) = O(log{p}) $

　　可见时间复杂度有了很大改善。直观地体会一下，例如指数p=18,446,744,073,709,551,616（2⁶⁴）时，只需执行64次乘法。(这里不考虑结果的溢出）

　　如何实现这个算法呢？朴素的做法是递归:

typedef long long ll;

ll fspow(ll a, ll p) {
    if(p==0)
        return 1;
    if(p&1) {
        ll t = fspow(a, p-1);
        return t*a;
    } else {
        ll t = fspow(a, p/2);
        return t*t;
    }
}

 

两行代码的非递归实现

　　递归算法缺陷在于空间复杂度为O(log p)，非递归算法避免了额外的存储空间。

for (s=1; p;p/=2,a=a*a)
        if (p&1) s=s*a;

　　其中，

输入变量：a为底数，p为指数

输出变量：s=a^p

　　再来看原理。思路仍是分治，这与我们之前的讨论是一致的，但不再递归地求解。

　　对于偶指数，每次划分指数折半（1/2），此时要确保等价，底数必须平方，即$ a^p = (a^2)^{frac{p}{2}} $。而指数p/2又可继续分为p/4 p/8……1，同时底数随之平方。当指数降为1时，底数即为幂的值。当指数不是偶数时，先拆出一个底数，则剩下的幂就是偶次幂，再按如上方法处理即可。

　　注意如下细节：

　　　　1. 当p为奇数时，$ lfloor frac{p}{2} floor = frac{p-1}{2} $。

　　　　2. 利用位运算避免%运算。

具体应用

1. RSA

　　幂运算在许多场合都有着重要的运用。例如RSA非对称加密算法中，明文通过如下公式进行加密：

C_i = M_i^e mod n

　　　　其中Ci为密文，Mi为与明文有关的数值，(e, n)为公钥对。

　　而密文通过如下公式解密：

M_i = C_i^d mod n

　　其中Mi为与明文有关的数值，(d, n)为私钥对。

 

　　可以看到，无论是加密还是解密，都要用到幂运算与模运算。这就与我们上文所述的快速幂运算很有关联。

 

　　不止RSA，在其它运用中也常常出现类型的形式，它们归结下来如下面这个式子：

a^p mod k

　　由于对幂取模，最终结果并不会大于等于k。再看看我们的快速幂算法，中间结果很可能远大于k，如果中间结果超出基本数据类型的范围，还需要通过高精度模拟运算。我们很快发现，当k可通过基本数据类型表示时，似乎没有必要采用高精度算法，上文的快速幂算法还可以改进。

　　随时取模性质指出，(a*b) mod k = (a mod k)*(b mod k) mod k。即两个数的乘积再取模，等于两个数分别取模再相乘，最后取模。利用这个性质，我们可以将a^p mod k等价变形为(a mod k)^p mod k，同时每次做乘法时，都先利用该性质缩小乘数，再相乘。

　　上述代码修改为

a %= k;
for (s=1; p;p/=2,a=(a*a)%k)
         if (p&1) s=(s*a)%k;

　　其中a，p，k，s的含义与上文一致。

　　这便是快速幂模k算法的非递归实现。

2. 矩阵幂加速

快速Fibonacci计算

　　Fibonacci数列的递推定义如下：

　　为了将递推式用线性代数的形式表达出来，设待定系数a,b,c,d。当n>=2时，建立如下等式：

　　解出各未知量，得$ a=1,b=1,c=1,d=0 $，即递推式等价于：

　　注意上式必须满足n>=2。结合n<2的基准情况，可以得到任意$ f_n $、$ f_{n-1} $的表达式如下：

　　即，只需计算出左矩阵的n次幂，即可求出Fibonacci的第n-1和第n项。若采用快速幂算法计算矩阵幂，我们会惊讶地发现求出该数列第n项只需O(log n)的时间复杂度。

　　如何实现矩阵乘法？只需封装矩阵ADT（抽象数据类型）。利用语言特性重载*=运算符。

typedef long long ll;

class mat{
    ll a,b,c,d;
    inline void operator*=(const mat &m) {
        mat n;
        n.a=a*m.a+b*m.c;
        n.b=a*m.b+b*m.d;
        n.c=c*m.a+d*m.c;
        n.d=c*m.b+d*m.d; 
        *this=n;
    }
};

ll fib(int n) {
    mat s={1,0,0,1};
    mat p={1,1,1,0};
    
    for(; n; n/=2,p*=p) {
        if(n&1) s*=p;
    }
    return s.a;
}

　　实际应用中由于数值较大，往往需要对结果取模。可利用随时取模性质，在计算矩阵乘法时对中间结果取模。