module autocorr {
  use helix
  use compress
contains
  function autocorr1( aa, s0, a0) result ( ss) {
    type( filter)                       :: aa, ss
    real, intent (in), optional         :: a0
    real, intent (out)                  :: s0
    integer                             :: i, j, k, n, na
    if (present (a0))
      s0 = a0*a0 + dot_product (aa%flt,aa%flt)
    else
      s0 = 1.    + dot_product (aa%flt,aa%flt)
    na = size (aa%lag)
    call allocatehelix (ss, na*(na+1)/2) 
    k = na
    if (present (a0))
      ss%flt (:na) = aa%flt * a0
    else
      ss%flt (:na) = aa%flt
    ss%lag (:na) = aa%lag
    do i = 1, na {
      do j = i+1, na {
        k = k + 1
        ss%flt (k) = aa%flt (j) * aa%flt (i)
        ss%lag (k) = aa%lag (j) - aa%lag (i)
        do n = 1, k-1 {
           if (ss%lag (n) == ss%lag (k) .and. ss%flt (n) /= 0.) {
              ss%flt (n) = ss%flt (n) + ss%flt (k)
              ss%flt (k) = 0.
           }
    }}}
    call resize (ss)
  }
}