P-recursive + LLL
例題
解法
P-recursive で解説した方法で漸化式を推定しようとすると、係数を一意に絞り込むのに 70 項ほど必要になります。一方素朴に思いつける bitDP だと時間的に 20 項ちょっとしか求めることができません。そこで、正しい漸化式の係数の絶対値は小さいはず！という大胆予想を立てて LLL アルゴリズムでノルムの小さい係数を復元します。実装は全て Mathematica 上で行いました。 code:prog.nb
seq = {1, 2, 2, 8, 28, 152, 952, 7208, 62296, 605864, 6522952, 76951496, 986411272, 13647501224, 202653304456, 3214170659528, 54223248489064, 969422707369448, 18308139748071496, 364185707232235784, 7610646121718276392, 166694564791285014056, 3818542633160305746952} (* bitDP で集めた項 *);
terms = 14 (* 何項間漸化式の存在を仮定するか *);
degree = 3 (* 係数多項式の次数を何次未満と仮定するか *);
mat = Table[Flatten@Table[(nn - i)^j seqnn - i, {i, 0, terms - 1}, {j, 0, degree - 1}], {nn, terms, Length@seq}];
{Dimensionsmat, MatrixRankmat} (* 式の数, 変数の数, 独立な式の数 の確認 *) dim = Dimensionsmat2 - MatrixRankmat; lat = SortBy[LatticeReduce[(SmithDecompositionmat3\Transpose)-dim ;; -1], N@Norm@# &]; ({N@Norm@#, #} & /@ lat1 ;; Min3, Length@lat) // Column (* 他の係数列の候補も見ておく *) coef = lat1 (* 2乗ノルム小さめの係数列を得る *);
eq = Sum[coef1 + i degree + j (nn - i)^j dpsubnn - i, {i, 0, terms - 1}, {j, 0, degree - 1}]; sol = Solve[eq == 0, dpsubnn]; CForm@FullSimplify@sol1, 1, 2 (* そのまま貼り付けられる形式で出力 *)
code:output.txt
{{10, 42}, 10}
{289.344,{-4,0,0,22,4,0,-38,-16,0,-3,-1,-2,56,79,9,-32,-74,-8,28,-90,-14,20,116,18,-136,64,20,8,-98,-22,66,-6,-18,39,45,16,12,17,3,-6,-8,-2}}
{417.539,{4,0,0,-22,-4,0,34,16,0,11,5,2,-69,-80,-9,68,69,5,-111,53,14,-3,-110,-5,-6,128,23,7,20,-57,12,-20,-56,-36,66,16,168,15,41,102,-174,-152}}
{420.058,{-2,0,0,1,2,0,8,2,0,-16,-12,-1,-8,2,-1,-99,17,0,123,3,0,20,7,63,-89,-31,34,-23,48,-155,104,-87,-221,-60,4,60,-30,-5,-17,-138,112,-46}}
(-2*(-12 + nn)*(-10 + nn)*dpsub(-13 + nn) + (240 + nn*(-55 + 3*nn))*dpsub(-12 + nn) + nn*((-307 + 16*nn)*dpsub(-11 + nn) + 6*(59 - 3*nn)*dpsub(-10 + nn) + 298*dpsub(-9 + nn) - 256*dpsub(-8 + nn) - 136*dpsub(-7 + nn) + 78*dpsub(-6 + nn) + 6*dpsub(-5 + nn) + 7*dpsub(-4 + nn) + nn*(-22*dpsub(-9 + nn) + 20*dpsub(-8 + nn) + 18*dpsub(-7 + nn) - 14*dpsub(-6 + nn) - 8*dpsub(-5 + nn) + 9*dpsub(-4 + nn) - 2*dpsub(-3 + nn)) + 11*dpsub(-3 + nn) - 16*dpsub(-2 + nn) + 4*dpsub(-1 + nn)) + 2*(740*dpsub(-11 + nn) - 837*dpsub(-10 + nn) - 446*dpsub(-9 + nn) + 316*dpsub(-8 + nn) + 45*dpsub(-7 + nn) + 32*dpsub(-6 + nn) + 69*dpsub(-5 + nn) - 58*dpsub(-4 + nn) - 9*dpsub(-3 + nn) - 3*dpsub(-2 + nn) + 9*dpsub(-1 + nn)))/4.
23 項しかないので、変数が 42 個もあるのに式は 10 本しかありません。
普通なら諦めるところですが、LLL を信じて係数を求めるとノルムが 289.344 しかないベクトルを得ることができました。
次善の 417.539 との差も大きいことからアタリを引けている感じがします。
試しに出力をそのまま C++ のコードに貼り付けて提出したところ、無事 AC することができました。
提出
他の使用例
DFS による全探索で項を集めますが、係数を一意に定めるには足りないので LLL を利用します。
BFS による全探索で項を集めますが、係数を一意に定めるには足りないので LLL を利用します。
P-recursive ではなく D-algebraic ですが同様の手法で母関数の満たす微分方程式が得られます。
後は Mathematica の RSolve で微分方程式を解いて母関数を求め、実装だけ FPS ライブラリで行えばいいです。
元ポスト
など