这样的十进数有多少个？（排列问题）

mathe

题目中还有一点比较有意思的结果，可以看到虽然题目中矩阵是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)$

jjnet

终于理解了到矩阵计算这一步， 但是下面一句
'经计算可以知道p(n)的特征多项式为"   不懂， 特征多项式跟递推公式和矩阵结果存在什么关系？ 是线性代数的书有讲吗？

jjnet

另外， 这个递推公式能转换成通项公式吗？

mathe

比如数值表示可以用Pari/Gp如下计算:
? b=[1,0,0,0,0,1]
%1 = [1, 0, 0, 0, 0, 1]
? v=[0,0,0,0,0,720]~
%2 = [0, 0, 0, 0, 0, 720]~
? m=[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]
%3 =
[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]

? 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]

[1 -1 1 -0.067178502879078694817658349328214971211 - 0.29922725942272670598722686443324258459*I -0.067178502879078694817658349328214971211 + 0.29922725942272670598722686443324258459*I]

[1 1 1 -0.23800383877159309021113243761996161228 + 0.053433439182629768926290511505936175822*I -0.23800383877159309021113243761996161228 - 0.053433439182629768926290511505936175822*I]

[1 1 1 1 1]

? 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]

[2.938735877055718770 E-39 - 1.7632415262334312620 E-38*I 1.4693679385278593850 E-38 - 1.7632415262334312620 E-38*I 2.938735877055718770 E-39 - 1.7632415262334312620 E-38*I -0.50000000000000000000000000000000000002 - 2.7838821814150109610597356494592747602*I -2.350988701644575016 E-38 - 2.938735877055718770 E-38*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 = [34.999999999999999999999999999999999999 + 1.0579449157400587572 E-36*I, 35.999999999999999999999999999999999998 - 1.0579449157400587572 E-36*I, 64.000000000000000000000000000000000000 + 3.134651602192766688 E-37*I, 292.50000000000000000000000000000000000 - 39.602969096903865607333658755210973202*I, 292.50000000000000000000000000000000000 + 39.602969096903865607333658755210973202*I]~
? 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'= [34.999999999999999999999999999999999999 + 1.0579449157400587572 E-36*I, 35.999999999999999999999999999999999998 - 1.0579449157400587572 E-36*I, 64.000000000000000000000000000000000000 + 3.134651602192766688 E-37*I, 292.50000000000000000000000000000000000 - 39.602969096903865607333658755210973202*I, 292.50000000000000000000000000000000000 + 39.602969096903865607333658755210973202*I]~$
而$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})))$

mathe

楼上计算过程说明准确计算过程需要将矩阵转化为Jordan标准型，但是pari/gp好像没有对应的函数。

hujunhua

mathe 发表于 2017-3-30 12:48
题目中的进一步要求把片断长度改为k会导致计算过程快速复杂化。
比如题目中k从3变化到4，可以有34种不同的模式，也就是题目需要处理一个34阶矩阵了.

k=4时，34种模式中，有10个镜像对称偶，剩余14个为镜像自对称。故递推矩阵可降为24阶。

手工排的，不知道有没有错误。如果没有错误的话，猜想k=5时，或许递推矩阵阶数可降为5！

mathe

1. #include <stdio.h>
   
2. #include <string.h>
   
3. #include <stdlib.h>
   
4. #define K 4
   
5. #define NUM ((K)-1)
   
6. #define MASK_FOR_LAST_N_BITS(n)   ((1<<(n))-1)
   
7. #define MASK_FOR_BEFORE_N_BITS(n) (~MASK_FOR_LAST_N_BITS(n))
   
8. #define IS_BIT_SET(x, b)   (((x)&(1<<(b)))!=0)
   
9. 
   
10. int match_count;
    
11. int total_count;
    
12. void output(char map[NUM])
    
13. {
    
14.     int unused=NUM;
    
15.     int i;
    
16.     char map2[NUM];
    
17.     int unmatched_found=0;
    
18.     memset(map2,-1,sizeof(map2));
    
19.     for(i=0;i<NUM;i++){
    
20.         printf("%c", i+'A');
    
21.     }
    
22.     printf("  ==>  ");
    
23.     for(i=0;i<NUM;i++){
    
24.         if(map[i]!=-1){
    
25.             map2[map[i]]=NUM-1-i;
    
26.             printf("%c", map[i]+'A');
    
27.         }else{
    
28.             printf("%c", unused++ + 'A');
    
29.         }
    
30.     }
    
31.     printf("  R  ");
    
32.     unused=NUM;
    
33.     for(i=NUM-1;i>=0;i--){
    
34.         if(map2[i]!=map[NUM-1-i])unmatched_found=1;
    
35.         if(map2[i]!=-1){
    
36.            printf("%c", map2[i]+'A');
    
37.         }else{
    
38.            printf("%c", unused++ + 'A');
    
39.         }
    
40.     }
    
41.     printf("\n");
    
42.     if(!unmatched_found){
    
43.        match_count++;
    
44.     }
    
45.     total_count++;
    
46. }
    
47. 
    
48. void recursive(int old_state, int search_bit, char map[NUM])
    
49. {
    
50.     if((old_state&MASK_FOR_BEFORE_N_BITS(search_bit))==0){//finished search
    
51.         output(map);
    
52.         return;
    
53.     }
    
54.     while(!IS_BIT_SET(old_state, search_bit))search_bit++;
    
55.     int i;
    
56.     for(i=0;i<NUM;i++){
    
57.         if(map[i]==-1){
    
58.             map[i]=search_bit;
    
59.             recursive(old_state, search_bit+1, map);
    
60.             map[i]=-1;
    
61.         }
    
62.     }
    
63. }
    
64. 
    
65. void search(int old_state)/*bit state to show whether a old state is used*/
    
66. {
    
67.     char map[NUM];
    
68.     memset(map, -1, sizeof(map));
    
69.     recursive(old_state, 0, map);
    
70. }
    
71. 
    
72. int main()
    
73. {
    
74.     int i;
    
75.     for(i=0; i<(1<<NUM);i++){
    
76.         search(i);
    
77.     }
    
78.     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);
    
79.     return 0;
    
80. }
    

复制代码

output:
ABC  ==>  DEF  R  DEF
ABC  ==>  ADE  R  DEC
ABC  ==>  DAE  R  DEB
ABC  ==>  DEA  R  DEA
ABC  ==>  BDE  R  DCE
ABC  ==>  DBE  R  DBE
ABC  ==>  DEB  R  DAE
ABC  ==>  ABD  R  DBC
ABC  ==>  ADB  R  DAC
ABC  ==>  BAD  R  DCB
ABC  ==>  DAB  R  DAB
ABC  ==>  BDA  R  DCA
ABC  ==>  DBA  R  DBA
ABC  ==>  CDE  R  CDE
ABC  ==>  DCE  R  BDE
ABC  ==>  DEC  R  ADE
ABC  ==>  ACD  R  BDC
ABC  ==>  ADC  R  ADC
ABC  ==>  CAD  R  CDB
ABC  ==>  DAC  R  ADB
ABC  ==>  CDA  R  CDA
ABC  ==>  DCA  R  BDA
ABC  ==>  BCD  R  BCD
ABC  ==>  BDC  R  ACD
ABC  ==>  CBD  R  CBD
ABC  ==>  DBC  R  ABD
ABC  ==>  CDB  R  CAD
ABC  ==>  DCB  R  BAD
ABC  ==>  ABC  R  ABC
ABC  ==>  ACB  R  BAC
ABC  ==>  BAC  R  ACB
ABC  ==>  CAB  R  CAB
ABC  ==>  BCA  R  BCA
ABC  ==>  CBA  R  CBA

mathe

K=5的状态算出来了，共209，去除对称性还有126个状态，当然最终递推数列的阶应该会更低一些

mathe

稍微修改上面代码，可以将两个矩阵输出成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分配的内存不够)

1. #include <stdio.h>
   
2. #include <string.h>
   
3. #include <stdlib.h>
   
4. #include <map>
   
5. #include <string>
   
6. using namespace std;
   
7. #define K 5
   
8. #define NUM ((K)-1)
   
9. #define MASK_FOR_LAST_N_BITS(n)   ((1<<(n))-1)
   
10. #define MASK_FOR_BEFORE_N_BITS(n) (~MASK_FOR_LAST_N_BITS(n))
    
11. #define IS_BIT_SET(x, b)   (((x)&(1<<(b)))!=0)
    
12. 
    
13. typedef map<string, string> SMap;
    
14. typedef map<string, int> IMap;
    
15. SMap string_map;
    
16. IMap imap1, imap2;
    
17. int ic1,ic2;
    
18. void decode_string(const string& s, int& present, int& last_used)
    
19. {
    
20.     present = 0;
    
21.     last_used=0;
    
22.     int i;
    
23.     for(i=0;i<NUM;i++){
    
24.        int d=s[i]-'A';
    
25.        if(d>=NUM)last_used=d-NUM+1;
    
26.        else present |= 1<<d;
    
27.     }
    
28. }
    
29. 
    
30. string shift_string(const string& s, char& new_char)
    
31. {
    
32.     char c=s[0];
    
33.     int first_new=0;
    
34.     if(c-'A'>=NUM){
    
35.         first_new=1;
    
36.     }
    
37.     int i;
    
38.     string out;
    
39.     for(i=1;i<NUM;i++){
    
40.         c=s[i];
    
41.         if(c-'A'>=NUM){
    
42.             c-=first_new;
    
43.         }
    
44.         out+=c;
    
45.     }
    
46.     if(first_new)new_char--;
    
47.     return out;
    
48. }
    
49. 
    
50. void gen_sum(void)
    
51. {
    
52.     printf("T=10;\n");
    
53.     printf("m1=matrix(%d,%d);\n",ic1,ic1);
    
54.     printf("m2=matrix(%d,%d);\n",ic2,ic2);
    
55.     IMap::iterator it;
    
56.     for(it=imap1.begin(); it!=imap1.end();++it){
    
57.         int from = it->second;
    
58.         string s=it->first;
    
59.         int present, last_used;
    
60.         decode_string(s, present, last_used);
    
61.         int i;
    
62.         char new_char='A'+NUM+last_used;
    
63.         string t = shift_string(s, new_char);
    
64.         for(i=0;i<NUM;i++){
    
65.           if(!IS_BIT_SET(present, i)){
    
66.              string h=t;
    
67.              h+=(char)('A'+i);
    
68.              IMap::iterator it2 = imap1.find(h);
    
69.              int to = it2->second;
    
70.              printf("m1[%d,%d]=1;\n", to+1, from+1);
    
71.           }
    
72.         }
    
73.         t+=new_char;
    
74.         IMap::iterator it2=imap1.find(t);
    
75.         int to = it2->second;
    
76.         printf("m1[%d,%d]=T-%d;\n", to+1, from+1, last_used+NUM);
    
77.     }
    
78.     for(it=imap1.begin(); it!=imap1.end();++it){
    
79.         string s = it->first;
    
80.         IMap::iterator it2 = imap2.find(s);
    
81.         if(it2==imap2.end()){
    
82.             SMap::iterator sit = string_map.find(s);
    
83.             string ns = sit->second;
    
84.             it2 = imap2.find(ns);
    
85.         }
    
86.         int from = it2->second;
    
87.         int present, last_used;
    
88.         decode_string(s, present, last_used);
    
89.         int i;
    
90.         char new_char='A'+NUM+last_used;
    
91.         string t = shift_string(s, new_char);
    
92.         for(i=0;i<NUM;i++){
    
93.            if(!IS_BIT_SET(present, i)){
    
94.               string h=t;
    
95.               h+=(char)('A'+i);
    
96.               IMap::iterator it2=imap2.find(h);
    
97.               if(it2 != imap2.end()){
    
98.                   int to=it2->second;
    
99.                   printf("m2[%d,%d]+=1;\n", to+1, from+1);
    
100.               }
     
101.            }
     
102.         }
     
103.         t+=new_char;
     
104.         it2=imap2.find(t);
     
105.         if(it2!=imap2.end()){
     
106.             int to = it2->second;
     
107.             printf("m2[%d,%d]+=T-%d;\n", to+1, from+1, last_used+NUM);
     
108.         }
     
109.     }
     
110. }
     
111. 
     
112. void output(char map[NUM])
     
113. {
     
114.     int unused=NUM;
     
115.     int i;
     
116.     char map2[NUM];
     
117.     int unmatched_found=0;
     
118.     memset(map2,-1,sizeof(map2));
     
119.     string s1,s2;
     
120.     for(i=0;i<NUM;i++){
     
121.         if(map[i]!=-1){
     
122.             map2[map[i]]=NUM-1-i;
     
123.             s1+=(char)(map[i]+'A');
     
124.         }else{
     
125.             s1+=(char)(unused++ + 'A');
     
126.         }
     
127.     }
     
128.     unused=NUM;
     
129.     for(i=NUM-1;i>=0;i--){
     
130.         if(map2[i]!=map[NUM-1-i])unmatched_found=1;
     
131.         if(map2[i]!=-1){
     
132.            s2+= (char)(map2[i]+'A');
     
133.         }else{
     
134.            s2+=(char)(unused++ + 'A');
     
135.         }
     
136.     }
     
137.     string_map.insert(make_pair(s1,s2));
     
138.     imap1.insert(make_pair(s1, ic1)); ic1++;
     
139.     if(s1<=s2){
     
140.        imap2.insert(make_pair(s1,ic2)); ic2++;
     
141.     }
     
142. }
     
143. 
     
144. void recursive(int old_state, int search_bit, char map[NUM])
     
145. {
     
146.     if((old_state&MASK_FOR_BEFORE_N_BITS(search_bit))==0){//finished search
     
147.         output(map);
     
148.         return;
     
149.     }
     
150.     while(!IS_BIT_SET(old_state, search_bit))search_bit++;
     
151.     int i;
     
152.     for(i=0;i<NUM;i++){
     
153.         if(map[i]==-1){
     
154.             map[i]=search_bit;
     
155.             recursive(old_state, search_bit+1, map);
     
156.             map[i]=-1;
     
157.         }
     
158.     }
     
159. }
     
160. 
     
161. void search(int old_state)/*bit state to show whether a old state is used*/
     
162. {
     
163.     char map[NUM];
     
164.     memset(map, -1, sizeof(map));
     
165.     recursive(old_state, 0, map);
     
166. }
     
167. 
     
168. int main()
     
169. {
     
170.     int i;
     
171.     for(i=0; i<(1<<NUM);i++){
     
172.         search(i);
     
173.     }
     
174.     gen_sum();
     
175.     return 0;
     
176. }

复制代码

mathe

现在我们以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)]$