2019.01.19【雅礼集训2019Day1T3】math（矩阵快速幂）

描述

给出 n， m， x，你需要求出下列式子的值:
$$\sum _{{\sum }_{i=1}^m{k}_i=n}\prod _{i=1}^m\sin (k_i\ast x)$$
其中ki为正整数,由于答案非常大，你只需要输出答案(保证不为 0)的正负(如果是负数输出负号，否则输出正号)和从左往右第一个非 0 数位上的数字即可。

输入

第一行一个整数 T 表示数据组数。

对于每组数据，每行有两个整数 m,n 和一个两位小数 x。第一行一个整数 T 表示数据组数。

对于每组数据，每行有两个整数 m,n 和一个两位小数 x。

输出

输出共 T 行，每行两个字符表示答案。

样例输入

2
3 5 0.01
3 6 0.02

样例输出

+2
+4

提示

对于 10%的数据， n,m≤5

对于 30%的数据， n≤500,m≤10

对于 50%的数据， n≤10000

对于 100%的数据， T≤10, n ≤ 1e9 ,m≤30

---

解析：

没想到吧，这道题有隐含的递推关系。

首先要知道$\sin$和$\cos$的$K$倍角公式：
$$\forall K\in {\mathbb{N}}_{\ast },x\in \mathbb{R},$$$$\sin (K\ast x)=\sin (x)\cos ((K-1)\ast x)+\cos (x)\sin ((K-1)\ast x)$$$$\cos (K\ast x)=\cos (x)\cos ((K-1)\ast x)-\sin (x)\sin ((K-1)\ast x)$$

证明可以看这里：https://blog.csdn.net/zxyoi_dreamer/article/details/85956203

然后我们可以开始强行推式子了。

推导：

首先刚才提到的$K$倍角公式是要用的。这里先写上：
$$\sin (K\ast x)=\sin (x)\cos ((K-1)\ast x)+\cos (x)\sin ((K-1)\ast x)$$

但是令人窒息的是我们还有一个$\cos ((K-1)\ast x)$需要处理。

作如下推导：
$$\begin{array}{rl}\sin (x)\cos (K\ast x)=& \sin (x)(\cos (x)\cos ((K-1)\ast x)-\sin (x)\sin ((K-1)\ast x))\\ =& \sin (x)\cos (x)\cos ((K-1)\ast x)+({\cos }^2(x)-1)\sin ((K-1)\ast x)\\ =& \cos (x)(\sin (x)\cos ((K-1)\ast x)+\cos (x)\sin ((K-1)\ast x))-\sin ((K-1)\ast x)\\ =& \cos (x)\sin (K\ast x)-sin((K-1)\ast x)\end{array}$$

所以容易得到$$\sin (K\ast x)=2\cos (x)\sin ((K-1)\ast x)-\sin ((K-2)\ast x)$$

我们现在已经确定了$x$，设$f(n,m)$表示$n,m$为对应值的时候的答案。

考虑现在我们需要知道$f(n_0,m_0)$的答案。

我们只针对$k_{m_0}$不同讨论不同的转移：
1.$k_{m_0}=1$，显然贡献为$f(n_0-1,m_0-1)\times \sin (x)\to f(n_0,m_0)$
2.$k_{m_0}>1$，此时所有序列的贡献为$2\ast \cos (x)f(n_0-1，m_0)-f(n_0-2,m_0)\to f(n_0,m_0)$

所以我们需要维护一个$(2m)\ast (2m)$的矩阵，为了同时维护前两次转移出的矩阵，直接按照递推关系构造矩阵快速幂就可以了。

---

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define re register
#define gc getchar
#define pc putchar
#define cs const

inline int getint(){
	re int num;
	re char c;
	while(!isdigit(c=gc()));num=c^48;
	while(isdigit(c=gc()))num=(num+(num<<2)<<1)+(c^48);
	return num;
}

inline double getdb(){
	re double x,y=1.0;
	re char c;
	re bool f=0;
	while(!isdigit(c=gc()))if(c=='-')f=1;x=c^48;
	while(isdigit(c=gc()))x=(x*10)+(c^48);
	if(c!='.')return f?-x:x;
	while(isdigit(c=gc()))x+=(y/=10)*(c^48);
	return f?-x:x;
}

cs int M=62;
int msiz;
struct matrix{
	double a[M][M];
	matrix(){memset(a,0,sizeof a);}
	inline cs double *cs operator[](cs int &offset)cs{return a[offset];}
	inline double *cs operator[](cs int &offset){return a[offset];}
	inline friend matrix operator*(cs matrix &A,cs matrix &B){
		matrix C;
		for(int re i=1;i<=msiz;++i)
		for(int re j=1;j<=msiz;++j)if(fabs(A[i][j])>1e-6)
		for(int re k=1;k<=msiz;++k)
		C[i][k]+=A[i][j]*B[j][k];
		return C;
	}
}A,B;

inline matrix quickpow(matrix a,int b,matrix res){
	for(;b;b>>=1,a=a*a)if(b&1)res=res*a;
	return res;
}

int T;
int n,m;
double x;
signed main(){
	T=getint();
	while(T--){
		m=getint();
		n=getint();msiz=m<<1;
		x=getdb();
		A=matrix();B=matrix();
		for(int re i=m+1;i<=m*2;++i){
			A[i][i-m]=1;
			A[i-m][i]=-1;
			A[i][i]=2*cos(x);
			if(i<m*2)A[i][i+1]=sin(x);
		}
		B[1][1]=sin(x);
		B[1][m+1]=sin(2*x);
		B[1][m+2]=sin(x)*sin(x);
		B=quickpow(A,n-1,B);
		double res=B[1][m];
		if(res<0)pc('-'),res=-res;
		else pc('+');
		while(res<1)res*=10;
		while(res>=10)res/=10;
		pc(48^(int)floor(res));pc('\n');
	}
	return 0;
}