fortran66のブログ

fortran について書きます。

【メモ帳】制限付き分割 p(n,5)

Andrews & Eriksson 「整数の分割」演習 88

だいぶ昔に読んでたアンドリューズとエリクソンの「整数の分割」の演習問題 88. が面倒くさかったので飛ばしていたのですが、夏休みなので暇つぶしに Maple/Mathematica の手を借りてやってみました。すぐ忘れるのでメモしておきます。

整数の分割

整数の分割

きんもーっ☆エリクソン
Integer Partitions

Integer Partitions

演習 88.

p(n, 5)=\{{1\over2880}(n+8)(n^3+22n^2+44n+248+180\lfloor {n\over2}\rfloor  \}
であることを証明せよ。

(ここで p(n,5) は整数nを1から5までの数に分割するときの組み合わせの個数、中括弧は四捨五入、\lfloor x\rfloor は床関数になってます。)

例 p(n, 4)、 4 までの制限つき分割
-n=1 p(1, 4)=1 
O
-n=2 p(2 ,4)=2
O, OO
O
-n=3 p(3, 4)=3
O, OO, OOO
O, O
O
-n=4 p(4, 4)=5
O, OO, OO, OOO, OOOO
O, O,  OO, O
O, O
O
-n=5 p(5, 4)=6
O, OO, OO, OOO, OOO, OOOO
O, O,  OO, O,   OO,  O
O, O,  O,  O
O, O
O   
導出


\begin{align}
p(n, 5)&=\prod_{k=1}^5{1\over(1-x^k)}\\
&={1\over120}{1\over(1-x)^5}+{1\over24}{1\over(1-x)^4}+{31\over288}{1\over(1-x)^3}+{3\over16}{1\over(1-x)^2}+{1\over64}{1\over(1+x)^2}\\
&+{20831\over86400}{1\over1-x}+{13\over128}{1\over1+x}+{1\over27}{x+2\over1+x+x^2}+{1\over25}{4+3x+2x^2+x^3\over1+x+x^2+x^3+x^4}\\
&={1\over120}{1\over(1-x)^5}+{1\over24}{1\over(1-x)^4}+{31\over288}{1\over(1-x)^3}+{3\over16}{1\over(1-x)^2}+{1\over64}{1\over(1+x)^2}\\
&+{1\over16}{1\over1-x}+{13\over128}({1\over1-x}+{1\over1+x})+{1\over27}({1\over1-x}+{x+2\over1+x+x^2})+
{1\over5}({1\over1-x}+{4+3x+2x^2+x^3\over1+x+x^2+x^3+x^4})\\
&={1\over120}{1\over(1-x)^5}+{1\over24}{1\over(1-x)^4}+{31\over288}{1\over(1-x)^3}+{3\over16}{1\over(1-x)^2}+{1\over64}{1\over(1+x)^2}\\
&+{1\over16}{1\over1-x}+{13\over64}{1\over1-x^2}+{1\over9}{1\over1-x^3}+{1\over5}{1\over1-x^5}\\
\end{align}

部分分数に展開の上で、いくつかの項をまとめた。(部分分数展開 Maple convert(f, parfrac), Mathematica Apart[f]) 

級数に展開すると
 
\begin{align}
p(n,5)&=\sum_n\left({1\over120}\binom{n+4}{4}+{1\over24}\binom{n+3}{3}+{31\over288}\binom{n+2}{2}+{3\over16}\binom{n+1}{1}+{(-1)^n\over64}\binom{n+1}{1}+{1\over16}    \right)x^n \\
&+\sum_n\left({13\over64}x^{2n}+{1\over9}x^{3n}+{1\over5}x^{5n}\right)\\
\end{align}
ここで、二段目の項は区間 [0, 14/45] に収まるので、とりあえず無視しておく。

1段目の項をさらに変形すると、

\begin{align}
&={1\over2880}
\left( n^4+30n^3+310n^2+1275n+1124 \right)+{45\over2880}(-1)^n\left((n+1)+7\right)\\
&-{45\over2880}(-1)^n7\\
&={1\over2880}
\left( n^4+30n^3+310n^2+1230n+764+45( (-1)^n (n+8)+(n+8) )\right)\\
&-{45\over2880}(-1)^n7\\
\end{align}
ここで、因数分解の係数を合わせるために(答えから逆算して)余分な項を付け足し、二段目でそれをキャンセルさせている。この項は先ほど無視することにした項と合わせて値は区間 [-7/64, 99/320] になるが、四捨五入に影響しないので無視できる。また前の方の項から値を借りてきて後ろの項を整えている。

ここで 、
(-1)^n+1=2\left(2\lfloor{n\over2}\rfloor-n+1\right)
であるので、これを用いて一段目の式を変形すると、

\begin{align}
&={1\over2880}\left(n^4+30n^3+310n^2+1230n+764+90(n+8)(2\lfloor{n\over2}\rfloor-n+1)\right)\\
&={1\over2880}\left(n^4+30n^3+220n^2+600n+1484+(n+8)180\lfloor{n\over2}\rfloor\right)\\
&={1\over2880}\left((n+8)(n^3+22n^2+44n+248+180\lfloor{n\over2}\rfloor)-500\right)\\
\end{align}
ここで、-{500\over2880}=-{25\over144} の項は因数分解の都合合わせで出る余りで、今までの半端な項と合わせてやって区間 [-163/576,391/2880] となって絶対値は 1/2 に満たないので四捨五入に影響しない。よって無視できる。

以上のことから

p(n,5)=\left\{  {1\over2880}(n+8)\left( n^3+22n^2+44n+248+180\lfloor{n\over2}\rfloor\right)  \right\}
であることが示せた。

無視できるカスの項は適当な計算なので間違ってたらゴメンw

その他

演習87. p(n, 4)=\left\{ {1\over144}(n+5)(n^2+n+22+18\lfloor{n\over2}\rfloor) \right\}
は、上記と同じ方法で求まる。

演習86. 
p(n, 4)=\left\{ {1\over36}\lfloor{n+4\over2}\rfloor^2(3\lfloor{n+9\over2}\rfloor-\lfloor{n+10\over2}\rfloor)\right\} は等式変形としては求められるが、何か求め方があってこの形になっているのではないかと思われるがよく分からんw

いじっているうちに、副産物として 
p(n, 4)=\left\{ {1\over36}\lfloor{n+7\over2}\rfloor^2(3\lfloor{n+2\over2}\rfloor-\lfloor{n+1\over2}\rfloor) 
\right\} が出てきた。こっちの方があぶれるカス項が小さい。

Fortran にて

4までの分割は、3までの分割の和で表わせるので、それをプログラムする。

4段のヤング図の4段目を0から1個づつ増やして行くと、残りの n-4i 個を3段のヤング図にするだけの可能性がある。その総和を求めれば p(n, 4)。なお p(n,3)=round((n+3)^2/12)。これは上記の論法を3段と2段でやれば求まる。2段は1段とやればよい。

この方法で演習86. が出来ないかと思ったが、うまくいかなかったw

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

ソース・プログラム

    module m_mod
        implicit none
    contains
        integer function p4(n) 
            integer, intent(in) :: n
            integer :: i
            p4 = 0
            do i = 0, n / 4
                p4 = p4 + p3(n - 4 * i) 
            end do
        end function p4
        
        integer function p3(n)
            integer, intent(in) :: n
            p3 = nint((n + 3)**2 / 12.0)        
        end function p3
    end module m_mod
        
    program partition4
      use m_mod
      implicit none
      integer :: i
      do i = 1, 10
          print *, i, p4(i)
      end do 
    end program partition4

実行結果

           1           1
           2           2
           3           3
           4           5
           5           6
           6           9
           7          11
           8          15
           9          18
          10          23
続行するには何かキーを押してください . . .