/*
矩阵乘法模板题
给定矩阵A,请快速计算出A^n(n个A相乘)的结果,输出的每个数都mod p。
由于矩阵乘法具有结合律,因此A^4 = A * A * A * A = (A*A) * (A*A) = A^2 * A^2。我们可以得到这样的结论:当n为偶数时,A^n = A^(n/2) * A^(n/2);当n为奇数时,A^n = A^(n/2) * A^(n/2) * A (其中n/2取整)。这就告诉我们,计算A^n也可以使用二分快速求幂的方法。例如,为了算出A^25的值,我们只需要递归地计算出A^12、A^6、A^3的值即可。根据这里的一些结果,我们可以在计算过程中不断取模,避免高精度运算。
递归实现POW函数
Matrix POW( Matrix t,int k )
{
if( k == 1 )
return t;
Matrix t1 = POW( t, k/2 );
t1 = t1*t1;
if( k & 1 )
return t1 * t;
else
return t1;
}
递归的容易理解,但时间花费较多。
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
typedef struct Node
{
int m[11][11];
}Matrix;
Matrix init, unit; //初始化输入矩阵,单位矩阵如果用递归写Pow函数可以不用单位矩阵
int n, K;
void Init( ) //初始化
{
scanf( "%d%d", &n, &K );
for( int i=0; i<n; ++ i )
for( int j=0; j<n; ++ j )
{
scanf( "%d", &init.m[i][j] );
unit.m[i][j]=( i == j );
}
}
Matrix Mul( Matrix a, Matrix b ) //矩阵乘法
{
Matrix c;
for( int i=0; i<n; ++ i )
{
for( int j=0; j<n; ++ j )
{
c.m[i][j]=0; //特别注意
for( int k=0; k<n; ++ k )
{
c.m[i][j] += a.m[i][k]*b.m[k][j];
}
c.m[i][j] %= 9973;
}
}
return c;
}
Matrix Pow( Matrix a, Matrix b, int k )
{
while( k>1 )
{
if( k&1 ) // k为奇数时
{
k --;
b=Mul( a, b );
}
else // k为偶数
{
k >>= 1;
a=Mul( a, a );
}
}
a=Mul( a, b );
return a;
}
int main( )
{
int T;
scanf( "%d", &T );
while( T -- )
{
Matrix x;
Init( );
x=Pow( init, unit, K );
int sum=0, i=0;
n--;
while( n >= 0 )
{
sum += x.m[n][n];
sum%=9973;
n --;
}
printf( "%d\n", sum%9973 );
}
}
A Simple Math Problem
/*
If x < 10 f(x) = x.
If x >= 10 f(x) = a0 * f(x-1) + a1 * f(x-2) + a2 * f(x-3) + …… + a9 * f(x-10);
|f(10) | |a0 a1 a2 ...a8 a9| |f(9)|
| f(9) | | 1 0 0 ... 0 0| |f(8)|
| .....| = |.. ... ... ... ..| | .. |
| f(2) | | 0 0 0 ... 0 0| |f(1)|
| f(1) | | 0 0 0 ... 1 0| |f(0)|
另A举证为10*10的举证,如上图。
可以推出:
(f(n),f(n-1),...,f(n-9))^(-1) = A^(n-9)*(f(9),f(8),...,f(0))^(-1)
*/
#include<iostream>
using namespace std;
__int64 k,m;
struct mat
{
int mar[10][10];
}a,b,tmp;
mat matrixmul(mat a,mat b)
{
int i,j,k;
for(i = 0;i < 10;i++)
for(j = 0;j < 10;j++)
{
tmp.mar[i][j] = 0;
for(k = 0;k < 10;k++)
tmp.mar[i][j] += (a.mar[i][k] * b.mar[k][j]) % m;
tmp.mar[i][j] %= m;
}
return tmp;
}
void matrix_binary()
{
while(k)
{
if(k & 1)
b = matrixmul(b,a);
a = matrixmul(a,a);
k = k >> 1;
}
}
int main()
{
int i;
while (scanf("%I64d%I64d",&k,&m) != EOF)
{
memset(a.mar,0,sizeof(a.mar));
for(i = 1;i < 10;i++)
a.mar[i][i-1] = 1;
memset(b.mar,0,sizeof(b.mar));
for(i = 0;i < 10;i++)
b.mar[i][i] = 1;
for(i = 0;i < 10;i++)
scanf("%d",&a.mar[0][i]);
if(k < 10)
{
printf("%d\n", k % m);
continue;
}
k -= 9;
matrix_binary();
int res = 0;
for(i = 0;i < 10;i++)
res += (b.mar[0][i] * (9-i)) % m;
printf("%d\n",res%m);
}
return 0;
}
Pendant
/*
本题从表面上看是排列组合题,但要推出公式还是有相当难度。所以想到用DP来做。方程可通过枚举几种情况来推出。
以F[i][j]表示长度为i的pendant,用了j种珍珠,所构成的方案数,则F[i][j]=F[i-1][j]*j+F[i-1][j-1]*(k-j+1)。结果就是F[1][k]+…+F[n][k]。但注意到N的范围很大,申请那么大的数组会MLE。
注意到当前状态只与上一个状态有关,那么可以使用循环数组来,并累加上每个F[I][K]的方法来做。但这样的复杂度为O(NK),会TLE。
优化的方法是使用矩阵来做。将F[i-1]到F[i]的转移用矩阵来描述,相当于一个k*k的线性变换矩阵。因此F[i]=A*F[i-1],这里A是转移矩阵,即F[i]=Ai-1*F[1],所以F[1]+…+F[n]=A0*F[1]+…+An-1*F[1]=(E+A+A2+…+An-1)*F[1]。
F[i-1]这个矩阵是
F[i-1][0] F[i-1][1] F[i-1][2] ...... F[i-1][k]
0 0 0 …… 0
. . . . .
. . . . .
0 0 0 0 0
(这是个K+1阶的阶阵)
A矩阵是
0 k 0 0 0 0
0 1 k-1 0 0 0
0 0 2 0 0 0
. . . . . .
. . . . . .
0 0 0 0 k-1 1
0 0 0 0 0 k
这两个矩阵为什么是这样的,用矩阵的乘法来乘一下,将结果与DP方程对比下就知道了。
接下来要解决的问题就是F[1]要怎样填写。每次再另写一个程序来填F[1]?大可不必。因为F[1]=F[0]*A,那么只要填F[0]就行了。将F[0]变为E(单位阵)就可以了。
最后的问题就是求矩阵的平方和了。这里的方法很多,到网上一找也很容易找到,就不多说了。下面用的是二分求和,二分求幂。
*/
#include<iostream>
void MIN(int& a,int b) { if( a>b ) a=b; }
void MAX(int& a,int b) { if( a<b ) a=b; }
#define CLR(NAME,VALUE) memset(NAME,VALUE,sizeof(NAME))
using namespace std;
const int M=1234567891;
int n,k;
int a[31][31],b[31][31],c[31][31];
int d[31][31],z[31][31];
void Mul(int x[31][31],int y[31][31]) { //矩阵乘法,结果放在X中
int i,j,jj;
for(i=0;i<=k;++i) {
for(j=0;j<=k;++j) {
z[i][j]=0;
for(jj=0;jj<=k;++jj) {
z[i][j]=((__int64)z[i][j]+((__int64)x[i][jj]*y[jj][j])%M)%M;
}
}
}
for(i=0;i<=k;++i) {
for(j=0;j<=k;++j) {
x[i][j]=z[i][j];
}
}
}
void Cal(int p) { //递归二分算和式:A^1+A^2+...+A^N
if( p==0 ) return;
Cal(p/2);
int i,j;
for(i=0;i<=k;++i) {
for(j=0;j<=k;++j) {
d[i][j]=b[i][j];
}
}
for(i=0;i<=k;++i) { //D=B+E
d[i][i]=(d[i][i]+1)%M;
}
Mul(c,d); //C=C*D=C*B+C;
Mul(b,b); //B=B*B;
if( p&1 ) {
Mul(b,a); //B=B*A
for(i=0;i<=k;++i) { //C=C+B
for(j=0;j<=k;++j) {
c[i][j]=((__int64)c[i][j]+b[i][j])%M;
}
}
}
}
int main() {
int ca,i;
scanf("%d",&ca);
while( ca-- ) {
scanf("%d%d",&n,&k);
CLR(a,0);
CLR(b,0);
CLR(c,0);
//将B矩阵变为单位阵
for(i=0;i<=k;++i) b[i][i]=1;
//F[i][j]=F[i-1][j]*j+F[i-1][j-1]*(k-j+1)
for(i=1;i<=k;++i) { //转移矩阵A
a[i-1][i]=k-i+1;
a[i][i]=i;
}
Cal(n);
printf("%d\n",c[0][k]);
}
return 0;
}
奥运
#include <iostream>
#include <cstring>
#include <cstdio>
#define MAX 35
using namespace std;
struct matrix
{
int num[MAX][MAX];
matrix()
{
memset(num,0,sizeof(num));
}
};
int n;
long long map[MAX][MAX];
matrix A[10005];
long long city[35];
matrix operator*(matrix &a,matrix &b)
{
matrix t;
int i,j,k;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
t.num[i][j]=0;
for(k=0;k<n;k++)
t.num[i][j]+=(a.num[i][k]*b.num[k][j]);
t.num[i][j]%=2008;
}
return t;
}
void init()
{
int i,j;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
A[0].num[i][j]=map[i][j];
}
}
}
int main()
{
int nn;
int Find(int u);
while(scanf("%d",&nn)!=EOF)
{
int i;
memset(map,0,sizeof(map));
memset(city,0,sizeof(city));
n=0;
for(i=0;i<nn;i++)
{
int u,v;
scanf("%d%d",&u,&v);
int lo1=Find(u);
int lo2=Find(v);
map[lo1][lo2]++;
}
int k;
init();
for(i=1;i<10005;i++)
A[i]=A[i-1]*A[0];
scanf("%d",&k);
while(k--)
{
int v1,v2,t1,t2;
scanf("%d%d%d%d",&v1,&v2,&t1,&t2);
if(t1>t2)
{
int t=t1;
t1=t2;
t2=t;
}
int l=n;
int l1=Find(v1);
int l2=Find(v2);
if(l1>=l||l2>=l)
{
printf("0\n");
n=l;
continue;
}
int sum=0;
for(i=t1-1;i<t2;i++)
sum=(sum%2008+A[i].num[l1][l2]%2008)%2008;
printf("%d\n",sum);
}
}
return 0;
}
int Find(int u)
{
int i;
for(i=0;i<n;i++)
{
if(city[i]==u)
break;
}
if(i>=n)
{
n++;
city[i]=u;
return i;
}
else
{
return i;
}
}
Kiki & Little Kiki 2
/*
题目大意:题目大意给定一系列灯的初始状态,0代表暗,1代表亮,每一秒所有的灯都有可能发生状态切换,
切换规则:当前灯的左边第一个灯是亮的,则当前的灯切换状态,如果当前灯的左边第一盏灯是暗的,则当前灯的状态无需变化!
注意:最左边的参考左右边那栈灯。
题目分析;
首先有两种情况:
左边第一盏灯亮的:
当前灯的动作: 1->0; 0->1;
左边第一盏灯暗的:
当前灯的动作:
1->1; 0->0;
我们可以看到的是, 可以用一个方程来模拟这个过程: F[i] = ( f [ i] +f [i+n-2]%n+1 )%2;
所以我们只要计算这个值就OK啦。
然后由于数据过大,开数组肯定会爆掉~
这里我们要想到的是 矩阵是一个模拟递归的高效算法
这里我们要构造一个 可以计算如上的方程的矩阵:
1 0 0...0 1
1 1 0...0 0
0 1 1..0 0
0 0 1..0 0
. . . 0 ....
. . .0.....
0 0 0..0 1
*/
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
const int N = 105;
const int mod = 2;
struct Mat
{
int num[N][N];
Mat()
{
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++)
num[i][j] = 0;
}
};
Mat mul(Mat a, Mat b, int n)
{
Mat r;
for(int i = 0; i < n; i++)
for(int k = 0; k < n; k++)
{
if(a.num[i][k] == 0)
continue;
for(int j = 0; j < n; j++)
{
if(b.num[k][j] == 0)
continue;
r.num[i][j] = (r.num[i][j] + a.num[i][k] * b.num[k][j]) % mod;
}
}
return r;
}
Mat mal(Mat init, Mat unit, int m, int n)
{
while(m)
{
if(m & 1)
{
unit = mul(init, unit, n);
m--;
}
else
{
init = mul(init, init, n);
m >>= 1;
}
}
return unit;
}
int main()
{
int m;
while(cin >> m)
{
char str[N];
cin >> str;
Mat init, unit;
int len = strlen(str);
for(int i = 0; i < len; i++)
{
unit.num[i][0] = str[i] - '0';
}
init.num[0][0] = 1; init.num[0][len - 1] = 1;
int k = 0;
for(int i = 1; i < len; i++)
{
for(int j = k; j < k + 2; j++)
init.num[i][j] = 1;
k++;
if(k == len - 1)
break;
}
Mat ans;
ans = mal(init, unit, m, len);
for(int i = 0 ; i < len; i++)
cout << ans.num[i][0];
cout << endl;
}
}
Tower
/*
由 a(n) = 2*a(2)*a(n-1) - a(n-2)
let p = 2*a(2);
==>a(n) =p*a(n-1) - a(n-2)
==>a(n)^2 = p^2*a(n-1)^2 - 2*p*a(n-1)*a(n-2) + a(n-2)^2................(1)
let s(n) =sum( a(i) );
s(n) = s(n-1) + a(n)^2 ................................................................... (2)
a(n)*a(n-1) = p*a(n-1)^2 - a(n-1)*a(n-2)...................(3)
so we have:
a(n-1)^2 0 , 1 , 0 , 0 a(n-2)^2
a(n)^2 = 1 , p^2, -2*p, 0 * a(n-1)^2
a(n)*a(n-1) 0 , p , -1 , 0 a(n-1)*a(n-2)
s(n-1) 0 , 1 , 0 , 1 s(n-2)
use matrix to do this.
Problem is sovle now.
*/
#include<cstdio>
#define FOR(i,a,b) for(int (i)=(a);(i)<=(b);(i)++)
#define nMax 10
#define LL long long
LL mod;
template<typename T>
struct Matrix{
T mtix[nMax][nMax];
int n;
Matrix(int n=1):n(n) {
FOR(i,0,n-1) FOR(j,0,n-1) mtix[i][j] = 0;
FOR(i,0,n-1) mtix[i][i] = 1;
}
friend Matrix operator * (const Matrix& a,const Matrix& b) {
Matrix c(a.n);
FOR(i,0,a.n-1) FOR(j,0,a.n-1) {
c.mtix[i][j] = 0;
FOR(k,0,a.n-1) if(a.mtix[i][k] && b.mtix[k][j]){
c.mtix[i][j] += a.mtix[i][k] * b.mtix[k][j] ;
c.mtix[i][j]%=mod;
}
}
return c;
}
friend Matrix Exp(Matrix a,int n){
Matrix c(a.n);
Matrix b = a;
while(n){
if(n&1) c=c*b;
n >>= 1;
b = b*b;
}
return c;
}
void init(T p) {
p%=mod;
FOR(i,0,n-1) FOR(j,0,n-1) mtix[i][j] = 0 ;
mtix[1][1] = p*p%mod;
mtix[0][1] = mtix[1][0] = mtix[3][1] = mtix[3][3] = 1;
mtix[1][2] = ((-2*p) % mod + mod) % mod;
mtix[2][1] = p;
mtix[2][2] = ((-1 ) % mod + mod) % mod;
}
};
LL b[5];
int main(){
//freopen("input.txt","r",stdin);
int n;
LL a;
int t;
scanf("%d",&t);
while(t--){
scanf("%I64d%d%I64d",&a,&n,&mod);
a%=mod;
if(mod == 1){
printf("0\n");
continue;
}
Matrix<LL> s(4);
s.init(a*2);
b[1]=a*a%mod,b[0]=1LL,b[2]=a,b[3]=1;
s = Exp( s ,n-1 );
LL ans = 0LL;
FOR(i,0,3) ans = (ans+s.mtix[3][i]*b[i]) % mod;
printf("%I64d\n",ans);
}
return 0;
}
Lucky Coins Sequence
//参考文献http://blog.sina.com.cn/s/blog_83ccc39d0101jv6d.html
//代码参考http://blog.csdn.net/acm_davidcn/article/details/5819193
//矩阵快速幂
#include<iostream>
#define MOD 10007
using namespace std;
int f[4][4], res[4][4];
void A(int a[4][4], int b[4][4])
{
int t[4][4];
memset(t, 0, sizeof(t));
int i, j, k;
for( i=1; i < 4; i++ )
for( j=1; j < 4; j++ )
for( k=1; k < 4; k++ )
t[i][j] += a[i][k] * b[k][j] % MOD;
for( i=1; i < 4; i++ )
for( j=1; j < 4; j++ )
a[i][j] = t[i][j];
}
int main()
{
int n;
int a[5] = {0, 0, 0, 2, 6};
while(scanf("%d", &n) != EOF )
{
if(n <= 4) { printf("%d\n", a[n]); continue;}
memset(f, 0, sizeof(f));
memset(res, 0, sizeof(res));
f[1][1] = f[1][2] = f[1][3] = f[2][1] = 1;
f[3][3] = 2;
res[1][1] = res[2][2] = res[3][3] = 1;
n -= 4;
while(n > 0)
{
if(n % 2 == 0)
n /=2, A(f, f);
else
n--, A(res, f);
}
printf("%d\n", (res[1][1]*6 + res[1][2]*2 + res[1][3]*8) % MOD);
}
return 0;
}
Buge's Fibonacci Number Problem
#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
const int mm=55;
__int64 x[mm][mm],y[mm][mm],s[mm][mm],A[mm],B[mm],F1[mm],F2[mm],tmp;
int f1,f2,a,b,K,n,m;
int i,j,k,t;
void mul(__int64 x[mm][mm],__int64 y[mm][mm])
{
int i,j,k;
for(i=0;i<K+2;++i)
for(j=0;j<K+2;++j)
for(s[i][j]=k=0;k<K+2;++k)
{
s[i][j]+=x[i][k]*y[k][j];
if(s[i][j]>m)s[i][j]%=m;
}
for(i=0;i<K+2;++i)
for(j=0;j<K+2;++j)
x[i][j]=s[i][j];
}
int main()
{
scanf("%d",&t);
while(t--)
{
scanf("%d%d%d%d%d%d%d",&f1,&f2,&a,&b,&K,&n,&m);
A[0]=B[0]=F1[0]=F2[0]=1;
for(i=1;i<K+2;++i)
{
A[i]=A[i-1]*a;
B[i]=B[i-1]*b;
F1[i]=F1[i-1]*f1;
F2[i]=F2[i-1]*f2;
if(A[i]>=m)A[i]%=m;
if(B[i]>=m)B[i]%=m;
if(F1[i]>=m)F1[i]%=m;
if(F2[i]>=m)F2[i]%=m;
}
memset(x,0,sizeof(x));
memset(y,0,sizeof(y));
for(i=0;i<K+2;++i)x[i][i]=1;
y[0][0]=y[0][1]=1;
for(i=1;i<K+2;++i)
for(j=1,tmp=1,k=K-i+2;k--;++j)
{
y[i][j]=(tmp%m)*A[k];
if(y[i][j]>=m)y[i][j]%=m;
y[i][j]*=B[j-1];
if(y[i][j]>=m)y[i][j]%=m;
tmp*=k;
tmp/=j;
}
--n;
while(n)
{
if(n&1)mul(x,y);
mul(y,y);
n>>=1;
}
tmp=F1[K]*x[0][0];
if(tmp>=m)tmp%=m;
for(i=1;i<K+2;++i)
{
tmp+=((F2[K-i+1]*F1[i-1])%m)*x[0][i];
if(tmp>=m)tmp%=m;
}
printf("%I64d\n",tmp);
}
return 0;
}