[Shootout-list] fannkuch.f90 update

[Simon Geard]

This is a multi-part message in MIME format.
--------------050903060300020003010902
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit

Actually it's only a comment that I've changed but it does make it a lot 
clearer!

Simon Geard

--------------050903060300020003010902
Content-Type: text/plain;
 name="fannkuch.f90"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="fannkuch.f90"

! Fannkuch access function implementation
! Simon Geard, 1/12/04
!
! Building info.
! ==============
!
! Linux  - using the Intel Fortran90 compiler:
!
!          ifort fannkuch.f90 -O3 -static-libcxa -o fannkuch
!
! WinXP  - Compaq Visual Fortran 6.6c
!
!          f90 fannkuch.f90 /link /libpath:"d:\Program Files\Microsoft Visual Studio\df98\lib"
!
! Cygwin - g95 compiler
!
!          g95 fannkuch.f90 -O3 -o fannkuch.exe
!
!!$"Take a permutation of {1,...,n}, for example: {4,2,1,5,3}. Take the first element, here 4, and reverse the order of the first 4 elements: {5,1,2,4,3}. Repeat this until the first element is a 1, so flipping won't change anything more: {3,4,2,1,5}, {2,4,3,1,5}, {4,2,3,1,5}, {1,3,2,4,5}. Count the number of flips, here 5. Do this for all n! permutations, and record the maximum number of flips needed for any permutation. The conjecture is that this maximum count is approximated by n*log(n) when n goes to infinity.
!!$
!!$FANNKUCH is an abbreviation for the German word Pfannkuchen, or pancakes, in analogy to flipping pancakes."
!!$
!!$Correct output N = 7 is:
!!$
!!$Pfannkuchen(7) = 16
!!$
!!$
!!$Correct output N = 8 is:
!!$
!!$Pfannkuchen(8) = 22
!!$
!!$
!!$Correct output N = 9 is:
!!$
!!$Pfannkuchen(9) = 30
!!$
!!$
!!$Correct output N = 10 is:
!!$
!!$Pfannkuchen(10) = 38

program fannkuch
  implicit none
  integer count, pcount, max_pc, i, k, num
  character(len=8) argv
  integer, dimension(:), allocatable :: data, refData

  call getarg(1,argv)
  read(argv,*) num
  allocate(data(num))
  allocate(refData(num))
  data = (/ (i,i=1,num) /)
  refData = data

  count = 0
  max_pc = factorial(num)
  do i=1,max_pc
     count = max(count,countFlips(data))
     call nextPerm(data)
  end do

  print *,'number of flips ',count
contains

  recursive integer function factorial(n) result(if)
    integer, intent(in) :: n
    if (n == 1) then
       if = 1
    else
       if = n*factorial(n-1)
    end if
  end function factorial

  ! Reverse an array
  subroutine reverse(data)
    integer, dimension(:), intent(inout) :: data
    integer, dimension(size(data)) :: work
    integer i
    do i=1,size(data)
       work(i) = data(size(data)-i+1)
    end do
    data = work
  end subroutine reverse

  ! Count the number of flips in a permutation
  integer function countFlips(data)
    integer, dimension(:), intent(in) :: data
    integer, dimension(size(data)) :: work
    work = data
    countFlips = 0
    do
       if (work(1) <= 1) exit
       call reverse(work(1:work(1)))
       countFlips = countFlips + 1
    end do
  end function countFlips

  ! Get the next permutation
  subroutine nextPerm(data)
    integer, dimension(:), intent(inout) :: data
    integer, dimension(:), allocatable, save :: dcheck
    integer, save :: N = 0
    integer, save :: s, p
    integer :: i, j, index

    if (N == 0) then
       N = size(data)
       s = (N+1)*N/2
       p = factorial(N)
       allocate(dcheck(N))
    end if

    do
       ! Add one to the current number in modulo N
       index = N
       do
          if (data(index) < N) then
             data(index) = data(index) + 1
             exit
          end if
          data(index) = 1
          index = index-1
       end do

       ! Need to ensure the new number has each digit once. I'm checking the sum and prduct first
       ! because I it is much (4x) faster. Ideally the number would be correct by construction and
       ! no such test would be necessary.
       if (sum(data) == s .and. product(data) == p) then
          ! Now ensure each digit is distinct
          dcheck = (/ (0,i=1,N) /)
          do i=1,N
             do j=1,N
                if (data(j) == i) then
                   dcheck(i) = 1
                   exit
                end if
             end do
          end do
          if (product(dcheck) == 1) then
             exit
          end if
       end if
    end do
  end subroutine nextPerm

end program fannkuch

--------------050903060300020003010902--