source code
! The Computer Language Benchmarks Game
! http://benchmarksgame.alioth.debian.org/
!
! contributed by Steve Decker
! compilation:
! g95 -O3 -funroll-loops -fomit-frame-pointer pidigits.f90
! ifort -O -ip pidigits.f90
module big_int_mod
implicit none
save
integer, parameter, private :: Pwr = 15, Base = 2**Pwr, maxDigs = 2558
type big_int
private
integer :: sigDigs
logical :: sign
integer, dimension(maxDigs) :: digits
end type big_int
interface assignment (=)
module procedure int_to_big_int
end interface
interface operator (*)
module procedure big_int_times_int
end interface
interface operator (+)
module procedure big_int_plus_big_int
end interface
interface operator (/)
module procedure big_int_div_big_int
end interface
contains
subroutine int_to_big_int(bi, n)
type(big_int), intent(inout) :: bi
integer, intent(in) :: n
integer :: i
if (n > 0) then
bi = big_int(1, .true., (/ n, (0, i = 2, maxDigs) /) )
else
bi = big_int(0, .true., 0)
end if
end subroutine int_to_big_int
function big_int_times_int(bi, n) result(c)
type(big_int), intent(in) :: bi
integer, intent(in) :: n
type(big_int) :: c
integer :: m, i, curDig, k, j, carry
c = big_int(0, .true., 0)
if (n == 0 .or. bi%sigDigs == 0) return
c%sign = n >= 0 .eqv. bi%sign
m = abs(n)
do i = 1, maxDigs
curDig = mod(m,Base)
k = 1
carry = 0
do j = i, i + bi%sigDigs + 1
c%digits(j) = c%digits(j) + curDig * bi%digits(k) + carry
carry = ibits(c%digits(j),Pwr,Pwr+1)
c%digits(j) = mod(c%digits(j),Base)
k = k + 1
end do
m = ibits(m,Pwr,Pwr+1)
if (m == 0) exit
end do
do j = i + bi%sigDigs, 1, -1
c%sigDigs = j
if (c%digits(j) /= 0) exit
end do
end function big_int_times_int
function big_int_plus_big_int(bi1, bi2) result(c)
type(big_int), target, intent(in) :: bi1, bi2
type(big_int) :: c
integer :: i, carry, n
type(big_int), pointer :: a, b
c = big_int(0, .true., 0)
if (bi1%sigDigs == 0) then
c = bi2
return
else if (bi2%sigDigs == 0) then
c = bi1
return
end if
if (bi1%sign .eqv. bi2%sign) then
c%sign = bi1%sign
carry = 0
n = max(bi1%sigDigs, bi2%sigDigs) + 1
do i = 1, n
c%digits(i) = bi1%digits(i) + bi2%digits(i) + carry
carry = ibits(c%digits(i),Pwr,Pwr+1)
c%digits(i) = mod(c%digits(i),Base)
end do
else
if (greater_in_mag(bi1, bi2)) then
a => bi1
b => bi2
else if (greater_in_mag(bi2, bi1)) then
a => bi2
b => bi1
else
return
end if
n = max(a%sigDigs, b%sigDigs)
c%sign = a%sign
do i = 1, n
if (a%digits(i) < b%digits(i)) then
a%digits(i) = a%digits(i) + Base
a%digits(i+1) = a%digits(i+1) - 1
end if
c%digits(i) = a%digits(i) - b%digits(i)
end do
end if
do i = n, 1, -1
c%sigDigs = i
if (c%digits(i) /= 0) exit
end do
end function big_int_plus_big_int
function big_int_div_big_int(a, b) result(c)
type(big_int), intent(in) :: a, b
integer :: c
integer :: k, carry, n, j
type(big_int) :: accumulator
c = 0
if (a%sigDigs == 0) return
accumulator = big_int(0, .true., 0)
do k = 0, Base-1
carry = 0
n = max(accumulator%sigDigs, b%sigDigs) + 1
do j = 1, n
accumulator%digits(j) = &
accumulator%digits(j) + b%digits(j) + carry
carry = ibits(accumulator%digits(j),Pwr,Pwr+1)
accumulator%digits(j) = mod(accumulator%digits(j),Base)
end do
do j = n, 1, -1
accumulator%sigDigs = j
if (accumulator%digits(j) /= 0) exit
end do
if (greater_in_mag(accumulator, a)) then
c = k
exit
end if
end do
end function big_int_div_big_int
logical function greater_in_mag(a, b)
type(big_int), intent(in) :: a, b
integer :: i
greater_in_mag = .false.
do i = max(a%sigDigs, b%sigDigs), 1, -1
if (a%digits(i) > b%digits(i)) then
greater_in_mag = .true.
exit
else if (a%digits(i) < b%digits(i)) then
exit
end if
end do
end function greater_in_mag
end module big_int_mod
module pi_mod
use big_int_mod
implicit none
contains
function lfts(k)
integer, intent(in) :: k
integer, dimension(2,2) :: lfts
lfts = reshape( (/ k, 0, 4*k + 2, 2*k + 1 /), (/ 2, 2 /) )
end function lfts
function comp1(a, b)
integer, dimension(2,2), intent(in) :: a
type(big_int), dimension(2,2), intent(in) :: b
type(big_int), dimension(2,2) :: comp1
comp1(1,1) = b(1,1)*a(1,1) + b(2,1)*a(1,2)
comp1(2,1) = b(1,1)*a(2,1) + b(2,1)*a(2,2)
comp1(1,2) = b(1,2)*a(1,1) + b(2,2)*a(1,2)
comp1(2,2) = b(1,2)*a(2,1) + b(2,2)*a(2,2)
end function comp1
function comp2(a, b)
type(big_int), dimension(2,2), intent(in) :: a
integer, dimension(2,2), intent(in) :: b
type(big_int), dimension(2,2) :: comp2
comp2(1,1) = a(1,1)*b(1,1) + a(1,2)*b(2,1)
comp2(2,1) = a(2,1)*b(1,1) + a(2,2)*b(2,1)
comp2(1,2) = a(1,1)*b(1,2) + a(1,2)*b(2,2)
comp2(2,2) = a(2,1)*b(1,2) + a(2,2)*b(2,2)
end function comp2
function prod(z, n)
type(big_int), dimension(2,2), intent(in) :: z
integer, intent(in) :: n
type(big_int), dimension(2,2) :: prod
prod = comp1(reshape( (/ 10, 0, -10*n, 1 /), (/ 2, 2 /) ), z)
end function prod
logical function safe(z, n)
type(big_int), dimension(2,2), intent(in) :: z
integer, intent(in) :: n
safe = n == (z(1,1) * 4 + z(1,2)) / (z(2,1) * 4 + z(2,2))
end function safe
integer function next(z)
type(big_int), dimension(2,2), intent(in) :: z
next = (z(1,1) * 3 + z(1,2)) / (z(2,1) * 3 + z(2,2))
end function next
end module pi_mod
program pidigits
use pi_mod
implicit none
character(len=12), parameter :: proForma = " " // achar(9) // ":"
type(big_int), dimension(2,2) :: z
integer :: num, y, i = 1, j = 1, k = 1
character(len=17) :: outLine = proForma
character(len=4) :: argv
call get_command_argument(1, argv)
read(argv, *) num
z(1,1) = 1; z(2,1) = 0; z(1,2) = 0; z(2,2) = 1
do
y = next(z)
if (safe(z, y)) then
write(unit=outLine(k:k), fmt="(i1)") y
if (k < 10 .and. i < num) then
k = k + 1
else
k = 1
write(unit=outLine(13:17), fmt="(i0)") i
write(*, "(a)") trim(outLine)
outLine = proForma
end if
if (i == num) exit
z = prod(z, y)
i = i + 1
else
z = comp2(z, lfts(j))
j = j + 1
end if
end do
end program pidigits