/*
 * Get_lsp_pol
 *
 *
 * Parameters:
 *    lsp               I: line spectral frequencies
 *    f                 O: polynomial F1(z) or F2(z)
 *
 * Function:
 *    Find the polynomial F1(z) or F2(z) from the LSPs.
 *
 *    F1(z) = product ( 1 - 2 lsp[i] z^-1 + z^-2 )
 *             i=0,2,4,6,8
 *    F2(z) = product   ( 1 - 2 lsp[i] z^-1 + z^-2 )
 *             i=1,3,5,7,9
 *
 *    where lsp[] is the LSP vector in the cosine domain.
 *
 *    The expansion is performed using the following recursion:
 *
 *    f[0] = 1
 *    b = -2.0 * lsp[0]
 *    f[1] = b
 *    for i=2 to 5 do
 *       b = -2.0 * lsp[2*i-2];
 *       f[i] = b*f[i-1] + 2.0*f[i-2];
 *       for j=i-1 down to 2 do
 *          f[j] = f[j] + b*f[j-1] + f[j-2];
 *       f[1] = f[1] + b;
 */

/**
 * Find the polynomial F1(z) or F2(z) from the lsps
 *
 * @param lsp               vector of lsps
 * @param f                 pointer to the polynomial F1(z) or F2(z)
 *
 * @return void
 */

static void lsp_poly(const int *lsp, int *f) {
    int b;
    int i, j;

    f[0] = 1;
    b = -(lsp[0]<<1);
    f[1] = b;
    for(i=2; i<6; i++) {
        b = lsp[2*i-2]<<1;
        f[i] = (f[i-2]<<1) - b*f[i-1];
        for(j=i-1; j>1; j--) {
            f[j] += f[j-2] - b*f[j-1];
        }
        f[1] -= b;
    }
}

/*

mine:
f[0] = 1;

f[1] =           - (lsp[0]<<1);

f[2] = (f[0]<<1) - (lsp[2]<<1)*f[1];
f[1] -=            (lsp[2]<<1);

f[3] = (f[1]<<1) - (lsp[4]<<1)*f[2];
f[2] += f[0]     - (lsp[4]<<1)*f[1];
f[1] -=            (lsp[4]<<1);

f[4] = (f[2]<<1) - (lsp[6]<<1)*f[3];
f[3] += f[1]     - (lsp[6]<<1)*f[2];
f[2] += f[2]     - (lsp[6]<<1)*f[1];
f[1] -=          - (lsp[6]<<1);


theirs:
f[2] = (f[0]<<1) - (( ((f[1]>>16)*lsp[2]) + (( (f[1] & 0xFFFE)*lsp[2] )>>16) )<<2);

*/

static void Get_lsp_pol(const int *lsp, int *f) {
   volatile int f0, f1, f2, f3, f4, f5;

   /* f[0] = 1.0; */
   f0 = 1<<24;

   /* f1 = *lsp * -1024; */
   f1 = -lsp[0] << 10;

   f2 = f0 << 1;
   f2 -= ( ( ( f1 >> 16 ) * lsp[2] ) + ( ( ( f1 & 0xFFFE ) * lsp[2] ) >> 16 ) ) << 2;
   f1 -= lsp[2] << 10;

   f3 = f1 << 1;
   f3 -= ( ( ( f2 >> 16 ) * lsp[4] ) + ( ( ( f2 & 0xFFFE ) * lsp[4] ) >> 16 ) ) << 2;
   f2 += f0;
   f2 -= ( ( ( f1 >> 16 ) * lsp[4] ) + ( ( ( f1 & 0xFFFE ) * lsp[4] ) >> 16 ) ) << 2;
   f1 -= lsp[4] << 10;

   f4 = f2 << 1;
   f4 -= ( ( ( f3 >> 16 ) * lsp[6] ) + ( ( ( f3 & 0xFFFE ) * lsp[6] ) >> 16 ) ) << 2;
   f3 += f1;
   f3 -= ( ( ( f2 >> 16 ) * lsp[6] ) + ( ( ( f2 & 0xFFFE ) * lsp[6] ) >> 16 ) ) << 2;
   f2 += f0;
   f2 -= ( ( ( f1 >> 16 ) * lsp[6] ) + ( ( ( f1 & 0xFFFE ) * lsp[6] ) >> 16 ) ) << 2;
   f1 -= lsp[6] << 10;

   f5 = f3 << 1;
   f5 -= ( ( ( f4 >> 16 ) * lsp[8] ) + ( ( ( f4 & 0xFFFE ) * lsp[8] ) >> 16 ) ) << 2;
   f4 += f2;
   f4 -= ( ( ( f3 >> 16 ) * lsp[8] ) + ( ( ( f3 & 0xFFFE ) * lsp[8] ) >> 16 ) ) << 2;
   f3 += f1;
   f3 -= ( ( ( f2 >> 16 ) * lsp[8] ) + ( ( ( f2 & 0xFFFE ) * lsp[8] ) >> 16 ) ) << 2;
   f2 += f0;
   f2 -= ( ( ( f1 >> 16 ) * lsp[8] ) + ( ( ( f1 & 0xFFFE ) * lsp[8] ) >> 16 ) ) << 2;
   f1 -= lsp[8] << 10;

   f[0] = f0;
   f[1] = f1;
   f[2] = f2;
   f[3] = f3;
   f[4] = f4;
   f[5] = f5;
   return;
}

int main() {
    const int lsp[10] = { 31216, 23530, 20571, 13690, -1869, -4576, -20408, -27784, -29826, 16777216 };
    int f1[6], f2[6];
    int i;

    Get_lsp_pol(lsp, f1);
    lsp_poly(lsp, f2);
    printf("\n\nTheirs :\n");
    for(i=0; i<6; i++) {
        printf("f[%d] = %d\n", i, f1[i]);
    }
    printf("\n\nMine :\n");
    for(i=0; i<6; i++) {
        printf("f[%d] = %d\n", i, f2[i]);
    }
    printf("\n\n");
}



And the output:

$ ./polytest 


Theirs :
f[0] = 16777216
f[1] = 323584
f[2] = -710000
f[3] = -5530272
f[4] = 5309548
f[5] = -1326264


Mine :
f[0] = 1
f[1] = 632
f[2] = -1119182183
f[3] = -1453367504
f[4] = 836835526
f[5] = 821401968



Obviously the first alteration I would make would be to 'correct' f[0] and f[1] but that doesn't appear to help too much as f[i] i>2 are completely wrong and looking at it I don't know how their expansion is supposed to be equivalent to the description. It looks like a load of rubbish to me...

Somehow I need to make:

f[2] = (f[0]<<1) - (lsp[2]<<1)*f[1];

equivalent to:

f[2] = (f[0]<<1) - (( ((f[1]>>16)*lsp[2]) + (( (f[1] & 0xFFFE)*lsp[2] )>>16) )<<2);

I don't know because I haven't seen it before but this almost looks like they're messing with bit/bytesex.