矩阵快速幂4

目录

Power of Matrix

Yet Another Number Sequence

---

问题十一：

Power of Matrix

UVA - 11149

Consider an n-by-n matrix A. We define Ak = A∗A∗...∗A (k times). Here, ∗ denotes the usual matrix multiplication. You are to write a program that computes the matrix A + A2 + A3 + ... + Ak.
Example Suppose A = 
0 2 0 0 0 2 0 0 0  . Then A2 = 
0 2 0 0 0 2 0 0 0   
0 2 0 0 0 2 0 0 0  = 
0 0 4 0 0 0 0 0 0  , thus:
A + A2 = 
0 2 0 0 0 2 0 0 0  + 
0 0 4 0 0 2 0 0 0  = 
0 2 4 0 0 2 0 0 0   Such computation has various applications. For instance, the above example actually counts all the paths in the following graph:
Input
Input consists of no more than 20 test cases. The first line for each case contains two positive integers n (≤ 40) and k (≤ 1000000). This is followed by n lines, each containing n non-negative integers, giving the matrix A. Input is terminated by a case where n = 0. This case need NOT be processed.
Output
For each case, your program should compute the matrix A+A2 +A3 +...+Ak. Since the values may be very large, you only need to print their last digit. Print a blank line after each case.
Sample Input
3 2 0 2 0 0 0 2 0 0 0 0 0
Sample Output
0 2 4 0 0 2 0 0 0

---

题意：给你n，k。有n组数据计算的值

分析：这道题两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：
　　A + A^2 + A^3 + A^4 + A^5 + A^6 =(A + A^2 + A^3) + A^3*(A + A^2 + A^3)
应用这个式子后，规模k减小了一半。我们二分求出A^3后再递归地计算A + A^2 + A^3，即可得到原问题的答案。

奇数: F[n]=F[n-1]+A^n        

偶数: F[n]=F[n/2]+F[n/2]*An/2

代码：

#include<stdio.h>
#include<string.h>
#define N 31
struct Matrix
{
    int a[N][N];
}res,origin,tmp,A,B,C,ans;
int n,m;
Matrix mul(Matrix x,Matrix y)   //矩阵乘法
{
    int i,j,k;
    memset(tmp.a,0,sizeof(tmp.a));
    for(i=1;i<=n;i++)
        for(j=1;j<=n;j++)
            for(k=1;k<=n;k++)
            {
                tmp.a[i][j]+=(x.a[i][k]*y.a[k][j]);
                tmp.a[i][j]%=m;
            }
    return tmp;
}
Matrix quickpow(Matrix origin,int k)   //矩阵快速幂
{
    int i;
    memset(res.a,0,sizeof(res.a));
    for(i=1;i<=n;i++)
        res.a[i][i]=1;
    while(k)
    {
        if(k&1)
            res=mul(res,origin);
        origin=mul(origin,origin);
        k>>=1;
    }
    return res;
}
Matrix sum(Matrix x,Matrix y)   //矩阵求和
{
    int i,j;
    for(i=1;i<=n;i++)
        for(j=1;j<=n;j++)
            tmp.a[i][j]=(x.a[i][j]+y.a[i][j])%m;
    return tmp;
}
Matrix bin(int k)
{
    if(k<=1)
        return A;
    else if(k%2)   //奇数：F[n]=F[n-1]+A^n   
    {
        B=bin(k-1);
        C=quickpow(A,k);   
        return sum(B,C);
    }
    else     //偶数： F[n]=F[n/2]+F[n/2]*A^(n/2)
    {
        B=bin(k/2);
        C=quickpow(A,k/2);
        C=mul(C,B);
        return sum(B,C);
    }
}
int main()
{
    int i,j,k;
    scanf("%d%d%d",&n,&k,&m);
    for(i=1;i<=n;i++)
        for(j=1;j<=n;j++)
            scanf("%d",&A.a[i][j]);
    ans=bin(k);   //二分
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++)
            if(j<n)
                printf("%d ",ans.a[i][j]);
            else
                printf("%d\n",ans.a[i][j]);
    }
    return 0;
}

问题十二：

Yet Another Number Sequence

CodeForces - 392C

Everyone knows what the Fibonacci sequence is. This sequence can be defined by the recurrence relation:

F1 = 1, F2 = 2, Fi = Fi - 1 + Fi - 2 (i > 2).

We'll define a new number sequence Ai(k) by the formula:

Ai(k) = Fi × ik (i ≥ 1).

In this problem, your task is to calculate the following sum: A1(k) + A2(k) + ... + An(k). The answer can be very large, so print it modulo 1000000007 (109 + 7).

Input

The first line contains two space-separated integers n, k (1 ≤ n ≤ 1017; 1 ≤ k ≤ 40).

Output

Print a single integer — the sum of the first n elements of the sequence Ai(k) modulo 1000000007 (109 + 7).

Examples

Input

1 1

Output

1

Input

4 1

Output

34

Input

5 2

Output

316

Input

7 4

Output

73825

---

题意：给出n，k。求的前n项和。

分析链接：优美的Fibonacci数列与矩阵

如何构造矩阵链接：根据递推公式构造系数矩阵用于快速幂

代码：

#include<iostream>
#include<cstdio>
#include<cstring>

using namespace std;
typedef long long LL;
const int mod = 1e9 + 7;
int len;
LL n, k, c[45][45];
struct matrix {
	LL m[100][100];
}mc,I;
matrix multi(matrix a, matrix b) {
	matrix ans;
	for (int i = 0; i < len; i++) {
		for (int j = 0; j < len; j++) {
			ans.m[i][j] = 0;
			for (int k = 0; k < len; k++) {
				ans.m[i][j] = (ans.m[i][j] + (a.m[i][k] * b.m[k][j]) % mod) % mod;
			}
		}
	}
	return ans;
}
matrix quick(matrix a, LL n) {
	matrix ans = I, t = a;
	while (n>0) {
		if (n & 1)
			ans = multi(ans, t);
		t = multi(t, t);
		n >>= 1;
	}
	return ans;
}
int main()
{
	c[0][0] = 1; c[1][0] = 1; c[1][1] = 1;
	for (int i = 2; i <= 40; i++) {
		for (int j = 0; j <= i; j++) {
			if (j == 0 || j == i)
				c[i][j] = 1;
			else
				c[i][j] = c[i - 1][j - 1] + c[i - 1][j];
		}
	}
	for (int i = 0; i <= 40; i++) {
		I.m[i][i] = 1;
	}
	while (cin >> n >> k) {
		memset(mc.m, 0, sizeof(mc.m));
		len = 2 * k + 3;
		mc.m[0][0] = 1;
		for (int j = 0; j <= k; j++) {
			mc.m[0][j + 1] = c[k][j];
			mc.m[0][j + k + 2] = c[k][j];
		}
		mc.m[1][1] = 1; mc.m[1][k + 2] = 1;
		for (int i = 1; i <= k; i++)
			for (int j = 0; j <= i; j++) {
				mc.m[i + 1][j + 1] = c[i][j];
				mc.m[i + 1][j + k + 2] = c[i][j];
			}
		mc.m[k + 2][1] = 1;
		for (int i = 1; i <= k; i++) {
			for (int j = 0; j <= i; j++) {
				mc.m[i + k + 2][j + 1] = c[i][j];
			}
		}
		mc = quick(mc, n - 1);
		LL ans = 0;
		for (int i = 0; i < len; i++) {
			ans = (ans + mc.m[0][i] % mod) % mod;
		}
		cout << ans << endl;
	}
	return 0;
}

大佬博客矩阵十题：https://www.cnblogs.com/frog112111/archive/2013/05/16/3082416.html