N棱柱M染色问题

P.S.1 以下默认对点染色，要求任意一对被边直接连接的点异色
P.S.2 备选色少于4色的情况略特殊，但先对更简单，这里不做讨论（但最后的计算程序包括了这种情况）
首先将棱柱沿着其中一条侧棱剪开，并强制规定左右两端染色情况相同，把棱柱染色转化为线性的平面染色

接下来考虑这种“铁轨”形状的染色
我们可以按照铁轨两端的染色把其分为5类
1.\({ }^{A}_{B} {}^{...}_{...} { }^{A}_{B}\)（用了2色，同色相对）（就是上文强制要求的那种情况）

2.\({ }^{A}_{B} {}^{...}_{...} { }^{B}_{A}\)（用了2色，异色相对）

3.\({ }^{A}_{B} {}^{...}_{...} { }^{A}_{C}\ (also\ { }^{A}_{B} {}^{...}_{...} { }^{C}_{B})\)（用了3色，同色相对）

4.\({ }^{A}_{B} {}^{...}_{...} { }^{C}_{A}\ (also\ { }^{A}_{B} {}^{...}_{...} { }^{B}_{C})\)（用了3色，异色相对）

5.\({ }^{A}_{B} {}^{...}_{...} { }^{C}_{D}\)（用了4色）
（建议思考一下为什么只需这5种情况就够了）
现在考虑递推式，也就是在最后新加入一列后的变化
令\(f_{i,j}\)表示目前共 \(i\) 列的第 \(j\) 种情况 方案数
记\(s_k\)为第\(k\)种情况的简称

1. \(f_{i,1}=0\ +f_{i-1,2}\ +0\ +2(m-2)f_{i-1,4}\ +(m-2)(m-3)f_{i-1,5}\)
   :
   由于染色条件限制\({ }^{A}_{B}\)不能直接跟在\({ }^{A}_{B}\ or\ { }^{A}_{C}\)后面，\(s_1/s_3\)对\(s_1\)的贡献是\(0\);
   \(s_2\)只能直接接上\({ }^{A}_{B}\)，贡献即是它本身
   \(s_4\)中的 \(C\)可以有除 \(A,B\) 外的 \(m-2\) 种颜色供选择，同时 与\(C\)相对的是\(A\)还是\(B\) 也是不同的情况，于是前面有个系数 \(2(m-2)\)
   \(s_5\)中的 \(C,D\)的颜色选择有 \(A^{2}_{m-2}=(m-2)(m-3)\) 种，于是前面有个系数\((m-2)(m-3)\)
   （这里只分析\(s_1\)递推式的推导，后面的留给读者自己思考吧）

2. \(f_{i,2}=f_{i-1,1}\ +0\ +2(m-2)f_{i-1,3}\ +0\ +(m-2)(m-3)f_{i-1,5}\)

3. \(f_{i,3}=0\ +f_{i-1,2}\ +(m-2)f_{i-1,3}\ +[(m-2)+(m-3)]f_{i-1,4}\ +(m-3)^2f_{i-1,5}\)

4. \(f_{i,4}=f_{i-1,1}\ +0\ +[(m-2)+(m-3)]f_{i-1,3}\ +(m-2)f_{i-1,4}\ +(m-3)^2f_{i-1,5}\)

5. \(f_{i,5}=f_{i-1,1}\ +f_{i-1,2}\ +2(m-3)f_{i-1,3}\ +2(m-3)f_{i-1,4}\ +[(m-4)^2+(m-3)]f_{i-1,5}\)

令向量\(V_i=[f_{i,1}\ ,\ f_{i,2}\ ,\ f_{i,3}\ ,\ f_{i,4}\ ,\ f_{i,5}]\)
初始\(V_1=[1,0,0,0,0]\)
设递推矩阵\(A\)
那么上述递推式可写成矩阵乘法的模式：

\[ A=\left[ \begin{matrix} 0 & 1 & 0 & 1 & 1 \\ 1 & 0 & 1 & 0 & 1 \\ 0 & 2(m-2) & (m-2) & [(m-2)+(m-3)] & 2(m-3) \\ 2(m-2) & 0 & [(m-2)+(m-3)] & (m-2) & 2(m-3) \\ (m-2)(m-3) & (m-2)(m-3) & (m-3)^2 & (m-3)^2 & [(m-4)^2+(m-3)] \end{matrix} \right] \]

\[V_{i+1}=V_i\times A \]

\[V_i=V_1\times A^{i-1} \]

最后我们要的\(N\)棱柱\(M\)染色方案数就是\(f_{n,1}\)，也就是\(V_n\)向量的第一维
然后我们可以求\(A\)矩阵的特征向量和特征值来简化计算过程，但是解5次方程有点困难
但我们有计算机这个好帮手:)
以下是我写的丑陋代码（其中第2,3种情况倒了个位置），可计算\(n \leq 10^{18},m \leq 10^9\)的情况，答案对\(10^9+7\)取模

#include
#include
#include
using namespace std;
#define O 1000000007
long long n;int m;
struct Mt{
	int A[5][5];
	inline Mt operator * (const Mt& U)const{
		long long tmpA[5][5]={{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0}};
		register int i,j;
		for(i=0;i<5;i++)for(j=0;j<5;j++){
			tmpA[i][0]+=1ll*A[i][j]*U.A[j][0];
			tmpA[i][1]+=1ll*A[i][j]*U.A[j][1];
			tmpA[i][2]+=1ll*A[i][j]*U.A[j][2];
			tmpA[i][3]+=1ll*A[i][j]*U.A[j][3];
			tmpA[i][4]+=1ll*A[i][j]*U.A[j][4];
		}
		Mt res;
		for(i=0;i<5;i++){
			res.A[i][0]=tmpA[i][0]%O;
			res.A[i][1]=tmpA[i][1]%O;
			res.A[i][2]=tmpA[i][2]%O;
			res.A[i][3]=tmpA[i][3]%O;
			res.A[i][4]=tmpA[i][4]%O;
		}
		return res;
	}
	inline void operator *= (const Mt& U){
		long long tmpA[5][5]={{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0}};
		register int i,j;
		for(i=0;i<5;i++)for(j=0;j<5;j++){
			tmpA[i][0]+=1ll*A[i][j]*U.A[j][0];
			tmpA[i][1]+=1ll*A[i][j]*U.A[j][1];
			tmpA[i][2]+=1ll*A[i][j]*U.A[j][2];
			tmpA[i][3]+=1ll*A[i][j]*U.A[j][3];
			tmpA[i][4]+=1ll*A[i][j]*U.A[j][4];
		}
		for(i=0;i<5;i++){
			A[i][0]=tmpA[i][0]%O;
			A[i][1]=tmpA[i][1]%O;
			A[i][2]=tmpA[i][2]%O;
			A[i][3]=tmpA[i][3]%O;
			A[i][4]=tmpA[i][4]%O;
		}
	}
}un,e;//矩阵定义
struct vec{
	int A[5];
	inline vec operator * (const Mt& U){
		long long tmpA[5]={0,0,0,0,0};
		register int i;
		for(i=0;i<5;i++){
			tmpA[0]+=1ll*A[i]*U.A[i][0];
			tmpA[1]+=1ll*A[i]*U.A[i][1];
			tmpA[2]+=1ll*A[i]*U.A[i][2];
			tmpA[3]+=1ll*A[i]*U.A[i][3];
			tmpA[4]+=1ll*A[i]*U.A[i][4];
		}
		vec res;
		res.A[0]=tmpA[0]%O;
		res.A[1]=tmpA[1]%O;
		res.A[2]=tmpA[2]%O;
		res.A[3]=tmpA[3]%O;
		res.A[4]=tmpA[4]%O;
		return res;
	}
}ans;//向量定义
void Init(){
	ans=(vec){ {0,0,1,1,1} };/*n=2 初始情况*/	e=(Mt){ {{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,1,0},{0,0,0,0,1}} };/*单位矩阵*/
	un=(Mt){ 
		{
			{0,0,1,1,1},
			{0,m-2,(int)(2ll*(m-2)%O),(int)((2ll*m-5)%O),(int)(2ll*(m-3)%O)},
			{1,1,0,0,1},
			{(int)(2ll*(m-2)%O),(int)((2ll*m-5)%O),0,m-2,(int)(2ll*(m-3)%O)},
			{(int)(1ll*(m-2)*(m-3)%O),(int)(1ll*(m-3)*(m-3)%O),(int)(1ll*(m-2)*(m-3)%O),(int)(1ll*(m-3)*(m-3)%O),(int)((1ll*(m-4)*(m-4)%O+(m-3))%O)}
		} 
	};//递推矩阵
}
void Make3(){
	ans.A[4]=0;un.A[4][0]=un.A[4][1]=un.A[4][2]=un.A[4][3]=un.A[4][4]=un.A[0][4]=un.A[1][4]=un.A[2][4]=un.A[3][4]=0;
}
Mt qkpw(Mt ds,long long ix){
	Mt re=e;while(ix){
		if(ix&1)re*=ds;
		ix>>=1;ds*=ds;
	}return re;
}//快速幂
void Solve(){
	ans=ans*qkpw(un,n-1);
	int Ans=ans.A[0];
	Ans=1ll*Ans*m%O;
	Ans=1ll*Ans*(m-1)%O;
	cout<>n>>m;//n棱柱m染色
        if(n==1)return cout<<(1ll*(m-1)*m%O)<<'\n',0;
        if(n==2&&m==2)return cout<<"2\n",0;
	//if(n<3)return cout<<"cannot form a prism!\n",0;
	switch(m){
		case 1:return cout<<0<<'\n',0;
		case 2:return cout<<0<<'\n',0;
		case 3:return Init(),Make3(),Solve(),0;
		default:return Init(),Solve(),0;
	}
	return 0;
}