Skip to content

Fix template generation overflow.

Luke Davis requested to merge fix/stop_LALInspiralIIR_overflow into spiir-O4-EW-development

This function is run when trying to generate templates for our banks:

int XLALInspiralGenerateIIRSet(REAL8Vector *phase, double epsilon, double alpha, double beta, double padding)
{  
    int j = phase->length-1;
    while (j > 3 ) {
        /* Get second derivative term */
        double phase_ddot = (phase->data[j-2] - 2.0 * phase->data[j-1] + phase->data[j]) / (2.0 * LAL_PI);
        double phase_tdot = (phase->data[j-3] - 3.0 * phase->data[j-2] + 3.0 * phase->data[j-1] - phase->data[j]) / (2.0 * LAL_PI);
        phase_ddot = fabs(phase_ddot);
        phase_tdot = fabs(phase_tdot);

        if (phase_ddot < 0 || phase_ddot > 8 * epsilon){
                j = j - 1;
                continue;
        }

        int jstep = (int) fabs(floor(sqrt(2.0 * epsilon / fabs(phase_ddot))+0.5));
        int jstepThird = (int) fabs(floor(pow(6.0 * epsilon / fabs(phase_tdot),1./3)+0.5));
        jstep = abs(jstep);
        jstepThird = abs(jstepThird);

When phase_ddot/phase_tdot is very small but still above 0, jstep/jstepThird can get to float.INFINITY, which when casted to int will become INT_MIN. When abs() is run on INT_MIN (-2^31), the result should be 2^31 but because INT_MAX is 2^31-1, this overflows right back to INT_MIN.

If we instead set the bound phase_ddot and phase_tdot to some value close to but not too close to zero, we can stop this from happening.

if (phase_ddot < 1e-8 || phase_ddot > 8 * epsilon){
    j = j - 1;
    continue;
}
  
if (phase_tdot < 1e-8){
    j = j - 1;
    continue;
}

We could instead copy the fixed function into this repo and explicitly use that instead during PyIIR in python/spiirbank/_spiir_decomp.c .

Merge request reports