c------------------------------------------------------------
c Chapter 12: Random partition of an n-set(p97)
c------------------------------------------------------------
c   Name of subroutine: RANEQU
c
c   Algorithm:Generating random equivalence relation on
c             {1,2,....,n}
c
c   input:    n=4 
c   number of calls: run=40
c   compiler: f77 ranper.f ranequ_2.f
c------------------------------------------------------------

      parameter(n=4,run=40)
      integer q(n),b(n)
      do 40 i=1,run
      call ranequ(n,l,q,b)
      write(*,*),(q(j),j=1,n)
   40 continue
      stop
      end

c-----Subroutine begins here--------------------------------
      subroutine ranequ(n,l,q,b)
      integer q(n)
      real b(n)
      data nlast/1/
      b(1)=1.
      if(n.le.nlast) go to 10
      nm1=n-1
      do 5  l=nlast,nm1
      sum=1./float(l)
      l1=l-1
      if(l1.eq.0)goto 5
      do 6  k=1,l1
    6 sum=(sum+b(k))/float(l-k)
    5 b(l+1)=(sum+b(l))/float(l+1)
      nlast=n
   10 m=n
      l=0
   20 z=m*b(m)*rand(iseed)
      k=0
      l=l+1
   30 q(m)=l
      m=m-1
      if(m.eq.0)goto 50
   40 z=z-b(m)
      k=k+1
      z=z*k
      if(z)20,30,30
   50 call ranper(n,q,.false.)
      return
      end

