# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "<stdin>"
module polydiv1
! Polynomial division (recursive filtering)
  use adj_mod
  implicit none
  real ,private, dimension (:), pointer :: aa
  contains
  subroutine polydiv1_init ( aa_in )
    real , dimension (:), pointer :: aa_in
    aa => aa_in
  end subroutine
  function polydiv1_lop ( adj, add, xx, yy) result(stat)
    integer :: stat
    logical,intent(in) :: adj,add
    real,dimension(:) :: xx,yy
    call adjnull (adj,add,xx,yy )
    call polydiv1_lop2(adj,add,xx,yy )
    stat=0
  end function
  subroutine polydiv1_lop2(adj,add,xx,yy)
    logical,intent(in) :: adj,add
    real, dimension (:) :: xx
    real, dimension (:) :: yy
    integer ia, ix, iy
    real tt
    if ( adj) then
      do ix= size(xx), 1, -1
        tt = yy( ix)
        do ia = 1, min( size(aa), size (xx) - ix)
          iy = ix + ia
          tt = tt - aa( ia) * xx( iy)
        end do
        xx( ix) = xx( ix) + tt
      end do
    else
      do iy= 1, size(xx)
        tt = xx( iy)
        do ia = 1, min( size(aa), iy-1)
          ix = iy - ia
          tt = tt - aa( ia) * yy( ix)
        end do
        yy( iy) = yy( iy) + tt
      end do
    end if
  end subroutine
  subroutine polydiv1_close()
  end subroutine
end module