2017-08-19

分割数 p(200) 直接計算

George E. Andrews, The Theory of Partitions (Encyclopedia of Mathematics and its Applications) の中に、Hardy-Ramanujan-Rademacher の p(n) の式の p(200) に対する具体的計算例が載っています。

追記H29/12/29　結果のズレに関して続編あり
分割数 p(200) 直接計算 (その２） - fortran66のブログ

amazon AA

Chap.5.1 p.70
f:id:fortran66:20170819014720p:plain

以前、計算した時のプログラムで再現してみたのですが、第三項が少し違うという結果が･･･

     1             3972998993185.896                        -0.000
     2                     36282.978                        -0.000
     3                       -87.584                         0.000
     4                         5.147                         0.000
     5                         1.424                         0.000
     6                         0.071                         0.000
     7                        -0.000                        -0.000
     8                         0.044                         0.000

p(200)             3972999029387.976                        -0.000
続行するには何かキーを押してください . . .

原因は今のところ不明です。
無限級数を第 50 項まで計算すると、整数化するので間違っているわけでもなさそうな･･･？
はぁ、やな感じ～。

     1             3972998993185.896                        -0.000
     2                     36282.978                        -0.000
     3                       -87.584                         0.000
     4                         5.147                         0.000
     5                         1.424                         0.000
     6                         0.071                         0.000
     7                        -0.000                        -0.000
     8                         0.044                         0.000
     9                         0.026                         0.000
    10                         0.009                         0.000
    11                         0.000                         0.000
    12                        -0.004                        -0.000
    13                        -0.000                        -0.000
    14                         0.000                         0.000
    15                        -0.003                        -0.000
    16                        -0.004                        -0.000
    17                         0.000                        -0.000
    18                        -0.002                         0.000
    19                        -0.000                        -0.000
    20                         0.001                        -0.000
    21                        -0.000                         0.000
    22                        -0.000                        -0.000
    23                         0.002                         0.000
    24                        -0.001                         0.000
    25                         0.002                        -0.000
    26                        -0.000                         0.000
    27                        -0.001                        -0.000
    28                         0.000                         0.000
    29                        -0.000                        -0.000
    30                        -0.001                        -0.000
    31                        -0.000                         0.000
    32                         0.000                        -0.000
    33                        -0.000                        -0.000
    34                         0.000                        -0.000
    35                         0.000                        -0.000
    36                         0.000                         0.000
    37                         0.001                        -0.000
    38                        -0.000                         0.000
    39                        -0.000                        -0.000
    40                        -0.000                        -0.000
    41                         0.001                         0.000
    42                        -0.000                         0.000
    43                        -0.001                        -0.000
    44                         0.000                        -0.000
    45                        -0.000                         0.000
    46                        -0.000                        -0.000
    47                        -0.001                        -0.000
    48                        -0.000                        -0.000
    49                        -0.000                        -0.000
    50                         0.000                         0.000

p(200)             3972999029388.000                        -0.000
続行するには何かキーを押してください . . .

fortran66.hatenablog.com

プログラム

以前のとちょっとだけ変えました（一部 Andrews の定義に書き直し）が結果は変わってません。

無限級数を第 8 項で打ち切っています。

    module m_part
      implicit none
      integer , parameter :: kd = kind(1.0q0)
      real(kd), parameter :: pi = 4 * atan(1.0_kd)
    contains  
      complex(kd) function cpartition(n)
        integer, intent(in) :: n
        integer, parameter :: kmax = 8
        real(kd) :: f, s, t, u
        complex(kd) :: c
        integer :: j, m, k
        cpartition = (0.0_kd, 0.0_kd)
        do k = 1, kmax
          f = sqrt(k / 2.0_kd) / pi 
          c = (0.0_kd, 0.0_kd)
          do m = 1, k
            if (igcd(m, k) == 1) then
              s = -real(2 * n * m, kd) / k + s_hk(m, k)
              c = c + exp(cmplx(0.0_kd, pi * s, kind = kd)) 
            end if
          end do
          cpartition = cpartition + f * c * da(k, n) 
          print '(i6, 2f30.3)', k, f * c * da(k, n) 
        end do  
      contains
        real(kd) function s_hk(h, k)
          integer, intent(in) :: h, k
          real(kd) :: t, u
          integer :: j
          s_hk = 0.0_kd
          do j = 1, k - 1
             t = real(m * j, kd) / k
             u = real(    j, kd) / k
             s_hk = s_hk + (u - floor(u) - 0.5_kd) * (t - floor(t) - 0.5_kd)                
          end do
        end function s_hk
      end function cpartition
  
      integer pure function igcd(ia, ib)
        integer, intent(in) :: ia, ib
        integer :: m(2)
        m = [ia, ib]
        do 
          m = [m(2), mod(m(1), m(2))]
          if (m(2) == 0) exit
        end do
        igcd = m(1)
      end function igcd  
      
      real(kd) pure function da(k, n) ! d/dx sinh( f * xs ) / xs ; xs = sqrt(x - 1/24)
        integer, intent(in) :: k, n 
        real(kd) :: xs, f
        xs = sqrt( real(n, kd) - 1.0_kd / 24.0_kd )
        f  = pi / k * sqrt(2.0_kd / 3.0_kd)
        da = ( f * cosh(f * xs) - sinh(f * xs) / xs ) / (2 * xs * xs) 
      end function da
    end module m_part
  
    program Pn
      use m_part
      implicit none
      integer :: i
      complex(kind = kd) :: z
      do i = 200, 200 
          z = cpartition(i)
        print *
        print '(a, i3, a, 2f30.3)', 'p(', i, ')', z
      end do  
    end program Pn

2017-08-17

メモ帳 Jacobiの三重積から

Jacobi の三重積
${\displaystyle\prod_{n=1}^\infty(1-x^{2n})(1+x^{2n-1}z)(1+x^{2n-1}z^{-1})=\sum_{n=-\infty}^{n=\infty}x^{n^2}z^n}$
で、 $x=x^{1/2}, z=x^{1/2}\zeta$ とおき、さらに $\zeta\rightarrow-1$ とすると、次の公式が得られます。
${\displaystyle\prod_{n=1}^\infty(1-x^n)^3=\sum_{m=0}^\infty(-1)^m(2m+1)x^{{1\over2}m(m+1)}}$
多分収束性の議論を気にしなければ、いきなり $z=-x^{1/2}$ でもいいんだと思いますｗ
これはオイラーの五角数定理のあとによく出てくる式です。

なんとなく、無限積と無限和でどちらが収束が早いのかと思って計算して見ました。

結果は、無限和の方が早く収束するようです。まだ詳しく見ていないのでメモ帳代わりに。

実行結果

           1 prod  0.512000000000000
           2 prod  0.452984832000000
           3 prod  0.442199937191510
           4 prod  0.440080771777257
           5 prod  0.439658429414744
           6 prod  0.439574020398704
           7 prod  0.439557140972379
           8 prod  0.439553765182178
           9 prod  0.439553090027941
          10 prod  0.439552954997245
          11 prod  0.439552927991112
          12 prod  0.439552922589886
          13 prod  0.439552921509640
          14 prod  0.439552921293591
          15 prod  0.439552921250382
          16 prod  0.439552921241740
          17 prod  0.439552921240011
          18 prod  0.439552921239665
          19 prod  0.439552921239596
          20 prod  0.439552921239583
          21 prod  0.439552921239580
          22 prod  0.439552921239579
          23 prod  0.439552921239579
          24 prod  0.439552921239579
          25 prod  0.439552921239579

           1 sum   0.400000000000000
           2 sum   0.440000000000000
           3 sum   0.439552000000000
           4 sum   0.439552921600000
           5 sum   0.439552921239552
           6 sum   0.439552921239579
           7 sum   0.439552921239579
           8 sum   0.439552921239579
           9 sum   0.439552921239579
          10 sum   0.439552921239579
          11 sum   0.439552921239579
          12 sum   0.439552921239579
          13 sum   0.439552921239579
          14 sum   0.439552921239579
          15 sum   0.439552921239579
          16 sum   0.439552921239579
          17 sum   0.439552921239579
          18 sum   0.439552921239579
          19 sum   0.439552921239579
          20 sum   0.439552921239579

  0.439552921239579
  0.439552921239579

ソース・プログラム

素朴に定義通り計算して見ました。

関数の中に PRINT文を置くというお行儀の悪いプログラムなので、よゐこのみんなは真似しないでください。

    program Jacobi3
      implicit none
      integer, parameter :: kd = kind(1.0d0), n = 100
      real(kd) :: x, p, s
      
      x =  0.2_kd
      p = jacobi_prod(x)
      print *
      s = jacobi_sum (x)
      
      print *
      print *, p
      print *, s
    contains
      real(kd) function jacobi_prod(x) result(p)
        real(kd), intent(in) :: x
        real(kd) :: q
        integer :: i
        p = 1.0_kd
        q = 1.0_kd
        do i = 1, 25
          q = q * x
          p = p * (1.0_kd - q)**3
          print *, i, 'prod', p
        end do   
      end function jacobi_prod

      real(kd) function jacobi_sum(x) result(s)
        real(kd), intent(in) :: x
        real(kd) :: t, u
        integer :: i
        s = 1.0_kd
        t = 1.0_kd
        u = 1.0_kd
        do i = 1, 20
          u = u * x  
          t = -t * u 
          s = s + t * (2 * i + 1)
          print *, i, 'sum ', s  
        end do    
      end function jacobi_sum
    end program Jacobi3

2017-08-17

【乞食速報】MIT PRESS タダ化（ただし古本に限る）

www.openculture.com
Internet Archive から読めるようになるようです。Internet Archive には昔の大型計算機の FORTRAN マニュアル類も沢山上がっており、たまに見ると面白いです。

MIT Press からは、Fortran95 Handbook、 MPI マニュアル類、High Performace Fortran など Fortran がらみの本も多く出ています。まだスキャン途上のようで、これから増えてゆくものと思われます。

その早速の成果が FORTRAN ぬり絵
www.reddit.com

なお reddit スレによると datamation 誌 1963 年九月号にもぬり絵があるとのことでスキャンがうｐされています。
f:id:fortran66:20170817012840p:plain

直リン http://thecomputerboys.com/wp-content/uploads/2010/12/programmers-coloring-book.pdf

類書のケンブリッジ大学出版の Illustrating FORTRAN のほうはまだ絶版になっていません。こちらはもっとまじめな内容になっています。

2017-08-16

メモ帳 S-function の正規化

Littlewood の S-function の standardization

夏休みの勉強ｗ

分割の共役

$[\lambda_1,\lambda_2,..., \lambda_p]$ 　
の共役は
$[p^{\lambda_p},(p-1)^{\lambda_{p-1}-\lambda_p},...,1^{\lambda_1-\lambda_2}]$

main = print (conj [4,3,3,2,2])

conj :: Partition -> Partition
conj [] = []
conj xs = conj' (xs ++ [0]) 

conj' :: Partition -> Partition
conj' [] = [] 
conj' [x] = [] 
conj' xs = (replicate (l2 - l1) p) ++ (conj' (init xs))
    where 
      p  = (length xs) - 1
      l1 = last xs
      l2 = last (init xs)

GHCi, version 8.0.1
   
[5,5,3,1]

□□□□
□□□
□□□
□□
□□

共役を取って

■■■■■
■■■■■
■■■
■

分割｛λ｝の正規化

分割｛λ｝を S-function とみて、｛λ｝が弱い減少列になっていない時に弱い減少列に直す。

規則１　尻の 0 は取り去る。
規則２　尻に負数が出たら 0。（ここではヌルとする）
規則３　隣り合う数の入れ替えは $\{...\lambda_i, \lambda_{i+1},...\}$ = $-\{...\lambda_{i+1}-1, \lambda_{i}+1,...\}$

行き当たりばったりで適当ｗ

main = print (standard (1, [2,1,3,6,1]))

type Partition = [Int]
type PartitionF = (Int, Partition) -- Partition with a (phase)factor

standard :: PartitionF -> PartitionF
standard (f, xs) = if stdq' xs then (f, xs)
                               else standard (std (f, xs))

std :: PartitionF -> PartitionF
std (f, [])= (f, [])
std (f, xs) = if (last xs) <  0 then (f, [])  
         else if (last xs) == 0 then std (f, init xs)
         else if (p1 xs)        then (f, [])
         else stdq (f, xs)

p1 :: Partition -> Bool
p1 xs = or [x == (y-1)|(x,y) <- ys ]
  where 
    ys = zip xs (tail xs)  
    
stdq :: PartitionF -> PartitionF
stdq (f, []) = (f, [])
stdq (f, xs) = if stdq' xs then (f, xs)
                           else std' (f, xs)
         
std' :: PartitionF -> PartitionF                           
std' (f, []) = (f, [])
std' (f, [x]) = (f, [x])
std' (f, x:y:xs) = if (x < y) then ex1 (f, x:y:xs)
                              else (f', x:xs')
    where 
      (f', xs') = std' (f, y:xs)

stdq' :: Partition -> Bool
stdq' xs = and [x >= y|(x,y) <- ys ]
  where 
    ys = zip xs (tail xs)     
    
ex1 :: PartitionF -> PartitionF
ex1 (f, x:y:xs) = (-f, (y-1):(x+1):xs)

GHCi, version 8.0.1
   
(1,[3,3,3,3,1])

2017-08-15

メモ帳　Modula-2 への不満リスト by Ogilvie

John W. L. Ogilvie, Modula-2 Programming (1985)

盆休みなのでチラチラと Ogilvie 氏の Modula-2 Programming を眺めていましたが、あんまり面白くないので斜め読みです。Modula-2 の module による visibility control のほかに、procedure type, type bound procedure, coroutine なども一通り説明されています。

pascal から受け継いだ variant record は、modula-2 の時点では残っていますが、 oberon では問題多しとして削除されています。別のところで Wirth は、使い方さえ間違わなければ便利だが危険性のある機能を入れると、抜け穴として常時危険に乱用されるとぼやいています。今は variant record は OO の polymorphism で解決するのだと思います。

python の妙なモジュール読み込み書式は、modula-2 そのままであることも知りました。元はヨーロッパ系の教育目的の言語だからでしょうか？

Modern Fortran は C 系統ではなく ALGOL 系統を取り込んだと言われていますが、基本的な発想が分かってきた気がします。

第２２章 A Discussion of Modular-2's Shortcomings

第２２章から不満リストをコピーしておきます。今の Fortran では解消しているものが多いですが、文字列処理まわりは進歩がゆっくりです。

f:id:fortran66:20170815162618p:plain

f:id:fortran66:20170815162636p:plain

f:id:fortran66:20170815162715p:plain

2017-08-14

Fortran のオブジェクト指向と N. Wirth

Fortran の OO と N. Wirth

以前紹介しました下記シンポジウムの M. Cohen 氏による講演資料にあるオブジェクト指向 (OO) を
Oberon から借りてきたという記述を見て、色々謎が解けた気がするのでメモしておきます。

fortran66.hatenablog.com

f:id:fortran66:20170814005725p:plain

Fortran ~ modernization or cumbersome ~
Project Editor, ISO/IEC Fortran standard / Numerical Algorithms Group Oxford/Tokyo
Malcolm Cohen
【発表資料】
Abstract: The Fortran standard has already undergone significant revision in this (21st) century, and further revision is underway. The larger changes in Fortran and 2008 will be described, with their effects on Fortran programming, compilers, and parallelization. Some of the forthcoming changes in Fortran 2015 will also be mentioned.

N. Wirth, MODULA-2 and object-oriented programming (1990)

MODULA-2 and object-oriented programming - ScienceDirect

この論文中の図１に、一般的 OO 用語と Wirth 流 OO 用語の対応表があり、それが Fortran における OO 用語におおむね対応しています。

f:id:fortran66:20170814010334p:plain

この論文は結構面白くて、序論ではコンピュータ界が流行り用語に弱いという今に続く悪習を腐しています。Wirth なりの OO の定義なども述べています。class を module 代わりに使っている（シングルトン？）ことも、概念の違うものを一つのもので表わすと混乱すると批判していて、至極もっともなことをすでに語っています。

よく読むと、oberon の type-bound procedure は、インスタンスごとに任意にくっつける一般的でないものらしく、今の Fortran のもの（つまり一般的なものと同じ）は Oberon-2 によるものになっているようです。

もう少し、このあたりを調べてみたいと思います。

[H29.8.24 追記]
N. Wirth, "From Modula to Oberon" (1989)

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.38.3932

に、より詳しい情報があります。

2017-08-11

メモ帳　LINK集

GNU Scientific Library Fortran 用インターフェース

github.com

Pyplot を Fortran から

ブログ記事
Fortran and Pyplot – Degenerate Conic

github
github.com

比べるのもおこがましいですがｗｗ、大昔に私が思い付きで試してみたようなことを、洗練されたインターフェースでまとめています。
fortran66.hatenablog.com

C++の場合、直接的な方法のライブラリがあるようです。
github.com

ブログは Fortran 専門で面白い記事があります。

フランス製の謎のソフトIRPF90

IRPF90: A Fortran Code Generator for High-Performance Computing
IRPF90 - Irpf90

github.com

www.anl.gov

Abstract: IRPF90 is a Fortran code generator thathelps the development of large Fortran codes. In Fortran programs, the programmer has to focus on the order of the instructions; before using a variable, the programmer has to be sure that it has already been computed in all possible situations. For large codes, this is a common source of error.
With IRPF90, most of the order of instructions is handled by the pre-processor, and a mechanism guarantees that every entity is built before being used. This mechanism relies on the {needs/needed by} relations between the entities, which are built automatically. The consequence is that the programmer does not need to know the production tree of each entity.
Codes written with IRPF90 execute usually faster than Fortran programs and are quicker to write and easier to maintain than standard Fortran programs.

"Modula-2 and Oberon" (2007) 講演の動画あり（有料）

N. Wirth
Modula-2 and Oberon

PDFは検索すると拾えます。