subroutine lazy(nmax,y,k,eps,it)
c?  noise reduction using constant fits.

c  Author: Thomas Schreiber, Niels Bohr Institute, Blegdamsvej 17,
c          DK-2100 Copenhagen O
c          schreib@complex.nbi.dk
c  should work under ULTRIX (DECstation) and VAX/VMS at least. 
c  On other machines the bitwise .and.127 could cause problems.
c  It should perform a modulo 128
c  Ask me in case of problems.

c  neighbor search described in P. Grassberger, Phys. Lett. A 148.3.

c  input parameters:
c    nmax:  length of the series
c    y:     the (scalar=real array) series
c    k:     (integer) embedding for past/future. 
c           total embedding dimension is 2*k+1 (past+present+future)
c    eps:   size of neighborhoods (changed on exit)
c  on exit:
c    y:     contains corrected series
c    it:    number of iterations made

      parameter(nx=100000,im=127,tiny=1e-6)
c    nx:    maximal length of series to expect

      real y(nx), z(nx)
      integer ptr(0:im,0:im),pntr(nx)

      it=0
1     it=it+1
      epsinv=1./eps

      do 1000 i=0,im
         do 1000 j=0,im
1000        ptr(i,j)=0
        
      do 2000 n=k+1,nmax-k   
c make linked list
         i1=int(y(n-k)*epsinv).and.im
         i2=int(y(n+k)*epsinv).and.im
         ipn=ptr(i1,i2)
         ptr(i1,i2)=n
2000     pntr(n)=ipn

      do 3000 n=k+1,nmax-k   
c find neighbors
         ne=0
         s=0
         i1=int(y(n-k)*epsinv).and.im
         i2=int(y(n+k)*epsinv).and.im

         do 3100 j2=i2-1,i2+1    
c scan neighboring boxes
            do 3100 j1=i1-1,i1+1
               l1=j1.and.im
               l2=j2.and.im
               ip=ptr(l1,l2)
               if(ip.ne.0) then
3                 if(abs(y(n)-y(ip)).ge.eps) goto 4
                  do 3110 m=k,1,-1
                     if(abs(y(n+m)-y(ip+m)).ge.eps) goto 4
3110                 if(abs(y(n-m)-y(ip-m)).ge.eps) goto 4
                  ne=ne+1
                  s=s+y(ip)
c compute the center of mass of the neighborhood = zero-th order fit.
c you could introduce weights at this point
4                 ip=pntr(ip)
                  if(ip.ne.0) goto 3
               endif
3100        continue
3000     z(n)=s/ne

      eps=0
      do 4000 n=k+1,nmax-k
         eps=eps+(z(n)-y(n))**2
4000     y(n)=z(n)
      eps=sqrt(eps/(nmax-2*k))
      if(eps.gt.tiny) goto 1
      end