Spigot 法

1991 年に Stanley Rabinowitz が、いくつかの数学定数を整数演算だけで 1 桁づつ求めてゆくアルゴリズムを論文にしましたが、その時アブストラクトに円周率 $\pi$ を小数点以下 1000 桁求める 14 行の FORTRAN プログラムを示しました。

ここではそれを modern fortran に書き換えます。と言ってもほとんどそのままですが･･･

オリジナル

FORTRAN 77　元のもの

小文字を使っていますが、77 末期にはありがちでした。L は小文字が見間違えやすいので大文字で書いているのでしょう。スカラー変数は暗黙の型宣言に従っています。桁上がりをする carry が karry と書かれているのはそのためと思われます。配列宣言は暗黙の型を破っていますが、こういうのはいちいち宣言部に戻って型を確かめるのが面倒なので嫌われていました。

      integer vect(3350), buffer(201)
      data vect/3350*2/, more/0/
      DO 2 n=1, 201
          karray = 0
      DO 3 L=3350, 1, -1
          num = 100000*vect(L) + karry*L
          karry = num/(2*L-1)
3         vect(L) = num - karry*(2*L-1)
      k = karry/100000
      buffer(n) = more+k
2     more = karry - k*100000
      write(*,100) buffer
100   format(1x ,I1,'.'/(1x,10I5.5))
      end

なおこのプログラムでは桁上がり処理を簡略化しているので 5 桁以上繰り上がりが出る場合このままではまずいようですが、$\pi$ 1000 桁まででは大丈夫な模様です。

3.
14159265358979323846264338327950288419716939937510
58209749445923078164062862089986280348253421170679
82148086513282306647093844609550582231725359408128
48111745028410270193852110555964462294895493038196
44288109756659334461284756482337867831652712019091
45648566923460348610454326648213393607260249141273
72458700660631558817488152092096282925409171536436
78925903600113305305488204665213841469519415116094
33057270365759591953092186117381932611793105118548
07446237996274956735188575272489122793818301194912
98336733624406566430860213949463952247371907021798
60943702770539217176293176752384674818467669405132
00056812714526356082778577134275778960917363717872
14684409012249534301465495853710507922796892589235
42019956112129021960864034418159813629774771309960
51870721134999999837297804995105973173281609631859
50244594553469083026425223082533446850352619311881
71010003137838752886587533208381420617177669147303
59825349042875546873115956286388235378759375195778
18577805321712268066130019278766111959092164201989

modern fortran 書き換え

変数は明示的に宣言することにします。
do loop での文番号は do ... end do 構文にします。

data 文での変数や配列の初期化は modern fortran では推奨されません。

Format 文は直接 READ/WRITE に書き込みます。これは 77 でもできたのですが、90 以降では文字列を示す引用符としてシングル・ダブルどちらのクォーテーション対でも許されるようになったのでやりやすくなりました。

    program pi
        implicit none
        integer :: vect(3350) = 2, buffer(201)
        integer :: carry, n, L, k, more = 0, num
        do n = 1, 201
            carry = 0
            do L = 3350, 1, -1
                num = 100000 * vect(L) + carry * L
                carry = num / (2*L - 1)
                vect(L) = num - carry * (2*L - 1)
            end do
            k = carry / 100000
            buffer(n) = more + k
            more = carry - k * 100000
        end do    
        write(*, "(1x, I1, '.'/(1x, 10I5.5))") buffer
    end program pi