快速幂运算的应用

来源：互联网发布：化工单元仿真软件编辑：程序博客网时间：2024/05/03 03:44

快速幂取余的方法，所得返回值为 a^b % p

int f(int a,int b,int p){    int temp = 1;    while(b)    {        if(b & 0x01)        {            temp = (temp * (a%p)) % p;        }        a = ( (a%p) * (a%p) ) % p;        b >>= 1;    }    return temp;}

例如：a^11，设b = 11，写成2进制为1011。

则a^11可以表示为 a^( 1*2^3 + 0*2^2 + 1*2^1 + 1*2^0 ) = a^(2^3) *a^(2^1) * a^(2^0)

二进制值为1的那项，在a的指数中就存在a^(2^i)，

所以可以把a^(2^0),a^(2^1),a^(2^2)......求出来，

如果在b的二进制表示中第i位为1，就把a^(2^i)乘到最终保存结果的变量中去。

而求a^(2^0),a^(2^1),a^(2^2)......a^(2^i)......时，

a^(2^ (i+1) ) = a^( 2^i * 2 ) = (a^(2i))^2

所以如果知道了a^(2^0) ，对它进行平方就可以得到a^(2^1)，以此类推，

所以在程序中我们就不断的将a平方，再将b向右移位。

Q1：Compute (a^(b^c)) % 317000011 ( 1 <= a,b,c <= 100000)

Tips:

⒈ 费马小定理: 如果p为素数，则a^(p-1) % p = 1;

⒉ 除法定理: 对任意整数a和任意正整数n，存在唯一的整数q和r，满足 0 <= r < n ，并且 a = qn + r;

解：设b^c = n(p-1) + d，则a^(b^c) = a^(n(p-1)) * a^(d),依据定理⒈ ，a^(b^c) = a^(d),

而d = b^c Mod (p-1)。所以 a^(b^c) % p = a^( b^c Mod (p-1) ) Mod p 。

#include <stdio.h>#define P 317000011unsigned long long f(unsigned long long a,unsigned long long b,unsigned long long p){    unsigned long long temp = 1uLL;    while(b)    {        if(b & 0x01uLL)        {            temp = ((temp%p) * (a%p)) % p;        }        a = ((a%p)*(a%p))%p;        b >>= 1uLL;    }    return temp;}int main(){    int t;    unsigned long long a,b,c;    scanf("%d",&t);    while(t--)    {        scanf("%llu %llu %llu",&a,&b,&c);        printf("%llu\n",f(a,f(b,c,P-1),P));    }    return 0;}

Q2：矩阵幂快速求Fibonacci数列的第 N 项

定义一个基础矩阵 base

┌ ┐ ┌ ┐

│ 1 1 │ = │F2 F1│

│ 1 0 │ │F1 F0│

└ ┘ └ ┘

┌ ┐

求base^n = │Fn+1 Fn │

│Fn Fn-1│

└ ┘

所以矩阵的快速幂乘方方法，可以高效的求的Fibonacci数列，O(log2n)。

int f(int n){    if(n == 0)  return 0;    else if(n == 1) return 1;    else if(n == 2) return 1;    else    {        int base[2][2] = {{1,1},{1,0}};        int b = n-1;        int res[2][2] = {{1,1},{1,1}};        int temp[2][2];        while(b)        {            if(b & 0x01)            {                temp[0][0] = res[0][0] * base[0][0] + res[0][1] * base[1][0];                temp[0][1] = res[0][0] * base[0][1] + res[0][1] * base[1][1];                temp[1][0] = res[1][0] * base[0][0] + res[1][1] * base[1][0];                temp[1][1] = res[1][0] * base[0][1] + res[1][1] * base[1][1];                res[0][0] = temp[0][0]; res[0][1] = temp[0][1];                res[1][0] = temp[1][0]; res[1][1] = temp[1][1];            }            temp[0][0] = base[0][0] * base[0][0] + base[0][1] * base[1][0];            temp[0][1] = base[0][0] * base[0][1] + base[0][1] * base[1][1];            temp[1][0] = base[1][0] * base[0][0] + base[1][1] * base[1][0];            temp[1][1] = base[1][0] * base[0][1] + base[1][1] * base[1][1];            base[0][0] = temp[0][0]; base[0][1] = temp[0][1];            base[1][0] = temp[1][0]; base[1][1] = temp[1][1];            b >>= 1;        }        return res[0][0];    }}