fortran66のブログ

fortran について書きます。

整数の分割 p(n) を求める4つの方法

メモ帳

4つの漸化式による方法。

オイラーによるもの? ヤング図長 k の分割

p(0, k) = 0
p(k, 0) = 0
p(n, k) = 0 n < k
p(n, n) = 1

p(n, k) = p(n - 1, k - 1) + p(n - k, k)

p(n) = sum(p(n, 1:n))

2 最小マス数が k 以上

p(0, k) = 0
p(k, 0) = 0
p(n, k) = 0 n < k
p(n, n) = 1

p(n, k) = p(n - k, k) + p(n, k + 1)

p(n) = p(n, 1)

p(n, 2) が攪乱順列に対応します。対称群の指標の位数に関係ある漸化式です。

3 約数和関数との畳み込み漸化式

p(n) = 1/n sum(p(k) * sigma(n - k)) k=0..n-1
sigma は約数の和
但し p(0) = 1, sigma(0) = 1 と置く。

4 オイラーの五角数による漸化式

p(n) = sum( (-1)^k * p( k(3k-1) / 2 ) + (-1)^k * p( k(3k+1) / 2 )
但し p(0) = 1 とする。
なお同じ漸化式で p(0) = n と置くと、約数和 sigma(n) が得られます。

戯れに p(0) = INT(sqrt(n)) と置いてみたところ第25項までは、インターネット整数列辞典の A237831 と一致する数列が得られました。この数列は、妙な分割に対応しているようです。その後、少しづつずれてゆきます。

メモ

前二者はヤング図から証明でき、後二者は母関数を求めると証明できます。

ソース・プログラム

素朴に前2者再帰的にプログラムしたので、時間がかかります。後二者も素朴にプログラムしましたが、ただの畳み込み式なので素早く終わります。

      module m_partition
        implicit none
        integer, parameter :: ki = 8
      contains
        recursive integer(ki) pure function p1(n, k) result(ires) ! recurrence relation 1
          integer(ki), intent(in) :: n, k
          if (n <= 0 .or. k <= 0 .or. n < k) then
            ires = 0
          else if (n == 1 .and. k == 1) then 
            ires = 1
          else
            ires = p1(n - 1, k - 1) + p1(n - k, k)
          end if          
        end function p1  
        
        recursive integer(ki) pure function p2(n, k) result(ires) ! recurrence relation 2
          integer(ki), intent(in) :: n, k
          if (n < k .or. k <= 0) then
            ires = 0
          else if (n == k) then
            ires = 1
          else
            ires = p2(n - k, k) + p2(n, k + 1)
          end if
        end function p2
        
        integer(ki) function p3(n)
          integer(ki), intent(in) :: n
          integer(ki) :: i, k, isigma(n), p(0:n)
          do i = 1, n  
            isigma(i) = sum([(k, k = 1, i)], mask = mod(i, [(k, k = 1, i)]) == 0)
          end do   
          !
          p(0) = 1
          do i = 1, n
            p(i) = sum( isigma(:i) * p(i - 1:0:-1) ) / i
          end do  
          p3 = p(n) 
        end function p3

        integer(ki) pure function p4(n)
          integer(ki), intent(in) :: n
          integer(ki) :: i, j, k, penm(n), penp(n), p(0:n)
          do i = 1, n
            penm(i) = i * (3 * i - 1) / 2 ! pentagonal numbers
            penp(i) = i * (3 * i + 1) / 2 
          end do  
          !
          do i = 1, n
            p(0) = 1 ! 1 partition, i sigma_1
            p(i) = 0 
            do j = 1, int( ( 1 + sqrt(1.0 + 24 * i)) / 6 )
              p(i) = p(i) - (-1)**j * p(i - penm(j)) 
            end do  
            do j = 1, int( (-1 + sqrt(1.0 + 24 * i)) / 6 )
              p(i) = p(i) - (-1)**j * p(i - penp(j)) 
            end do
          end do
          p4 = p(n)
        end function p4
      end module m_partition   
  
      program test
        use m_partition
        implicit none
        integer(ki) :: i, j
        do i = 1, 300
         ! do j = i, 1, -1
         !   print *, i, j, p1(i, j)
         ! end do 
         !print *, 'p(n)=', sum([(p1(i, j), j = 1, i)]), p2(i, 1)
         !print *, 'p(n)=', i, sum([(p1(i, j), j = 1, i)])
         ! print *, 'p(n)=',i,  p2(i, 1)
         ! print '(a, i3, a, i18)', 'p(', i, ')=', sum([(p1(i, j), j = 1, i)])
         ! print '(a, i3, a, i18)', 'p(', i, ')=', p2(i, 1) 
          print '(a, i3, a, i18)', 'p(', i, ')=', p3(i) 
          print '(a, i3, a, i18)', 'p(', i, ')=', p4(i) 
        end do
      end program test

実行結果

p(  1)=                 1
p(  2)=                 2
p(  3)=                 3
p(  4)=                 5
p(  5)=                 7
p(  6)=                11
p(  7)=                15
p(  8)=                22
p(  9)=                30
p( 10)=                42
p( 11)=                56
p( 12)=                77
p( 13)=               101
p( 14)=               135
p( 15)=               176
p( 16)=               231
p( 17)=               297
p( 18)=               385
p( 19)=               490
p( 20)=               627
p( 21)=               792
p( 22)=              1002
p( 23)=              1255
p( 24)=              1575
p( 25)=              1958
p( 26)=              2436
p( 27)=              3010
p( 28)=              3718
p( 29)=              4565
p( 30)=              5604
p( 31)=              6842
p( 32)=              8349
p( 33)=             10143
p( 34)=             12310
p( 35)=             14883
p( 36)=             17977
p( 37)=             21637
p( 38)=             26015
p( 39)=             31185
p( 40)=             37338
p( 41)=             44583
p( 42)=             53174
p( 43)=             63261
p( 44)=             75175
p( 45)=             89134
p( 46)=            105558
p( 47)=            124754
p( 48)=            147273
p( 49)=            173525
p( 50)=            204226
p( 51)=            239943
p( 52)=            281589
p( 53)=            329931
p( 54)=            386155
p( 55)=            451276
p( 56)=            526823
p( 57)=            614154
p( 58)=            715220
p( 59)=            831820
p( 60)=            966467
p( 61)=           1121505
p( 62)=           1300156
p( 63)=           1505499
p( 64)=           1741630
p( 65)=           2012558
p( 66)=           2323520
p( 67)=           2679689
p( 68)=           3087735
p( 69)=           3554345
p( 70)=           4087968
p( 71)=           4697205
p( 72)=           5392783
p( 73)=           6185689
p( 74)=           7089500
p( 75)=           8118264
p( 76)=           9289091
p( 77)=          10619863
p( 78)=          12132164
p( 79)=          13848650
p( 80)=          15796476
p( 81)=          18004327
p( 82)=          20506255
p( 83)=          23338469
p( 84)=          26543660
p( 85)=          30167357
p( 86)=          34262962
p( 87)=          38887673
p( 88)=          44108109
p( 89)=          49995925
p( 90)=          56634173
p( 91)=          64112359
p( 92)=          72533807
p( 93)=          82010177
p( 94)=          92669720
p( 95)=         104651419
p( 96)=         118114304
p( 97)=         133230930
p( 98)=         150198136
p( 99)=         169229875
p(100)=         190569292
p(101)=         214481126
p(102)=         241265379
p(103)=         271248950
p(104)=         304801365
p(105)=         342325709
p(106)=         384276336
p(107)=         431149389
p(108)=         483502844
p(109)=         541946240
p(110)=         607163746
p(111)=         679903203
p(112)=         761002156
p(113)=         851376628
p(114)=         952050665
p(115)=        1064144451
p(116)=        1188908248
p(117)=        1327710076
p(118)=        1482074143
p(119)=        1653668665
p(120)=        1844349560
p(121)=        2056148051
p(122)=        2291320912
p(123)=        2552338241
p(124)=        2841940500
p(125)=        3163127352
p(126)=        3519222692
p(127)=        3913864295
p(128)=        4351078600
p(129)=        4835271870
p(130)=        5371315400
p(131)=        5964539504
p(132)=        6620830889
p(133)=        7346629512
p(134)=        8149040695
p(135)=        9035836076
p(136)=       10015581680
p(137)=       11097645016
p(138)=       12292341831
p(139)=       13610949895
p(140)=       15065878135
p(141)=       16670689208
p(142)=       18440293320
p(143)=       20390982757
p(144)=       22540654445
p(145)=       24908858009
p(146)=       27517052599
p(147)=       30388671978
p(148)=       33549419497
p(149)=       37027355200
p(150)=       40853235313
p(151)=       45060624582
p(152)=       49686288421
p(153)=       54770336324
p(154)=       60356673280
p(155)=       66493182097
p(156)=       73232243759
p(157)=       80630964769
p(158)=       88751778802
p(159)=       97662728555
p(160)=      107438159466
p(161)=      118159068427
p(162)=      129913904637
p(163)=      142798995930
p(164)=      156919475295
p(165)=      172389800255
p(166)=      189334822579
p(167)=      207890420102
p(168)=      228204732751
p(169)=      250438925115
p(170)=      274768617130
p(171)=      301384802048
p(172)=      330495499613
p(173)=      362326859895
p(174)=      397125074750
p(175)=      435157697830
p(176)=      476715857290
p(177)=      522115831195
p(178)=      571701605655
p(179)=      625846753120
p(180)=      684957390936
p(181)=      749474411781
p(182)=      819876908323
p(183)=      896684817527
p(184)=      980462880430
p(185)=     1071823774337
p(186)=     1171432692373
p(187)=     1280011042268
p(188)=     1398341745571
p(189)=     1527273599625
p(190)=     1667727404093
p(191)=     1820701100652
p(192)=     1987276856363
p(193)=     2168627105469
p(194)=     2366022741845
p(195)=     2580840212973
p(196)=     2814570987591
p(197)=     3068829878530
p(198)=     3345365983698
p(199)=     3646072432125
p(200)=     3972999029388
p(201)=     4328363658647
p(202)=     4714566886083
p(203)=     5134205287973
p(204)=     5590088317495
p(205)=     6085253859260
p(206)=     6622987708040
p(207)=     7206841706490
p(208)=     7840656226137
p(209)=     8528581302375
p(210)=     9275102575355
p(211)=    10085065885767
p(212)=    10963707205259
p(213)=    11916681236278
p(214)=    12950095925895
p(215)=    14070545699287
p(216)=    15285151248481
p(217)=    16601598107914
p(218)=    18028182516671
p(219)=    19573856161145
p(220)=    21248279009367
p(221)=    23061871173849
p(222)=    25025873760111
p(223)=    27152408925615
p(224)=    29454549941750
p(225)=    31946390696157
p(226)=    34643126322519
p(227)=    37561133582570
p(228)=    40718063627362
p(229)=    44132934884255
p(230)=    47826239745920
p(231)=    51820051838712
p(232)=    56138148670947
p(233)=    60806135438329
p(234)=    65851585970275
p(235)=    71304185514919
p(236)=    77195892663512
p(237)=    83561103925871
p(238)=    90436839668817
p(239)=    97862933703585
p(240)=   105882246722733
p(241)=   114540884553038
p(242)=   123888443077259
p(243)=   133978259344888
p(244)=   144867692496445
p(245)=   156618412527946
p(246)=   169296722391554
p(247)=   182973889854026
p(248)=   197726516681672
p(249)=   213636919820625
p(250)=   230793554364681
p(251)=   249291451168559
p(252)=   269232701252579
p(253)=   290726957916112
p(254)=   313891991306665
p(255)=   338854264248680
p(256)=   365749566870782
p(257)=   394723676655357
p(258)=   425933084409356
p(259)=   459545750448675
p(260)=   495741934760846
p(261)=   534715062908609
p(262)=   576672674947168
p(263)=   621837416509615
p(264)=   670448123060170
p(265)=   722760953690372
p(266)=   779050629562167
p(267)=   839611730366814
p(268)=   904760108316360
p(269)=   974834369944625
p(270)=  1050197489931117
p(271)=  1131238503938606
p(272)=  1218374349844333
p(273)=  1312051800816215
p(274)=  1412749565173450
p(275)=  1520980492851175
p(276)=  1637293969337171
p(277)=  1762278433057269
p(278)=  1896564103591584
p(279)=  2040825852575075
p(280)=  2195786311682516
p(281)=  2362219145337711
p(282)=  2540952590045698
p(283)=  2732873183547535
p(284)=  2938929793929555
p(285)=  3160137867148997
p(286)=  3397584011986773
p(287)=  3652430836071053
p(288)=  3925922161489422
p(289)=  4219388528587095
p(290)=  4534253126900886
p(291)=  4872038056472084
p(292)=  5234371069753672
p(293)=  5622992691950605
p(294)=  6039763882095515
p(295)=  6486674127079088
p(296)=  6965850144195831
p(297)=  7479565078510584
p(298)=  8030248384943040
p(299)=  8620496275465025
p(300)=  9253082936723602
続行するには何かキーを押してください . . .