11149 - Power of Matrix

[Darko]

This question will probably be silly, but here it goes...

I've never done the matrix exponentiation, so I didn't even try this problem. But what I did was try if what works for integers worked for that matrix in the example. Well, it did work, but I am not sure I could really generalize (because matrix multiplication is not commutative), although it doesn't seem that crazy.

In that example, you can get the answer using the following equation:
x^(k+1) - 1 = (x - 1) * (x^k + ... + x + 1)

Now, it works for 1x1 matrices (heh) and I just tried it for kicks, but it did work for that matrix ({{0,2,0},{0,0,2},{0,0,0}}). To be honest, I was surprised I got the right result. And I am not sure how to deal with a more complicated inverse, this one was nice.

So - can the equation above be generalized for n>1?

[Observer]

So - can the equation above be generalized for n>1?

The answer is YES.

Recall that for the formula

Code: Select all

                    x^(k+1) - 1
x^k + ... + x + 1 = -----------
                       x - 1

It does not work for x = 1 (in which case denominator = 0).

Now back to the matrix case. The formula below works (you can check it by direct computation):

Code: Select all

(A - I)(A^k + ... + I) = A^(k+1) - I
 

Therefore if A-I is non-singular (i.e. det(A-I) != 0 <=> inverse exists), then

Code: Select all

(A^k + ... + I) = (A - I)^(-1) * (A^(k+1) - I),
                              where (A - I)^(-1) means the inverse of A-I.

But if A-I is singular, then you'll have to think of other methods~

[cytmike]

removed

[Darko]

Duh, of course it works :)
(because matrix multiplication is associative)

I still have to deal with that inverse, but I will figure it out eventually.

[EonStrife]

The formula I am using is almost similar :
Sum = A * ((A^n)-I) * (A-I)^(-1)

As for the matrix inversion, I have done that (but not sure if It can handle that singular thingy). The problem is, it involves division and floating point, so the final result sometimes off -1 or 1 from desired result (Hello, WA)

[ibroker]

I use operator overloading calculating matrix.
And I divide equation like this.
A^1+A^2+...A^(2k) = A+A^2+...A^k + A^k(A+A^2..._A^k)
It throws me TLE.
What can I do to cut running time?

Code: Select all

[/code]

[Jan]

Code: Select all

a = calc(k/2);

You are calling it for every recursive call from divide. Try to store the values. Suppose you have stored A^k, then you can find A^(2k) and A^(2k+1) easily and efficiently (Just use the A^k you have stored already). I used this approach and got accepted.

[ibroker]

Thanks ! I got Accepted

[Sanny]

Are the following i/o correct?

Input:

Code: Select all

1 5
5

5 512
1 2 3 4 5
6 7 8 9 10
11 12 13 14 15
16 17 18 19 20
21 22 23 24 25

4 15
1 2 5 9
8 5 6 7
3 4 8 9
1 3 4 6

3 4095
1 2 3
4 5 6
7 8 9

0 0

Output:

Code: Select all

5

6 2 8 4 0
6 7 8 9 0
6 2 8 4 0
6 7 8 9 0
6 2 8 4 0

8 0 9 4
3 8 9 9
9 5 1 6
0 7 4 5

5 8 1
8 3 8
1 8 5

Thanks in advance.

[coolguy]

yes they are correct

[Sanny]

Some more i/o to verify:

Input:

Code: Select all

6 524289
1 2 3 4 5 6
6 5 4 3 2 1
1 3 5 2 4 6
1 3 6 2 5 4
6 2 5 1 3 4
3 4 5 6 2 1 

10 524287
45 46 132 75 45 62 21 19 38 67
45 46 132 75 45 62 21 19 38 67
45 46 132 75 45 62 21 19 38 67
45 46 132 75 45 62 21 19 38 67
45 46 132 75 45 62 21 19 38 67
45 46 132 75 45 62 21 19 38 67
45 46 132 75 45 62 21 19 38 67
45 46 132 75 45 62 21 19 38 67
45 46 132 75 45 62 21 19 38 67
45 46 132 75 45 62 21 19 38 67

2 524287
1 2
3 4

Output:

Code: Select all

6 0 2 1 9 1
1 3 6 4 6 9
6 6 0 7 8 2
7 9 7 0 7 9
2 3 2 7 5 0
2 3 7 5 4 8

5 6 2 5 5 2 1 9 8 7
5 6 2 5 5 2 1 9 8 7
5 6 2 5 5 2 1 9 8 7
5 6 2 5 5 2 1 9 8 7
5 6 2 5 5 2 1 9 8 7
5 6 2 5 5 2 1 9 8 7
5 6 2 5 5 2 1 9 8 7
5 6 2 5 5 2 1 9 8 7
5 6 2 5 5 2 1 9 8 7
5 6 2 5 5 2 1 9 8 7

9 0
5 4

[sohel]

These are also correct !!

[Sanny]

Ok now it is time to post the code i think. I've spent many hours with it. This program gets WA in ~.027 seconds.

Code: Select all

#pragma warning(disable:4786)
#include<vector>
#include<cstdio>
#include<map>
#include<algorithm>
#include<cassert>
using namespace std;

#define rep(i,n) for(i=0;i<(n);i++)

struct matrix
{
	int a[42][42];
}A,I;
int p;
map<int,matrix>mp;//stores different A^n
vector<int>v,vs;
map<int,matrix>S;//S[n] = A + A^2 + ... + A^n

matrix Sum(matrix a,matrix b)
{
	int i,j;
	matrix c;
	rep(i,p) rep(j,p) c.a[i][j] = (a.a[i][j] + b.a[i][j])%10;
	return c;
}

matrix Multiply(matrix a,matrix b)
{
	int i,j,k;
	matrix c;
	rep(i,p) rep(j,p) {
		c.a[i][j] = 0;
		rep(k,p) c.a[i][j] = (c.a[i][j] + a.a[i][k] * b.a[k][j])%10;
	}
	return c;
}

matrix Pow(int n)//finds A^n
{
	if(n==0) return I;
	if(n==1) return A;
	if(mp.find(n) != mp.end()) return mp[n];

	if(n&1) return mp[n] = Multiply(A,Pow(n-1));
	else {
		matrix n2;
		n2 = Pow(n/2);
		return mp[n] = Multiply(n2,n2);
	}
}


//S[n] = S[n/2] + A^(n/2) * S[n/2]; when n is even
//S[n] = A + A * S[n-1]; when n is odd
matrix find_sum(int n)
{
	if(n==0) return I;
	if(n==1) return A;
	if(S.find(n) != S.end()) return S[n];

	if(n&1) return S[n] = Sum(A,Multiply(A,find_sum(n-1)));
	else {
		matrix n2;
		n2 = find_sum(n/2);
		return S[n] = Sum(n2,Multiply(Pow(n/2),n2));
	}
}

//this is used to find which S[n] & A[n] we need in the calculation
void Dummy(int n)
{
	vs.push_back(n);//keeps the S[n]s
	if(n<2) return;
	if(n&1) Dummy(n-1);
	else {
		v.push_back(n/2);//keeps the A[n]s
		Dummy(n/2);
	}
}

int main()
{
	int i,j,n;
	matrix r;
	//p is the size of matrix & n is the power
	while(scanf(" %d %d",&p,&n)==2 && p) {
		assert(p>0 && n>0);
		
		//identity matrix
		rep(i,p) rep(j,p) I.a[i][j] = 0;
		rep(i,p) I.a[i][i] = 1;

		rep(i,p) rep(j,p) {
			scanf("%d",&A.a[i][j]);
			A.a[i][j]%=10;
		}

		mp.clear(); v.clear(); S.clear(); vs.clear();
		mp[0] = I, mp[1] = A;
		S[0] = I, S[1] = A;
		Dummy(n);
		sort(v.begin(),v.end());
		rep(i,v.size()) mp[v[i]] = Pow(v[i]);
		sort(vs.begin(),vs.end());
		rep(i,vs.size()) S[vs[i]] = find_sum(vs[i]);

		r = find_sum(n);
		for(i=0;i<p;i++,printf("\n")) for(printf("%d",r.a[i][0]),j=1;j<p;j++) printf(" %d",r.a[i][j]);
		printf("\n");
	}
	return 0;
}

[sohel]

I just commented out these five lines

Code: Select all

      Dummy(n); 
      sort(v.begin(),v.end()); 
      rep(i,v.size()) mp[v[i]] = Pow(v[i]); 
      sort(vs.begin(),vs.end()); 
      rep(i,vs.size()) S[vs[i]] = find_sum(vs[i]); 

And it gets AC.

But nothing seems to be wrong with the Dummy()... strange...

[Sanny]

When I remove all the 5 lines, I get AC in 1.473 sec.
If I remove last 2 lines, I get WA in 0.795 sec.
And if I remove 3rd & 4th line only, I get WA in 0.021 sec.

Can anybody explain these?