module cross_wilson{ # Wilson's factorization of cross-correlation
  use helicon
  use polydiv
  integer, parameter,         private :: n = 1024        
  real,    dimension( 2*n-1), private :: cross, ab, cc    
  real,    dimension( n),     private :: b, c            
contains
  subroutine cross_factor( niter, s0, ss, a0, aa, b0, bb) { 
    integer,                intent( in)     :: niter      # Newton iterations
    real,                   intent( in)     :: s0         # in: cross-correlat.
    type( filter),          intent( in)     :: ss         #        
    real,                   intent( in out) :: a0, b0     # out: min. phase
    type( filter),          intent( out)    :: aa, bb     #
    real                                    :: eps
    integer                                 :: i, stat
    cross = 0.; cross( n) = s0; b( 1) = 1.                # initialize 
    cross( n+ss%lag) = ss%flt                             # cross-correlation
    do i = 1, niter {
       call   polydiv_init( 2*n-1, aa)                    # divide   polynoms
       stat = polydiv_lop( .false., .false., cross, ab)   # ab = S/A
#       ab = ab/(a0*b0)                                   # scale
       call   polydiv_init( 2*n-1, bb)                    # divide   polynoms
       stat = polydiv_lop(  .true., .false., cc,    ab)   # cc = S/(AB')
       b( 2:n) = cc( n+1:2*n-1)/cc( n)                    # b = + side 
#       b( 1)   = 0.5*(cc( n) + 1.)                       #     of (1 + cc)
       eps = maxval(abs(b(2:n)))
       call   helicon_init( aa)                           # mutliply polynoms
       stat = helicon_lop( .false., .false., b, c)        # c = A b
#       c = c*a0; a0 = c( 1); aa%flt = c(1+aa%lag)/a0     # scale
        aa%flt = c(1+aa%lag)
       b( 2:n) = cc( n-1:1:-1)/cc( n)                     # b = - side 
       eps = max(maxval(abs(b(2:n))),eps)
       write( 0,*) i, eps                                 # "L1 norm"
       if( eps < epsilon(a0)) exit                        # convergence
       call   helicon_init( bb)                           # mutliply polynoms
       stat = helicon_lop( .false., .false., b, c)        # c = B b
#       c = c*b0; b0 = c( 1); bb%flt = c(1+bb%lag)/b0     # scale
	bb%flt = c(1+bb%lag)
    }
    call polydiv_close()
  }
}