module wilson { 
# Wilson's factorization
  use helicon
  use polydiv
  integer,                          private :: n
  real, dimension( :), allocatable, private :: auto, bb, cc, b, c    
contains
  subroutine wilson_init( nmax) {
    integer, intent (in) :: nmax;    n = nmax
    allocate ( auto( 2*n-1), bb( 2*n-1), cc( 2*n-1), b(n), c(n))
    }
  subroutine wilson_factor( niter, s0, ss, a0, aa, verb) { 
    integer,       intent( in)  :: niter # Newton iterations
    real,          intent( in)  :: s0    # autocorrelation zero lag
    type( filter), intent( in)  :: ss    # autocorrelation, other lags
    real,          intent( out) :: a0    # factor, zero lag
    type( filter)               :: aa    # factor, other lags
    logical,       intent( in)  :: verb
    optional                    :: verb
    real                        :: eps
    integer                     :: i, stat


    auto = 0.;  auto( n) = s0;  b( 1) =1.       # initialize 
    auto( n+ss%lag) =  ss%flt                   # autocorrelation
    auto( n-ss%lag) =  ss%flt                   # symmetrize input auto.
    call helicon_init( aa)                      # multiply polynoms
    call polydiv_init( 2*n-1, aa)               # divide   polynoms
    do i = 1, niter {
       stat= polydiv_lop(.false.,.false., auto, bb)  # bb = S/A
       stat= polydiv_lop( .true.,.false., cc,   bb)  # cc = S/(AA')
       b( 2:n) = 0.5*( cc( n+1:2*n-1  ) + 
                       cc( n-1:1    :-1)) / cc( n)     # b = plusside(1+cc) 
       eps = maxval( abs( b(2:n)));   # "L1 norm"
       if (present (verb)) { if (verb) write (0,*) i, eps}
       stat= helicon_lop( .false., .false., b, c)      # c = A b
       aa%flt = c( 1+aa%lag)                           # put on helix
       if( eps < epsilon( a0))  break                  # convergence
       }
    a0 = sqrt( cc( n))
    }
  subroutine wilson_close () {
    deallocate( auto, bb, cc, b, c)
    call polydiv_close()  	
    }
}