1. 程式人生 > >Codechef:Sine Partition Function/PARSIN(矩陣快速冪)

Codechef:Sine Partition Function/PARSIN(矩陣快速冪)

傳送門

題解:
f i , j = k

= 0 i f i , k sin
( ( j k ) x ) f_{i,j}=\sum_{k=0}^i f_{i,k} \sin((j-k)x)

不妨設 c i , j = k = 0 j f i , k cos ( ( j k ) x ) c_{i,j}=\sum_{k=0}^jf_{i,k}\cos((j-k)x)

可以推出:
f i , j = f i , j 1 cos ( x ) + c i 1 , j 1 sin ( x ) f_{i,j}=f_{i,j-1}\cos(x)+c_{i-1,j-1}\sin(x)
c i , j = c i , j 1 cos ( x ) f i , j 1 sin ( x ) + f i 1 , j c_{i,j}=c_{i,j-1}\cos(x)-f_{i,j-1}\sin(x)+f_{i-1,j}

直接矩陣快速冪就行了。

#include <bits/stdc++.h>
using namespace std;
typedef double DB;

const int M=66;

int n,m,tot; DB x,sn,cs;
struct matrix {
	DB a[M][M];
	matrix(int t=0) {
		memset(a,0,sizeof(a));
		for(int i=0;i<=tot;i++) a[i][i]=t;
 	}
	friend inline matrix operator *(const matrix &a,const matrix &b) {
		matrix c(0);
		for(int i=0;i<=tot;i++)
			for(int k=0;k<=tot;k++) if(a.a[i][k])
				for(int j=0;j<=tot;j++) if(b.a[k][j])
					c.a[i][j]+=a.a[i][k]*b.a[k][j];
		return c;
	}	
	friend inline matrix operator ^(matrix a,int b) {
		matrix c(1);
		for(;b;b>>=1,a=a*a) if(b&1) c=c*a;
		return c;
	}
} A,tr;

inline void solve() {
	scanf("%d%d%lf",&m,&n,&x);
	sn=sin(x); cs=cos(x);
	A=matrix(0);
	A.a[1][0]=1; A.a[1][m+2]=sn;
	tot=2*m+2;A.a[1][tot]=1; tr=matrix(0);
	for(int i=0;i<=m;i++) tr.a[i+m+1][i]=1;
	tr.a[tot][tot]=1;
	for(int i=1;i<=m;i++) {
		tr.a[i+m+1][i+m+1]=2*cs;
		tr.a[i][i+m+1]=-1;
		tr.a[i+m][i+m+1]=sn;
	}
	A=A*(tr^(n-1));
	printf("%.10f\n",A.a[1][2*m+1]);
}

int main() {
	int T; scanf("%d",&T);
	while(T--) solve();
}