module polydiv1 { # Polynomial division (recursive filtering) real, dimension (:), pointer :: aa #% _init ( aa) #% _lop ( xx, yy) integer ia, ix, iy real tt if( adj) do ix= size(xx), 1, -1 { tt = yy( ix) do ia = 1, min( size(aa), size (xx) - ix) { iy = ix + ia tt -= aa( ia) * xx( iy) } xx( ix) = xx( ix) + tt } else do iy= 1, size(xx) { tt = xx( iy) do ia = 1, min( size(aa), iy-1) { ix = iy - ia tt -= aa( ia) * yy( ix) } yy( iy) = yy( iy) + tt } }