! $Id: bbp.f90 834 2006-06-13 08:49:40Z olau $
!
! Copyright (c) Johannes Nendwich, http://ams.astro.univie.ac.at/~nendwich/
!

program bbp

implicit none

integer               :: i, k, l, n
integer, dimension(4) :: arr1, arr2, arr4
real   , dimension(4) :: arr3
real                  :: log2

log2 = log(2.0)

write (*, '(a)', advance='no') '3.'
do n = 0, 5791
   arr3 = [0.,0.,0.,0.]     ! ~ arr2 / arr4
   do k = 0, n
      l = n - k
      arr1 = [4,4,4,4]       ! 16^(2^i) mod arr4
      arr2 = [1,1,1,1]       ! 
      arr4 = [1,4,5,6] + 8*k ! 8k+j
      if (l/=0) then
         do i = 0, int(log(l+0.5)/log2)
            arr1 = mod(arr1**2,arr4)
            if (btest(l,i)) arr2 = mod(arr2*arr1,arr4)
         end do ! i
      end if
      arr3 = mod(arr3+real(arr2)/arr4,1.0)
   end do ! k
   write (*, '(z1)', advance='no') int(16*modulo(dot_product([4.,-2.,-1.,-1.],arr3),1.0))
end do
write (*,*)

end program bbp
