:	This is a shell archive.
:	Remove everything above this line and
:	run the following text with /bin/sh to create:
:	NJAS-CONTENTS
:	backtr.f
:	bincoef.f
:	det1.f
:	hpsort.f
:	lexsub.f
:	nexcom.f
:	nexequ.f
:	nexksb.f
:	nexpar.f
:	nexper.f
:	nexsub.f
:	nexvec.f
:	nxksrd.f
:	pawser.f
:	rancom.f
:	ranequ.f
:	ranksb.f
:	ranpar.f
:	ranper.f
:	ransub.f
:	renumb.f
: This archive created: Tue Jun 13 17:35:48 1995
cat << 'SHAR_EOF' > NJAS-CONTENTS
From njas@research.att.com Tue May 23 10:15:06 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id KAA12113 for <skiena@dimacs.rutgers.edu>; Tue, 23 May 1995 10:15:05 -0400
From: njas@research.att.com
Message-Id: <199505231415.KAA12113@dimacs.rutgers.edu>
Date: Tue, 23 May 95 09:31 EDT
To: skiena@dimacs.rutgers.edu
Subject: Re: Wilf's programs?
Status: RO

i found 4 fioles wilf1 .. wilf4 which contain this :

wilf1:      subroutine renumb(m,n,sig,tau,a)
wilf1:      subroutine spanfo(n,e,endpt,k,x,nv,y)
wilf1:      subroutine poly(n,a,x0,option,val,b)
wilf1:      subroutine chromp(n,e,endpt,a,b,c,stack,nstk)
wilf2:      subroutine powser(a,b,c,n,alpha,option,d,f)
wilf2:      subroutine netflo(n,e,endpt,source,sink,floval,cut,cap,vert,aux)
wilf2:      subroutine perman(n,a,in,x,perm)
wilf2:      subroutine invert (n,a,ainv)
wilf2:      subroutine triang(n,zeta,sig)
wilf2:      subroutine mobius(n,h,mu,sigma,sig1)
wilf2:      subroutine backtr(l,a,index,k,m,stack,nstk)
wilf2:      subroutine colvrt(n,a,k,m,stack,nstk,lambda,adj,col)
wilf2:      subroutine eulcrc(e,a,k,m,stack,nstk,option,endpt,z1,ed)
wilf2:      subroutine hamcrc(n,a,k,m,stack,nstk,adj,vert,option)
wilf2:      subroutine spntre(e,n,a,k,m,stack,nstk,endpt,end,x,nv,y)
wilf2:      subroutine rantre(n,end,a,b,m)
wilf2:      subroutine ranrut(nn,out,stack,t)
wilf2:      subroutine signum(sigma,n,sign,cycles)
wilf2:      subroutine hpsort(n,b)
wilf2:      subroutine minspt(dist,n,endpt,u,y)
wilf3:      subroutine nexsub(n,in,mtc,ncard,j)
wilf3:      subroutine ransub(n,a)
wilf3:      subroutine nexksb(n,k,a,mtc)
wilf3:      subroutine nxksrd(n,k,a,mtc,in,out)
wilf3:      subroutine ranksb(n,k,a)
wilf3:      subroutine nexcom(n,k,r,mtc)
wilf3:      subroutine rancom(n,k,r)
wilf3:      subroutine nexper(n,a,mtc)
wilf3:      subroutine ranper(n,a)
wilf3:      subroutine nexpar(n,r,m,d,mtc)
wilf3:      subroutine ranpar(n,k,mult,p)
wilf3:      subroutine nexequ(n,nc,p,q,mtc)
wilf3:      subroutine ranequ(n,l,q,a,b,c)
wilf3:      subroutine renumb(m,n,sig,tau,a)
wilf3:      subroutine spanfo(n,e,endpt,k,x,nv,y)
wilf3:      subroutine poly(n,a,x0,option,val,b)
wilf3:      subroutine chromp(n,e,endpt,a,b,c,stack,nstk)
wilf4:      subroutine powser(a,b,c,n,alpha,option,d,f)
wilf4:      subroutine netflo(n,e,endpt,source,sink,floval,cut,cap,vert,aux)
wilf4:      subroutine perman(n,a,in,x,perm)
wilf4:      subroutine invert (n,a,ainv)
wilf4:      subroutine triang(n,zeta,sig)
wilf4:      subroutine mobius(n,h,mu,sigma,sig1)
wilf4:      subroutine backtr(l,a,index,k,m,stack,nstk)
wilf4:      subroutine colvrt(n,a,k,m,stack,nstk,lambda,adj,col)
wilf4:      subroutine eulcrc(e,a,k,m,stack,nstk,option,endpt,z1,ed)
wilf4:      subroutine hamcrc(n,a,k,m,stack,nstk,adj,vert,option)
wilf4:      subroutine spntre(e,n,a,k,m,stack,nstk,endpt,end,x,nv,y)
wilf4:      subroutine rantre(n,end,a,b,m)
wilf4:      subroutine ranrut(nn,out,stack,t)
wilf4:      subroutine signum(sigma,n,sign,cycles)
wilf4:      subroutine hpsort(n,b)
wilf4:      subroutine minspt(dist,n,endpt,u,y)



is this any good to you?

i also have some other comb. programs
in that directory, namely:

2surj3.f   canfm1d    det3d      harw7      nexpar     print.7    sort3
2surj3d    canfm2.f   det4.f     hpsort     nexper     rancom     stein1
3surj1.f   canfm3.f   det4.out1  in         nexper1.f  ranequ     temp1
3surj1d    ceil.f     det4d      in4        nexsub.f   rank1.f    timen
3surj2.f   cliq2.f    det5.f     inver1     nexvec     rank2      unpac2
3surj3.f   cliq3.f    det5d      invert.f   nxksrd     rank3.f    unpac3
3surj3d    colvrt1.f  erb.f      invert3.f  oldvor5    rank4.f    veit1.f
3surj4.f   combn2     floor.f    jin1       out        rank4d     veit2
3surj5.f   count      gcd1       jin2       out5       rankp.f    veit3.f
README     countm     gcd2.f     lexsub     pac2       rankpd     veit4
a.out      denn1      gcd3       mask       pard2.f    ranksb     veit5.f
backtr     denn2      gcd4       maxod1.f   perm1.f    rannosc    wilf1
binco1     denn3      had1       memo       perm2.f    rannosf    wilf2
binry1.f   denn4      had2       memo2      prime.f    ranpar     wilf3
binry2.f   des.c      harw1      memoj      print.1    ranper     wilf4
binry3.f   des2.c     harw2      mifix      print.2    ransub     wu1.f
binry3d    det1.f     harw3      nbinco2    print.3    rdmx1      wu1d
biplane.16 det1d      harw4      nexcom     print.4    reeds.c
canfm0     det2.f     harw5      nexequ     print.5    root2
canfm1.f   det3.f     harw6      nexksb.f   print.6    shift2


this is my main collection of purely comb.
stuff yoiu see.  i mention this in case
there are subsrouytines from NW that i have
extracted from the wilf files


let me know what you woiuld like

neil

SHAR_EOF
cat << 'SHAR_EOF' > backtr.f

ccccccc p.245
      subroutine backtr(l,a,index,k,m,stack,nstk)
c Nijenhuis Wilf p 245. Checked 8.31.90
      implicit integer(a-z)
      dimension a(l),stack(nstk)
10    if(index.ne.0) go to 50
20    k=1
      m=0
30    index=2
      return
50    nc=stack(m)
      m=m-1
60    if(nc.ne.0) go to 100
70    k=k-1
80    if(k.ne.0) go to 50
90    index=3
      return
100   a(k)=stack(m)
      stack(m)=nc-1
110   if(k.ne.l) go to 120
      index=1
      return
120   k=k+1
      go to 30
      end

SHAR_EOF
cat << 'SHAR_EOF' > bincoef.f

      double precision function bincoef(x,k)
c double precision binomial coefficients
      double precision b,x,y
      b=0
      if (k .lt. 0) go to 100
      if (k .ge. 0) b=1
      if (k .eq. 0) go to 100
      do 10 i=1,k
      y=x-i+1
      y=y/i
   10 b=b*y
  100 bincoef=b
      return
      end



SHAR_EOF
cat << 'SHAR_EOF' > det1.f

c det1.f. Computes the determinant of a real matrix a of size n X n.
c seems to work correctly 4/16/86
c To use type f77 det1.f -lport3    then a.out
c or run on cray which has more precision
c  usl "cray!fortran det1.f \
c  11<in \
c  12>out \
c  L=PORT3 \
c  time=02 grade=2 \
c  //+notify// "
c reads data from file in
c reads n in free format
c then on next line the format number to use 
c 1 = 72i1 format
c 2 = 36i2 format
c 4 = 15i4
c 7 = free format floating point
c 8 = free format integers
c then the rows.
c reference L Kaufman's linear equations package, GELU example.
c the determinant of a is given by detman*beta**idetex
c where beta is the base of the machine
c and detman is between 1/beta and 1
c main
      dimension a(100,100),it(100)
      maxn=100
      open(unit=11,file='in')
      rewind(11)
      open(unit=12,file='out')
      rewind(12)
      ibeta=i1mach(10)
    1 read(11,*)n
   99 format(i4)
      if(n.le.maxn)goto 20
      write(12,98)
   98 format(" n too large")
      stop
   20 continue
      write(12,31) n
   31 format(" det1.f. n=",i6)
c read format number
      read(11,*)iform
c read matrix
      write(12,33)
   33 format(" matrix is")
      do 21 i=1,n
      goto (201,202,203,204,205,206,207,208,209),iform
  201 read(11,301)(it(j),j=1,n)
  301 format(72i1)
      write(12,501)(it(j),j=1,n)
  501 format(1h ,72i1)
      goto 401
  202 read(11,302)(it(j),j=1,n)
  302 format(36i2)
      write(12,502)(it(j),j=1,n)
  502 format(1h ,36i2)
      goto 401
  203 continue
  204 read(11,304)(it(j),j=1,n)
  304 format(15i4)
      write(12,504)(it(j),j=1,n)
  504 format(1h ,15i4)
      goto 401
  205 continue
  206 continue
  207 continue
      read(11,*)(a(i,j),j=1,n)
      write(12,*)(a(i,j),j=1,n)
      goto 401
  208 read(11,*)(it(j),j=1,n)
      write(12,504)(it(j),j=1,n)
      goto 401
c toeplitz
 209  continue
 401  continue
      if(iform.eq.7)goto 21
      do 30 j=1,n
   30 a(i,j)=float(it(j))
   21 continue
c now have matrix
c temp
c      do 774 i=1,n
c      write(12,775)(a(i,j),j=1,n)
c 774  continue
 775  format(1h ,12f6.3)
      call det(n,a,maxn,detman,idetex,dett)
      write(12,40)detman,ibeta,idetex
   40 format(" det = ",e14.8,"*",i4,"**",i6)
      write(12,44)dett
   44 format(" in other words, det = ",e14.8)
      stop
      end
      subroutine det(n,a,ia,detman,idetex,dett)
      integer n,ia,idetex
      integer e,ipoint,istkgt,i1mach,isign,i
      integer in(1000)
      real a(ia,n),m,onovbe,mkfl
      double precision d(500)
      common/cstak/d
      equivalence(d(1),in(1))
c allocate space from the stack for the pivot array
      ipoint=istkgt(n,2)
      call gelu(n,a,ia,in(ipoint),1.0e-5)
c the det is the product of the diag entries and the last elt
c of the interchange array. We try to comput this product in a way
c that will avoid underflow etc
      beta=float(i1mach(10))
      onovbe=1.0/beta
      isign=ipoint+n-1
      detman=in(isign)*onovbe
      idetex=1
      do 10 i=1,n
      call umkfl(a(i,i),e,m)
      detman=detman*m
      idetex=idetex+e
      if(abs(detman).ge.onovbe)goto 10
      idetex=idetex-1
      detman=detman*beta
   10 continue
      dett = mkfl(idetex,detman)
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > hpsort.f

ccccccccc P.140
      subroutine hpsort(n,b)
c N & W p 140.
      integer b(500),bstar
      n1=n
      l=1+n/2
   11 l=l-1
      bstar=b(l)
      goto 30
   25 bstar=b(n1)
      b(n1)=b(1)
   29 n1=n1-1
   30 l1=l
   31 m=2*l1
      if(m-n1)32,33,37
   32 if(b(m+1).ge.b(m)) m=m+1
   33 if(bstar.ge.b(m))goto 37
      b(l1)=b(m)
      l1=m
      goto 31
   37 b(l1)=bstar
      if(l.gt.1)goto 11
      if(n1.ge.2)goto 25
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > lexsub.f

cccccc P.19
      subroutine lexsub(n,k,in,jmp,ndim)
c  ref NW p19. Subsets in lex order with extra cheese.
      dimension in(ndim)
      logical jmp
   10 if(k.ne.0)goto 40
   20 if(jmp)return
   30 is=0
  100 if(.not.jmp)k=k+1
  110 in(k)=is+1
      return
   40 if(in(k).eq.n)goto 50
   80 is=in(k)
   90 goto 100
   50 k=k-1
   60 if(k.eq.0)return
      is=in(k)
      goto 110
      end

SHAR_EOF
cat << 'SHAR_EOF' > nexcom.f

cccccc P.49
      subroutine nexcom(n,k,r,mtc)
c Next composition of n into k parts. NW 50.
      integer r(k),t,h
      logical mtc
   10 if(mtc)goto 20
      r(1)=n
      t=n
      h=0
      if(k.eq.1)goto 15
      do 11 i=2,k
   11 r(i)=0
   15 mtc=r(k).ne.n
      return
   20 if(t.gt.1) h=0
   30 h=h+1
      t=r(h)
      r(h)=0
      r(1)=t-1
      r(h+1)=r(h+1)+1
      goto 15
      end

SHAR_EOF
cat << 'SHAR_EOF' > nexequ.f

ccccc P.90
      subroutine nexequ(n,nc,p,q,mtc)
c next partition of an n-set. Ref NW 91.
      logical mtc
      integer p(n),q(n)
      if(mtc)goto 20
   10 nc=1
      do 11 i=1,n
   11 q(i)=1
      p(1)=n
   60 mtc=nc.ne.n
      return
   20 m=n
   30 l=q(m)
   90 if(p(l).ne.1) go to 40
   95 q(m)=1
      m=m-1
      go to 30
   40 nc=nc+m-n
      p(1)=p(1)+n-m
  110 if(l.ne.nc) go to 50
  120 nc=nc+1
      p(nc)=0
   50 q(m)=l+1
      p(l)=p(l)-1
      p(l+1)=p(l+1)+1
      go to 60
      end

SHAR_EOF
cat << 'SHAR_EOF' > nexksb.f

c next k-subset Ref NW p 19?
c Use of nexsub
      logical mtc
      dimension in(5)
      n=5
	k=2
      mtc=.false.
   10 call nexksb(n,k,in,mtc)
      write(*,20) (in(i),i=1,n)
   20 format(1h ,39i2)
      if(mtc) goto 10
      stop
      end

cccccc P.32
      subroutine nexksb(n,k,a,mtc)
      integer a(k),h
      logical mtc
      if(k.gt.0)goto 30
      do 1 i=1,n
 1    a(i)=0
      mtc=.false.
      return
   30 if(mtc) go to 40
   20 m2=0
      h=k
      go to 50
   40 if(m2.lt.n-h) h=0
      h=h+1
      m2=a(k+1-h)
   50 do 51 j=1,h
   51 a(k+j-h)=m2+j
      mtc=a(1).ne.n-k+1
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > nexpar.f

ccccc P.69
      subroutine nexpar(n,r,m,d,mtc)
c  next partition of n. Ref NW 69.
      implicit integer(a-z)
      logical mtc
      dimension m(n),r(n)
      data nlast/0/
   10 if(n.eq.nlast) go to 20
      nlast=n
   30 s=n
      d=0
   50 d=d+1
      r(d)=s
      m(d)=1
   40 mtc=m(d).ne.n
      return
   20 if(.not.mtc) go to 30
      sum=1
      if(r(d).gt.1) go to 60
      sum=m(d)+1
      d=d-1
   60 f=r(d)-1
      if(m(d).eq.1) go to 70
      m(d)=m(d)-1
      d=d+1
   70 r(d)=f
      m(d)=1+sum/f
      s=mod(sum,f)
      if(s) 40,40,50
      end

SHAR_EOF
cat << 'SHAR_EOF' > nexper.f

      integer a(40)
      logical mtc,even
      n=4
      mtc=.false.
 10   continue
      call nexper(n,a,mtc,even)
      write(06,100)(a(i),i=1,n)
 100  format(16i4)
      if(mtc)goto 10
      write(06,*)"all done"
      stop
      end

cccccc P. 59
      subroutine nexper(n,a,mtc,even)
c next permutation of {1,...,n}. Ref NW p 59.
      integer a(n),s,d
      logical mtc,even
      if(mtc)goto 10
      nm3=n-3
      do 1 i=1,n
    1 a(i)=i
      mtc=.true.
    5 even=.true.
      if(n.eq.1)goto 8
    6 if(a(n).ne.1.or.a(1).ne.2+mod(n,2))return
      if(n.le.3)goto 8
      do 7 i=1,nm3
      if(a(i+1).ne.a(i)+1)return
    7 continue
    8 mtc=.false.
      return
   10 if(n.eq.1)goto 27
      if(.not.even)goto 20
      ia=a(1)
      a(1)=a(2)
      a(2)=ia
      even=.false.
      goto 6
   20 s=0
      do 26 i1=2,n
   25 ia=a(i1)
      i=i1-1
      d=0
      do 30 j=1,i
   30 if(a(j).gt.ia) d=d+1
      s=d+s
      if(d.ne.i*mod(s,2)) goto 35
   26 continue
   27 a(1)=0
      goto 8
   35 m=mod(s+1,2)*(n+1)
      do 40 j=1,i
      if(isign(1,a(j)-ia).eq.isign(1,a(j)-m))goto 40
      m=a(j)
      l=j
   40 continue
      a(l)=ia
      a(i1)=m
      even=.true.
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > nexsub.f
From njas@research.att.com Fri May 26 08:09:14 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29905 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:09:13 -0400
From: njas@research.att.com
Message-Id: <199505261209.IAA29905@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:08 EDT
To: skiena@dimacs.rutgers.edu
Status: RO

c next subset Ref NW p 19.
c Use of nexsub
      logical mtc
      dimension in(5)
      n=5
      mtc=.false.
   10 call nexsub(n,in,mtc,ncard,j)
      print 20,(in(i),i=1,n),j,ncard
   20 format(1h ,39i2)
      if(mtc) goto 10
      stop
      end

ccccccc P.18
      subroutine nexsub(n,in,mtc,ncard,j)
      logical mtc
      dimension in(n)
      if(mtc)goto 20
      do 11 i=1,n
   11 in(i)=0
      ncard=0
      mtc=.true.
      return
   20 j=1
      if(mod(ncard,2).eq.0)goto 30
   40 j=j+1
      if(in(j-1).eq.0)goto 40
c if(j.gt.n)j=n
   30 in(j)=1-in(j)
      ncard=ncard+2*in(j)-1
      mtc=ncard.ne.in(n)
      return
      end

From njas@research.att.com Fri May 26 08:10:20 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29966 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:19 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29966@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO

      subroutine nexsub(n,in,mtc,ncard,j)
      logical mtc
      dimension in(n)
      data nlast/0/
10    if(n.eq.nlast) go to 15
20    m=0
      mtc=.true.
      do 21  i=1,n
21    in(i)=0
      ncard=0
      nlast=n
      return
15    if(.not.mtc) go to 20
30    m=m+1
      m1=m
      j=0
39    j=j+1
40    if(mod(m1,2).eq.1) go to 60
50    m1=m1/2
      go to 39
60    l=in(j)
      in(j)=1-l
      ncard=ncard+1-2*l
      mtc=ncard.ne.1.or.in(n).eq.0
      return
      end
      subroutine ransub(n,a)
      integer a(n)
      do 10  i=1,n
10    a(i)=2.*uni(1)
      return
      end

cccccccc P.32
      subroutine nexksb(n,k,a,mtc)
      integer a(k),h
      logical mtc
      data nlast/0/,klast/0/
10    if(k.ne.klast.or.n.ne.nlast)go to 20
30    if(mtc) go to 40
20    m2=0
      h=k
      nlast=n
      klast=k
      mtc=.true.
      go to 50
40    do 41  h=1,k
      i=k+1-h
      m2=a(i)
      if(m2.ne.n+1-h) go to 50
41    continue
50    do 51 j=1,h
      i=k+j-h
51    a(i)=m2+j
      mtc=(a(1).ne.n-k+1)
      return
      end

cccccc P.33
      subroutine nxksrd(n,k,a,mtc,in,out)
      integer a(k),out
      logical mtc
      data nlast,klast/0,0/
5     if(n.eq.nlast.and.k.eq.klast) go to 20
      nlast=n
      klast=k
6     mtc=.true.
10    do 11  i=1,k
11    a(i)=i
      in=0
      out=0
      go to 120
20    if(.not.mtc)  go to 6
25    if(mod(k,2).eq.0) go to 29
26    if(k.eq.1.or.a(2)-a(1).gt.1)  go to 100
27    if(a(1).gt.1) go to 110
28    km1=k-1
      do 40  n1=1,km1
      if(a(n1+1).gt.a(n1)+1)  go to 41
40    continue
      n1=k
41    if(mod(n1+k,2).eq.0)  go to 60
30    if(n1.eq.km1.or.(n1.lt.km1.and.a(n1+2).gt.a(n1+1)+1)) go to 90
      go to 80
29    if(a(1)-1) 29,28,70
60    out=a(n1-1)
      in=out+2
      a(n1-1)=a(n1)
      a(n1)=out+2
      go to 120
70    in=a(1)-1
      out=a(1)
      a(1)=in
      go to 120
80    in=a(n1)+1
      out=a(n1+2)
      a(n1+2)=a(n1+1)
      a(n1+1)=in
      go to 120
90    in=a(n1+1)+1
      out=a(n1)
      a(n1)=in-1
      a(n1+1)=in
      go to 120
100   in=a(1)+1
      out=a(1)
      a(1)=in
      go to 120
110   in=1
      out=a(2)
      a(2)=a(1)
      a(1)=1
120   mtc=a(k).ne.n.or.(k.ne.1.and.a(k-1).ne.k-1)
      return
      end

ccccc P.43
      subroutine ranksb(n,k,a)
      integer r,a(k)
      l2=k+1
      r=l2
      do 15 i=1,k
15    a(i)=0
      go to 20
45    a(i)=a(i)+r
      i=r
50    a(i)=l2*m
20    m=1.0+uni(1)*float(n)
      i=1+mod(m-1,k)
      if (a(i).eq.0) go to 50
30    if (m.eq.a(i)/l2) go to 20
      link=mod(a(i),l2)
      if (link.eq.0) go to 40
      i=link
      go to 30
40    r=r-1
      if (r.eq.0) go to 55
      if (a(r)) 45,45,40
55    do 60  i=1,k
60    a(i)=a(i)/l2
      return
      end

ccccc P.49
      subroutine nexcom(n,k,r,mtc)
      integer r(k),t
      logical mtc
      data klast,nlast/0,0/
10    if (n.eq.nlast.and.k.eq.klast) go to 60
20    nlast=n
      klast=k
      do 21  i=1,k
21    r(i)=0
      r(1)=n
30    mtc=(r(k).ne.n)
      return
60    if(.not.mtc) go to 20
70    do 71  i=1,k
      if(r(i).ne.0) go to 100
71    continue
100   t=r(i)
      r(i)=0
      r(1)=t-1
      r(i+1)=r(i+1)+1
      go to 30
      end

cccccccc P.52
      subroutine rancom(n,k,r)
      integer r(k)
      call ranksb(n+k-1,k-1,r)
      call hpsort(k-1,r)
      r(k)=n+k
      l=0
      do 10 i=1,k
      m=r(i)
      r(i)=m-l-1
10    l=m
      return
      end

cccccccc P.59
      subroutine nexper(n,a,mtc)
      integer a(n),b,h1,h,t,v
      logical mtc
      data nlast/0/
10    if(n.eq.nlast) go to 20
30    nlast=n
      m=1
      v=1
      nf=1
      do 31  j=1,n
      nf=nf*j
31    a(j)=j
40    mtc=(m.ne.nf)
      return
20    if(.not.mtc) go to 30
      go to (70,80),v
70    t=a(2)
      a(2)=a(1)
      a(1)=t
      v=2
      m=m+1
      go to 40
80    h=3
      m1=m/2
90    b=mod(m1,h)
100   if(b.ne.0) go to 120
110   m1=m1/h
      h=h+1
      go to 90
120   m1=n
      h1=h-1
      do 160  j=1,h1
130   m2=a(j)-a(h)
      if(m2.lt.0)  m2=m2+n
140   if(m2.ge.m1) go to 160
150   m1=m2
      j1=j
160   continue
180   t=a(h)
      a(h)=a(j1)
      a(j1)=t
      v=1
      m=m+1
      return
      end

ccccc P.63
      subroutine ranper(n,a)
      integer a(n)
      do 10  i=1,n
10    a(i)=i
20    do 40  m=1,n
30    l=float(m)+uni(1)*float(n+1-m)
      l1=a(l)
      a(l)=a(m)
40    a(m)=l1
      return
      end

ccccc P.69
      subroutine nexpar(n,r,m,d,mtc)
      integer d,f,r(n),s,sum
      logical mtc
      dimension m(n)
      data nlast/0/
10    if(n.eq.nlast) go to 20
      nlast=n
30    s=n
      d=0
50    d=d+1
      r(d)=s
      m(d)=1
40    mtc=m(d).ne.n
      return
20    if(.not.mtc) go to 30
      sum=1
      if(r(d).gt.1) go to 60
      sum=m(d)+1
      d=d-1
60    f=r(d)-1
      if(m(d).eq.1) go to 70
      m(d)=m(d)-1
      d=d+1
70    r(d)=f
      m(d)=1+sum/f
      s=mod(sum,f)
      if(s) 40,40,50
      end

cccccc P.75
      subroutine ranpar(n,k,mult,p)
      integer p(n),d,mult(n)
      data nlast/0/
10    if(n.le.nlast) go to 30
20    p(1)=1
      m=nlast+1
      nlast=n
      if(n.eq.1) go to 30
      do 21  i=m,n
      isum=0
26    do 22  d=1,i
      is=0
      i1=i
24    i1=i1-d
      if(i1) 22,25,23
23    is=is+p(i1)
      go to 24
25    is=is+1
22    isum=isum+is*d
21    p(i)=isum/i
30    m=n
      k=0
      do 31  i=1,n
31    mult(i)=0
40    z=uni(i)*float(m*p(m))
      d=0
110   d=d+1
60    i1=m
      j=0
150   j=j+1
70    i1=i1-d
80    if(i1)  110,90,120
120   z=z-float(d*p(i1))
130   if(z) 145,145,150
90    z=z-float(d)
100   if(z) 145,145,110
145   mult(d)=mult(d)+j
      k=k+j
160   m=i1
170   if(m.ne.0) go to 40
      return
      end

ccccccc P. 90
      subroutine nexequ(n,nc,p,q,mtc)
      logical mtc
      integer p(n),q(n)
      data nlast/0/
10    if(n.eq.nlast) go to 20
30    nlast=n
      nc=1
      do 35  i=1,n
35    q(i)=1
      p(1)=n
40    mtc=(nc.ne.n)
      return
20    if(.not.mtc)  go to 30
70    m=n
80    l=q(m)
90    if(p(l).ne.1) go to 100
95    q(m)=1
      m=m-1
      go to 80
100   nc=nc+m-n
      p(1)=p(1)+n-m
110   if(l.ne.nc) go to 130
120   nc=nc+1
      p(nc)=0
130   q(m)=l+1
      p(l)=p(l)-1
      p(l+1)=p(l+1)+1
      go to 40
      end


ccccc P.97
      subroutine ranequ(n,l,q,a,b,c)
      integer q(n),a(n),c(n)
      dimension b(n)
      data nlast/1/
      b(1)=1
      if(n.le.nlast) go to 10
3     m=nlast
      nlast=n
      nm1=n-1
      do 5  l=m,nm1
      sum=1./float(l)
      l1=l-1
      do 6  k=1,l1
6     sum=(sum+b(k))/float(l-k)
5     b(l+1)=(sum+b(l))/float(l+1)
10    do 11  i=1,n
11    q(i)=0
      l=0
      m=n
20    z1=uni(1)
      k=m-1
      t=1.0/float(m)
      m1=m-1
60    if(k.eq.0) go to 70
30    z=t*b(k)/b(m)
40    if(z1.lt.z) go to 80
50    k=k-1
      t=t/float(m-1-k)
      z1=z1-z
      go to 60
80    l1=n
81    if(q(l1).eq.0) go to 82
      l1=l1-1
      go to 81
82    l=l+1
      q(l1)=l
      call ranksb(m1,k,a)
      m2=1
      i=1
90    if(q(i).eq.0) go to 110
100   if(i.eq.m) go to 130
120   i=i+1
      go to 90
110   c(m2)=i
      m2=m2+1
      q(i)=l
      go to 100
130   do 131  i=1,k
      j=a(i)
      j=c(j)
131   q(j)=0
      m=k
      go to 20
70    l=l+1
      do 71  i=1,n
71    if(q(i).eq.0)  q(i)=l
      return
      end

ccccccc P.155
      subroutine renumb(m,n,sig,tau,a)
      integer sig(m),tau(n),a(m,n),t1,t2
      do 5  i=1,m
      i1=sig(i)
6     if(i1.le.i) go to 5
      i2=sig(i1)
      sig(i1)=-i2
      i1=i2
      go to 6
5     sig(i)=-sig(i)
      if(tau(1).lt.0) go to 9
      do 7  j=1,n
      j1=tau(j)
8     if(j1.le.j) go to 7
      j2=tau(j1)
      tau(j1)=-j2
      j1=j2
      go to 8
7     tau(j)=-tau(j)
9     do 10  i=1,m
      i1=-sig(i)
      if(i1.lt.0) go to 10
      lc=0
20    i1=sig(i1)
      lc=lc+1
      if(i1.gt.0) go to 20
      i1=i
      do 30  j=1,n
      if(tau(j).gt.0) go to 30
      j2=j
      k=lc
40    j1=j2
      t1=a(i1,j1)
50    i1=iabs(sig(i1))
      j1=iabs(tau(j1))
      t2=a(i1,j1)
      a(i1,j1)=t1
      t1=t2
      if(j1.ne.j2) go to 50
      k=k-1
      if(i1.ne.i) go to 50
      j2=iabs(tau(j2))
55    if(k.ne.0) go to 40
30    continue
10    continue
      do 60  i=1,m
60    sig(i)=iabs(sig(i))
      if(tau(1).gt.0) return
      do 70  j=1,n
70    tau(j)=iabs(tau(j))
      return
      end

ccccccc P.167
      subroutine spanfo(n,e,endpt,k,x,nv,y)
      integer e1,e,endpt,s,t1,t2,t,v1,v2,x,y,z
      dimension endpt(2,e),x(n),nv(n),y(n),s(2)
      data s(1),s(2),m/1,2,2/
      do 10  i=1,n
      x(i)=-i
      nv(i)=1
10    y(i)=0
      j=1
      e1=e
20    v1=endpt(1,j)
      v2=endpt(2,j)
25    t1=x(v1)
      if(t1.lt.0)  t1=v1
      t2=x(v2)
      if(t2.lt.0)  t2=v2
      if(t1.ne.t2) go to 40
      if(j.lt.e1)  go to 30
      e1=e1-1
      go to 60
30    endpt(1,j)=endpt(1,e1)
      endpt(2,j)=endpt(2,e1)
      endpt(1,e1)=v1
      endpt(2,e1)=v2
      e1=e1-1
      go to 20
40    if(nv(t1).le.nv(t2))  go to 50
      t=t1
      t1=t2
      t2=t
      i3=-x(t2)
50    y(i3)=t1
      x(t2)=x(t1)
      i=t1
55    x(i)=t2
      i=y(i)
      if(i.ne.0) go to 55
      nv(t2)=nv(t2)+nv(t1)
      nv(t1)=0
      j=j+1
      if(j.le.e1.and.j.lt.n)  go to 20
60    k=0
      do 70  i=1,n
      if(nv(i).eq.0)  go to 70
      k=k+1
      nv(k)=nv(i)
      y(i)=k
70    continue
      do 80  i=1,n
      t=x(i)
      if(t.lt.0)  t=i
80    x(i)=y(t)
      if(k.eq.1)  return
90    i2=nv(1)
      nv(1)=1
      do 100  l=2,k
      i1=nv(l)
      nv(l)=nv(l-1)+i2-1
100   i2=i1
      do 110  i=1,e1
      i3=endpt(1,i)
      z=x(i3)
      y(i)=nv(z)
110   nv(z)=nv(z)+1
      call renumb(m,e1,s,y,endpt)
      i1=1
      do 120  l=1,k
      i2=nv(l)
      nv(l)=i2-i1+1
120   i1=i2
      return
      end

ccccccc P.175
      subroutine poly(n,a,x0,option,val,b)
      integer a,b,option,v,val,x0,z
      dimension a(n),b(n)
      val=a(n)
      if(n.eq.1) return
      n1=n-1
      if(option.eq.0) go to 20
      if(option.eq.(-1)) go to 26
      do 10 i=1,n
10    b(i)=a(i)
      if(option) 50,30,30
20    do 25 i=1,n1
      i1=n-i
25    val=val*x0+a(i1)
      return
26    do 27 i=1,n1
      i1=n-i
27    val=val*(x0-n1+i)+a(i1)
      return
30    max=min0(n1,option)
      do 35 j=1,max
      m=n1
      v=val
37    v=b(m)+v*x0
      b(m)=v
      m=m-1
      if(m.ge.j) go to 37
35    continue
      val=b(1)
      return
50    if(n.eq.2) return
      n2=n-2
      if(option.eq.(-3)) go to 70
      do 55  j=1,n2
      v=val
      m=n1
60    v=b(m)+j*v
      b(m)=v
      m=m-1
      if(m.gt.j) go to 60
55    continue
      return
70    do 75  j=1,n2
      z=n1-j
      m=z+1
80    b(m)=b(m)-z*b(m+1)
      m=m+1
      if(m.le.n1) go to 80
75    continue
      return
      end


cccccc P.183
      subroutine chromp(n,e,endpt,a,b,c,stack,nstk)
      integer a,ai,b,bj,c,e1,e,ed,en1,en2,en,end,endpt
      integer oth,r,stack,tr,v
      dimension endpt(2,e),stack(2,nstk),a(n),b(n),c(n)
10    call spanfo(n,e,endpt,k,a,b,c)
      if (k.gt.1) stop
12    do 13  i=1,n
      a(i)=0
13    b(i)=1
20    n1=n-1
      do 90  j=1,n1
30    do 40  l=1,2
50    i0=endpt(l,j)
      j0=0
60    j1=a(i0)
      a(i0)=j0
70    if(j1.eq.0) go to 40
80    j0=j1
      bj=b(j0)
      b(j0)=-bj
      i3=(3-bj)/2
      i0=endpt(i3,j0)
      go to 60
40    continue
      i3=endpt(2,j)
90    a(i3)=j
      r=endpt(1,n1)
      do 140  j=1,n1
      do 140  l=1,2
      i3=(3+b(j)*(2*l-3))/2
140   stack(l,j)=endpt(i3,j)
      do 170  i=1,n
      b(i)=0
      if(a(i).eq.0) go to 170
      k=0
      i1=i
150   ai=a(i1)
      if(ai.eq.0) go to 155
      a(i1)=0
      k=k+1
      c(k)=ai
      i1=stack(1,ai)
      go to 150
155   i1=i
      j0=j0+k
      do 160  j=1,k
      do 160  l=1,2
      i3=j0-j+1
      i4=c(j)
160   endpt(l,i3)=stack(l,i4)
170   continue
200   n1=n
      e1=e
      is=0
210   if(e1+1.ne.n1) go to 280
220   b(n1)=b(n1)+1
230   if(is.eq.0) go to 600
240   n1=stack(1,is)
      e1=stack(2,is)
      j0=is-e1
      do 241  j=1,e1
      do 242  l=1,2
242   endpt(l,j)=stack(l,j0)
241   j0=j0+1
250   if(e1.ne.n1) go to 260
270   b(n1)=b(n1)+1
      is=is-e1-1
      go to 310
280   if(e1.ne.n1) go to 300
290   b(n1)=b(n1)+1
      go to 310
300   do 301  j=1,e1
      is=is+1
      do 301  l=1,2
301   stack(l,is)=endpt(l,j)
      stack(1,is)=n1
      stack(2,is)=e1-1
      go to 310
260   is=is-1
      stack(1,is)=n1
      stack(2,is)=e1-1
310   do 311  i=1,n
311   a(i)=0
320   i1=endpt(1,e1)
      i2=endpt(2,e1)
      e1=e1-1
      tr=n1-1
      j=0
      if(i1.eq.r) go to 327
      if(i2.eq.r) go to 325
      do 321  j=1,tr
      en2=endpt(2,j)
322   if(en2.eq.i2) go to 325
323   if(en2.eq.i1) go to 326
321   continue
325   i0=i1
      i1=i2
      i2=i0
326   if(j.le.0) go to 327
      i3=endpt(1,j)
      a(i3)=1
327   ed=j
      j1=j+1
330   do 340  j=j1,tr
350   en1=endpt(1,j)
      en2=endpt(2,j)
360   if(en2.ne.i2) go to 400
370   end=en1
      go to 340
400   if(en1.ne.i2) go to 420
      en1=i1
      go to 430
420   if(en1.ne.i1) go to 440
430   a(en2)=1
440   ed=ed+1
      endpt(1,ed)=en1
      endpt(2,ed)=en2
340   continue
450   n1=tr
      endpt(1,n1)=i1
      endpt(2,n1)=end
460   do 470  j=n1,e1
480   do 485  l=1,2
490   en=endpt(l,j)
500   if(en.eq.i1) go to 530
510   if(en.eq.i2) go to 520
485   continue
555   oth=endpt(1,j)
      go to 560
520   en=i1
530   i3=3-l
      oth=endpt(i3,j)
540   if(a(oth).eq.1) go to 470
550   a(oth)=1
560   ed=ed+1
      endpt(1,ed)=en
      endpt(2,ed)=oth
470   continue
570   e1=ed
      go to 210
600   do 601  i=1,n
601   a(i)=(1-2*mod(n-i,2))*b(i)
      call poly(n,a,0,-2,v,c)
      call poly(n,a,-1,n,v,a)
      do 602  i=1,n
602   a(i)=iabs(a(i))
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > nexvec.f
From njas@research.att.com Fri May 26 08:09:15 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29910 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:09:15 -0400
From: njas@research.att.com
Message-Id: <199505261209.IAA29910@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO

      dimension nvec(20)
      logical mtc
      do 1 n=1,5
      print,n
      mtc=.false.
 2    call nexvec(n,nvec,mtc)
      if(.not.mtc)goto 1
      print 4,(nvec(i),i=1,n)
 4    format(5h      ,30i1)
      goto 2
 1    continue
      stop; end

      subroutine nexvec(n,nvec,mtc)
c returns next binary vector of length n in nvec
      dimension nvec(n)
      logical mtc
      if(mtc)goto 3
      mtc=.true.
      do 4 i=1,n
 4    nvec(i)=0
      return
 3    np=1
 2    if(nvec(np).eq.1)goto 1
      nvec(np)=1
      return
 1    nvec(np)=0
      np=np+1
      if(np.le.n)goto 2
      mtc=.false.
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > nxksrd.f
From njas@research.att.com Fri May 26 08:10:14 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29946 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:13 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29946@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO

 integer a(20),out
 logical mtc
 10  mtc=.false.
   read,n,k
 if(n.lt.0)stop
 1 call nxksrd(n,k,a,mtc,in,out)
 print 2,(a(i),i=1,k),in,out
 2 format(1h ,25i5)
 if(mtc)goto 1
 goto 10
 stop
 end

cccccc P.33
      subroutine nxksrd(n,k,a,mtc,in,out)
c next k-subset of n-set, in revolving door order. Ref NW p34.
      integer a(k),out
      logical mtc
      if(mtc) goto 10
      do 1 i=1,k
    1 a(i)=i
      mtc=k.ne.n
      return
   10 j=0
   20 if(mod(k,2).ne.0)goto 100
   30 j=j+1
c 40 if(j.le.k)goto 60
c 50 a(k)=k
c in=k
c out=n
c return
   60 if(a(j).eq.j)goto 100
   70 out=a(j)
      in=out-1
      a(j)=in
   80 if(j.eq.1)goto 200
   90 in=j-1
      a(j-1)=in
      goto 200
  100 j=j+1
  130 m=n
  110 if(j.lt.k)m=a(j+1)-1
  140 if(m.eq.a(j)) goto 30
  150 in=a(j)+1
      a(j)=in
      out=in-1
      if(j.eq.1)goto 200
      a(j-1)=out
      out=j-1
  200 if(k.eq.1)goto 201
      mtc=a(k-1).eq.k-1
  201 mtc=(.not.mtc).or.a(k).ne.n
      return
c200 return
      end

SHAR_EOF
cat << 'SHAR_EOF' > pawser.f
From njas@research.att.com Fri May 26 08:10:18 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29961 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:17 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29961@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO

cccccc P.191
      subroutine powser(a,b,c,n,alpha,option,d,f)
      integer q,option
      dimension a(n),b(n),c(n),d(n),f(n)
      go to (90,100,3,7),option
3     n2=n
      j3=1
      go to 10
7     do 8  i=2,n
8     c(i)=0
      c(1)=a(1)/b(1)
      n2=2
      j3=2
10    do 11  l=1,n2
11    d(l)=b(1)*c(l)
20    do 21  q=1,n2
      if(c(q).ne.0.0) go to 30
21    continue
30    n1=n2/q
40    do 41  m=2,n1
      if(b(m).eq.0.0) go to 41
      m1=m*q
      f(m1)=c(q)**m
      j=m1
50    d(j)=d(j)+b(m)*f(j)
      if(j.eq.n2) go to 41
      j=j+1
      s=0
      j1=j-1
60    do 61  mu=m1,j1
      i1=m*(j+q-mu)-mu
      i2=j+q-mu
61    s=s+float(i1)*f(mu)*c(i2)
      f(j)=s/(c(q)*float(j-m1))
      go to 50
41    continue
      go to (70,80),j3
70    do 71  i=1,n
71    a(i)=d(i)
      return
80    c(n2)=(a(n2)-d(n2))/b(1)
81    if(n2.eq.n)   return
82    n2=n2+1
      go to 10
90    a(1)=alpha*c(1)
      do 91  j=2,n
      j1=j-1
      s=0
      do 95  mu=1,j1
      i1=j-mu
95    s=s+a(mu)*c(i1)*(alpha*float(i1)-float(mu))
91    a(j)=alpha*c(j)+s/float(j)
      return
100   a(1)=c(1)
      do 105  j=2,n
      j1=j-1
      s=0
      do 110  mu=1,j1
      i1=j-mu
110   s=s+a(mu)*c(i1)*float(i1)
105   a(j)=c(j)+s/float(j)
      return
      end

cccccc P.209
      subroutine netflo(n,e,endpt,source,sink,floval,cut,cap,vert,aux)
      integer aux,c,cap,cut,delta,e,endpt,floval
      integer p,q,rd,sink,source,vert,wr
      dimension endpt(4,e),cut(n),cap(n,n),vert(n,n),aux(n)
c
c***   initialization
c
      do 1 i=1,n
      do 1 j=1,n
1     cap(i,j)=0
      do 2 i=1,e
      i1=endpt(1,i)
      i2=endpt(2,i)
2     cap(i1,i2)=endpt(3,i)
      do 3 i=1,n
      k=0
      do 4 j=1,n
      if(cap(i,j)+cap(j,i).eq.0) go to 4
      k=k+1
      vert(i,k)=j
4     continue
3     vert(i,n)=k
      nmin=-n-1
      floval=0
c
c***   scanning and labeling
c
10    lblsnk=n
      do 14 i=1,n
14    cut(i)=nmin
      rd=0
      wr=0
      p=source
      label=-1
15    m=vert(p,n)
      i=1
20    if(i.gt.m) go to 25
      q=vert(p,i)
      if(cap(p,q).eq.0) go to 23
      if(q.eq.sink) lblsnk=-label
      if(cut(q)-label) 21,22,23
21    cut(q)=label
      wr=wr+1
      aux(wr)=q
22    i=i+1
      go to 20
23    vert(p,i)=vert(p,m)
      vert(p,m)=q
      m=m-1
      go to 20
25    cut(p)=m
      rd=rd+1
      if(rd.gt.wr) go to 50
      p=aux(rd)
      if(cut(p)+lblsnk.eq.0) go to 30
      label=cut(p)-1
      go to 15
c
c***   construction of path from source to sink
c
30    q=source
      k=0
35    k=k+1
      aux(k)=q
      if(k.gt.lblsnk) go to 42
40    p=aux(k)
41    m=cut(p)
      if(m.eq.0) go to 43
      q=vert(p,m)
      go to 35
43    k=k-1
      if(k.eq.0) go to 10
      p=aux(k)
      cut(p)=cut(p)-1
      go to 41
42    if(q.ne.sink) go to 43
      delta=cap(p,q)
      do 45 i=2,k
      i1=aux(i-1)
      i2=aux(i)
45    delta=min0(delta,cap(i1,i2))
46    k=k-1
      if(k.eq.0) go to 47
      p=aux(k)
      c=cap(p,q)-delta
      if(c.gt.0) go to 48
      cut(p)=cut(p)-1
      k0=k
48    cap(p,q)=c
      cap(q,p)=cap(q,p)+delta
      q=p
      go to 46
47    floval=floval+delta
      k=k0
      go to 40
c
c***   exit procedure
c
50    do 51  i=1,n
51    cut(i)=min0(1,cut(i)-nmin)
      do 52  i=1,e
      i1=endpt(1,i)
      i2=endpt(2,i)
52    endpt(4,i)=endpt(3,i)-cap(i1,i2)
      return
      end



ccccccc P.224
      subroutine perman(n,a,in,x,perm)
      double precision a,p,perm,prod,sgn,sum,x1,x,z
      logical mtc
      dimension a(n,n),in(n),x(n)
10    p=0
      n1=n-1
      do 11  i=1,n
      sum=0
      do 15  j=1,n
15    sum=sum+a(i,j)
11    x(i)=a(i,n)-sum/2.d0
      sgn=-1
20    sgn=-sgn
      prod=sgn
30    call nexsub(n1,in,mtc,ncard,j)
      if(ncard.eq.0) go to 38
      z=2*in(j)-1
      do 35  i=1,n
35    x(i)=x(i)+z*a(i,j)
38    do 39  i=1,n
39    prod=prod*x(i)
      p=p+prod
      if(mtc) go to 20
40    x1=2*mod(n,2)-1
      perm=2.0*x1*p
      return
      end



cccccc P.227
      subroutine invert (n,a,ainv)
      integer a,ainv,sum
      dimension a(n,n),ainv(n,n)
      j=n
10    i=n
20    sum=0
      if (i.eq.j) sum=1
      k=i+1
25    if (k.gt.j) go to 30
      sum=sum-a(i,k)*ainv(k,j)
      k=k+1
      go to 25
30    ainv (i,j)=sum
      i=i-1
      if (i.gt.0) go to 20
      j=j-1
      if (j.gt.0) go to 10
      return
      end



ccccccc P.230
      subroutine triang(n,zeta,sig)
      integer q,r,sig,t,zeta
      dimension sig(n), zeta(n,n)
10    m=0
      l=0
      do 11 i=1,n
11    sig(i)=0
20    m=m+1
30    if (sig(m).eq.0) go to 40
130   if (m.eq.n) return
      go to 20
40    t=m+1
50    r=t
60    if (r.gt.n) go to 100
70    if (sig(r).ne.0.or.zeta(r,m).eq.0) go to 90
80    sig(r)=m
      m=r
      go to 50
90    r=r+1
      go to 60
100   l=l+1
      q=sig(m)
      sig(m)=l
110   if (q.eq.0) go to 130
      r=m+1
120   m=q
      go to 60
      end



cccccc P.237
      subroutine mobius(n,h,mu,sigma,sig1)
      integer h,sig1,sigma
      dimension h(n,n),sigma(n),mu(n,n),sig1(n)
      call triang(n,h,sigma)
      do 1 i=1,n
      do 1 j=1,n
1     mu(i,j)=h(i,j)
      call renumb(n,n,sigma,sigma,mu)
      n1=n-1
      do 11  i=1,n1
      j1=i+1
      do 11  j=j1,n
11    mu(i,j)=-mu(i,j)
      call invert(n,mu,mu)
      do 12  i=1,n
      do 12  j=i,n
12    if(mu(i,j).ne.0)mu(i,j)=1
      call invert(n,mu,mu)
      do 20  i=1,n
      i1=sigma(i)
20    sig1(i1)=i
      call renumb(n,n,sig1,sig1,mu)
      return
      end


ccccccc P.245
      subroutine backtr(l,a,index,k,m,stack,nstk)
      integer a,stack
      dimension a(l),stack(nstk)
10    if(index.ne.0) go to 50
20    k=1
      m=0
30    index=2
      return
50    nc=stack(m)
      m=m-1
60    if(nc.ne.0) go to 100
70    k=k-1
80    if(k.ne.0) go to 50
90    index=3
      return
100   a(k)=stack(m)
      stack(m)=nc-1
110   if(k.ne.l) go to 120
      index=1
      return
120   k=k+1
      go to 30
      end


cccccccc P.247
      subroutine colvrt(n,a,k,m,stack,nstk,lambda,adj,col)
      integer a,stack
      logical adj,col
      dimension a(n),stack(nstk),adj(n,n),col(n)
      if(k.gt.1) go to 10
      stack(1)=1
      stack(2)=1
      m=2
      return
10    k1=k-1
      do 20  i=1,lambda
20    col(i)=.true.
      do 30  i=1,k1
      i1=a(i)
30    if(adj(i,k)) col(i1)=.false.
      m1=m
      do 40  i=1,lambda
      if(.not.col(i)) go to 40
      m1=m1+1
      stack(m1)=i
40    continue
      stack(m1+1)=m1-m
      m=m1+1
      return
      end



ccccccc P.250
      subroutine eulcrc(e,a,k,m,stack,nstk,option,endpt,z1,ed)
      integer a,e,endpt,option,stack,t,z1
      logical ed(e)
      dimension a(e),stack(nstk),endpt(2,e),z1(e)
10    if(k.ne.1) go to 30
20    z1(1)=endpt(2,1)
      stack(1)=1
      stack(2)=1
      m=2
      return
30    if(k.eq.2) go to 60
40    i1=a(k-1)
      z1(k-1)=endpt(1,i1)+endpt(2,i1)-z1(k-2)
60    t=z1(k-1)
      if(option.eq.2) go to 80
61    do 62  i=1,e
62    ed(i)=t.eq.endpt(1,t).or.t.eq.endpt(2,i)
64    k1=k-1
65    do 66  i=1,k1
      i1=a(i)
66    ed(i1)=.false.
70    m1=m
      do 71  i=1,e
      if(.not.ed(i)) go to 71
      m1=m1+1
      stack(m1)=i
71    continue
      stack(m1+1)=m1-m
      m=m1+1
      return
80    do 81  i=1,e
81    ed(i)=t.eq.endpt(1,i)
      go to 64
      end



cccccc P.257
      subroutine hamcrc(n,a,k,m,stack,nstk,adj,vert,option)
      integer a1,a,option,stack
      logical adj(n,n),vert(n)
      dimension a(n),stack(nstk)
10    if(k.ne.1) go to 30
20    stack(1)=1
      stack(2)=1
      m=2
      return
30    k1=k-1
      a1=a(k1)
      do 31  i=1,n
31    vert(i)=adj(a1,i)
      do 32  i=1,k1
      m1=a(i)
32    vert(m1)=.false.
      m1=m
      if(k.eq.n) go to 50
40    do 41  i=1,n
      if(.not.vert(i)) go to 41
      m1=m1+1
      stack(m1)=i
41    continue
44    stack(m1+1)=m1-m
      m=m1+1
      return
50    do 51  i=1,n
      if(.not.vert(i)) go to 51
      if(option.eq.2) go to 52
      if(i.gt.a(2)) go to 44
52    if(.not.adj(i,1)) go to 44
      m=m+2
      stack(m-1)=i
      stack(m)=1
      return
51    continue
      go to 44
      end



ccccccccc P.263
      subroutine spntre(e,n,a,k,m,stack,nstk,endpt,end,x,nv,y)
      integer a,comp,e,end,endpt,stack,x,y
      dimension a(n),stack(nstk),endpt(2,e),end(2,n),nv(n),y(n),x(n)
10    if(k.ne.1) go to 30
20    n2=e-n+1
      do 21  i=1,n2
21    stack(i)=i
      m=n2+1
      stack(m)=n2
      return
30    k1=k-1
      do 31  i=1,k1
      i3=a(i)
      end(1,i)=endpt(1,i3)
31    end(2,i)=endpt(2,i3)
      n3=n+1
      call spanfo(n3,k1,end,comp,x,nv,y)
      i1=a(k1)+1
      i2=e-n+k
      m1=m
32    do 35  i=i1,i2
      i3=endpt(1,i)
      i4=endpt(2,i)
      if(x(i3).eq.x(i4)) go to 35
      m1=m1+1
      stack(m1)=i
35    continue
      stack(m1+1)=m1-m
      m=m1+1
      return
      end


      subroutine rantre(n,end,a,b,m)
      integer a,end
      logical b(n)
      dimension end(2,n),a(n),m(n)
      n2=n-2
10    do 11  i=1,n
      a(i)=uni(1)*float(n)
      a(i)=a(i)+1
      b(i)=.true.
11    m(i)=0
      m1=1
      do 20  i=1,n2
      i1=a(i)
20    m(i1)=m(i1)+1
30    l=a(m1)
40    do 50  i=1,n
      if(b(i).and.(m(i).eq.0)) go to 60
50    continue
60    end(1,m1)=i
      end(2,m1)=l
70    if(m1.eq.n2) go to 90
80    b(i)=.false.
      m(l)=m(l)-1
      m1=m1+1
      go to 30
90    end(2,n-1)=0
      do 130  i=1,n
100   if(.not.b(i)) go to 130
110   end(1,n-1)=end(2,n-1)
      end(2,n-1)=i
130   continue
      return
      end


ccccccc P.279
      subroutine ranrut(nn,out,stack,t)
      integer d,out,stack,sum,t,td,z
      dimension out(2,nn),stack(2,nn),t(nn)
      data nlast/1/
      t(1)=1
1     if(nn.le.nlast) go to 10
      sum=0
      do 2  d=1,nlast
      i=nlast+1
      td=t(d)*d
      do 3  j=1,nlast
      i=i-d
      if(i.le.0) go to 2
3     sum=sum+t(i)*td
2     continue
      nlast=nlast+1
      t(nlast)=sum/(nlast-1)
      go to 1
10    n=nn
      is1=0
      is2=0
12    if(n.le.2) go to 70
20    z=float((n-1)*t(n))*uni(1)
      d=0
30    d=d+1
      td=d*t(d)
      m=n
      j=0
40    j=j+1
      m=m-d
      if(m.lt.1) go to 30
50    z=z-t(m)*td
      if(z.ge.0) go to 40
60    is1=is1+1
      stack(1,is1)=j
      stack(2,is1)=d
      n=m
      go to 12
70    if(n.le.1) go to 71
      is2=is2+1
      out(1,is2)=1
      out(2,is2)=2
71    is2=is2+1
      out(1,is2)=n
80    n=stack(2,is1)
      if(n.eq.0) go to 90
      stack(2,is1)=0
      go to 12
90    j=stack(1,is1)
      is1=is1-1
100   n2=out(1,is2)
      is2=is2-n2
      n1=out(1,is2)
      inc1=n1
      inc2=n1
      l1=is2+1
      l2=is2+n2-1
      do 110  m=1,j
      out(1,is2)=1
      out(2,is2)=1+inc1
      is2=is2+1
      if(n2.eq.1) go to 130
      do 120  k=l1,l2
      out(1,is2)=out(1,k)+inc2
      out(2,is2)=out(2,k)+inc2
120   is2=is2+1
130   inc1=inc1+n2
110   inc2=inc1-n1
140   out(1,is2)=n1+j*n2
150   if(is1.eq.0) return
      go to 80
      end


      subroutine signum(sigma,n,sign,cycles)
      integer sigma(n),sign,cycles
      cycles=0
      do 5  i=1,n
5     sigma(i)=-sigma(i)
      do 10  i=1,n
      if(sigma(i).gt.0) go to 10
      cycles=cycles+1
      m=i
15    sigma(m)=-sigma(m)
      m=sigma(m)
      if(m.ne.i) go to 15
10    continue
      sign=1-2*mod(n+cycles,2)
      return
      end


ccccc P.140
      subroutine hpsort(n,b)
      integer b(n),bstar
      n1=n
      l=1+n/2
11    l=l-1
      bstar=b(l)
      go to 30
25    bstar=b(n1)
      b(n1)=b(1)
29    n1=n1-1
30    l1=l
31    m=2*l1
      if(m-n1) 32,33,37
32    if(b(m+1).ge.b(m)) m=m+1
33    if(bstar.ge.b(m)) go to 37
      b(l1)=b(m)
      l1=m
      go to 31
37    b(l1)=bstar
      if(l.gt.1) go to 11
      if(n1.ge.2) go to 25
      return
      end


ccccc P.285
      subroutine minspt(dist,n,endpt,u,y)
      integer u(n),endpt(2,n)
      dimension y(n),dist(n,n)
10    l=0
20    do 21  i=2,n
      u(i)=1
21    y(i)=dist(1,i)
30    dmin=1.e37
40    do 41  i=2,n
50    if(y(i).eq.0.0) go to 41
60    if(y(i).ge.dmin) go to 41
70    dmin=y(i)
      imin=i
41    continue
80    l=l+1
      endpt(1,l)=imin
      endpt(2,l)=u(imin)
90    if(l.eq.n-1) return
100   y(imin)=0
110   do 111  i=2,n
      if(y(i).eq.0.0) go to 111
      d1=dist(i,imin)
112   if(y(i).le.d1) go to 111
113   u(i)=imin
      y(i)=d1
111   continue
      go to 30
      end

SHAR_EOF
cat << 'SHAR_EOF' > rancom.f
From njas@research.att.com Fri May 26 08:10:08 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29921 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:08 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29921@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO


ccccccc P.52
      subroutine rancom(n,k,r)
c random k-composition of n-set. Ref NW 53.
      integer r(k)
      call ranksb(n+k-1,k-1,r)
      r(k)=n+k
      l=0
      do 10 i=1,k
      m=r(i)
      r(i)=m-l-1
10    l=m
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > ranequ.f
From njas@research.att.com Fri May 26 08:10:09 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29926 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:09 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29926@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO


cccccccccc P.97
      subroutine ranequ(n,l,q,b)
c random partition of n-set. NW 97.
      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)*uni(0)
      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

SHAR_EOF
cat << 'SHAR_EOF' > ranksb.f
From njas@research.att.com Fri May 26 08:10:10 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29931 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:10 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29931@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO

c ranksb
 integer a(100)
 1 read,n,k
 if(n.le.0)stop
 2  read,j
 if(j.le.0)goto 1
 call ranksb(n,k,a)
 print 30,(a(i),i=1,k)
 30 format(1h ,10i6)
 goto 2
 end


cccccccc P.43
      subroutine ranksb(n,k,a)
c random k-subset of n-set. Ref NW p 43.
      integer a(k),x,r,ds,p,s,c
      c=k
      do 1 i=1,k
    1 a(i)=(i-1)*n/k
   10 x=1+n*uni(1)
      l=1+(x*k-1)/n
      if(x.le.a(l))goto 10
      a(l)=a(l)+1
      c=c-1
      if(c.ne.0)goto 10
      p=0
      s=k
      do 20 i=1,k
      m=a(i)
      a(i)=0
      if(m.eq.(i-1)*n/k) goto 20
      p=p+1
      a(p)=m
   20 continue
   30 l=1+(a(p)*k-1)/n
      ds=a(p)-(l-1)*n/k
      a(p)=0
      a(s)=l
      s=s-ds
      p=p-1
      if(p.gt.0)goto 30
      l=k
   40 if(a(l).eq.0)goto 50
      r=l
      m0=1+(a(l)-1)*n/k
      m=a(l)*n/k-m0+1
   50 x=m0+m*uni(1)
      i=l
   60 i=i+1
      if(i.le.r)goto 80
   70 a(i-1)=x
      m=m-1
      l=l-1
      if(l.eq.0)return
      goto 40
   80 if(x.lt.a(i)) goto 70
      x=x+1
      a(i-1)=a(i)
      goto 60
      end

SHAR_EOF
cat << 'SHAR_EOF' > ranpar.f
From njas@research.att.com Fri May 26 08:10:12 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29936 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:11 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29936@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO



cccccc P.75
      subroutine ranpar(n,k,mult,p)
c random partition of n. Ref NW 76.
      integer p(n),d,mult(n)
      data nlast/0/
10    if(n.le.nlast) go to 30
20    p(1)=1
      m=nlast+1
      nlast=n
      if(n.eq.1) go to 30
      do 21  i=m,n
      isum=0
26    do 22  d=1,i
      is=0
      i1=i
24    i1=i1-d
      if(i1) 22,25,23
23    is=is+p(i1)
      go to 24
25    is=is+1
22    isum=isum+is*d
21    p(i)=isum/i
30    m=n
      k=0
      do 31  i=1,n
31    mult(i)=0
40    z=uni(i)*float(m*p(m))
      d=0
110   d=d+1
60    i1=m
      j=0
150   j=j+1
70    i1=i1-d
80    if(i1)  110,90,120
120   z=z-float(d*p(i1))
130   if(z) 145,145,150
90    z=z-float(d)
100   if(z) 145,145,110
145   mult(d)=mult(d)+j
      k=k+j
160   m=i1
170   if(m.ne.0) go to 40
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > ranper.f
From njas@research.att.com Fri May 26 08:10:13 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29941 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:12 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29941@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO

ccccc P.62
      subroutine ranper(n,a,setup)
c random perm of n letters. NW p63.
      integer a(n)
      logical setup
      if(.not.setup)goto 20
      do 10  i=1,n
   10 a(i)=i
   20 do 40  m=1,n
   30 l=m+uni(0)*(n+1-m)
      l1=a(l)
      a(l)=a(m)
   40 a(m)=l1
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > ransub.f
From njas@research.att.com Fri May 26 08:10:15 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29951 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:15 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29951@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO


ccccc P.23
      subroutine ransub(n,a)
c generate random subset of an n-set. Ref NW p 24.
      integer a(n)
      do 10  i=1,n
10    a(i)=2.*uni(1)
      return
      end

SHAR_EOF
cat << 'SHAR_EOF' > renumb.f
From njas@research.att.com Fri May 26 08:10:16 1995
Received: from research.att.com (research.att.com [192.20.225.3]) by dimacs.rutgers.edu (8.6.12+bestmx+oldruq+newsunq+grosshack/8.6.12) with SMTP id IAA29956 for <skiena@dimacs.rutgers.edu>; Fri, 26 May 1995 08:10:16 -0400
From: njas@research.att.com
Message-Id: <199505261210.IAA29956@dimacs.rutgers.edu>
Date: Fri, 26 May 95 08:09 EDT
To: skiena@dimacs.rutgers.edu
Status: RO

ccccc P.155
      subroutine renumb(m,n,sig,tau,a)
      integer sig(m),tau(n),a(m,n),t1,t2
      do 5  i=1,m
      i1=sig(i)
6     if(i1.le.i) go to 5
      i2=sig(i1)
      sig(i1)=-i2
      i1=i2
      go to 6
5     sig(i)=-sig(i)
      if(tau(1).lt.0) go to 9
      do 7  j=1,n
      j1=tau(j)
8     if(j1.le.j) go to 7
      j2=tau(j1)
      tau(j1)=-j2
      j1=j2
      go to 8
7     tau(j)=-tau(j)
9     do 10  i=1,m
      i1=-sig(i)
      if(i1.lt.0) go to 10
      lc=0
20    i1=sig(i1)
      lc=lc+1
      if(i1.gt.0) go to 20
      i1=i
      do 30  j=1,n
      if(tau(j).gt.0) go to 30
      j2=j
      k=lc
40    j1=j2
      t1=a(i1,j1)
50    i1=iabs(sig(i1))
      j1=iabs(tau(j1))
      t2=a(i1,j1)
      a(i1,j1)=t1
      t1=t2
      if(j1.ne.j2) go to 50
      k=k-1
      if(i1.ne.i) go to 50
      j2=iabs(tau(j2))
55    if(k.ne.0) go to 40
30    continue
10    continue
      do 60  i=1,m
60    sig(i)=iabs(sig(i))
      if(tau(1).gt.0) return
      do 70  j=1,n
70    tau(j)=iabs(tau(j))
      return
      end


ccccc P.167
      subroutine spanfo(n,e,endpt,k,x,nv,y)
      integer e1,e,endpt,s,t1,t2,t,v1,v2,x,y,z
      dimension endpt(2,e),x(n),nv(n),y(n),s(2)
      data s(1),s(2),m/1,2,2/
      do 10  i=1,n
      x(i)=-i
      nv(i)=1
10    y(i)=0
      j=1
      e1=e
20    v1=endpt(1,j)
      v2=endpt(2,j)
25    t1=x(v1)
      if(t1.lt.0)  t1=v1
      t2=x(v2)
      if(t2.lt.0)  t2=v2
      if(t1.ne.t2) go to 40
      if(j.lt.e1)  go to 30
      e1=e1-1
      go to 60
30    endpt(1,j)=endpt(1,e1)
      endpt(2,j)=endpt(2,e1)
      endpt(1,e1)=v1
      endpt(2,e1)=v2
      e1=e1-1
      go to 20
40    if(nv(t1).le.nv(t2))  go to 50
      t=t1
      t1=t2
      t2=t
      i3=-x(t2)
50    y(i3)=t1
      x(t2)=x(t1)
      i=t1
55    x(i)=t2
      i=y(i)
      if(i.ne.0) go to 55
      nv(t2)=nv(t2)+nv(t1)
      nv(t1)=0
      j=j+1
      if(j.le.e1.and.j.lt.n)  go to 20
60    k=0
      do 70  i=1,n
      if(nv(i).eq.0)  go to 70
      k=k+1
      nv(k)=nv(i)
      y(i)=k
70    continue
      do 80  i=1,n
      t=x(i)
      if(t.lt.0)  t=i
80    x(i)=y(t)
      if(k.eq.1)  return
90    i2=nv(1)
      nv(1)=1
      do 100  l=2,k
      i1=nv(l)
      nv(l)=nv(l-1)+i2-1
100   i2=i1
      do 110  i=1,e1
      i3=endpt(1,i)
      z=x(i3)
      y(i)=nv(z)
110   nv(z)=nv(z)+1
      call renumb(m,e1,s,y,endpt)
      i1=1
      do 120  l=1,k
      i2=nv(l)
      nv(l)=i2-i1+1
120   i1=i2
      return
      end

cccc p.175
      subroutine poly(n,a,x0,option,val,b)
      integer a,b,option,v,val,x0,z
      dimension a(n),b(n)
      val=a(n)
      if(n.eq.1) return
      n1=n-1
      if(option.eq.0) go to 20
      if(option.eq.(-1)) go to 26
      do 10 i=1,n
10    b(i)=a(i)
      if(option) 50,30,30
20    do 25 i=1,n1
      i1=n-i
25    val=val*x0+a(i1)
      return
26    do 27 i=1,n1
      i1=n-i
27    val=val*(x0-n1+i)+a(i1)
      return
30    max=min0(n1,option)
      do 35 j=1,max
      m=n1
      v=val
37    v=b(m)+v*x0
      b(m)=v
      m=m-1
      if(m.ge.j) go to 37
35    continue
      val=b(1)
      return
50    if(n.eq.2) return
      n2=n-2
      if(option.eq.(-3)) go to 70
      do 55  j=1,n2
      v=val
      m=n1
60    v=b(m)+j*v
      b(m)=v
      m=m-1
      if(m.gt.j) go to 60
55    continue
      return
70    do 75  j=1,n2
      z=n1-j
      m=z+1
80    b(m)=b(m)-z*b(m+1)
      m=m+1
      if(m.le.n1) go to 80
75    continue
      return
      end
      subroutine chromp(n,e,endpt,a,b,c,stack,nstk)
      integer a,ai,b,bj,c,e1,e,ed,en1,en2,en,end,endpt
      integer oth,r,stack,tr,v
      dimension endpt(2,e),stack(2,nstk),a(n),b(n),c(n)
10    call spanfo(n,e,endpt,k,a,b,c)
      if (k.gt.1) stop
12    do 13  i=1,n
      a(i)=0
13    b(i)=1
20    n1=n-1
      do 90  j=1,n1
30    do 40  l=1,2
50    i0=endpt(l,j)
      j0=0
60    j1=a(i0)
      a(i0)=j0
70    if(j1.eq.0) go to 40
80    j0=j1
      bj=b(j0)
      b(j0)=-bj
      i3=(3-bj)/2
      i0=endpt(i3,j0)
      go to 60
40    continue
      i3=endpt(2,j)
90    a(i3)=j
      r=endpt(1,n1)
      do 140  j=1,n1
      do 140  l=1,2
      i3=(3+b(j)*(2*l-3))/2
140   stack(l,j)=endpt(i3,j)
      do 170  i=1,n
      b(i)=0
      if(a(i).eq.0) go to 170
      k=0
      i1=i
150   ai=a(i1)
      if(ai.eq.0) go to 155
      a(i1)=0
      k=k+1
      c(k)=ai
      i1=stack(1,ai)
      go to 150
155   i1=i
      j0=j0+k
      do 160  j=1,k
      do 160  l=1,2
      i3=j0-j+1
      i4=c(j)
160   endpt(l,i3)=stack(l,i4)
170   continue
200   n1=n
      e1=e
      is=0
210   if(e1+1.ne.n1) go to 280
220   b(n1)=b(n1)+1
230   if(is.eq.0) go to 600
240   n1=stack(1,is)
      e1=stack(2,is)
      j0=is-e1
      do 241  j=1,e1
      do 242  l=1,2
242   endpt(l,j)=stack(l,j0)
241   j0=j0+1
250   if(e1.ne.n1) go to 260
270   b(n1)=b(n1)+1
      is=is-e1-1
      go to 310
280   if(e1.ne.n1) go to 300
290   b(n1)=b(n1)+1
      go to 310
300   do 301  j=1,e1
      is=is+1
      do 301  l=1,2
301   stack(l,is)=endpt(l,j)
      stack(1,is)=n1
      stack(2,is)=e1-1
      go to 310
260   is=is-1
      stack(1,is)=n1
      stack(2,is)=e1-1
310   do 311  i=1,n
311   a(i)=0
320   i1=endpt(1,e1)
      i2=endpt(2,e1)
      e1=e1-1
      tr=n1-1
      j=0
      if(i1.eq.r) go to 327
      if(i2.eq.r) go to 325
      do 321  j=1,tr
      en2=endpt(2,j)
322   if(en2.eq.i2) go to 325
323   if(en2.eq.i1) go to 326
321   continue
325   i0=i1
      i1=i2
      i2=i0
326   if(j.le.0) go to 327
      i3=endpt(1,j)
      a(i3)=1
327   ed=j
      j1=j+1
330   do 340  j=j1,tr
350   en1=endpt(1,j)
      en2=endpt(2,j)
360   if(en2.ne.i2) go to 400
370   end=en1
      go to 340
400   if(en1.ne.i2) go to 420
      en1=i1
      go to 430
420   if(en1.ne.i1) go to 440
430   a(en2)=1
440   ed=ed+1
      endpt(1,ed)=en1
      endpt(2,ed)=en2
340   continue
450   n1=tr
      endpt(1,n1)=i1
      endpt(2,n1)=end
460   do 470  j=n1,e1
480   do 485  l=1,2
490   en=endpt(l,j)
500   if(en.eq.i1) go to 530
510   if(en.eq.i2) go to 520
485   continue
555   oth=endpt(1,j)
      go to 560
520   en=i1
530   i3=3-l
      oth=endpt(i3,j)
540   if(a(oth).eq.1) go to 470
550   a(oth)=1
560   ed=ed+1
      endpt(1,ed)=en
      endpt(2,ed)=oth
470   continue
570   e1=ed
      go to 210
600   do 601  i=1,n
601   a(i)=(1-2*mod(n-i,2))*b(i)
      call poly(n,a,0,-2,v,c)
      call poly(n,a,-1,n,v,a)
      do 602  i=1,n
602   a(i)=iabs(a(i))
      return
      end

SHAR_EOF
:	End of shell archive
exit 0
