Code for crosscorrelation factorization

The code below will factor the polynomial ( 2/Z + 7 + 3Z) into its two factors (3+1/Z)(2+Z) except for a scale factor which is indeterminate. Choosing to set u₀=0 leads to the interesting factorization ( 2/Z + 7 + 3Z)/6 = (1+1/(3Z)) (1+Z/2) where the scale is chosen so the coefficient of Z⁰ is 1 (a convenience when you have factored a band matrix into two parts, and plan polynomial division for back substitution). Being based on fast FT, the work space size, n, must be a power of 2. Periodicity of FT requires that one side of the crosscorrelation is at the beginning of rr(:) while the other side is at the end. A test case for ( 2/Z + 7 + 3Z) is:

         rr( n)  =  2.        # first negative lag
         rr( 1)  =  7.        # zero lag
         rr( 2)  =  3.        # first positive lag