module npolydiv {			  # Helix polynomial division
use nhelix
integer                           :: nd
type( nfilter)                    :: aa
real, dimension (nd), allocatable :: tt
integer, dimension(:), pointer :: lag
real,    dimension(:), pointer :: flt
#%  _init ( nd, aa)
#%  _lop  ( xx, yy)
integer                        :: ia, ix, iy, ip
tt = 0.



allocate(flt(size(aa%hlx(1)%lag)))
allocate(lag(size(aa%hlx(1)%flt)))
if( adj) {
        tt = yy
	do iy= nd, 1, -1 { 
#		if(associated(aa%mis)){if(aa%mis(iy)){ xx(iy)=0.; next}}
		if(associated(aa%mis)){if(aa%mis(iy)){ tt(iy)=0.; next}}
		ip = aa%pch( iy)
		lag = aa%hlx( ip)%lag
		flt = aa%hlx( ip)%flt
		do ia = 1, size( lag){
			ix = iy - lag( ia);  
   		if( ix < 1)  cycle
#			if(ix <0 || ix >size(xx)) write(0,*) "1x error",ix,size(xx)
			tt( ix) -=  flt( ia) * tt( iy)
			} 
		}
	xx += tt
 } else { 
        tt = xx
	do iy= 1, nd { 
#		if(associated(aa%mis)){if(aa%mis(iy)){ yy(iy)=0.; next}}
		if(associated(aa%mis)){if(aa%mis(iy)){ tt(iy)=0.; next}}
		ip = aa%pch( iy)
		lag = aa%hlx( ip)%lag; 
flt = aa%hlx( ip)%flt
		do ia = 1, size( lag) {
			ix = iy - lag( ia);      if( ix < 1)  cycle
#			if(ix <0 || ix >size(xx)) write(0,*) "2x error",ix,size(xx),size(yy),lag(ia)
			tt( iy) -=  flt( ia) * tt( ix)
			} 
		}
	yy += tt
        }
deallocate(lag,flt)
}