      PROGRAM MAIN
      implicit double precision(a-h,o-z)
      parameter nfncs=12
      parameter bad=2.459e-14
      double precision av(nfncs),sq(nfncs)
      double precision avb(nfncs),sqb(nfncs)
      double precision x,tax
      character*30 fnstring(nfncs)
      data fnstring/'x','t','u=x+t','v=x-t','[(u-1)(u+1)(v-1)(v+1)]^2',
     +  '(1+COSpiu)(1+COSpiv)','x[(u-1)(u+1)(v-1)(v+1)]^2',
     +  't[(u-1)(u+1)(v-1)(v+1)]^2','x(1+COSpiu)(1+COSpiv)',
     +  't(1+COSpiu)(1+COSpiv)','(1+COS5piu)(1+COS5piv)',
     +  '(1-u)[(u-1)(u+1)(v-1)(v+1)]^2'/
      do j=1,nfncs
        av(j)=0.
        avb(j)=0.
        sq(j)=0.
        sqb(j)=0.
      end do
      n=0
      nb=0
      istop=1
      i=0
10    i=i+1
      istaxhigh=0
      do j=1,nfncs
15    x=pobox(i,j,tax)
      if (j.gt.1) goto 18
      if (tax.gt.bad) then
        print *
        print *,'        tax is an exceptionally large error'
        print *
        istaxhigh=1
      end if
      nb=nb+1
      if (istaxhigh.eq.0) n=n+1
18    print *,' RUN #',i
      print *,' FUNCTION ',fnstring(j)
      if (n.eq.0) goto 30
      if (n.eq.1) goto 20
      sq(j)=sq(j)*(n-2)
      sq(j)= sq(j) + (av(j)**2)*(n-1)
20    av(j)=((n-1)*av(j)+x)/n
      if (n.eq.1) goto 30
      sq(j)=sq(j)-n*(av(j)**2)
      sq(j)=sq(j)+x**2
      sq(j)=sq(j)/(n-1)
30    if (tax.eq.0.0) goto 50
      if (nb.eq.1) goto 40
      sqb(j)=sqb(j)*(nb-2)
      sqb(j)= sqb(j) + (avb(j)**2)*(nb-1)
40    avb(j)=((nb-1)*avb(j)+x)/nb
      if (nb.eq.1) goto 50
      sqb(j)=sqb(j)-nb*(avb(j)**2)
      sqb(j)=sqb(j)+x**2
      sqb(j)=sqb(j)/(nb-1)
50    print *,' action =',x
      if (n.eq.0) goto 60
      print *,n,' good data points give:'
      print *,' running average=',av(j)
      print *,' running standard deviation=',sq(j)**.5
      print *,' running standard error of the mean=',(sq(j)/n)**.5
      print *
60    if (tax.eq.0.0) goto 70
      print *,nb,' total data points give:'
      print *,' running average=',avb(j)
      print *,' running standard deviation=',sqb(j)**.5
      print *,' running standard error of the mean=',(sqb(j)/nb)**.5
      print *
70    end do
      if (i.lt.istop) goto 10
      end
 
      DOUBLE PRECISION FUNCTION POBOX(little_dumb,jfn,tax)
      parameter nfncs=12
      parameter nexp=200
      parameter mex=nexp+60
      parameter ndim=2
      parameter fr=1.
      dimension l(mex,mex)
      double precision f(mex,nfncs),c(mex,4),cin(mex,mex)
      double precision dummy(mex,mex),indx(mex)
      dimension icheck(mex,mex)
      double precision pi,tax,sum,s,action(nfncs)
c     dimension lilac(42)
c     dimension lilac(2)
      character*3 meth
      common pi
c     data lilac/111222,111333/
c     data lilac/111222,111333,111444,111555,111666,111777,111888,111999,
c     +         222111,222333,222444,222555,222666,222777,222888,222999,
c     +         333111,333222,333444,333555,333666,333777,333888,333999,
c     +         444111,444222,444333,444555,444666,444777,444888,444999,
c     +         555111,555222,555333,555444,555666,555777,555888,555999,
c     +         314159,271828/
c     ___________________________
c     meth='sym' gives symmetric
c          green's function
c     meth='ret' gives retarded
c          green's function
c     ---------------------------
      meth='sym'
c     meth='ret'
      if (jfn.gt.1) goto 200
      pi=2*dasin(1.0d0)
      iseed=itime(1)
c     ****************************************
c     to make portable, can use ran1 or ran3 in
c     place of ran in all calls for random #'s
c     also the call for time in ITIME must be
c     modified to correct form for specific
c     machine being used
c     ****************************************
      iseed=222666
c     iseed=lilac(little_dumb)
      print *,' RUN #',little_dumb
      print *,' iseed=',iseed
      x=ran3(-iseed)
10    n=npt(nexp,iseed)
      n=6
      if (n.gt.mex) goto 10
      if (little_dumb.eq.1) goto 111
      do i=1,n
        do j=i,n
                cin(j,i)=0.
                dummy(j,i)=0.
                icheck(j,i)=0
                l(j,i)=0
                cin(i,j)=0.
                dummy(i,j)=0.
                icheck(i,j)=0
                l(i,j)=0
        end do
      end do
c     ***********************************************
c     This allows interactive control over
c       # of points and iseed
c     print *,' input # of points:'
c     read *,n
c     print *,' type in random seed(6 digits or less):'
c     read *,iseed
c     ***********************************************
111   print *,' number of points =',n
      call sprinkle(mex,n,ndim,c,f,iseed,nfncs)
      print *,' coordinates :'
      print *,((c(i,j),j=1,ndim),i=1,n)
      print *,' field values :'
      print *,(f(i,9),i=1,n)
      call setc(mex,n,c,ndim,l)
      print *,' causal matrix :'
      print *,((l(i,j),j=1,n),i=1,n)
      call relab(mex,n,l,f,nfncs)
      print *,' relabelled causal matrix :'
      print *,((l(i,j),j=1,n),i=1,n)
      print *,' relabelled field values :'
      print *,(f(i,9),i=1,n)
      if (ndim.eq.4) call link(mex,l,n)
      if (meth.eq.'ret') call gret(mex,n,l,dummy,icheck,fr)
      if (meth.eq.'sym') then
        call gsym(mex,n,l,dummy,icheck)
        call discard1(mex,n,dummy,icheck,f,nfncs)
        call discard2(mex,n,dummy,icheck,f,nfncs)
        print *,' number of points (after discards) =',n
      end if
      print *,' matrix to be inverted :'
      print *,((dummy(i,j),j=1,n),i=1,n)
c     print *,((icheck(i,j),j=1,n),i=1,n)
      call inv(mex,n,dummy,cin,indx)
      print *,' inverse of causal matrix :'
      print *,((cin(i,j),j=1,n),i=1,n)
      tax=0.0
c     if (meth.eq.'sym') then
      do i=1,n
        do j=1,n
                sum=0.0
                do k=1,n
                        sum=sum+icheck(i,k)*cin(k,j)
                end do
c               print *,' identity matrix elements :'
c               print *,i,j,sum
                if ((i.eq.j).and.(dabs(sum-1.).gt.tax)) tax=dabs(sum-1.)
                if ((i.ne.j).and.(dabs(sum).gt.tax)) tax=dabs(sum)
        end do
      end do
      print *,' tax=',tax
c     end if
      do j=1,nfncs
      s=act(mex,n,cin,f,j,nfncs)
      if (meth.eq.'sym') s=4*s
      if (meth.eq.'ret') s=2*s
      action(j)=s
      end do
200   pobox=action(jfn)
      end
 
      SUBROUTINE GSYM(mex,n,l,dummy,icheck)
      dimension l(mex,mex),icheck(mex,mex)
      double precision dummy(mex,mex)
      diag=0.0
      idiag=0
      do j=2,n
        do i=1,j-1
                dummy(i,j)=l(i,j)
                dummy(j,i)=l(i,j)
                icheck(i,j)=l(i,j)
                icheck(j,i)=l(i,j)
        end do
        dummy(j,j)=diag
        icheck(j,j)=idiag
      end do
      dummy(1,1)=diag
      icheck(1,1)=idiag
      end
 
      SUBROUTINE GRET(mex,n,l,dummy,icheck,fr)
      dimension l(mex,mex),icheck(mex,mex)
      double precision dummy(mex,mex)
      do j=2,n
        do i=1,j-1
                dummy(i,j)=l(i,j)
                dummy(j,i)=0.0
                icheck(i,j)=l(i,j)
                icheck(j,i)=0
        end do
        dummy(j,j)=fr
      end do
      dummy(1,1)=fr
      end
 
      DOUBLE PRECISION FUNCTION ACT(mex,n,c,f,jl,nfncs)
      double precision c(mex,mex),f(mex,nfncs),pi,sum
      common pi
      sum=0
      do i=1,n
        do j=1,n
                        sum=sum+f(i,jl)*c(j,i)*f(j,jl)
        end do
      end do
      act=sum
      end
      DOUBLE PRECISION FUNCTION FLD(x,t,j,nfncs)
      double precision x,t,pi,u,v
      common pi
      u=x+t
      v=x-t
      if (j.eq.1) fld=x
      if (j.eq.2) fld=t
      if (j.eq.3) fld=u
      if (j.eq.4) fld=v
      if (j.eq.5) then
        fld=(u-1)*(u+1)*(v-1)*(v+1)
        fld=fld**2
      end if
      if (j.eq.6) fld=(1+dcos(pi*u))*(1+dcos(pi*v))
      if (j.eq.7) then
        fld=(u-1)*(u+1)*(v-1)*(v+1)
        fld=fld**2
        fld=x*fld
      end if
      if (j.eq.8) then
        fld=(u-1)*(u+1)*(v-1)*(v+1)
        fld=fld**2
        fld=t*fld
      end if
      if (j.eq.9) fld=x*(1+dcos(pi*u))*(1+dcos(pi*v))
      if (j.eq.10) fld=t*(1+dcos(pi*u))*(1+dcos(pi*v))
      if (j.eq.11) fld=(1+dcos(5.*pi*u))*(1+dcos(5.*pi*v))
      if (j.eq.12) then
        fld=(u-1)*(u+1)*(v-1)*(v+1)
        fld=fld**2
        fld=(1-u)*fld
      end if
      if ((j.lt.1).or.(j.gt.nfncs)) then
        print *,' function # out of range'
        call halt
      end if
      end
      SUBROUTINE INV(np,n,a,y,indx)
      dimension indx(np)
      double precision a(np,np),y(np,np)
      do i=1,n
        do j=1,n
                y(j,i)=0.
        end do
        y(i,i)=1.
      end do
      call ludcmp(a,n,np,indx,d)
      do j=1,n
        call lubksb(a,n,np,indx,y(1,j))
      end do
      end
      SUBROUTINE discard1(mex,n,du,ich,f,nfncs)
      dimension ich(mex,mex)
      double precision du(mex,mex),f(mex,nfncs)
      if (n.eq.1) then
        print *,' matrix of rank 1 and so no deletion possible'
        print *,' exiting discard1 immediately'
        goto 40
      end if
      if (n.lt.1) then
        print *,' matrix of rank < 1 not possible !'
        call halt
      end if
      i=0
10    i=i+1
      j=0
20    j=j+1
      if (ich(i,j).eq.1) goto 30
      if (j.lt.n) goto 20
      call del(mex,n,ich,du,i,f,nfncs)
      print *,' a row of zeros has been deleted'
      i=i-1
30    if (i.lt.n) goto 10
40    end
      SUBROUTINE discard2(mex,n,du,ich,f,nfncs)
      dimension ich(mex,mex)
      double precision du(mex,mex),f(mex,nfncs)
      if (n.eq.1) then
        print *,' matrix of rank 1 and so no deletion possible'
        print *,' exiting discard 2 immediately'
        goto 50
      end if
      if (n.lt.1) then
        print *,' matrix of rank < 1 not possible !'
        call halt
      end if
      i=0
10    i=i+1
      j=i
20    j=j+1
      k=0
30    k=k+1
      if (ich(i,k).ne.ich(j,k)) goto 40
      if (k.lt.n) goto 30
      call del(mex,n,ich,du,j,f,nfncs)
      print *,' non-Hegelian point has been removed'
      j=j-1
40    if (j.lt.n) goto 20
      if (i.lt.n-1) goto 10
50    end
      SUBROUTINE del(mex,n,ich,du,i,f,nfncs)
      dimension ich(mex,mex)
      double precision du(mex,mex),f(mex,nfncs)
      if (n.le.1) then
        print *,' matrix of rank < 2 and so will not delete'
        call halt
      end if
      if ((i.lt.1).or.(i.gt.n)) then
        print *,' trying to delete point which does not exist !'
        call halt
      end if
      if (i.eq.n) goto 10
        do k=i+1,n
                do m=1,n
                        du(k-1,m)=du(k,m)
                        ich(k-1,m)=ich(k,m)
                end do
                do jc=1,nfncs
                        f(k-1,jc)=f(k,jc)
                end do
        end do
        do m=1,n-1
                do k=i+1,n
                        du(m,k-1)=du(m,k)
                        ich(m,k-1)=ich(m,k)
                end do
        end do
10      n=n-1
      end
      SUBROUTINE HALT
      print *,' execution of program halted'
      stop
      end
      FUNCTION ITIME(n)
      character *23 date
      character a,b,c,d,e,f
      call lib$date_time(date)
      a=date(14:14)
      b=date(17:17)
      c=date(19:19)
      d=date(20:20)
      e=date(22:22)
      f=date(23:23)
      if ((ichar(e).eq.0).and.(ichar(f).eq.0)) goto 10
      itime=(ichar(e)-48)*10000+(ichar(d)-48)*1000+(ichar(c)-48)*100
      itime=itime+(ichar(b)-48)*10+ichar(a)-48
      itime=itime+(ichar(f)-48)*100000
      x=ran3(-itime)
10    x=ran3(1)
      if (x.lt.0.1) goto 10
      itime=int(x*1000000.)
      end
      FUNCTION RAN1(idum)
c     to start must set idum=any neg. integer with first call
      dimension r(97)
      parameter m1=259200,ia1=7141,ic1=54773,rm1=1./m1
      parameter m2=134456,ia2=8121,ic2=28411,rm2=1./m2
      parameter m3=243000,ia3=4561,ic3=51349
      data iff /0/
      if ((idum.lt.0).or.(iff.eq.0)) then
        iff=1
        ix1=mod(ic1-idum,m1)
        ix1=mod(ia1*ix1+ic1,m1)
        ix2=mod(ix1,m2)
        ix1=mod(ia1*ix1+ic1,m1)
        ix3=mod(ix1,m3)
        do j=1,97
                ix1=mod(ia1*ix1+ic1,m1)
                ix2=mod(ia2*ix2+ic2,m2)
                r(j)=(float(ix1)+float(ix2)*rm2)*rm1
        end do
        idum=1
      end if
      ix1=mod(ia1*ix1+ic1,m1)
      ix2=mod(ia2*ix2+ic2,m2)
      ix3=mod(ia3*ix3+ic3,m3)
      j=1+(97*ix3)/m3
      if ((j.gt.97).or.(j.lt.1)) pause
      ran1=r(j)
      r(j)=(float(ix1)+float(ix2)*rm2)*rm1
      return
      end
      FUNCTION RAN3(idum)
c     initially call with a negative integer
      implicit real*4(m)
      parameter mbig=4000000.,mseed=1618033.,mz=0.,fac=1./mbig
      dimension ma(55)
      data iff /0/
      if ((idum.lt.0).or.(iff.eq.0)) then
        iff=1
        mj=mseed-iabs(idum)
        mj=mod(mj,mbig)
        ma(55)=mj
        mk=1
        do i=1,54
                ii=mod(21*i,55)
                ma(ii)=mk
                mk=mj-mk
                if (mk.lt.mz) mk=mk+mbig
                mj=ma(ii)
        end do
        do k=1,4
                do i=1,55
                        ma(i)=ma(i)-ma(1+mod(i+30,55))
                        if (ma(i).lt.mz) ma(i)=ma(i)+mbig
                end do
        end do
        inext=0
        inextp=31
        idum=1
      end if
      inext=inext+1
      if (inext.eq.56) inext=1
      inextp=inextp+1
      if (inextp.eq.56) inextp=1
      mj=ma(inext)-ma(inextp)
      if (mj.lt.mz) mj=mj+mbig
      ma(inext)=mj
      ran3=mj*fac
      return
      end
      SUBROUTINE link(mex,k,num)
      dimension k(mex,mex),kh(1024,1024)
      do i=1,num
        do j=1,num
                do l=1,num
                        if ((k(i,l).eq.1).and.(k(l,j).eq.1)) goto 10
                end do
                goto 20
10              kh(i,j)=0
20      end do
      end do
      end
      FUNCTION NPT(nexp,iseed)
      x=ran3(iseed)
      n=nexp
      p=0
      j=0
100   p=p+pdv(j,n)
      if (x.lt.p) goto 200
      j=j+1
      goto 100
200   npt=j
      end
      SUBROUTINE RELAB(mex,n,l,f,nfs)
      dimension l(mex,mex)
      double precision x,f(mex,nfs)
      do i=n,2,-1
50      do j=1,i-1
                if (l(i,j).eq.0) goto 100
                do k=1,n
                        m=l(i,k)
                        l(i,k)=l(j,k)
                        l(j,k)=m
                end do
                do k=1,n
                        m=l(k,i)
                        l(k,i)=l(k,j)
                        l(k,j)=m
                end do
                do k=1,nfs
                        x=f(i,k)
                        f(i,k)=f(j,k)
                        f(j,k)=x
                end do
                goto 50
100     end do
      end do
      end
      SUBROUTINE SETC(mex,n,c,ndim,l)
      dimension l(mex,mex)
      double precision c(mex,4)
      double precision dx,dy,dz,dt,ds2
      if (ndim.eq.4) then
      do i=1,n-1
        do j=i+1,n
                dx=c(i,1)*dsin(c(i,2))*dcos(c(i,3))
                dx=dx-c(j,1)*dsin(c(j,2))*dcos(c(j,3))
                dy=c(i,1)*dsin(c(i,2))*dsin(c(i,3))
                dy=dy-c(j,1)*dsin(c(j,2))*dsin(c(j,3))
                dz=c(i,1)*dcos(c(i,2))-c(j,1)*dcos(c(j,2))
                dt=c(i,4)-c(j,4)
                ds2=dx**2+dy**2+dz**2-dt**2
                if (ds2.gt.0) goto 100
                if (dt.gt.0) l(j,i)=1
                if (dt.lt.0) l(i,j)=1
100     end do
      end do
      end if
      if (ndim.eq.2) then
      do i=1,n-1
        do j=i+1,n
                dx=c(i,1)-c(j,1)
                dt=c(i,2)-c(j,2)
                ds2=dx**2-dt**2
                if (ds2.gt.0) goto 200
                if (dt.gt.0) l(j,i)=1
                if (dt.lt.0) l(i,j)=1
200     end do
      end do
      end if
      end
      SUBROUTINE SPRINKLE(mex,n,ndim,coor,field,iseed,nf)
      double precision pi,x1,x2,x3,x4,x,t,coor(mex,4),field(mex,nf)
      common pi
      print *,' at start of routine sprinkle iseed is',iseed
      iseed=-iseed
      do i=1,n
        x1=ran3(iseed)
        x2=ran3(iseed)
        if (ndim.eq.4) then
                x3=ran3(iseed)
                x4=ran3(iseed)
        end if
        if (ndim.eq.2) then
                coor(i,1)=x1-x2
                coor(i,2)=x1+x2-1
        end if
        if (ndim.eq.4) then
                coor(i,4)=2*x4-1
                coor(i,1)=x1**(1./3)*coor(i,4)
                coor(i,2)=dacos(1-2*x2)
                coor(i,3)=2*pi*x3
        end if
        x=coor(i,1)
        t=coor(i,ndim)
        do j=1,nf
                field(i,j)=fld(x,t,j,nf)
        end do
      end do
      end
      FUNCTION PDV(i,n)
      f=0
      do j=1,i
        f=f+alog(real(j))
      end do
      pdv=i*alog(real(n))-f-real(n)
      pdv=exp(pdv)
      end
      SUBROUTINE LUBKSB(a,n,np,indx,b)
      dimension indx(n)
      double precision sum,a(np,np),b(n)
      ii=0
      do i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0) then
                do j=ii,i-1
                        sum=sum-a(i,j)*b(j)
                end do
        else if (sum.ne.0.) then
                ii=i
        end if
        b(i)=sum
      end do
      do i=n,1,-1
        sum=b(i)
        if (i.lt.n) then
                do j=i+1,n
                        sum=sum-a(i,j)*b(j)
                end do
        end if
c       print *,sum,a(i,i)
        b(i)=sum/a(i,i)
      end do
      end
      SUBROUTINE LUDCMP(a,n,np,indx,d)
      parameter nmax=1000
      parameter tiny=1.0e-10
      dimension indx(n)
      double precision aamax,sum,dum,a(np,np),vv(nmax)
      d=1.
      do i=1,n
        aamax=0.
        do j=1,n
                if (dabs(a(i,j)).gt.aamax) aamax=dabs(a(i,j))
        end do
        if (aamax.eq.0) pause 'Singular matrix'
        vv(i)=1./aamax
      end do
      do j=1,n
        if (j.gt.1) then
                do i=1,j-1
                        sum=a(i,j)
                        if (i.gt.1) then
                                do k=1,i-1
                                        sum=sum-a(i,k)*a(k,j)
                                end do
                                a(i,j)=sum
                        end if
                end do
        end if
        aamax=0.
        do i=j,n
                sum=a(i,j)
                if (j.gt.1) then
                        do k=1,j-1
                                sum=sum-a(i,k)*a(k,j)
                        end do
                        a(i,j)=sum
                end if
                dum=vv(i)*dabs(sum)
                if (dum.ge.aamax) then
                        imax=i
                        aamax=dum
                end if
        end do
        if (j.ne.imax) then
                do k=1,n
                        dum=a(imax,k)
                        a(imax,k)=a(j,k)
                        a(j,k)=dum
                end do
                d=-d
                vv(imax)=vv(j)
        end if
        indx(j)=imax
        if (j.ne.n) then
                if (a(j,j).eq.0.) a(j,j)=tiny
                dum=1./a(j,j)
                do i=j+1,n
                        a(i,j)=a(i,j)*dum
                end do
        end if
      end do
      if (a(n,n).eq.0.) a(n,n)=tiny
      end
 
