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