题目中还有一点比较有意思的结果，可以看到虽然题目中矩阵是7阶的，但是这个数列却是一个5阶线性递推数列，
于是可以猜测其中应该还有一些隐含的约束在里面。
如果我们把3#中每个模式反序观之，可以看出4和7是镜像偶，相应的计数函数可以互相代换，省掉一个，其它5个模式是镜像自对称的，计数函数也各不相等。
于是对应的矩阵刚好可以降低到6阶，7用4代换后的矩阵如下：
$M_2=[(6,0,0,7,0,7),(0,0,0,0,1,0),(0,0,0,1,0,0),(0,0,8,0,7,0),(1,0,0,0,0,1),(0,8,0,7,0,0)]$
经计算不出所料$M$和$M_2$的最小多项式都是$x^5−7x^4−x^3−57x^2+64=(x-8)(x-1)(x+1)(x^2+x+8)$,但是出乎意料的是$M_2$的特征多项式不是$M$的特征多项式的因子
其中$M$的特征多项式是$(x-8)(x-1)(x+1)(x^2+x+8)^2$,而$M_2$的特征多项式是$(x-8)(x-1)(x+1)^2(x^2+x+8)$
ps: 最出人意料的是只有1对计数相等的模式，而不是两对从而降为5阶矩阵。
另外一般情况，我们把10个数改为T个数，那么递推数列的特征多项式就是$(x-T+2)(x-1)(x+1)(x^2+x+T-2)$
终于理解了到矩阵计算这一步， 但是下面一句
'经计算可以知道p(n)的特征多项式为" 不懂， 特征多项式跟递推公式和矩阵结果存在什么关系？ 是线性代数的书有讲吗？
另外， 这个递推公式能转换成通项公式吗？
比如数值表示可以用Pari/Gp如下计算:
? b=
%1 =
? v=~
%2 = ~
? m=
%3 =
? mateigen(m)
%4 =
[-2 0 7 -0.73224568138195777351247600767754318618 + 0.63585792627329425022285708692064049226*I -0.73224568138195777351247600767754318618 - 0.63585792627329425022285708692064049226*I]
[-1 1 1/8 -0.0037188099808061420345489443378119001906 - 0.086161420681990502393643449803322083513*I -0.0037188099808061420345489443378119001906 + 0.086161420681990502393643449803322083513*I]
[-1 -1 1/8 0.10832533589251439539347408829174664107 - 0.0046754259284801047810504197567694153845*I 0.10832533589251439539347408829174664107 + 0.0046754259284801047810504197567694153845*I]
? P=%4;
? P^-1*m*P
%6 =
[-1.0000000000000000000000000000000000000 + 0.E-37*I -1.4693679385278593850 E-39 + 0.E-37*I -3.673419846319648463 E-38 + 0.E-38*I 1.4693679385278593850 E-39 + 8.265194654219209041 E-40*I -7.346839692639296925 E-39 - 2.3877229001077715006 E-39*I]
[-5.142787784847507848 E-39 + 0.E-37*I 1.0000000000000000000000000000000000000 + 0.E-37*I -6.538687326448974264 E-38 - 3.159141067834897678 E-38*I 2.5713938924237539236 E-39 + 7.255004196481305714 E-39*I 9.550891600431086002 E-39 - 7.255004196481305714 E-39*I]
[-8.816207631167156310 E-39 + 4.353682780823287065 E-40*I -8.816207631167156310 E-39 + 1.0884206952058217666 E-39*I 7.9999999999999999999999999999999999999 - 9.360417978770067193 E-39*I -2.938735877055718770 E-39 + 7.074734518837841484 E-40*I 0.E-38 - 1.8639204405399697752 E-39*I]
[-1.7632415262334312620 E-38 + 5.289724578700293786 E-38*I 1.7632415262334312620 E-38 + 5.289724578700293786 E-38*I 5.877471754111437540 E-39 + 5.289724578700293786 E-38*I -1.4693679385278593850 E-38 - 5.877471754111437540 E-39*I -0.50000000000000000000000000000000000000 + 2.7838821814150109610597356494592747603*I]
? P^-1*v
%7 = ~
? b*P
%8 = [-1, 1, 8, 0.26775431861804222648752399232245681382 + 0.63585792627329425022285708692064049226*I, 0.26775431861804222648752399232245681382 - 0.63585792627329425022285708692064049226*I]
?
由此可以看出，由于数列通项为$b*M^{n-3}*v' = b*P*(P^{-1}MP)^{n-3}*P^{-1}v'$
其中$b*P=[-1, 1, 8, 0.26775431861804222648752399232245681382 + 0.63585792627329425022285708692064049226*I, 0.26775431861804222648752399232245681382 - 0.63585792627329425022285708692064049226*I]$
$P^{-1}v'= ~$
而$P^{-1}MP$是对角阵$diag{-1，1，8， -0.5-2.7838821814150109610597356494592747602*I，-0.5+2.7838821814150109610597356494592747602*I}$
表示$a(n)=-35*(-1)^(n-3)+36+512*8^(n-3)+(103.5+175.38457742914569054676334591593430989i)*(-0.5-2.7838821814150109610597356494592747602i)^(n-3)+(103.5-175.38457742914569054676334591593430989i)*(-0.5+2.7838821814150109610597356494592747602i)^(n-3)$
只是很奇怪，这里P不是方阵,Pari-gp是如何算出$P^{-1}$的而且最后通项好像也没有错
如果我们设$r_1,r_2$分别是方程$x^2+x+8=0$的根，那么通解应该是
$a(n)=-35*(-1)^{n-3}+36+512*8^{n-3}+(135+63*r_2)r_1^{n-3}+(135+63*r_1)r_2^{n-3}$
上面表达式还可以写成
$a(n)=36-35*(-1)^{n-3}+512*8^{n-3}+135*(r_1^{n-3}+r_2^{n-3})+504*(r_1^{n-4}+r_2^{n-4})=36-35*(-1)^{n-3}+512*8^{n-3}+135b(n-3)+504b(n-4)$
其中$b(n)$是二阶递推数列$b(0)=2,b(1)=-1,b(n+2)=-b(n+1)-8b(n)$
或者
$a(n)=36+35*(-1)^n+8^n+270*\sqrt{8}^{n-3}*cos((n-3)(\pi-arctan(\sqrt{31})))+1008*\sqrt{8}^{n-4}*cos((n-4)(\pi-arctan(\sqrt{31})))$
楼上计算过程说明准确计算过程需要将矩阵转化为Jordan标准型，但是pari/gp好像没有对应的函数。
mathe 发表于 2017-3-30 12:48
题目中的进一步要求把片断长度改为k会导致计算过程快速复杂化。
比如题目中k从3变化到4，可以有34种不同的模式，也就是题目需要处理一个34阶矩阵了.
k=4时，34种模式中，有10个镜像对称偶，剩余14个为镜像自对称。故递推矩阵可降为24阶。
手工排的，不知道有没有错误。如果没有错误的话，猜想k=5时，或许递推矩阵阶数可降为5！:Q:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#define K 4
#define NUM ((K)-1)
#define MASK_FOR_LAST_N_BITS(n) ((1<<(n))-1)
#define MASK_FOR_BEFORE_N_BITS(n) (~MASK_FOR_LAST_N_BITS(n))
#define IS_BIT_SET(x, b) (((x)&(1<<(b)))!=0)
int match_count;
int total_count;
void output(char map)
{
int unused=NUM;
int i;
char map2;
int unmatched_found=0;
memset(map2,-1,sizeof(map2));
for(i=0;i<NUM;i++){
printf("%c", i+'A');
}
printf("==>");
for(i=0;i<NUM;i++){
if(map!=-1){
map2]=NUM-1-i;
printf("%c", map+'A');
}else{
printf("%c", unused++ + 'A');
}
}
printf("R");
unused=NUM;
for(i=NUM-1;i>=0;i--){
if(map2!=map)unmatched_found=1;
if(map2!=-1){
printf("%c", map2+'A');
}else{
printf("%c", unused++ + 'A');
}
}
printf("\n");
if(!unmatched_found){
match_count++;
}
total_count++;
}
void recursive(int old_state, int search_bit, char map)
{
if((old_state&MASK_FOR_BEFORE_N_BITS(search_bit))==0){//finished search
output(map);
return;
}
while(!IS_BIT_SET(old_state, search_bit))search_bit++;
int i;
for(i=0;i<NUM;i++){
if(map==-1){
map=search_bit;
recursive(old_state, search_bit+1, map);
map=-1;
}
}
}
void search(int old_state)/*bit state to show whether a old state is used*/
{
char map;
memset(map, -1, sizeof(map));
recursive(old_state, 0, map);
}
int main()
{
int i;
for(i=0; i<(1<<NUM);i++){
search(i);
}
printf("Total %d states, %d mapped to itself and total %d states after removing duplicated\n", total_count, match_count, match_count+(total_count-match_count)/2);
return 0;
}
output:
ABC==>DEFRDEF
ABC==>ADERDEC
ABC==>DAERDEB
ABC==>DEARDEA
ABC==>BDERDCE
ABC==>DBERDBE
ABC==>DEBRDAE
ABC==>ABDRDBC
ABC==>ADBRDAC
ABC==>BADRDCB
ABC==>DABRDAB
ABC==>BDARDCA
ABC==>DBARDBA
ABC==>CDERCDE
ABC==>DCERBDE
ABC==>DECRADE
ABC==>ACDRBDC
ABC==>ADCRADC
ABC==>CADRCDB
ABC==>DACRADB
ABC==>CDARCDA
ABC==>DCARBDA
ABC==>BCDRBCD
ABC==>BDCRACD
ABC==>CBDRCBD
ABC==>DBCRABD
ABC==>CDBRCAD
ABC==>DCBRBAD
ABC==>ABCRABC
ABC==>ACBRBAC
ABC==>BACRACB
ABC==>CABRCAB
ABC==>BCARBCA
ABC==>CBARCBA
K=5的状态算出来了，共209，去除对称性还有126个状态，当然最终递推数列的阶应该会更低一些
稍微修改上面代码，可以将两个矩阵输出成gp格式，然后用gp求两矩阵特征多项式并且求它们公因子
可以得到K=4时，数列特征多项式为14次（暗示最终结果等于自对称的数目?）多项式$(x-7)(x+1)^2(x^2-x+1)(x^3-x^2+x-7)(x^3+x^2-x-7)(x^3+x^2+7*x+49)$
同样，如果将10个数改为T个数，特征多项式为$(x-T+3)(x+1)^2(x^2-x+1)(x^3-x^2+x-T+3)(x^3+x^2-x-T+3)(x^3+x^2+(T-3)*x+(T-3)^2)$
而对于K=5, 计算特征多项式在我的笔记本上垮了(应该是gp分配的内存不够)
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <map>
#include <string>
using namespace std;
#define K 5
#define NUM ((K)-1)
#define MASK_FOR_LAST_N_BITS(n) ((1<<(n))-1)
#define MASK_FOR_BEFORE_N_BITS(n) (~MASK_FOR_LAST_N_BITS(n))
#define IS_BIT_SET(x, b) (((x)&(1<<(b)))!=0)
typedef map<string, string> SMap;
typedef map<string, int> IMap;
SMap string_map;
IMap imap1, imap2;
int ic1,ic2;
void decode_string(const string& s, int& present, int& last_used)
{
present = 0;
last_used=0;
int i;
for(i=0;i<NUM;i++){
int d=s-'A';
if(d>=NUM)last_used=d-NUM+1;
else present |= 1<<d;
}
}
string shift_string(const string& s, char& new_char)
{
char c=s;
int first_new=0;
if(c-'A'>=NUM){
first_new=1;
}
int i;
string out;
for(i=1;i<NUM;i++){
c=s;
if(c-'A'>=NUM){
c-=first_new;
}
out+=c;
}
if(first_new)new_char--;
return out;
}
void gen_sum(void)
{
printf("T=10;\n");
printf("m1=matrix(%d,%d);\n",ic1,ic1);
printf("m2=matrix(%d,%d);\n",ic2,ic2);
IMap::iterator it;
for(it=imap1.begin(); it!=imap1.end();++it){
int from = it->second;
string s=it->first;
int present, last_used;
decode_string(s, present, last_used);
int i;
char new_char='A'+NUM+last_used;
string t = shift_string(s, new_char);
for(i=0;i<NUM;i++){
if(!IS_BIT_SET(present, i)){
string h=t;
h+=(char)('A'+i);
IMap::iterator it2 = imap1.find(h);
int to = it2->second;
printf("m1[%d,%d]=1;\n", to+1, from+1);
}
}
t+=new_char;
IMap::iterator it2=imap1.find(t);
int to = it2->second;
printf("m1[%d,%d]=T-%d;\n", to+1, from+1, last_used+NUM);
}
for(it=imap1.begin(); it!=imap1.end();++it){
string s = it->first;
IMap::iterator it2 = imap2.find(s);
if(it2==imap2.end()){
SMap::iterator sit = string_map.find(s);
string ns = sit->second;
it2 = imap2.find(ns);
}
int from = it2->second;
int present, last_used;
decode_string(s, present, last_used);
int i;
char new_char='A'+NUM+last_used;
string t = shift_string(s, new_char);
for(i=0;i<NUM;i++){
if(!IS_BIT_SET(present, i)){
string h=t;
h+=(char)('A'+i);
IMap::iterator it2=imap2.find(h);
if(it2 != imap2.end()){
int to=it2->second;
printf("m2[%d,%d]+=1;\n", to+1, from+1);
}
}
}
t+=new_char;
it2=imap2.find(t);
if(it2!=imap2.end()){
int to = it2->second;
printf("m2[%d,%d]+=T-%d;\n", to+1, from+1, last_used+NUM);
}
}
}
void output(char map)
{
int unused=NUM;
int i;
char map2;
int unmatched_found=0;
memset(map2,-1,sizeof(map2));
string s1,s2;
for(i=0;i<NUM;i++){
if(map!=-1){
map2]=NUM-1-i;
s1+=(char)(map+'A');
}else{
s1+=(char)(unused++ + 'A');
}
}
unused=NUM;
for(i=NUM-1;i>=0;i--){
if(map2!=map)unmatched_found=1;
if(map2!=-1){
s2+= (char)(map2+'A');
}else{
s2+=(char)(unused++ + 'A');
}
}
string_map.insert(make_pair(s1,s2));
imap1.insert(make_pair(s1, ic1)); ic1++;
if(s1<=s2){
imap2.insert(make_pair(s1,ic2)); ic2++;
}
}
void recursive(int old_state, int search_bit, char map)
{
if((old_state&MASK_FOR_BEFORE_N_BITS(search_bit))==0){//finished search
output(map);
return;
}
while(!IS_BIT_SET(old_state, search_bit))search_bit++;
int i;
for(i=0;i<NUM;i++){
if(map==-1){
map=search_bit;
recursive(old_state, search_bit+1, map);
map=-1;
}
}
}
void search(int old_state)/*bit state to show whether a old state is used*/
{
char map;
memset(map, -1, sizeof(map));
recursive(old_state, 0, map);
}
int main()
{
int i;
for(i=0; i<(1<<NUM);i++){
search(i);
}
gen_sum();
return 0;
}
现在我们以K=3为例子分析为什么最终阶数等于自对称状态数目。因为AC和CB相等，除了抛弃状态CB外，利用AC和CB递推式中右边相等，我们还可以将AC表示为其它各项(不包括CB)的线性组合AC(n)=8BA(n)+7CA(n)-CD(n)，于是我们还可以将递推式中所有AC(n)继续用这个线性组合替换，从而可以继续抛弃AC
利用这种方法对于K=3,我们可以得出通项公式
$a(n)=[(1,0,0,0,1)][(-1,0,56,49,7),(0,0,0,1,0),(-1,0,8,7,0),(1,0,0,0,1),(-7,8,56,49,0)]^{n-3}[(0),(0),(0),(0),(720)]$