矩陣快速冪在ACM中的應用

矩陣快速冪在ACM中的應用

16計算機2黃睿博

首發於我的博客http://www.cnblogs.com/BobHuang/

做爲一個acmer，矩陣在這個算法競賽中仍是蠻多的，一個優秀的算法能夠影響到一個程序的運行速度的快慢，在算法競賽中經常採用快速冪算法，由於有些遞推式及有些問題均可以化爲矩陣，因此矩陣快速冪也就有了意義。

首先引入快速冪：

咱們要求a^b，那麼其實b是能夠拆成二進制的，該二進制數第i位的權爲2^(i-1)，

例如當b=11時  11的二進制是1011

11 = 2³×1 + 2²×0 + 2¹×1 + 2º×1所以，咱們將a¹¹轉化爲算 由於存在這一相等關係，因此將O(n)複雜度的算法降到了O(logn)，代碼以下

int poww(int x, int n)
{
    int ans=1;
    while(n)
    {
        if(n&1)ans=ans*x;
        x=x*x;
        n>>=1;
    }
    return ans;
}

這裏面用到了一些奇怪的運算，位運算，關於這個詳見個人另外一篇博客（>>至關於/2，&1至關於判斷末位是否爲1，比較接近計算機底層，取模還有除法運算常數比較大）

個人矩陣快速冪模板

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=105;
int G;
struct MX
{
    ll v[N][N];
    void O()
    {
        memset(v,0,sizeof v);
    }
    void E()
    {
        memset(v,0,sizeof v);
        for(int i=0; i<G; i++)v[i][i]=1;
    }
    void P()
    {
        for(int i=0; i<G; i++)
            for(int j=0; j<G; j++)printf(j==G-1?"%d\n":"%d ",v[i][j]);
    }
    MX operator+(const MX &b) const
    {
        MX c;
        c.O();
        for(int i=0; i<G; i++)
            for(int j=0; j<G; j++)c.v[i][j]=v[i][j]+b.v[i][j];
        return c;
    }
    MX operator*(const MX &b)const
    {
        MX c;
        c.O();
        for(int k=0; k<G; k++)
            for(int i=0; i<G; i++)
                if(v[i][k])for(int j=0; j<G; j++)c.v[i][j]+=v[i][k]*b.v[k][j];
        return c;
    }
    MX operator^(int p)const
    {
        MX y,x;
        y.E(),memcpy(x.v,v,sizeof(v));
        for(; p; x=x*x,p>>=1)if(p&1)y=y*x;
        return y;
    }
} a,ans;

 

目前c/c++最大支持__int128，也就是最大的有符號整數爲2^127-1，數量級大概是1e38，不過經常使用的仍是long long，是長整形，2^63-1。因此經常使用一個MOD質數的方法，這樣獲得的答案更加豐富。

引入完畢，接下來進入正題，無非是將以上的乘法轉換爲矩陣乘法。

咱們從最簡單的斐波那契數列入手，求菲波那切數列的第n項，n很大，因此MOD1e9+7

斐波的爲前兩項之和，即發f[0]=1,f[1]=1,f[i] = f[i-1]+f[i-2](i>1)

比較簡單的能夠直接獲得遞推式

 

，進一步遞推

因爲矩陣乘法知足結合律，所它也能夠用快速冪算法求解。

代碼以下

struct Matrix
{
    long long a[2][2];
    Matrix()
    {
        memset(a,0,sizeof(a));
    }
    Matrix operator*(const Matrix  y)
    {
        Matrix  ans;
        for(int i=0; i<=1; i++)
            for(int j=0; j<=1; j++)
                for(int k=0; k<=1; k++)
                    ans.a[i][j]+=a[i][k]*y.a[k][j];
        for(int i=0; i<=1; i++)
            for(int j=0; j<=1; j++)
                ans.a[i][j]%=M;
        return ans;
    }
    void operator=(const Matrix b)
    {
        for(int i=0; i<=1; i++)
            for(int j=0; j<=1; j++)
                a[i][j]=b.a[i][j];
    }
};
long long solve(long long x)
{
    Matrix ans,t;
    ans.a[0][0]=ans.a[1][1]=1;
    t.a[0][0]=t.a[1][0]=t.a[0][1]=1;
    while(x)
    {
        if(x&1)ans=ans*t;
        t=t*t;
        x>>=1;
    }
    return ans.a[0][0];
}

 強行遞推式

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int mod = 1000000009;
unordered_map<ll,ll> M;
ll fic(ll x)
{
    if(M.count(x))return M[x];
    ll a=(x+1)/2,b=x/2+1;
    return M[x]=((fic(a)*fic(b))%mod+(fic(a-1)*fic(b-1)))%mod;
}
int main()
{
    ll n;
    M[0]=0,M[1]=M[2]=1;
    cin>>n;
    cout<<fic(n);
    return 0;
}

 

若是通過推理獲得一個遞推式呢

若是對矩陣乘法仍是不太懂的話，能夠先看下這個知乎回答

例子HDU2842

Dumbear likes to play the Chinese Rings (Baguenaudier). It’s a game played with nine rings on a bar. The rules of this game are very simple: At first, the nine rings are all on the bar.

The first ring can be taken off or taken on with one step.

If the first k rings are all off and the (k + 1)th ring is on, then the (k + 2)th ring can be taken off or taken on with one step. (0 ≤ k ≤ 7)

 

Now consider a game with N (N ≤ 1,000,000,000) rings on a bar, Dumbear wants to make all the rings off the bar with least steps. But Dumbear is very dumb, so he wants you to help him.

九連環：若是要拆第n個環，那麼第n-1個環就必須在竿上，

前n-2個環都必須已經被拆下，題目是如今是一個n連環，求將n連環所有

拆掉須要的最少的步驟%200907.

若是前k個環被拆掉,第k+1個還被掛着,那麼第k+2個就能夠拿下或者裝上

f[n] = 2 * f[n - 2] + f[n - 1] + 1

 

 

就是f(n)=f[n-1]+2*f[n-2]+1,f[n-1]=1*f[n-1],1=1;

就是能夠從前兩個狀態推到當前狀態，這個關係矩陣的係數仍是比較好找的

下一個例子

n個六邊形排成一行，相鄰兩個六邊形共用一條邊，以下圖所示：

 

 

 

記這個圖形的生成樹個數爲t(n)（因爲每條邊都是不一樣的，不存在同構的問題）。那麼t(1)=6，t(2)=35……

給出n，求 mod 1000000007

第i個矩形右邊那條邊會不會去掉的方案數

dp[i][1]=dp[i-1][0]+dp[i-1][1],

dp[i][0]=4*dp[i-1][1]+5*dp[i-1][0]

i個六邊形總個數就是dp[i][0] + dp[i][1];

矩陣就是

 

分析下dp[i+1][2]，這一項已經作了前綴和，因此和以前的總數是不同的

dp[i+1][2]=dp[i+1][1]+dp[i+1][0]+dp[i][2]=dp[i][0]+dp[i][1]+4dp[i][0]+5dp[i-1][1]+dp[i][2]

=5dp[i][0]+6dp[i][1]+dp[i][2]

因此這個就很好解決了

假設a[i]爲當前項的前綴和，很容易發現

a[i] =6*a[i-1] - a[i-2] +5

這個用矩陣怎麼填呢

 

能夠這樣填，還有看到他們最後把a[0][0]-1的，可是相信懂了矩陣乘法的你懂這戲額都是在幹嗎

 

（由於初值的設置不一樣）

 最後一個例子

一個n位數，它的每位都是奇數，且數字1和3出現偶數次，這樣的n位數有多少個。好比說n=1，只有3個，它們分別是5，7和9。讓你求下知足這些條件的數的個數MOD9973，對於給定的n都會包含有四種狀態，1和3的個數都是奇數；1是奇數，3是偶數；1是偶數，3是奇數；1是偶數，3是偶數。

1是偶數，3是偶數是我想要的狀態

在相互轉換中實際上是能夠直接寫出係數的

 

a[0][3]是咱們要的狀態即爲所求

代碼以下

#include <stdio.h>
#include <string.h>
const int MD=9973;
struct matrix
{
    int mat[5][5];
};
matrix matmul(matrix a,matrix b,int n)
{
    int i,j,k;
    matrix c;
    memset(c.mat,0,sizeof(c.mat));
    for(i=0; i<n; i++)
    {
        for(j=0; j<n; j++)
        {
            for(k=0; k<n; k++)
            {
                c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%MD;
            }

        }
    }
    return c;
}
matrix matpow(matrix a,int k,int n)
{
    matrix b;
    int i;
    memset(b.mat,0,sizeof(b.mat));
    for(i=0; i<n; i++) b.mat[i][i]=1;
    while(k)
    {
        if(k&1) b=matmul(a,b,n);
        a=matmul(a,a,n);
        k>>=1;
    }
    return b;
}
int main()
{
    int k;
    while(~scanf("%d",&k))
    {
        matrix a,b,ans;
        memset(a.mat,0,sizeof(a.mat));
        memset(b.mat,0,sizeof(b.mat));
        a.mat[0][3]=1;
        for(int i=0; i<4; i++)
        {
            for(int j=0; j<4; j++)
            {
                if(i==j) b.mat[i][j]=3;
                else if(i+j==3) b.mat[i][j]=0;
                else b.mat[i][j]=1;
            }
        }
        ans=matmul(a,matpow(b,k,4),4);
        printf("%d\n",ans.mat[0][3]);
    }
    return 0;
}