static inline int silk_is_lpc_stable(const int16_t lpc[16], int order)

{

    int k, j, DC_resp = 0;

    int32_t lpc32[2][16];       // Q24

    int totalinvgain = 1 << 30; // 1.0 in Q30

    int32_t *row = lpc32[0], *prevrow;



    /* initialize the first row for the Levinson recursion */

    for (k = 0; k < order; k++) {

        DC_resp += lpc[k];

        row[k] = lpc[k] * 4096;

    }



    if (DC_resp >= 4096)

        return 0;



    /* check if prediction gain pushes any coefficients too far */

    for (k = order - 1; 1; k--) {

        int rc;      // Q31; reflection coefficient

        int gaindiv; // Q30; inverse of the gain (the divisor)

        int gain;    // gain for this reflection coefficient

        int fbits;   // fractional bits used for the gain

        int error;   // Q29; estimate of the error of our partial estimate of 1/gaindiv



        if (FFABS(row[k]) > 16773022)

            return 0;



        rc      = -(row[k] * 128);

        gaindiv = (1 << 30) - MULH(rc, rc);



        totalinvgain = MULH(totalinvgain, gaindiv) << 2;

        if (k == 0)

            return (totalinvgain >= 107374);



        /* approximate 1.0/gaindiv */

        fbits = opus_ilog(gaindiv);

        gain  = ((1 << 29) - 1) / (gaindiv >> (fbits + 1 - 16)); // Q<fbits-16>

        error = (1 << 29) - MULL(gaindiv << (15 + 16 - fbits), gain, 16);

        gain  = ((gain << 16) + (error * gain >> 13));



        /* switch to the next row of the LPC coefficients */

        prevrow = row;

        row = lpc32[k & 1];



        for (j = 0; j < k; j++) {

            int x = prevrow[j] - ROUND_MULL(prevrow[k - j - 1], rc, 31);

            row[j] = ROUND_MULL(x, gain, fbits);

        }

    }

}
