C Arrange n+1 weights of size 2...n+1 at equal distances
C on a circle such that their center of gravity is at (0,0)
C The weigth 1 is at (1,0)
C Author: Hugo Pfoertner http://www.pfoertner.org/
C Version history:
C Sep 20 2005 Variable length of fixed front part of
C             permuted vector
C Sep 19 2005 Read initial permutation from input
C             and keep 1st and 2nd entry fixed to enable
C             parallel runs for n>=16
C Sep 15 2005 Initial version
      parameter ( nmax=20 )
      implicit doubleprecision (d)
      integer aini(nmax)
      doubleprecision a(nmax)
      doubleprecision x(nmax), y(nmax), pi, rmin
      pi = atan2 ( 0.0D0, -1.0D0 )
2     continue
      nfix = 1
      write (*,*) ' Number of weights, fixed part:'
      read (*,*) n, nfix
      write (*,*) ' Initial arrangement'
      read (*,*) (aini(k),k=1,n)
C Create angular positions
      dphi = (pi+pi) / dble(n)
      x(1) = 1.0D0
      y(1) = 0.0D0
      do 10 i = 2, n
      x(i) = cos(dble(i-1)*dphi)
      y(i) = sin(dble(i-1)*dphi)
      write (*,*) i,x(i), y(i)
10    continue
C preset permutation array
      do 20 I = 1, n
      a(i) = dble(aini(i))
20    continue
      nbalan = 0
C      i2prev = 0
C Start of loop over all permutations
      rmin = 100000.0D0
100   continue
C Progress indicator
C      if ( a(2) .ne. i2prev ) then
C        write (*,*) ' a(2)=', a(2)
C        i2prev = a(2)
C      endif
C Accumulators for center of gravity
C First mass remains fixed at position (1,0)
      dcgx = 1.0D0
      dcgy = 0.0D0
      do 30 i = 2, n
      dcgx = dcgx + a(i) * x(i)
      dcgy = dcgy + a(i) * y(i)
30    continue
      dist = dcgx*dcgx + dcgy*dcgy
      delr = dist-rmin
      if ( delr .le. 1.0D-8 ) then
        if ( dist .lt. rmin-1.0D-8 ) then
          rmin = dist
          nbalan = 1
          write (*,1000) sqrt(rmin)/dble(n), nbalan,
     &                  (nint(a(k)),k=1,n)
1000      format ( F12.9, I5, 2x, 18I3 )
        else
C      if ( abs(dcgx) .lt. 1.0E-4 .and. abs(dcgy) .lt. 1.0E-4 ) then
          nbalan = nbalan + 1
          write (*,1000) sqrt(rmin)/dble(n), nbalan,
     &                   (nint(a(k)),k=1,n)
        endif
      endif
C Permute only the upper n-nfix positions
      call lpg ( n-nfix, a(nfix+1), next )
C Check if permutation loop is finished
      if ( next .ne. 0 ) goto 100
      write (*,*) ' Enter 0 to stop'
      read (*,*)  i
      if ( i .ne. 0 ) goto 2
      end