Hi there,

I am coding a library to perform a statistical method. I have been using a -O2
optimisation level but I have just realised that switching it on and off
yields different results.

I have two functions cal_probability_var() and cal_heterozygosity(), the first
calls the second (see below). If I compile my code without optimisation I get
the expected result but if optimisations are allowed I get the wrong result.
My best guess is that the problem is that with optimisations the function
cal_heterozygosity() is inlined within cal_probability_var() (using gdb I
confirmed it because gdb behaves as expected for inlined functions).

To my surprise, if I swap the operators of the multiplication from:
  ln_2z =  log(2.0) * cal_heterozygosity(dist_matrix, table_type);
to:
  ln_2z = cal_heterozygosity(dist_matrix, table_type) * log(2.0);
then the problem disappears.

My main concerned is whether the way I have coded cal_heterozygosity()
conflicts with GCC or whether there is a bug in GCC. If it is indeed my
mistake I would like to know what is wrong with the coding as I may be using
the same standards to code other functions in my library. Or, it is possible
to get access to the source code once it has been optimised?.

static long double
cal_probability_var                (enumDistMatrix         *dist_matrix,
                                    enumTableType           table_type)
{
    gulong                     *table;
    guint                       i, j;
    long double                 ln_fij = 0.0;
    long double                 ln_2z;

    table = get_table(dist_matrix, table_type);

    ln_2z =  log(2.0) * cal_heterozygosity(dist_matrix, table_type);
//    ln_2z = cal_heterozygosity(dist_matrix, table_type) * log(2.0);

    for (i = 0 ; i < dist_matrix->k ; i++) {
        for (j = 0 ; j <= i ; j++) {
            ln_fij += fac_get_factorial(dist_matrix->factorials,
                                        table[OFFSET(i, j, dist_matrix->k)]);
        }
    }

    return (ln_2z - ln_fij);
}


static gulong
cal_heterozygosity                 (enumDistMatrix         *dist_matrix,
                                    enumTableType           table_type)
{
    gulong                     *table;
    guint                       i;
    gulong                      total = dist_matrix->n;

    table = get_table(dist_matrix, table_type);

    for (i = 0 ; i < dist_matrix->k ; i++) {
        total -= table[OFFSET(i, i, dist_matrix->k)];
    }

    return total;
}


Thanks

Hazael


Reply via email to