The Computer Language
Benchmarks Game

pidigits Fortran Intel program

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
    

notes, command-line, and program output

NOTES:
32-bit Ubuntu one core
Intel(R) Fortran Compiler XE for applications running on IA-32, Version 15.0.2.164 Build 20150121
Copyright (C) 1985-2015 Intel Corporation.  All rights reserved.
FOR NON-COMMERCIAL USE ONLY


Wed, 25 Mar 2015 01:15:11 GMT

MAKE:
/usr/local/src/intel/bin/ifort -O3 -xHost -ipo -lgmp pidigits.f90 -o pidigits.ifc_run
rm pidigits.f90
0.75s to complete and log all make actions

COMMAND LINE:
./pidigits.ifc_run 2000

PROGRAM FAILED 


PROGRAM OUTPUT:
3141592653	:10
5897932384	:20
6264338327	:30
9502884197	:40
1693993751	:50
0582097494	:60
4592307816	:70
4062862089	:80
9862803482	:90
5342117067	:100
9821480865	:110
1328230664	:120
7093844609	:130
5505822317	:140
2535940812	:150
8481117450	:160
2841027019	:170
3852110555	:180
9644622948	:190
9549303819	:200
6442881097	:210
5665933446	:220
1284756482	:230
3378678316	:240
5271201909	:250
1456485669	:260
2346034861	:270
0454326648	:280
2133936072	:290
6024914127	:300
3724587006	:310
6063155881	:320
7488152092	:330
0962829254	:340
0917153643	:350
6789259036	:360
0011330530	:370
5488204665	:380
2138414695	:390
1941511609	:400
4330572703	:410
6575959195	:420
3092186117	:430
3819326117	:440
9310511854	:450
8074462379	:460
9627495673	:470
5188575272	:480
4891227938	:490
1830119491	:500
2983367336	:510
2440656643	:520
0860213949	:530
4639522473	:540
7190702179	:550
8609437027	:560
7053921717	:570
6293176752	:580
3846748184	:590
6766940513	:600
2000568127	:610
1452635608	:620
2778577134	:630
2757789609	:640
1736371787	:650
2146844090	:660
1224953430	:670
1465495853	:680
7105079227	:690
9689258923	:700
5420199561	:710
1212902196	:720
0864034418	:730
1598136297	:740
7477130996	:750
0518707211	:760
3499999983	:770
7297804995	:780
1059731732	:790
8160963185	:800
9502445945	:810
5346908302	:820
6425223082	:830
5334468503	:840
5261931188	:850
1710100031	:860
3783875288	:870
6587533208	:880
3814206171	:890
7766914730	:900
3598253490	:910
4287554687	:920
3115956286	:930
3882353787	:940
5937519577	:950
8185778053	:960
2171226806	:970
6130019278	:980
7661119590	:990
9216420198	:1000
9380952572	:1010
0106548586	:1020

forrtl: severe (174): SIGSEGV, segmentation fault occurred

Stack trace terminated abnormally.