ACM | 算法 | 快速冪

時間  2019-12-09
標籤 acm 算法 快速 欄目 應用數學 简体版
原文   原文鏈接

目錄php

快速冪ios

快速冪取模c++

矩陣快速冪算法

矩陣快速冪取模函數

HDU1005練習this

快速冪

​ 冪運算：\(x ^ n\)spa

​ 根據其通常定義咱們能夠簡單實現其非負整數狀況下的函數code

定義法：ip

int Pow (int x, int n) {
    int result = 1;
    
    while(n--) {
        result *= x;
    }
    
    return result;
}

​ 不難看出此時算法的時間複雜度是\(O(n)\)，一旦n取較大數值，計算時間就會大大增長，極其容易出現超時的狀況。ci

快速冪：

​ 首先要在此列舉兩個前提：

  1. 計算機是經過二進制進行存儲數據的，對二進制數據能夠直接進行操做。

  2. \(2^n+2^n=2*2^n=2^{n + 1}\)

​ 對於\(x ^ n\)，其中n能夠表示爲m位的二進制形式，即\(n=n_12^0+n_22^1+n_32^3+\cdots+n_m2^{m-1}\)

​ 那麼\(x ^ n=x ^ {n_12^0+n_22^1+n_32^3+\cdots+n_m2^{m-1}}\)

​ 即\(x^n=x ^ {n_12^0}x^{n_22^1}x^{n_32^3}\cdots x^{n_m2^{m-1}}\)

​ 根據前提1，計算機能夠直接對二進制格式存儲的數據進行讀取，那麼咱們就能夠對\(n\)一個個讀取再對其進行累乘。

​ 當取到第\(a\)位時，其乘數有通項公式：\(x^{n_a2^{a-1}}\)

​ 經過標準庫math，用代碼實現：

int Pow (int x, int n) {

    int result = 1;

    int a = 1;
    
    while(n) {
        if(n & 1) result *= round( pow(x, 1 << (a - 1)) );//round的做用在於double轉int時防止丟失精度，對1進行位運算是一種計算2的n次冪的快速方法
        n >>= 1;
        
        a++;
    }
    
    return result;
}

​ 但實際上這個調用了標準庫的代碼並無實現快速冪，由於仍然是採用pow()進行的運算

此處由 \(2^n+2^n=2\times2^n=2^{n + 1}\)

即\((x ^ {2 ^ {n}} )^2=x ^ {2\times 2 ^ {n}} =x ^ {2 ^ {n + 1}}\)

所以咱們能夠經過對前項進行二次冪運算獲得後項

先獲得首項\(f(1)=x^{2^{1-1}}=x\)

即令int f = x

具體實現：

int Pow (int x, int n) {

    int result = 1;

    int f = x;
    
    while(n) {
        if(n & 1) result *= f;

        f *= f;
        
        n >>= 1;
    }
    
    return result;
}

不難發現此時算法的時間複雜度由其次數的二進制位數而決定，即\(O(m)\)也就是\(O(log_2n)\)

另外由於此算法經常用於計算大數，因此int類型最好都換成long long類型，防止溢出。

快速冪取模

​ 對\(x^n\)取模運算：\(x^n\mod p\)，當n極大時，一樣其運算量也會增大

​ 這就須要咱們尋求一種快速解決的方法，此時能夠運用咱們上文所提到的快速冪算法

​ 引論1：\((np+a)\mod p=a\mod p\quad (n\in\mathbb{Z} )\)

​ 證實以下：設\(f(n,p,a)=(np+a)\mod p\quad (n\in\mathbb{Z} )\)

​ 則由定義有\((np+a)\mod p\\=\frac{np+d}{p}-([\frac{np+d}{p}+1]-1)\\ =d-p([\frac{d}{p}+1]-1)\)

​ 顯而易見，\((np+a)\mod p=a\)與\(n\)無關

​ 令\(n=0\)得\((np+a)\mod p=a\mod p\quad (n\in\mathbb{Z} )\\ Q.E.D\)

​ 引論2：\((mn)\mod p=((m\mod p)(n\mod p))\mod p\)

​ 證實以下：令\(\left\{ \begin{array}{c} m=(ip+c)& (i\in\mathbb{Z} )\\ n=(jp+d) & (j\in\mathbb{Z} )\end{array} \right.\)

​ 則有\(\left\{ \begin{array}{c} m\mod p=c\\ n\mod p=d \end{array} \right.\)

​ 原式\((m*n)％p\\=((ip+c)(jp+d))\mod p\\=(jip^2+jpc+ipd+cd)\mod p\\=(jip^2+jpc+ipd+cd)\mod p\\=((jip+jc+id)p+cd)\mod p\\由於(jip+jc+id)\in\mathbb{Z}，由引論1得：\\=(cd)\mod p\\=((m\mod p)(n\mod p))\mod p\)

​ 即\((mn)\mod p=((m\mod p)(n\mod p))\mod p\\ Q.E.D\)

​

​ 所以對於\(x^n\mod p\)

​ 能夠寫成\((x ^ {n_12^0}x^{n_22^1}x^{n_32^3}\cdots x^{n_m2^{m-1}})\mod p\\ =((x ^ {n_12^0}\mod p)(x^{n_22^1}\mod p)(x^{n_32^3}\mod p)\cdots (x^{n_m2^{m-1}}\mod p))\mod p\)

​ 而且由以前的推定\((x ^ {2 ^ {n}} )^2=x ^ {2\times 2 ^ {n}} =x ^ {2 ^ {n + 1}}\)

​ 有\((x ^ {2 ^ {n}} \mod p)^2\mod p =(x ^ {2 ^ {n}})^2\mod p=x ^ {2 ^ {n + 1}}\mod p\)

​ 代碼實現：

int Mod (int x, int n, int p) {

    int result = 1;

    int f = x % p;
    
    while(n) {
        if(n & 1) result = (result * f)%p;

        f = (f * f)%p;
        
        n >>= 1;
    }
    
    return result;
}

矩陣快速冪

​ 當\(X\)爲任意方塊矩陣，即\(X=\left| \begin{matrix} x_{11} & \cdots & x_{1n}\\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{nn}\\ \end{matrix} \right|\)時

​ 同理\(X^a\)，有\(X^a=X ^ {a_12^0}X^{a_22^1}X^{a_32^3}\cdots X^{a_m2^{m-1}}\)

​ 一樣的，任意方塊矩陣也適用於快速冪

​ 代碼實現：

#include <iostream>
#include <cstring>

using namespace std;

#define R 2

struct Matrix{
    int data[R][R];
};

Matrix multi(Matrix a,Matrix b,int rec){

    Matrix result;

    memset(&(result.data), 0, sizeof(result.data));

    for(int i = 0; i < rec; i++)
        for(int j = 0; j < rec; j++)
            for(int k = 0; k < rec; k++)
             result.data[i][j] += a.data[i][k] * b.data[k][j];
           
    
    return result;
}

Matrix poww (Matrix x, int n) {

    Matrix result;

    memset(&(result.data), 0,sizeof( result.data));

    for(int i = 0; i < R; i++) result.data[i][i] = 1;

    Matrix f = x;
    
    while(n) {
        if(n & 1) result = multi(result, f, R);

        f = multi(f, f, R);

        n >>= 1;
    }
    
    return result;
}

void MatrixPrint(Matrix target) {

    for(int i = 0; i < R; i++) {

        for(int j = 0; j < R; j++) cout << target.data[i][j]<< " ";

        cout << endl;
    
    }

}

int main() {

    Matrix result;

    memset(&(result.data), 0,sizeof( result.data));

    result.data[0][0] = 1;

    result.data[0][1] = 2;

    result.data[1][0] = 3;

    result.data[1][1] = 4;

    MatrixPrint(result);
    
    cout << endl;

    result =  poww(result, 3);

    MatrixPrint(result);

    return 0;
}

輸出結果：

1 2
3 4

37 54
81 118

矩陣快速冪取模

​ 運算定義：當\(X\)爲任意方塊矩陣，即\(X=\left| \begin{matrix} x_{11} & \cdots & x_{1n}\\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{nn}\\ \end{matrix} \right|\)時

​ 則\(X\mod p=\left| \begin{matrix} x_{11} \mod p & \cdots & x_{1n}\mod p\\ \vdots & \ddots & \vdots \\ x_{n1}\mod p & \cdots & x_{nn}\mod p\\ \end{matrix} \right|\)

​

​ 引論1：\((m+n)\mod p=((m\mod p)+(n\mod p))\mod p\)

​ 證實以下：令\(\left\{ \begin{array}{c} m=(ip+c)& (i\in\mathbb{Z} )\\ n=(jp+d) & (j\in\mathbb{Z} )\end{array} \right.\)

​ 則有\(\left\{ \begin{array}{c} m \mod p=c\\ n \mod p=d \end{array} \right.\)

​ 原式\((m+n)\mod p\\=((ip+c)+(jp+d))\mod p\\=((i+j)p+c+d)\mod p\\由於(i+j)\in\mathbb{Z}，由上文引論得：\\=(c+d)\mod p\\=((m\mod p)+(n\mod p))\mod p\)

​ 即\((m+n)\mod p=((m\mod p)+(n\mod p))\mod p\\ Q.E.D\)

​ 引論2：有方塊矩陣\(X=\left| \begin{matrix} x_{11} & \cdots & x_{1n}\\ \vdots & \ddots & \vdots \\ x_{nn} & \cdots & x_{nn}\\ \end{matrix} \right|\)，\(Y=\left| \begin{matrix} y_{11} & \cdots & y_{1m}\\ \vdots & \ddots & \vdots \\ y_{n1} & \cdots & y_{nm}\\ \end{matrix} \right|\)

​ 則有\(XY\mod p=((X\mod p)·(Y\mod p))\mod p\)

​ 證實以下：

​ 令\(X=\left| \begin{matrix} x_{11} & \cdots & x_{1n}\\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{nn}\\ \end{matrix} \right|\)，\(Y=\left| \begin{matrix} y_{11} & \cdots & y_{1n}\\ \vdots & \ddots & \vdots \\ y_{n1} & \cdots & y_{nn}\\ \end{matrix} \right|\)

​ 則

\(X·Y\mod p=\left| \begin{matrix} (x_{11}y_{11} + \cdots +x_{1n}y_{n1})\mod p & \cdots & (x_{11}y_{1n} + \cdots +x_{1n}y_{nn})\mod p\\ \vdots & \ddots & \vdots \\ (x_{n1}y_{11} + \cdots +x_{nn}y_{n1})\mod p & \cdots & (x_{n1}y_{1n} + \cdots +x_{nn}y_{nn})\mod p\\ \end{matrix} \right|\)

\(=\left| \begin{matrix} ((x_{11}y_{11})\mod p + \cdots +(x_{1n}y_{n1})\mod p )\mod p & \cdots & ((x_{11}y_{1n})\mod p + \cdots +(x_{1n}y_{nn})\mod p\mod p\\ \vdots & \ddots & \vdots \\ ((x_{n1}y_{11})\mod p + \cdots +(x_{nn}y_{n1})\mod p)\mod p & \cdots & ((x_{n1}y_{1n})\mod p + \cdots +(x_{nn}y_{nn})\mod p)\mod p\\ \end{matrix} \right|\)

\(=\left| \begin{matrix} (((x_{11}\mod p)(y_{11}\mod p))\mod p + \cdots +((x_{1n}\mod p)(y_{n1}\mod p))\mod p )\mod p & \cdots & (((x_{11}\mod p)(y_{1n}\mod p))\mod p + \cdots +(((x_{1n}\mod p)(y_{nn}\mod p))\mod p)\mod p\\ \vdots & \ddots & \vdots \\ (((x_{n1}\mod p)(y_{11}\mod p))\mod p + \cdots +((x_{nn}\mod p)(y_{n1}\mod p))\mod p)\mod p & \cdots & (((x_{n1}\mod p)(y_{1n}\mod p))\mod p + \cdots +((x_{nn}\mod p)(y_{nn}\mod p))\mod p)\mod p\\ \end{matrix} \right|\)

\(=\left| \begin{matrix} ((x_{11}\mod p)(y_{11}\mod p) + \cdots +(x_{1n}\mod p)(y_{n1}\mod p))\mod p & \cdots & ((x_{11}\mod p)(y_{1n}\mod p) + \cdots +((x_{1n}\mod p)(y_{nn}\mod p))\mod p\\ \vdots & \ddots & \vdots \\ ((x_{n1}\mod p)(y_{11}\mod p) + \cdots +(x_{nn}\mod p)(y_{n1}\mod p))\mod p & \cdots & ((x_{n1}\mod p)(y_{1n}\mod p) + \cdots +(x_{nn}\mod p)(y_{nn}\mod p))\mod p\\ \end{matrix} \right|\)

\(=\left( \left| \begin{matrix} x_{11}\mod p & \cdots & x_{1n}\mod p\\ \vdots & \ddots & \vdots \\ x_{n1}\mod p & \cdots & x_{nn}\mod p\\ \end{matrix} \right|\left| \begin{matrix} y_{11}\mod p & \cdots & y_{1n}\mod p\\ \vdots & \ddots & \vdots \\ y_{n1}\mod p & \cdots & y_{nn}\mod p\\ \end{matrix} \right| \right)\mod p\)

\(=((X\mod p)·(Y\mod p))\mod p \\\text{即}XY\mod p=((X\mod p)·(Y\mod p))\mod p\\Q.E.D\)

說明任意矩陣也符合快速冪取模的算法

具體實現：

#include <iostream>
#include <cstring>

using namespace std;

#define R 2

struct Matrix{
    int data[R][R];
};

Matrix multi(Matrix a,Matrix b,int rec,int p){

    Matrix result;

    memset(&(result.data), 0,sizeof( result.data));

    for(int i = 0; i < rec; i++)
        for(int j = 0; j < rec; j++) {

            for(int k = 0; k < rec; k++)
             result.data[i][j] += a.data[i][k] * b.data[k][j];

             result.data[i][j] %= p;

        }
           
    
    return result;
}

Matrix poww (Matrix x, int n, int p) {

    Matrix result;

    memset(&(result.data), 0,sizeof( result.data));

    for(int i = 0; i < R; i++) result.data[i][i] = 1;

    Matrix f = x;
    
    while(n) {
        if(n & 1) result = multi(result, f, R, p);

        f = multi(f, f, R, p);

        n >>= 1;
    }
    
    return result;
}

void MatrixPrint(Matrix target) {

    for(int i = 0; i < R; i++) {

        for(int j = 0; j < R; j++) cout << target.data[i][j]<< " ";

        cout << endl;
    
    }

}

int main() {

    Matrix result;

    memset(&(result.data), 0,sizeof( result.data));

    result.data[0][0] = 1;

    result.data[0][1] = 2;

    result.data[1][0] = 3;

    result.data[1][1] = 4;

    MatrixPrint(result);
    
    cout << endl;

    result =  poww(result, 3, 3);

    MatrixPrint(result);

    return 0;
}

輸出結果：

1 2
3 4

1 0
0 1

練習（HDU1005）：

Number Sequence

Problem Description

A number sequence is defined as follows:

\(f(1) = 1, f(2) = 1, f(n) = (A * f(n - 1) + B * f(n - 2)) \mod 7.\)

Given A, B, and n, you are to calculate the value of f(n).

Input

The input consists of multiple test cases. Each test case contains 3 integers A, B and n on a single line (1 <= A, B <= 1000, 1 <= n <= 100,000,000). Three zeros signal the end of input and this test case is not to be processed.

Output

For each test case, print the value of f(n) on a single line.

Sample Input

1 1 3
1 2 10
0 0 0

Sample Output

2
5

Author

CHEN, Shunbao

Source

ZJCPC2004

​ 由題意不難發現，題目輸入三個數A，B，n，要求咱們求出廣義的斐波那契數列（generalized Fibonacci sequence）第n項的取模

​ 即\(f(n) = (A * f(n - 1) + B * f(n - 2)) \mod 7 = (Aa_{n-1}+Ba_{n-2})= a_n \mod 7\)

​ 所以經過先後項關係，能夠構建矩陣\(\left| \begin{matrix} a_{n} & 0\\ a_{n-1} & 0\\ \end{matrix} \right|\mod 7=\left( \left| \begin{matrix} A & B\\ 1 & 0\\ \end{matrix} \right|\left| \begin{matrix} a_{n-1} & 0\\ a_{n-2} & 0\\ \end{matrix} \right| \right)\mod 7\)

\(=\left( \left| \begin{matrix} A & B\\ 1 & 0\\ \end{matrix} \right|^{n-2}\left| \begin{matrix} a_2 & 0\\ a_1 & 0\\ \end{matrix} \right| \right)\mod 7\)

\(=\left(\left( \left| \begin{matrix} A & B\\ 1 & 0\\ \end{matrix} \right|^{n-2}\mod 7\right)\left(\left| \begin{matrix} a_2 & 0\\ a_1 & 0\\ \end{matrix} \right|\mod7\right) \right)\mod 7\)

如此，就轉化爲一個矩陣快速冪的問題

具體實現（謝謝提醒！原來個人solution忽略了\(n = 1\)和\(n = 2\)的狀況，致使無限循環）：

#include <iostream>
#include <cstring>

using namespace std;

#define R 2

struct M{
    int data[R][R];
};

M multi(M a,M b,int rec,int p){

    M result;

    memset(&(result.data), 0,sizeof( result.data));

    for(int i = 0; i < rec; i++)
        for(int j = 0; j < rec; j++) {

            for(int k = 0; k < rec; k++)
             result.data[i][j] += a.data[i][k] * b.data[k][j];

             result.data[i][j] %= p;

        }
           
    
    return result;
}

M poww (M x, int n, int p) {

    M result;

    M f = x;

    memset(&(result.data), 0,sizeof( result.data));

    for(int i = 0; i < R; i++) result.data[i][i] = 1;
    
    while(n) {

        if(n & 1) result = multi(result, f, R, p);

        f = multi(f, f, R, p);

        n >>= 1;
        
    }
    
    return result;
}

M solve(int a, int b, int n, int p) {

    M result;

    M trans;

    memset(&(result.data), 0,sizeof( result.data));

    memset(&(trans.data), 0,sizeof( trans.data));

    result.data[0][0] = 1;

    result.data[1][0] = 1;

    trans.data[0][0] = a;

    trans.data[0][1] = b;

    trans.data[1][0] = 1;

    trans =  poww(trans, n, p);
    
    result = multi(trans, result, R, p);

    return result;

}

int main() {

    int a, b, n;


    while(true) {

        cin >> a >> b >> n;

        if(!a && !b && !n) break;

        if(n == 2 || n == 1 ){

            cout << 1 << endl;

            continue;
        }

        cout << solve(a, b, n - 2, 7).data[0][0] << endl;

    }


    return 0;
}
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程