拡張ユークリッド互除法を使って特性多項式を求める

定数係数線形漸化式を持つような数列の先頭 n 項が与えられるので漸化式を求めよ、という問題を解く話。特性多項式を求める方針で攻められる。

線形漸化式を持つタイプの数列は母関数が有理式になるので、mod x^2m 上で有理式近似を掛けるときれいに求まるという感じらしい。有理式近似は Pade approximation を参考にするといいかも。

母関数と x^2m に対して拡張ユークリッド互除法を走らせると、その過程でかならず現れるらしいが、本当にそうなのか証明はしていない。まあそんな気はする。

import sympy

x = sympy.Symbol('x')
mod = 10**9+7

def extended_euclidean_algorithm(a, b):
	init = b
	x0 = sympy.Poly(0, x, modulus=mod)
	x1 = sympy.Poly(1, x, modulus=mod)

	while b != 0:
		c = sympy.Poly(b.coeffs()[-1], x, modulus=mod)
		print('P=', sympy.div(x1, c)[0])
		print('Q=', sympy.div(b, c)[0])
		print('DOF=', x1.degree() + b.degree() + 1)
		print()

		q, r = sympy.div(a, b, x)
		a, b = b, r
		x0, x1 = x1, x0 - q * x1

a = [1, 1, 1, 2, 100, 1]
c = [4, 3, 2, 1, 3, 10]
for i in range(13):
	a.append(sum(x * y for x, y in zip(a[-len(c)-1:-1], c)))
a.reverse()

print(len(a))

p = sympy.Poly(x**len(a), modulus=mod)
q = sympy.Poly(a, x, modulus=mod)

extended_euclidean_algorithm(p, q)
19
P= Poly(1, x, modulus=1000000007)
Q= Poly(433300268*x**18 + 91176767*x**17 + 40346412*x**16 + 7817089*x**15 + 3778212*x**14 + 660061*x**13 + 356105*x**12 + 54558*x**11 + 33824*x**10 + 4354*x**9 + 3224*x**8 + 321*x**7 + 311*x**6 + x**5 + 100*
x**4 + 2*x**3 + x**2 + x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(-458720627*x + 1, x, modulus=1000000007)
Q= Poly(-189109440*x**17 + 98052629*x**16 + 242370173*x**15 + 410119453*x**14 + 292925697*x**13 + 120563428*x**12 + 233675522*x**11 - 269562155*x**10 + 84713259*x**9 - 249317014*x**8 + 337886325*x**7 - 45872
0316*x**6 + 127937623*x**5 + 82558853*x**4 - 458720625*x**3 - 458720626*x**2 - 458720626*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(110737070*x**2 - 429080800*x + 1, x, modulus=1000000007)
Q= Poly(385824477*x**16 + 12285214*x**15 - 472131743*x**14 - 223531587*x**13 - 219136967*x**12 - 79623645*x**11 - 201445119*x**10 + 190114060*x**9 - 295704085*x**8 - 333390478*x**7 - 355373566*x**6 + 3133944
42*x**5 + 252575577*x**4 - 318343728*x**3 - 318343729*x**2 - 429080799*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(331854849*x**3 + 256761159*x**2 + 313273835*x + 1, x, modulus=1000000007)
Q= Poly(158575077*x**15 + 111745351*x**14 + 412440497*x**13 + 179510050*x**12 + 12279572*x**11 + 317678473*x**10 - 377969900*x**9 - 254521710*x**8 - 129591852*x**7 - 346900445*x**6 + 172760444*x**5 + 2151637
71*x**4 - 98110162*x**3 - 429965012*x**2 + 313273836*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(-158195669*x**4 + 162488741*x**3 + 440244834*x**2 + 399082250*x + 1, x, modulus=1000000007)
Q= Poly(266864312*x**14 + 296985854*x**13 + 423900252*x**12 + 466865813*x**11 - 86494178*x**10 + 335562504*x**9 + 364468946*x**8 + 487306687*x**7 - 409652541*x**6 - 206992546*x**5 + 242702499*x**4 + 1815820*
x**3 - 160672922*x**2 + 399082251*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(353121972*x**5 - 153947681*x**4 + 488417650*x**3 + 216965224*x**2 - 190009502*x + 1, x, modulus=1000000007)
Q= Poly(271541463*x**13 - 311030879*x**12 - 152497574*x**11 + 452003256*x**10 + 111347971*x**9 + 283031296*x**8 + 11002103*x**7 - 317477361*x**6 + 120572316*x**5 + 171416289*x**4 - 484626633*x**3 + 26955723*
x**2 - 190009501*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(265498670*x**6 + 133561350*x**5 - 221331431*x**4 + 200136041*x**3 - 261171489*x**2 + 46417327*x + 1, x, modulus=1000000007)
Q= Poly(-164450549*x**12 - 126662090*x**11 + 177221842*x**10 - 278472701*x**9 + 275247039*x**8 + 144618549*x**7 - 492730416*x**6 + 231755655*x**5 - 189532125*x**4 - 14618119*x**3 - 214754161*x**2 + 46417328*
x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(-4*x**7 - 3*x**6 - 2*x**5 - x**4 - 3*x**3 - 10*x**2 + 1, x, modulus=1000000007)
Q= Poly(-701*x**6 - 25*x**5 + 86*x**4 - 11*x**3 - 9*x**2 + x + 1, x, modulus=1000000007)
DOF= 14

P= Poly(-432334901*x**13 + 451937767*x**12 - 84849138*x**11 + 414755438*x**10 + 164413392*x**9 + 479304138*x**8 + 188841490*x**7 + 465095370*x**6 - 88085838*x**5 + 15015864*x**4 - 357520392*x**3 + 89112240*x
**2 + 411791140*x + 1, x, modulus=1000000007)
Q= Poly(-73252172*x**5 - 429809915*x**4 + 143382990*x**3 - 499096626*x**2 + 411791141*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(-485549536*x**14 - 307852562*x**13 + 286646299*x**12 - 454066241*x**11 + 28187677*x**10 - 210588267*x**9 + 267204505*x**8 - 430863595*x**7 - 275825196*x**6 + 24500831*x**5 - 83206444*x**4 + 262317670
*x**3 + 454010118*x**2 + 408883680*x + 1, x, modulus=1000000007)
Q= Poly(450888797*x**4 + 125211463*x**3 - 137106208*x**2 + 408883681*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(181445891*x**15 + 376104137*x**14 + 28983539*x**13 + 430277410*x**12 + 306655481*x**11 - 244200915*x**10 + 363824597*x**9 + 370987859*x**8 - 37946091*x**7 - 195739975*x**6 + 337143087*x**5 + 11009798
1*x**4 - 168587021*x**3 - 99935396*x**2 + 79212168*x + 1, x, modulus=1000000007)
Q= Poly(-189310247*x**3 - 20723227*x**2 + 79212169*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(-159146717*x**16 - 415262651*x**15 - 475589347*x**14 + 17046883*x**13 - 359952422*x**12 + 391261832*x**11 + 1114016*x**10 - 81974509*x**9 + 2314290*x**8 + 101306447*x**7 - 360072227*x**6 - 78717762*x
**5 - 104442693*x**4 + 52213782*x**3 - 156656379*x**2 + 104442595*x + 1, x, modulus=1000000007)
Q= Poly(-52213783*x**2 + 104442596*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(79300793*x**17 - 111374585*x**16 + 301307826*x**15 - 391072155*x**14 + 407405482*x**13 + 135449279*x**12 + 433963753*x**11 - 455024730*x**10 - 186233746*x**9 - 257421677*x**8 - 72664918*x**7 - 276188
924*x**6 - 402079573*x**5 + 26782586*x**4 - x**3 + 26782683*x**2 - 26782684*x + 1, x, modulus=1000000007)
Q= Poly(-26782683*x + 1, x, modulus=1000000007)
DOF= 19

P= Poly(14595107*x**18 - 221455372*x**17 - 4362905*x**16 + 23247227*x**15 - 9359715*x**14 + 2569025*x**13 - 417225*x**12 - 143680*x**11 + 60732*x**10 - 26418*x**9 + 6317*x**8 + 495*x**7 - 408*x**6 + 197*x**5
 - 97*x**4 - x**3 - x + 1, x, modulus=1000000007)
Q= Poly(1, x, modulus=1000000007)
DOF= 19

Note

  • DOF=degree of freedom
  • P=denominator
  • Q=numerator
  • 一個だけ自由度が 14 になってるけど、これは mod F[x]/<x^19> で考えても F[x] で考えても結果が同じ(精度が足りてる)ことを意味していて(多分)、これが欲しい有理式になっている。
  • 多項式演算を書くのが面倒だったので sympy を使ったけど、そんなにテクニカルな機能を使っているわけではないので C++ で書き直すのは難しくはないと思う。愚直な除算を使っても O(n^2) になるし、そもそも剰余は要らなくて減算だけ実装するだけで良い。