fortran66のブログ

fortran について書きます。

メモ帳 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}でもいいんだと思いますw
これは オイラーの五角数定理のあとによく出てくる式です。

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

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

実行結果

           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