I tested MUX 2.1 on a concrete machine and found the same behavior, so
the Hercules emulator is vindicated. It looks like a compiler bug.

Without the calls to irrelevant(), AddDoubles(2, {1.6, .1}) returns 0.2.
With the calls to irrelevant, it returns 1.7.

static void irrelevant(void)
{
    static volatile int x;
    x++;
}

static double AddWithError(double& err, double a, double b)
{
    irrelevant();
    double sum = a+b;
    irrelevant();
    err = b-(sum-a);
    irrelevant();
    return sum;
}

// Double compensation method. Extended by Priest from Knuth and Kahan.
//
// Error of sum is less than 2*epsilon or 1 ulp except for very large n.
// Return the result that yields the shortest number of base-10 digits.
//
static double AddDoubles(int n, double pd[])
{
    qsort(pd, n, sizeof(double), f_comp_abs);
    irrelevant();
    double sum = pd[0];
    irrelevant();
    double sum_err = 0.0;
    irrelevant();
    int i;
    for (i = 1; i < n; i++)
    {
        irrelevant();
        double addend_err;
        irrelevant();
        double addend = AddWithError(addend_err, sum_err, pd[i]);
        irrelevant();
        double sum1_err;
        irrelevant();
        double sum1 = AddWithError(sum1_err, sum, addend);
        irrelevant();
        sum = AddWithError(sum_err, sum1, addend_err + sum1_err);
        irrelevant();
    }
    return NearestPretty(sum);
}

Stephen Dennis

Reply via email to