連分数展開でπその他

Numerical Recipes 2nd ed. §5.2の計算法を使って連分数展開の式を求めました。

Fortran2003では、ALLOCATABLEな配列に、最初に確保した大きさとは異なる配列代入を行うと、自動的に新しい大きさにreallocateされます。これはバグの元になりそうな機能ですが、array constructorと組み合わせて使うと、定数の数をいちいち数えなくてもいので、なかなか便利です。

Intel Fortran v.11では、コンパイラ・オプションとして/assume:realloc_lhsを加えないと、この機能は有効になりません。

最近、数学系のサイトでいくつか評判になっているものがあるようです。
・オンライン整数列辞典
・Inverse symbolic calculator
・Wolfram alpha

オンライン整数列辞典は、整数列を入れてやると、漸化式等の情報を返してくれます。
Inverse symbolic calculatorは、無理数っぽい数値などを入れてやると、表式などを教えてくれます。Mapleのエンジンを使っているようです。
Wolfram alphaはMathematica的な処理をした検索情報を返してくれます。連分数展開の係数の一部はここからコピペしてきました。

■実行結果

■ソース・プログラム

オプション　/assume:realloc_lhs が必要。

MODULE m_conti
  IMPLICIT NONE
  INTEGER, PARAMETER :: kd = KIND(1.0q0), nmax = 35
 
 CONTAINS

  REAL(kd) FUNCTION continued_fraction(a, b)
    IMPLICIT NONE
    REAL(kd), INTENT(IN) :: a(:), b(0:)
    REAL(kd) :: c0, c1, d0, d1, e, f
    INTEGER :: i
 
    f  = b(0)
    c0 = f
    d0 = 0.0_kd
    DO i = 1, SIZE(b) - 1
     d1 = b(i) + a(i) * d0
     IF (d1 == 0.0_kd) d1 = TINY(0.0_kd)
     c1 = b(i) + a(i) / c0
     IF (c1 == 0.0_kd) c1 = TINY(0.0_kd)
     d0 = 1.0_kd / d1
     c0 = c1
     e = c0 * d0
     f = f * e
     IF (ABS(e - 1.0_kd) < EPSILON(1.0_kd)) EXIT
    END DO
    continued_fraction = f

   RETURN
 END FUNCTION continued_fraction

END MODULE m_conti

!==============================================================

PROGRAM conti
 USE m_conti
 IMPLICIT NONE
 REAL(kd), ALLOCATABLE :: a(:), b(:), x
 REAL(kd) :: pi
 INTEGER :: i
 !
 ! pi
 !
 ALLOCATE( b(0:0) ) ! Set LBOUND of b to 0
 b = [REAL(kd):: 3, 7, 15, 1, 292, 1,  1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1,&
                84, 2,  1, 1,  15, 3, 13, 1, 4, 2, 6, 6, 99, 1, 2, 2, 6, 3, 5, 1 ]
 a = [( 1.0_kd, i = 1, UBOUND(b,1) )]
 pi = continued_fraction(a, b)
 PRINT *, '          ', pi
 PRINT *, '      pi =', 4.0_kd * ATAN(1.0_kd)
 !
 ! Tangent(pi/4)
 !
 x = pi / 4
 a = [x, (-x**2, i = 2, UBOUND(a, 1))]
 b = [TINY(0.0_kd), (REAL(2 * i - 1, kd), i = 1, UBOUND(b, 1))]
 PRINT *, '           ', continued_fraction(a, b)
 PRINT *, 'tan(pi/4)= ', TAN(x)
 !
 ! Euler-Mascheroni Constant
 !
 b = [TINY(0.0_kd), 1, 1, 2, 1, 2, 1, 4, 3, 13, 5, 1,  1, 8, 1, 2,  4, 1,   1, 40, 1, &
                11, 3, 7, 1, 7, 1, 1, 5, 1, 49, 4, 1, 65, 1, 4, 7, 11, 1, 399, 2]
 a = [( 1.0_kd, i = 1, UBOUND(b,1) )]
 PRINT *, '           ', continued_fraction(a, b)
 PRINT *, '   gamma = ', 0.5772156649015328606065120900824024310421593359399235988057_kd
 !               
 ! exp(1.0)
 !
 b = [2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1, 14,&
      1, 1,16, 1, 1,18, 1, 1,20, 1, 1,22, 1, 1, 24, 1, 1, 26, 1, 1, 28 ]
 a = [( 1.0_kd, i = 1, UBOUND(b,1) )]
 PRINT *, '          ', continued_fraction(a, b)
 PRINT *, '      e  =', EXP(1.0_kd)
 !
 !
 STOP
END PROGRAM conti