Bellard's Formulae

bellard.f90

      real(kind=16) function term (p)
          implicit none
          integer p
          real(kind=16) t1
          real(kind=16) t2
          real(kind=16) t3
          real(kind=16) t6
          real(kind=16) n

          n = dble(p)
          t1 = (-1) ** n
          t2 = 10 * n
          t3 = 2 ** t2
          t6 = 4 * n
          term = dble(t1 / t3 * (-32 / (t6 + 1) - 1 / (t6 + 3) + 256 / (1 &
          &+ t2) - 64 / (t2 + 3) - 4 / (t2 + 5) - 4 / (t2 + 7) + 1 / (t2 + 9) &
          &)) / 0.64D2
          return

          return
      end

      program  main
          implicit none
          integer :: n
          integer :: m
          integer :: max_m
          real(kind=16) :: temp_term
          real(kind=16) :: temp_pi
          real(kind=16) :: term
          temp_pi = 0
          do n=0,20
            temp_term = 0.0
            max_m = 13
            if (n.eq.0) max_m = 13 + 3

            do m = 0,12
                !write(0,*) n*13 + m
                temp_term = temp_term + term(n*13 + m)
            end do
            temp_pi = temp_pi + temp_term
            write(0,'(F192.190)') temp_term
          end do


      end program