1. 程式人生 > >矩陣快速冪專題(一)

矩陣快速冪專題(一)

最近閒來無事,準備集中精力刷一波數論與圖論。矩陣快速冪是數論裡面的重要組成部分,值得我好好學習一下。因為題目比較多,分析也比較多,所以將此專題分成幾個部分。做完這一專題,可能會暫時轉向圖論部分,然後等我組合數學學得差不多了,再回過頭來繼續做數論題。

矩陣快速冪演算法的核心思想是將問題建模轉化為數學模型(有一些簡單題目是裸的矩陣模型,但是大部分難題就是難在要構造矩陣,用矩陣方法解決問題),推倒遞推式,構造計算矩陣,用快速冪的思想求解矩陣A的n次方取mod,從而得到矩陣裡面你需要的資料。

矩陣快速冪問題有某些特徵,是類似快速冪的形式。首先,從題目中是一定可以得到一個遞推式的,而且這個遞推式不是簡單的An與An-1的關係,而是關係式中帶有An-2或之後的項,還有的情況甚至帶有常數項和多個數列的組合(加減乘數);其次,資料量比較大,不可能用打表方法做出來;最後,有些題目可能與矩陣一點沒有關係,但是通過規律分析,能夠轉化為已知遞推式求解某一項(取模)的。

這裡給出模版:  我覺得我的矩陣模版不是特別好,但是也不是特別麻煩,所以請讀者自行選擇是否採用

#define maxn 10+5;
typedef long long ll;
const ll mod = 1e9+7;
#define clr(x,y) memset(x,y,sizeof(x))
struct matrix 
{
	int n;
	ll maze[maxn][maxn];
	void init(int n)
	{
		this->n=n;
		clr(maze,0);
	}
	matrix operator *(const matrix& rhs)
	{
		matrix ans;
		ans.init(n);
		for(int i=0;i<n;i++)
			for(int j=0;j<n;j++)
				for(int k=0;k<n;k++)
					ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
		return ans;
	}
};
matrix qlow(matrix a,int n)
{
	matrix ans;
	ans.init(a.n);
	for(int i=0;i<a.n;i++)ans.maze[i][i]=1;
	while(n)
	{
		if(n&1)ans=ans*a;
		a=a*a;
		n>>=1;
	}
	return ans;
}


第一題 hdu-1757 

分析:題目直接給出了遞推式(帶有f(x-1)到f(x-10)),而且資料量非常大需要取模(這裡並不要管a0~a9的值,若不為01那麼只需多考慮一步取模),這是簡單的矩陣快速冪裸題。

#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 6+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
//const ll mod = 1e9+7;

int a[10];
ll k,m;
ll qmul(ll a,ll b)
{
     ll ans=0;
    while(b)
    {
        if(b&1)ans=(ans+a)%m;
        a=(a+a)%m;
        b>>=1;
    }
    return ans;
}
struct matrix
{
    int n;
    ll maze[maxn][maxn];
    void init(int n)
    {
        this->n=n;
        clr(maze,0);
    }
    matrix operator * (matrix & rhs)
    {
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                for(int k=0;k<n;k++)
                    ans.maze[i][j]=(ans.maze[i][j]+qmul(maze[i][k],rhs.maze[k][j]))%m;
        return ans;
    }
};
matrix qlow(matrix a,ll n)
{
    matrix ans;
    ans.init(a.n);
    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
    while(n)
    {
        if(n&1)ans=ans*a;
        a=a*a;
        n>>=1;
    }
    return ans;
}
int main()
{
    //freopen("d:\\acm\\in.in","r",stdin);
    while(~scanf("%lld %lld",&k,&m))
    {
        for(int i=0;i<10;i++)
            scanf("%d",&a[i]);
        if(k<10)
        {
            printf("%lld\n",k%m);
            continue;
        }
        matrix ans;
        ans.init(10);
        for(int i=0;i<10;i++)ans.maze[0][i]=i;
        matrix ant;
        ant.init(10);
        for(int i=0;i<9;i++)ant.maze[i+1][i]=1;
        for(int i=0;i<10;i++)ant.maze[i][9]=a[9-i];
        matrix tmp;
        tmp=qlow(ant,k-9);
        ans=ans*tmp;
        printf("%lld\n",ans.maze[0][9]);
    }
    return 0;
}
第二題 hdu-1575

分析:這道題就是非常直白的矩陣快速冪,將主對角線上的數加起來取模-。-

#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 6+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 9973;//1e9+7;

ll qmul(ll a,ll b)
{
    ll ans=0;
    while(b)
    {
        if(b&1)ans=(ans+a)%mod;
        a=(a+a)%mod;
        b>>=1;
    }
    return ans;
}
struct matrix
{
    int n;
    ll maze[maxn][maxn];
    void init(int n)
    {
        this->n=n;
        clr(maze,0);
    }
    matrix operator * (matrix& rhs)
    {
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                for(int k=0;k<n;k++)
                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
        return ans;
    }
};
matrix qlow(matrix a,int n)
{
    matrix ans;
    ans.init(a.n);
    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
    while(n)
    {
        if(n&1)ans=(ans*a);
        a=(a*a);
        n>>=1;
    }
    return ans;
}
int main()
{
    //freopen("d:\\acm\\in.in","r",stdin);
    int t;
    scanf("%d",&t);
    while(t--)
    {
        int n,k;
        scanf("%d %d",&n,&k);
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                scanf("%d",&ans.maze[i][j]);
        ans=qlow(ans,k);
        ll anw=0;
        for(int i=0;i<n;i++)
            anw=(anw+ans.maze[i][i])%mod;
        printf("%lld\n",anw);
    }
    return 0;
}
第三題 hdu-2604

分析:我這裡只提出我的方法,其他網上大神的方法請到別的部落格上查詢。

對於L較大的情況(L大於2,至於為什麼是2,你沒有看見fmf和fff都是三個嗎?),我考慮對於在E-queue裡面的所有合法佇列(-。-,其實和佇列沒有任何關係)。他們都是以mm,fm,mf,ff之一結尾的,我就分別將以他們結尾的並且總長度為n的個數記為mm[n],fm[n],mf[n],ff[n]。

考慮mm[n]後面加m是mm[n+1]的一部分,加f是mf[n+1]的一部分;

fm[n]加f屬於O-queue,捨去,加m是mm[n+1]的一部分;

依次類推。。。。。。

得到遞推式mm[n+1]=mm[n]+fm[n];    mf[n+1]=mm[n+1];   fm[n+1]=ff[n]+mf[n];    ff[n+1]=mf[n+1];

構造矩陣    


#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 0+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1e9+7;
int m;
struct matrix
{
    int n;
    ll maze[maxn][maxn];
    void init(int n)
    {
        this->n=n;
        clr(maze,0);
    }
    matrix operator * (const matrix& rhs)
    {
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                for(int k=0;k<n;k++)
                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%m;
        return ans;
    }
};
matrix qlow(matrix a,int n)
{
    matrix ans;
    ans.init(a.n);
    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
    while(n)
    {
        if(n&1)ans=ans*a;
        a=a*a;
        n>>=1;
    }
    return ans;
}
int main()
{
    //freopen("d:\\acm\\in.in","r",stdin);
    int l;
    while(~scanf("%d %d",&l,&m))
    {
        if(l<=2)
        {
            if(l==0)puts("0");
            else if(l==1)printf("%d\n",2%m);
            else printf("%d\n",4%m);
            continue;
        }
        matrix ans;
        ans.init(4);
        for(int i=0;i<4;i++)ans.maze[0][i]=1;
        matrix ant;
        ant.init(4);
        ant.maze[0][1]=ant.maze[1][3]=ant.maze[2][0]=ant.maze[2][1]=ant.maze[3][2]=ant.maze[3][3]=1;
        ant=qlow(ant,l-2);
        ans=ans*ant;
        int anw=0;
        for(int i=0;i<4;i++)anw=(anw+ans.maze[0][i])%m;
        printf("%d\n",anw);
    }
    return 0;
}


第四題  hdu-2256

分析:一開始會以為與矩陣無關,但是看見根號還要取模就會發現滿滿都是套路。-。-雖然根號2與根號3是兩個根號,但是還有個2次方,先將2次方乘進去,得到5+2sqrt(6),發現根號減少到一個,而且對於5+2sqrt(6)的n次方,一定是一個常數項和一個帶根號6的項,而且互相之間有遞推關係,這就把問題引向矩陣上了。


但是這裡不能直接取模,因為Bn還有一部分是sqrt(6),直接取模顯然是錯的


這樣構造的矩陣和最後要求的資料都有了 -。-

#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 0+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1024;//1e9+7;

struct matrix 
{
    int n;
    ll maze[2][2];
    void init(int n)
    {
        this->n=n;
        clr(maze,0);
    }
    matrix operator * (const matrix& rhs)
    {
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                for(int k=0;k<n;k++)
                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod; 
        return ans;
    }
};
matrix qlow(matrix a,int n)
{
    matrix ans;
    ans.init(a.n);
    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
    while(n)
    {
        if(n&1)ans=ans*a;
        a=a*a;
        n>>=1;
    }
    return ans;
}
int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        int n;
        scanf("%d",&n);
        matrix ans;
        ans.init(2);
        ans.maze[0][0]=1;
        matrix ant;
        ant.init(2);
        ant.maze[0][0]=ant.maze[1][1]=5;
        ant.maze[0][1]=2;
        ant.maze[1][0]=12;
        ant=qlow(ant,n);
        ans=ans*ant;
        printf("%d\n",(2*ans.maze[0][0]-1)%mod);
    }
    return 0;
}


第五題 codeforces-185A

分析:做完前面的題,這道題反而簡單了。不妨考慮當n年後(操作n次),正三角的個數是An,倒三角的個數是Bn。那麼繼續再操作一次,An個正三角變成了3An個正三角和An個倒三角,Bn個倒三角變成Bn個正三角和3Bn個倒三角。所以有遞推關係  An+1=3An+Bn    Bn+1=An+3Bn

構造矩陣   


#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 0+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1e9+7;

ll mul(ll a,ll b)
{
	ll ans=0;
	while(b)
	{
		if(b&1)ans=(ans+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return ans;
}
struct matrix 
{
	int n;
	ll maze[maxn][maxn];
	void init(int n)
	{
		this->n=n;
		clr(maze,0);
	}
	matrix operator * (const matrix& rhs)
	{
		matrix ans;
		ans.init(n);
		for(int i=0;i<n;i++)
			for(int j=0;j<n;j++)
				for(int k=0;k<n;k++)
					ans.maze[i][j]=(ans.maze[i][j]+mul(maze[i][k],rhs.maze[k][j]))%mod;
		return ans;
	}
};
matrix qlow(matrix a,ll n)
{
	matrix ans;
	ans.init(a.n);
	for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
	while(n)
	{
		if(n&1)ans=ans*a;
		a=a*a;
		n>>=1;
	}
	return ans;
}
int main()
{
    //freopen("d:\\acm\\in.in","r",stdin);
	ll n;
	scanf("%I64d",&n);
	matrix ans;
	ans.init(2);
	ans.maze[0][0]=1;
	matrix ant;
	ant.init(2);
	ant.maze[0][0]=ant.maze[1][1]=3;
	ant.maze[1][0]=ant.maze[0][1]=1;
	ant=qlow(ant,n);
	ans=ans*ant;
	printf("%I64d\n",ans.maze[0][0]);
    return 0;
}
第六題 hdu-2842
分析:題目講的好像是九連環(好像又不是,英語不好),-。-    可憐我並沒有玩過

這道題非常像漢諾塔,但是題目只是求最小步數,而不需要求出具體的步驟。我假設對於長度為i的中國環,我最少需要f(i)步(1<=i<=n)。那麼對於長度為n+1的中國環,我們先考慮將最後一個拿下來(沒辦法,奇葩的規則,只有拿下前n-1個且保留第n個,才能拿下第n+1個)。那麼先將前n-1個取下來,需要f(n-1)步;再將第n+1個取下來,需要1步;再把前n-1個放回去,需要f(n-1)步(這裡需要解釋一下,前n-1個必須放回去,不然第n個無法取下,而且取下需要f(n-1)步,那麼放回去也需要f(n-1),直接相反的步驟再來一遍,不明白的多讀幾遍題目);最後,將前n個都取下,需要f(n)步。所以得到遞推式 f(n+1)=f(n)+2f(n-1)+1。對於遞推式中有常數項,只需要給出一行或列留給常數項。

構造矩陣


#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 0+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 200907;//1e9+7;

struct matrix 
{
    int n;
    ll maze[maxn][maxn];
    void init(int n)
    {
        this->n=n;
        clr(maze,0);
    }
    matrix operator * (const matrix& rhs)
    {
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                for(int k=0;k<n;k++)
                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
        return ans;
    }
};
matrix qlow(matrix a,int n)
{
    matrix ans;
    ans.init(a.n);
    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
    while(n)
    {
        if(n&1)ans=ans*a;
        a=a*a;
        n>>=1;
    }
    return ans;
}
int main()
{
    //freopen("d:\\acm\\in.in","r",stdin);
    int n;
    matrix ant;
    ant.init(3);
    ant.maze[0][0]=ant.maze[0][1]=ant.maze[2][0]=ant.maze[2][2]=1;
    ant.maze[1][0]=2;
    while(scanf("%d",&n),n)
    {
        if(n<=2)
        {
            if(n==1)printf("%d\n",1);
            else printf("%d\n",2);
            continue;
        }
        matrix ans;
        ans.init(3);
        ans.maze[0][0]=2;
        ans.maze[0][1]=ans.maze[0][2]=1;
        ans=ans*qlow(ant,n-2);
        printf("%I64d\n",ans.maze[0][0]);
    }
    return 0;
}


第七題 hdu-2276

分析:對於01串,如果前一位的數位為1,那麼這一位就翻轉(0變1,1變0);前一位為0,那麼這一位就不變。這就像前一位與這一位相加,並且mod2,如果前一位是0,那麼這一位加上前一位再mod2數值不會改變;如果前一位是1,那麼這一位就相當於+1mod2,正好和翻轉的定義相同。這樣我們就可以構造出矩陣


#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 100+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1e9+7;

char str[maxn];
struct matrix 
{
    int n;
    ll maze[maxn][maxn];
    void init(int n)
    {
        this->n=n;
        clr(maze,0);
    }
    matrix operator * (const matrix & rhs)
    {
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                for(int k=0;k<n;k++)
                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%2;
        return ans;
    }
};
matrix qlow(matrix a,int n)
{
    matrix ans;
    ans.init(a.n);
    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;
    while(n)
    {
        if(n&1)ans=ans*a;
        a=a*a;
        n>>=1;
    }
    return ans;
}
int main()
{
    //freopen("d:\\acm\\in.in","r",stdin);
    int m;
    while(~scanf("%d",&m))
    {
        scanf("%s",str);
        int n=strlen(str);
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            if(str[i]=='1')ans.maze[0][i]=1;
        matrix ant;
        ant.init(n);
        for(int i=0;i<n;i++)
        {
            ant.maze[i][i]=1;
            if(i==0)ant.maze[n-1][0]=1;
            else ant.maze[i-1][i]=1;
        }
        ant=qlow(ant,m);
        ans=ans*ant;
        for(int i=0;i<n;i++)
            printf("%d",ans.maze[0][i]);
        puts("");
    }
    return 0;
}


第八題 hdu-2254

分析:學過離散數學的都知道,如果用鄰接矩陣來存圖的話,矩陣A的n次方里面的(A^n)[i][j]表示從i到j且路徑長度為n的不同路徑條數。這一定理正好適用這一題目。我用鄰接矩陣來建圖(單向邊,一條邊對應一個++操作)得到矩陣A,題目要求的就是{A^(t1)}[v1][v2]+{A^(t1+1)}[v1][v2]+{A^(t1+2)}[v1][v2]+......+{A^(t2)}[v1][v2]的和取模。但是如果一個一個算的話特別麻煩,還有超時的可能,所以我們進一步改進思路。

對於,我們不妨設Z(n)=A^1+A^2+...+A^n

那麼可以構造遞推關係 A^(n+1)=A^n * A      Z(n+1)=Z(n)+ A^(n+1) =Z(n)+A^n * A 

構造矩陣


但是這道題目略坑,有幾個注意點:

1、城市的編號不是從0到n-1,而是隨便的一個數字,需要離散化否則不能存相關資訊

2、城市數不超過30,也就是說我的方法開矩陣不超過60,但是我殘念的一開始以為最多可能有20000個不同城市    血崩!

3、圖中可能有重邊,所以別用=1,要用++操作

4、詢問中v1,v2可能在前面的城市編號集中沒有出現,那麼此時答案為0

5、t1可能比t2大,這種情況你就交換下t1,t2

這道題要做對也是真心不容易!

#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 60+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 2008;//1e9+7;

ll gra[10000+5][2];
ll num[20000+5];
int sum;
int bs(ll k)
{
    int m,l=0,r=sum-1;
    while(l<=r)
    {
        m=(l+r)>>1;
        if(k==num[m])return m;
        else if(k>num[m])l=m+1;
        else r=m-1;
    }
    return -1;
}
struct matrix 
{
    int n;
    ll maze[maxn][maxn];
    void init(int n)
    {
        this->n=n;
        clr(maze,0);
    }
    matrix operator * (const matrix & rhs)
    {
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                for(int k=0;k<n;k++)
                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
        return ans;
    }
};
matrix qlow(matrix a,int n)
{
    matrix ans;
    ans.init(a.n);
    for(int i=0;i<a.n;i++)ans.maze[i][i]=1;
    while(n)
    {
        if(n&1)ans=ans*a;
        a=a*a;
        n>>=1;
    }
    return ans;
}
int main()
{
    //freopen("d:\\acm\\in.in","r",stdin);
    int n;
    while(~scanf("%d",&n))
    {
        sum=n*2;
        for(int i=0;i<n;i++)
        {
            scanf("%lld %lld",&gra[i][0],&gra[i][1]);
            num[i<<1]=gra[i][0];
            num[i<<1|1]=gra[i][1];
        }
        sort(num,num+sum);
        sum=unique(num,num+sum)-num;
        matrix ans;
        ans.init(sum<<1);
        matrix ant;
        ant.init(sum<<1);
        int a,b;
        for(int i=0;i<n;i++)
        {
            a=bs(gra[i][0]);
            b=bs(gra[i][1]);
            ans.maze[a][b]++;
            ans.maze[a][b+sum]++;
        }
        for(int i=0;i<sum;i++)
            ans.maze[i+sum][i+sum]=1;
        int q;
        scanf("%d",&q);
        while(q--)
        {
            int t1,t2;
            ll v1,v2;
            scanf("%lld %lld %d %d",&v1,&v2,&t1,&t2);
            v1=bs(v1);
            v2=bs(v2);
            if(v1==-1||v2==-1)
            {
                puts("0");
                continue;
            }
            if(t1>t2)swap(t1,t2);
            t1--;
            int d,d1,d2;
            if(t1<=0)d1=0;
            else 
            {
                ant=qlow(ans,t1);
                d1=ant.maze[v1][v2+sum];
            }
            if(t2<=0)d2=0;
            else 
            {
                ant=qlow(ans,t2);
                d2=ant.maze[v1][v2+sum];
            }
            d=d2-d1;
            while(d<0)d+=mod;
            printf("%d\n",d);
        }
    }
    return 0;
}


第九題 hdu-3117

分析:其實這道題不應該放在矩陣快速冪專題的,因為這道題關於矩陣的部分很簡單-。-,難的是題目的數學部分-。-

求解斐波那契數列的值,對於數值位不超過8位的直接輸出數值;但是對於大於8位的,輸出值的前四位和後四位(格式自己看-。-)。後四位當然非常簡單了,f(n)=f(n-1)+f(n-2) ,mod 10000,我就不再贅述了

但是對於前四位還是挺麻煩的,這裡要用到對數的優良性質。logA(b^c)=c*logA(c)   logA(b*c)=logA(b)+logA(c)

對於一個大數(不如就選10234432),怎麼知道它的前四位呢

log10(10234432)=log10(1.0234432*10^7)=log10(1.0234432)+7

7是log10(10234432)的整數部分,log10(1.0234432)是log10(10234432)的小數部分

可以明顯看出pow(10,log10(1.0234432))=1.0234432   從而可以簡單地得到了前幾個位數

有人會問為什麼要這麼麻煩,不如以字串的形式讀入然後輸出前四位不就好了嗎 ?

對!沒錯,對於確定的大數,確實是可以的,但是對b^c來說卻不能使用,但是取對數的方法卻可以

對於這道題  我們考慮斐波那契數列的通項(有人問,比賽的時候你也知道通項?我的回答是可以,給我類似斐波那契的遞推式我就可以求出其通項,這個在數學競賽中有學到,請讀者自行百度)

F(n)=(1/sqrt(5))*{[(1+sqrt(5))/2]^n-[(1-sqrt(5))/2]^n}=(1/sqrt(5))*{[(1+sqrt(5))/2]^n}*{1-[(1-sqrt(5))/(1+sqrt(5))]^n}

設p=(1-sqrt(5))/2

s=log10(F(n))=-1/2*log10(5)+n*log10(p)+Δ(Δ非常小,可以忽略不計,Δ小於log10(1)

s去掉它的整數部分,然後(int) [ pow(10,s-(int)s)*1000 ] 就得到了結果

我只是簡單搬運而且=、=

#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 60+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 10000;//1e9+7;

int fa[60];
struct matrix 
{
    int n;
    ll maze[5][5];
    void init(int n)
    {
        this->n=n;
        clr(maze,0);
    }
    matrix operator * (const matrix & rhs)
    {
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                for(int k=0;k<n;k++)
                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
        return ans;
    }
};
matrix qlow(matrix a,int n)
{
    matrix ans;
    ans.init(a.n);
    for(int i=0;i<a.n;i++)ans.maze[i][i]=1;
    while(n)
    {
        if(n&1)ans=ans*a;
        a=a*a;
        n>>=1;
    }
    return ans;
}
int main()
{
    //freopen("d:\\acm\\in.in","r",stdin);
    fa[0]=0;
    fa[1]=1;
    for(int i=2;i<40;i++)
        fa[i]=fa[i-1]+fa[i-2];
    int n;
    while(~scanf("%d",&n))
    {
        if(n<40)
        {
            printf("%d\n",fa[n]);
            continue;
        }
        double p=(sqrt(5.0)+1.0)/2.0;
        double s=1.0*n*log10(p)-log10(5.0)/2;
        s=s-(int)s;
        int anw=(int)(pow(10.0,s)*1000);
        printf("%d...",anw);
        matrix ant;
        ant.init(2);
        ant.maze[0][0]=ant.maze[0][1]=ant.maze[1][0]=1;
        ant=qlow(ant,n-1);
        printf("%04d\n",ant.maze[0][0]);
    }
    return 0;
}

第十題 hdu-4686

分析:這是一道顯然的矩陣題(刷完前面那麼多道,如果還沒有一眼看出來的話,建議停下來重新理解演算法-。-)

感覺難度不是很大,但是這道題有多個遞推關係,而且還涉及到常數項和兩個不相關的遞推關係乘積,稍微有點麻煩。但是冷靜下來,對每一個遞推關係進行分析和對你所要求的項進行分解。

由AoD(n)=a(0)b(0)+a(1)b(1)+...+a(n-1)b(n-1)   可得    AoD(n+1)=AoD(n)+a(n)b(n)

那麼a(n)b(n)怎麼求呢   現在只知道  a(n)=a(n-1)*ax+ay    b(n)=b(n-1)*bx+by

那麼可不可以利用已經得到的a(n-1)和b(n-1)和a(n-1)b(n-1)   

可以得到 a(n)b(n)=ax*bx*a(n-1)*b(n-1)+ax*by*a(n-1)+ay*bx*b(n-1)+ay*by

那麼問題又來了  a(n)和 b(n)又怎麼單獨求呢     這個題目已經給出  a(n)=a(n-1)*ax+ay    b(n)=b(n-1)*bx+by

分析清楚了所有的遞推關係後,可以來構造矩陣了


注意點:a0,ax,ay,b0,bx,by都可能大於1e7+9,所以需要先取一次模

#include<set>
#include<map>
#include<cmath>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cctype>
#define maxn 1+5
#define clr(x,y) memset(x,y,sizeof(x))
using namespace std;
const int inf = 0x3f3f3f3f;
typedef long long ll;
const double pi = acos( -1 );
const ll mod = 1e9+7;

ll mul(ll a,ll b)
{
    ll ans=0;
    while(b)
    {
        if(b&1)ans=(ans+a)%mod;
        a=(a+a)%mod;
        b>>=1;
    }
    return ans;
}
struct matrix 
{
    int n;
    ll maze[maxn][maxn];
    void init(int n)
    {
        this->n=n;
        clr(maze,0);
    }
    matrix operator * (const matrix & rhs)
    {
        matrix ans;
        ans.init(n);
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                for(int k=0;k<n;k++)
                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;
        return ans;
    }
};
matrix qlow(matrix a,ll n)
{
    matrix ans;
    ans.init(a.n);
    for(int i=0;i<a.n;i++)ans.maze[i][i]=1;
    while(n)
    {
        if(n&1)ans=ans*a;
        a=a*a;
        n>>=1;
    }
    return ans;
}
int main()
{
    //freopen("d:\\acm\\in.in","r",stdin);
    ll n;
    while(~scanf("%lld",&n))
    {
        ll a0,b0,ax,bx,ay,by;
        scanf("%lld %lld %lld %lld %lld %lld",&a0,&ax,&ay,&b0,&bx,&by);
        a0%=mod;
        b0%=mod;
        ax%=mod;
        bx%=mod;
        ay%=mod;
        by%=mod;
        if(n==0)
        {
            puts("0");
            continue;
        }
        matrix ans;
        ans.init(5);
        ans.maze[0][0]=1;
        ans.maze[0][1]=a0;
        ans.maze[0][2]=b0;
        ans.maze[0][3]=a0*b0%mod;
        matrix ant;
        ant.init(5);
        ant.maze[0][0]=1;
        ant.maze[0][1]=ay;
        ant.maze[1][1]=ax;
        ant.maze[0][2]=by;
        ant.maze[2][2]=bx;
        ant.maze[0][3]=ay*by%mod;
        ant.maze[1][3]=ax*by%mod;
        ant.maze[2][3]=ay*bx%mod;
        ant.maze[3][3]=ax*bx%mod;
        ant.maze[3][4]=ant.maze[4][4]=1;
        ant=qlow(ant,n);
        ans=ans*ant;
        printf("%lld\n",ans.maze[0][4]);
    }
    return 0;
}


這裡先小結一下,矩陣快速冪型別的題目非常像科學歸納法(第一科學歸納法、第二科學歸納法還有其他型別的),主要是求得當前項與前幾項或者常數項的遞推關係,從而構造矩陣求解。