synthetic division algorithm to deflate conjugate pairs

Question

I am implementing a synthetic division algorithm in C as follows:

int deflate( double r, double im, double* poly, int n ) {
    int retval;
    int i;
    if( im == 0 ) {
            if( n < 1 ) {
                    retval = 1;
            } else {
                    double* out = ( double* )malloc( ( n )*sizeof( double ) );
                    out[ n - 1 ] = poly[ n ];
                    for( i = n - 2; i >= 0; i-- ) {
                            out[ i ] = out[ i + 1 ]*r + poly[ i + 1 ];
                    }

                    for( i = 0; i < n; i++ ) {
                            poly[ i ] = out[ i ];
                    }
                    poly[ n ] = 0;
                    free( out );
                    retval = 0;
            }
    } else {
            if( n < 2 ) {
                    retval = 1;
            } else {
                    /*code to handle complex numbers here*/
                    retval = 0;
            }
    }
    return retval;
}

I am trying to think of an efficient way to implement this for a non-zero imaginary component. Specifically, I would like to deflate both complex conjugate roots in one pass, without having to work with a complex coefficient polynomial. Can anyone think of a way to do this?

Answer 1

I've managed to figure out a way to do this. The solution is a simple matter of understanding that deflating a conjugate pair is just deflating two roots at once, both of which form a second polynomial, which is the divisor of the operation:

double div1 = -2*r;
double div0 = r*r + im*im;

Now use these two values to do synthetic division in the same way as you would if you were doing long division between a polynomial and a quadratic whose leading coefficient is 1 (which is essentially what synthetic division is in the case of a real root, except the divisor is linear whether than quadratic):

int i;
out[ n - 2 ] = poly[ n ];
out[ n - 3 ] = -out[ n - 2 ]*div1 + poly[ n - 1 ];
for( i = n - 4; i >= 0; i-- ) {
    out[ i ] = -out[ i + 2 ]*div0 + -out[ i + 1 ]*div1 + poly[ i + 2 ];
}
for( i = 0; i < n; i++ ) {
    poly[ i ] = out[ i ];
}
poly[ n - 1 ] = 0.0;
poly[ n ] = 0.0;
free( out );
retval = 0;

This entirely bypasses the need to form a complex coefficient polynomial, and deflates both complex roots in a single pass.