0419 外观数列
* * * *
拉格朗日计划
* * * *
外观数列

外观数列的首项是1,之后每一项$a_n$都籍由观察$a_{n-1}$而得。

将$a_{n-1}$分割成若干段,每一段中各数位上的数字都相同,而相邻两段数位上的数字不同,而$a_n$就由各段的长度和该段数位上的数字拼接而成。

它的前几项分别是:1, 11, 21, 1211, 111221, 312211, 13112221, 1113213211, ...

其中1是1个1,因此第二项是 11。

11是2个1, 因此第三项是21。

21是1个2和1个1,因此第四项是1211。

1211是1个1、1个2和2个1,因此第五项是111221。

111221是3个1、2个2、1个1,因此第六项是312211。

以此类推。

分别用$A(n), B(n), C(n)$表示第n项元素中数字1、2、3的数量,可以验证$A(40)=31254, B(40)=20259, C(40)=11625$。

$A(10^12), B(10^12), C(10^12)$模$2^{30}$各自的余数,并将这三个余数用逗号隔开作为本题的答案。

(例如对$n=40$,答案的形式应为31254,20259,11625。)

本题难度:



解答

根据外观数列的所谓Cosmological定理(可以参考比如此页面此页面),从第8项开始该数列中的每项都总是可以分解成92个生成元的相互拼接,生成元之间的相互转化已在引文的页面中给出。

先初始化每种状态中数字1、2、3的数量和生成元的转化表,每一次转换后数字数量的变化可以用一个$92\times 92$的矩阵表示,接下来用一个函数计算给定数量变化矩阵V的10次方,递归调用12次该函数即得结果$998567458,1046245404,43363922$。

注:为便于格式化输出,以下代码为Python 3。

f=lambda s:[s.count("1"),s.count("2"),s.count("3")]

c=[]
c.append([0,0,0]) #0 Unused
c.append(f("22")) #1 H
c.append(f("13112221133211322112211213322112")) #2 He
c.append(f("312211322212221121123222112")) #3 Li
c.append(f("111312211312113221133211322112211213322112")) #4 Be
c.append(f("1321132122211322212221121123222112")) #5 B
c.append(f("3113112211322112211213322112")) #6 C
c.append(f("111312212221121123222112")) #7 N
c.append(f("132112211213322112")) #8 O
c.append(f("31121123222112")) #9 F
c.append(f("111213322112")) #10 Ne
c.append(f("123222112")) #11 Na
c.append(f("3113322112")) #12 Mg
c.append(f("1113222112")) #13 Al
c.append(f("1322112")) #14 Si
c.append(f("311311222112")) #15 P
c.append(f("1113122112")) #16 S
c.append(f("132112")) #17 Cl
c.append(f("3112")) #18 Ar
c.append(f("1112")) #19 K
c.append(f("12")) #20 Ca
c.append(f("3113112221133112")) #21 Sc
c.append(f("11131221131112")) #22 Ti
c.append(f("13211312")) #23 V
c.append(f("31132")) #24 Cr
c.append(f("111311222112")) #25 Mn
c.append(f("13122112")) #26 Fe
c.append(f("32112")) #27 Co
c.append(f("11133112")) #28 Ni
c.append(f("131112")) #29 Cu
c.append(f("312")) #30 Zn
c.append(f("13221133122211332")) #31 Ga
c.append(f("31131122211311122113222")) #32 Ge
c.append(f("11131221131211322113322112")) #33 As
c.append(f("13211321222113222112")) #34 Se
c.append(f("3113112211322112")) #35 Br
c.append(f("11131221222112")) #36 Kr
c.append(f("1321122112")) #37 Rb
c.append(f("3112112")) #38 Sr
c.append(f("1112133")) #39 Y
c.append(f("12322211331222113112211")) #40 Zr
c.append(f("1113122113322113111221131221")) #41 Nb
c.append(f("13211322211312113211")) #42 Mo
c.append(f("311322113212221")) #43 Tc
c.append(f("132211331222113112211")) #44 Ru
c.append(f("311311222113111221131221")) #45 Rh
c.append(f("111312211312113211")) #46 Pd
c.append(f("132113212221")) #47 Ag
c.append(f("3113112211")) #48 Cd
c.append(f("11131221")) #49 In
c.append(f("13211")) #50 Sn
c.append(f("3112221")) #51 Sb
c.append(f("1322113312211")) #52 Te
c.append(f("311311222113111221")) #53 I
c.append(f("11131221131211")) #54 Xe
c.append(f("13211321")) #55 Cs
c.append(f("311311")) #56 Ba
c.append(f("11131")) #57 La
c.append(f("1321133112")) #58 Ce
c.append(f("31131112")) #59 Pr
c.append(f("111312")) #60 Nd
c.append(f("132")) #61 Pm
c.append(f("311332")) #62 Sm
c.append(f("1113222")) #63 Eu
c.append(f("13221133112")) #64 Gd
c.append(f("3113112221131112")) #65 Tb
c.append(f("111312211312")) #66 Dy
c.append(f("1321132")) #67 Ho
c.append(f("311311222")) #68 Er
c.append(f("11131221133112")) #69 Tm
c.append(f("1321131112")) #70 Yb
c.append(f("311312")) #71 Lu
c.append(f("11132")) #72 Hf
c.append(f("13112221133211322112211213322113")) #73 Ta
c.append(f("312211322212221121123222113")) #74 W
c.append(f("111312211312113221133211322112211213322113")) #75 Re
c.append(f("1321132122211322212221121123222113")) #76 Os
c.append(f("3113112211322112211213322113")) #77 Ir
c.append(f("111312212221121123222113")) #78 Pt
c.append(f("132112211213322113")) #79 Au
c.append(f("31121123222113")) #80 Hg
c.append(f("111213322113")) #81 Tl
c.append(f("123222113")) #82 Pb
c.append(f("3113322113")) #83 Bi
c.append(f("1113222113")) #84 Po
c.append(f("1322113")) #85 At
c.append(f("311311222113")) #86 Rn
c.append(f("1113122113")) #87 Fr
c.append(f("132113")) #88 Ra
c.append(f("3113")) #89 Ac
c.append(f("1113")) #90 Th
c.append(f("13")) #91 Pa
c.append(f("3")) #92 U

def g(v):
  z=[0]*93
  for i in v:
    z[i]+=1
  return z

t=[]
t.append([0]*93) #0 Unused
t.append(g([1])) #1 H
t.append(g([1,72,91,20,3])) #2 He
t.append(g([2])) #3 Li
t.append(g([3,32,20])) #4 Be
t.append(g([4])) #5 B
t.append(g([5])) #6 C
t.append(g([6])) #7 N
t.append(g([7])) #8 O
t.append(g([8])) #9 F
t.append(g([9])) #10 Ne
t.append(g([10])) #11 Na
t.append(g([11,61])) #12 Mg
t.append(g([12])) #13 Al
t.append(g([13])) #14 Si
t.append(g([14,67])) #15 P
t.append(g([15])) #16 S
t.append(g([16])) #17 Cl
t.append(g([17])) #18 Ar
t.append(g([18])) #19 K
t.append(g([19])) #20 Ca
t.append(g([20,67,91,1,27])) #21 Sc
t.append(g([21])) #22 Ti
t.append(g([22])) #23 V
t.append(g([23])) #24 Cr
t.append(g([24,14])) #25 Mn
t.append(g([25])) #26 Fe
t.append(g([26])) #27 Co
t.append(g([27,30])) #28 Ni
t.append(g([28])) #29 Cu
t.append(g([29])) #30 Zn
t.append(g([30,63,20,89,1,20])) #31 Ga
t.append(g([31,67])) #32 Ge
t.append(g([32,11])) #33 As
t.append(g([33])) #34 Se
t.append(g([34])) #35 Br
t.append(g([35])) #36 Kr
t.append(g([36])) #37 Rb
t.append(g([37])) #38 Sr
t.append(g([38,92])) #39 Y
t.append(g([39,1,20,43])) #40 Zr
t.append(g([40,68])) #41 Nb
t.append(g([41])) #42 Mo
t.append(g([42])) #43 Tc
t.append(g([43,63,20])) #44 Ru
t.append(g([44,67])) #45 Rh
t.append(g([45])) #46 Pd
t.append(g([46])) #47 Ag
t.append(g([47])) #48 Cd
t.append(g([48])) #49 In
t.append(g([49])) #50 Sn
t.append(g([50,61])) #51 Sb
t.append(g([51,63,20])) #52 Te
t.append(g([52,67])) #53 I
t.append(g([53])) #54 Xe
t.append(g([54])) #55 Cs
t.append(g([55])) #56 Ba
t.append(g([56])) #57 La
t.append(g([57,1,20,27])) #58 Ce
t.append(g([58])) #59 Pr
t.append(g([59])) #60 Nd
t.append(g([60])) #61 Pm
t.append(g([61,20,30])) #62 Sm
t.append(g([62])) #63 Eu
t.append(g([63,20,27])) #64 Gd
t.append(g([64,67])) #65 Tb
t.append(g([65])) #66 Dy
t.append(g([66])) #67 Ho
t.append(g([67,61])) #68 Er
t.append(g([68,20,27])) #69 Tm
t.append(g([69])) #70 Yb
t.append(g([70])) #71 Lu
t.append(g([71])) #72 Hf
t.append(g([72,91,1,20,74])) #73 Ta
t.append(g([73])) #74 W
t.append(g([74,32,20])) #75 Re
t.append(g([75])) #76 Os
t.append(g([76])) #77 Ir
t.append(g([77])) #78 Pt
t.append(g([78])) #79 Au
t.append(g([79])) #80 Hg
t.append(g([80])) #81 Tl
t.append(g([81])) #82 Pb
t.append(g([82,61])) #83 Bi
t.append(g([83])) #84 Po
t.append(g([84])) #85 At
t.append(g([85,67])) #86 Rn
t.append(g([86])) #87 Fr
t.append(g([87])) #88 Ra
t.append(g([88])) #89 Ac
t.append(g([89])) #90 Th
t.append(g([90])) #91 Pa
t.append(g([91])) #92 U

mod=2**30

def step(v,t):
  res=[0]*93
  for j in range(1,93):
    for k in range(1,93):
      res[k]=(res[k]+v[j]*t[j][k])%mod
  return res

tList=[]
v=t
for i in range(11):
  x=v[:]
  for i in range(9):
    y=[0]
    for j in range(1,93):
      z=[0]*93
      for k in range(1,93):
        for l in range(1,93):
          z[l]=(z[l]+x[j][k]*v[k][l])%mod
      y.append(z)
    x=y
  tList.append(x)
  v=x

#starting state - step 8 = 72 & 50
x=step(step(g([72,50]),t),t)
for i in tList:
  for j in range(9):
    x=step(x,i)

print(f"{sum(c[j][0]*x[j] for j in range(1,93))%mod},{sum(c[j][1]*x[j] for j in range(1,93))%mod%mod},{sum(c[j][2]*x[j] for j in range(1,93))%mod%mod}")