glaurung

Fixture 181

compensated summation

C · 5 functions · 4 lanes · 18 of 20 function-lanes behave identically

2 of 4 lanes have a function that returns a different result after decompilation: clang-O0 (4/5), clang-O2 (4/5).

Kahan/Neumaier compensated summation at binary64.

COVERAGE TARGET: Subsd. The corpus reaches addsd, mulsd and divsd through 172 and 175, but nothing subtracts at binary64 — and subtraction is the one scalar operation whose lowering cannot be checked by symmetry, since a - b and b - a differ. Compensated summation is where a real program subtracts doubles: the compensation term IS (sum - old) - value, and the whole point of the algorithm is that this quantity is not zero.

It is also a recovery test with teeth. The compensation is exactly the rounding error the naive sum discards, so a recovery that reassociates the arithmetic, keeps an intermediate at the wrong width, or swaps a subtraction's operands returns the NAIVE sum — a number that is close enough to look right and is bit-exactly wrong, which is what the differential compares.

Inputs come in as float and widen, so the terms are exactly representable and the only inexactness in the answer is the one the algorithm compensates.

tests/decompiler_fixtures/src/181_compensated_summation.c source 
#include <stdint.h>

/* Kahan/Neumaier compensated summation at binary64.
 *
 * COVERAGE TARGET: `Subsd`. The corpus reaches `addsd`, `mulsd` and `divsd`
 * through 172 and 175, but nothing subtracts at binary64 — and subtraction is
 * the one scalar operation whose lowering cannot be checked by symmetry, since
 * `a - b` and `b - a` differ. Compensated summation is where a real program
 * subtracts doubles: the compensation term IS `(sum - old) - value`, and the
 * whole point of the algorithm is that this quantity is not zero.
 *
 * It is also a recovery test with teeth. The compensation is exactly the
 * rounding error the naive sum discards, so a recovery that reassociates the
 * arithmetic, keeps an intermediate at the wrong width, or swaps a subtraction's
 * operands returns the NAIVE sum — a number that is close enough to look right
 * and is bit-exactly wrong, which is what the differential compares.
 *
 * Inputs come in as `float` and widen, so the terms are exactly representable
 * and the only inexactness in the answer is the one the algorithm compensates. */

#define FP181_TERM_LIMIT 32

/* The naive sum, for the difference to be measured against. */
__attribute__((noinline)) double naive_sum_f64(const float *terms,
                                               int32_t count) {
    double total = 0.0;
    int32_t index;
    if (terms == 0 || count < 0 || count > FP181_TERM_LIMIT) {
        return 0.0;
    }
    for (index = 0; index < count; ++index) {
        total += (double)terms[index];
    }
    return total;
}

/* Kahan summation. `compensation` is carried by SUBTRACTION at binary64 and is
 * the only reason this differs from `naive_sum_f64`. */
__attribute__((noinline)) double kahan_sum_f64(const float *terms,
                                               int32_t count) {
    double total = 0.0;
    double compensation = 0.0;
    int32_t index;
    if (terms == 0 || count < 0 || count > FP181_TERM_LIMIT) {
        return 0.0;
    }
    for (index = 0; index < count; ++index) {
        double adjusted = (double)terms[index] - compensation;
        double next = total + adjusted;
        /* `(next - total)` is the part of `adjusted` that survived rounding;
         * subtracting `adjusted` leaves exactly the part that did not. Both
         * subtractions are binary64 and neither may be reassociated away. */
        compensation = (next - total) - adjusted;
        total = next;
    }
    return total;
}

/* The compensation itself, as a value: nonzero exactly when the naive sum lost
 * a digit. Returned as an integer so the verdict does not depend on how the
 * harness marshals a tiny double. */
__attribute__((noinline)) int32_t summation_disagrees(const float *terms,
                                                      int32_t count) {
    double naive = naive_sum_f64(terms, count);
    double kahan = kahan_sum_f64(terms, count);
    return (naive == kahan) ? 0 : 1;
}

/* A single compensated step, isolated. Three binary64 subtractions with no loop
 * around them, so a failure here is a lowering bug and not a loop-shape one. */
__attribute__((noinline)) double compensation_of_step(double total,
                                                      double value) {
    double next = total + value;
    return (next - total) - value;
}

/* Difference of products — the determinant shape, at binary64, where operand
 * order is observable. `a*d - b*c` and `b*c - a*d` are negatives of each other,
 * so a swapped `subsd` is caught by any input whose result is not zero. */
__attribute__((noinline)) double difference_of_products(double a, double b,
                                                        double c, double d) {
    return a * d - b * c;
}

Recovered C

Generated by glaurung decompile --style decbench at b47f6b43. baseline.json records the result after recompiling the C and calling it beside the original with seeded inputs.

clang -O0

4/5
compensation_of_step pass 8 lines 
// glaurung: compensation_of_step @ 0x12e0
double compensation_of_step(double arg0, double arg1) {
    double next;
    // x86-64 prologue: save rbp
    next = (arg0 + arg1);
    // x86-64 epilogue: restore rbp
    return ((next - arg0) - arg1);
}
difference_of_products fail 6 lines 
// glaurung: difference_of_products @ 0x1310
double difference_of_products(double arg0, double arg1, double arg2, double arg3) {
    // x86-64 prologue: save rbp
    // x86-64 epilogue: restore rbp
    return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(((arg1 * arg2) ^ (-0x7fffffffffffffffLL - 1LL))) }).value;
}
kahan_sum_f64 pass 36 lines 
// glaurung: kahan_sum_f64 @ 0x11b0
double kahan_sum_f64(const float * arg0, int32_t arg1) {
    double total;
    double compensation;
    int index;
    double adjusted;
    double next;
    double local_8;
    // x86-64 prologue: save rbp
    total = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    compensation = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    if ((arg0 == 0)) {
        local_8 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
        // x86-64 epilogue: restore rbp
        return local_8;
    }
    if (((long)(arg1) < 0)) {
        local_8 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
        // x86-64 epilogue: restore rbp
        return local_8;
    }
    if (((((unsigned long)((unsigned int)(arg1)) == 32) | ((long)(arg1) < 32)) == 0)) {
        local_8 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
        // x86-64 epilogue: restore rbp
        return local_8;
    }
    for (index = 0; (index < arg1); index++) {
        adjusted = ((double)(*(float *)(((long)arg0 + ((long)(index) * 4)))) - compensation);
        next = (total + adjusted);
        compensation = ((next - total) - adjusted);
        total = next;
    }
    local_8 = total;
    // x86-64 epilogue: restore rbp
    return local_8;
}
naive_sum_f64 pass 29 lines 
// glaurung: naive_sum_f64 @ 0x1120
double naive_sum_f64(const float * arg0, int32_t arg1) {
    double total;
    int index;
    double local_8;
    // x86-64 prologue: save rbp
    total = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    if ((arg0 == 0)) {
        local_8 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
        // x86-64 epilogue: restore rbp
        return local_8;
    }
    if (((long)(arg1) < 0)) {
        local_8 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
        // x86-64 epilogue: restore rbp
        return local_8;
    }
    if (((((unsigned long)((unsigned int)(arg1)) == 32) | ((long)(arg1) < 32)) == 0)) {
        local_8 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
        // x86-64 epilogue: restore rbp
        return local_8;
    }
    for (index = 0; (index < arg1); index++) {
        total = ((double)(*(float *)(((long)arg0 + ((long)(index) * 4)))) + total);
    }
    local_8 = total;
    // x86-64 epilogue: restore rbp
    return local_8;
}
summation_disagrees pass 16 lines 
// glaurung: summation_disagrees @ 0x1280
int32_t summation_disagrees(const float * arg0, int32_t arg1) {
    extern double kahan_sum_f64(int *, int);
    extern double naive_sum_f64(int *, int);
    double naive;
    double kahan;
    double var0;
    double var2;
    // x86-64 prologue: save rbp, frame 32 bytes
    var0 = naive_sum_f64((int *)(arg0), (unsigned long)((unsigned int)(arg1)));
    naive = var0;
    var2 = kahan_sum_f64((int *)(arg0), (unsigned long)((unsigned int)(arg1)));
    kahan = var2;
    // x86-64 epilogue: restore rbp
    return ((((naive == kahan) || ((naive != naive) || (kahan != kahan))) && (((naive != naive) | (kahan != kahan)) == 0)) ? 0 : 1);
}

clang -O2

4/5
compensation_of_step pass 8 lines 
// glaurung: compensation_of_step @ 0x12c0
double compensation_of_step(double arg0, double arg1) {
    double next;
    double var15;
    var15 = (((arg0 + arg1) - arg0) - arg1);
    arg0 = var15;
    return var15;
}
difference_of_products pass 4 lines 
// glaurung: difference_of_products @ 0x12e0
double difference_of_products(double arg0, double arg1, double arg2, double arg3) {
    return ((arg0 * arg3) - (arg1 * arg2));
}
kahan_sum_f64 pass 53 lines 
// glaurung: kahan_sum_f64 @ 0x11d0
double kahan_sum_f64(const float * arg0, int32_t arg1) {
    int index;
    double adjusted;
    double compensation;
    double next;
    double total;
    long ret;
    double var16;
    double var34;
    double var42;
    long var5;
    double var64;
    long var8;
    long var80;
    double var81;
    double var82;
    double var9;
    ret = 0;
    if (((unsigned long)((unsigned long)((unsigned int)((arg1 - 1)))) <= (unsigned long)(31))) {
        if ((arg0 == 0)) {
            return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
        }
        var5 = (unsigned long)((unsigned int)(arg1));
        if (((unsigned long)((unsigned int)(arg1)) != 1)) {
            var8 = (unsigned long)((unsigned int)(((unsigned long)((unsigned int)(var5)) & -2)));
            var9 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
            index = 0;
            var16 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
            do {
                var34 = ((double)(*(float *)(((long)arg0 + index * 4))) - var16);
                var42 = (var9 + var34);
                var64 = ((double)(*(float *)(((long)arg0 + index * 4 + 0x4))) - ((var42 - var9) - var34));
                ret = ((union { unsigned long long bits; double value; }){ .value = (var42 + var64) }).bits;
                var16 = ((((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value - var42) - var64);
                index = (index + 2);
                var9 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
                var80 = (unsigned long)((unsigned int)(index));
                var81 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
                var82 = var16;
            } while ((var8 != index));
        } else {
            ret = 0;
            var82 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
            var80 = 0;
            var81 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
        }
        if (((unsigned long)((unsigned char)((var5 & 1))) != 0)) {
            ret = ((union { unsigned long long bits; double value; }){ .value = (var81 + ((double)(*(float *)(((long)arg0 + var80 * 4))) - var82)) }).bits;
        }
    }
    return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
}
naive_sum_f64 pass 51 lines 
// glaurung: naive_sum_f64 @ 0x1120
double naive_sum_f64(const float * arg0, int32_t arg1) {
    int index;
    double total;
    long ret;
    long var11;
    double var12;
    long var5;
    long var67;
    double var68;
    long var75;
    long var77;
    double var78;
    long var9;
    ret = 0;
    if (((unsigned long)((unsigned long)((unsigned int)((arg1 - 1)))) <= (unsigned long)(31))) {
        if ((arg0 == 0)) {
            return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
        }
        var5 = (unsigned long)((unsigned int)(arg1));
        var9 = (unsigned long)((unsigned int)(((unsigned long)((unsigned int)(arg1)) & 3)));
        if (((unsigned long)(3) <= (unsigned long)(((unsigned long)((unsigned int)(arg1)) - 1)))) {
            var11 = (unsigned long)((unsigned int)((var5 & -4)));
            var12 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
            index = 0;
            do {
                ret = ((union { unsigned long long bits; double value; }){ .value = ((double)(*(float *)(((long)arg0 + index * 4 + 0xc))) + ((double)(*(float *)(((long)arg0 + index * 4 + 0x8))) + ((double)(*(float *)(((long)arg0 + index * 4 + 0x4))) + ((double)(*(float *)(((long)arg0 + index * 4))) + var12)))) }).bits;
                index = (index + 4);
                var12 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
                var67 = (unsigned long)((unsigned int)(index));
                var68 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
            } while ((var11 != index));
        } else {
            ret = 0;
            var67 = 0;
            var68 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
        }
        if ((var9 == 0)) {
            return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
        }
        var75 = (long)(((long)arg0 + (var67 * 4)));
        var77 = 0;
        var78 = var68;
        do {
            ret = ((union { unsigned long long bits; double value; }){ .value = (var78 + (double)(*(float *)((var75 + var77 * 4)))) }).bits;
            var77 = (var77 + 1);
            var78 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
        } while ((var9 != var77));
    }
    return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
}
summation_disagrees fail 20 lines 
// glaurung: summation_disagrees @ 0x1280
__attribute__((no_stack_protector)) int32_t summation_disagrees(const float * arg0, int32_t arg1) {
    extern double kahan_sum_f64(int *, int);
    extern double naive_sum_f64(int *, int);
    double kahan;
    unsigned char local_18[24];
    long var0;
    int * var1;
    double var2;
    double var4;
    // x86-64 prologue: save callee registers, frame 24 bytes
    var0 = (unsigned long)((unsigned int)(arg1));
    var1 = (int *)arg0;
    var2 = naive_sum_f64((int *)(arg0), arg1);
    *(double *)(&local_18[0]) = var2;
    var4 = kahan_sum_f64(var1, (unsigned long)((unsigned int)(var0)));
    /* asm: cmpsd */
    // x86-64 epilogue: restore callee registers
    return (unsigned int)((var4 & 1));
}

gcc -O0

5/5
compensation_of_step pass 11 lines 
// glaurung: compensation_of_step @ 0x12d8
double compensation_of_step(double arg0, double arg1) {
    double next;
    double var16;
    // x86-64 prologue: save rbp
    next = (arg0 + arg1);
    var16 = ((next - arg0) - arg1);
    arg0 = var16;
    // x86-64 epilogue: restore rbp
    return var16;
}
difference_of_products pass 9 lines 
// glaurung: difference_of_products @ 0x1314
double difference_of_products(double arg0, double arg1, double arg2, double arg3) {
    double var16;
    // x86-64 prologue: save rbp
    var16 = ((arg0 * arg3) - (arg1 * arg2));
    arg0 = var16;
    // x86-64 epilogue: restore rbp
    return var16;
}
kahan_sum_f64 pass 31 lines 
// glaurung: kahan_sum_f64 @ 0x11ba
double kahan_sum_f64(const float * arg0, int32_t arg1) {
    double total;
    double compensation;
    int index;
    double adjusted;
    double next;
    // x86-64 prologue: save rbp
    total = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    compensation = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    if ((arg0 == 0)) {
        // x86-64 epilogue: restore rbp
        return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    }
    if (((long)(arg1) < 0)) {
        // x86-64 epilogue: restore rbp
        return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    }
    if (((((unsigned long)((unsigned int)(arg1)) == 32) | ((long)(arg1) < 32)) == 0)) {
        // x86-64 epilogue: restore rbp
        return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    }
    for (index = 0; (index < arg1); index++) {
        adjusted = ((double)(*(float *)(((long)arg0 + (0 + ((long)((int)((unsigned long)((unsigned int)(index)))) * 4))))) - compensation);
        next = (total + adjusted);
        compensation = ((next - total) - adjusted);
        total = next;
    }
    // x86-64 epilogue: restore rbp
    return total;
}
naive_sum_f64 pass 24 lines 
// glaurung: naive_sum_f64 @ 0x1139
double naive_sum_f64(const float * arg0, int32_t arg1) {
    double total;
    int index;
    // x86-64 prologue: save rbp
    total = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    if ((arg0 == 0)) {
        // x86-64 epilogue: restore rbp
        return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    }
    if (((long)(arg1) < 0)) {
        // x86-64 epilogue: restore rbp
        return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    }
    if (((((unsigned long)((unsigned int)(arg1)) == 32) | ((long)(arg1) < 32)) == 0)) {
        // x86-64 epilogue: restore rbp
        return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    }
    for (index = 0; (index < arg1); index++) {
        total = ((double)(*(float *)(((long)arg0 + (0 + ((long)((int)((unsigned long)((unsigned int)(index)))) * 4))))) + total);
    }
    // x86-64 epilogue: restore rbp
    return total;
}
summation_disagrees pass 16 lines 
// glaurung: summation_disagrees @ 0x126d
int32_t summation_disagrees(const float * arg0, int32_t arg1) {
    extern double kahan_sum_f64(int *, int);
    extern double naive_sum_f64(int *, int);
    double naive;
    double kahan;
    double var2;
    double var7;
    // x86-64 prologue: save rbp, frame 32 bytes
    var2 = naive_sum_f64((int *)(arg0), (unsigned long)((unsigned int)(arg1)));
    naive = var2;
    var7 = kahan_sum_f64((int *)(arg0), (unsigned long)((unsigned int)(arg1)));
    kahan = var7;
    // x86-64 epilogue: restore rbp
    return (unsigned int)((unsigned char)((((((naive == kahan) | ((naive != naive) | (kahan != kahan))) == 0) ? 1 : (((naive != naive) | (kahan != kahan)) & 255)) & 255)));
}

gcc -O2

5/5
compensation_of_step pass 5 lines 
// glaurung: compensation_of_step @ 0x1230
double compensation_of_step(double arg0, double arg1) {
    double next;
    return (((arg0 + arg1) - arg0) - arg1);
}
difference_of_products pass 4 lines 
// glaurung: difference_of_products @ 0x1250
double difference_of_products(double arg0, double arg1, double arg2, double arg3) {
    return ((arg0 * arg3) - (arg1 * arg2));
}
kahan_sum_f64 pass 37 lines 
// glaurung: kahan_sum_f64 @ 0x1180
double kahan_sum_f64(const float * arg0, int32_t arg1) {
    double adjusted;
    double compensation;
    int index;
    double next;
    double total;
    long ret;
    double var0;
    long var10;
    long var11;
    double var16;
    long var17;
    double var27;
    var0 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
    if ((arg0 == 0)) {
        ret = ((union { unsigned long long bits; double value; }){ .value = var0 }).bits;
        return var0;
    }
    var10 = (unsigned long)((unsigned int)((arg1 - 1)));
    if (((unsigned long)(31) < (unsigned long)((unsigned long)((unsigned int)(var10))))) {
        ret = ((union { unsigned long long bits; double value; }){ .value = var0 }).bits;
        return var0;
    }
    var11 = (long)((((long)arg0 + (var10 * 4)) + 4));
    var16 = var0;
    var17 = (long)arg0;
    do {
        var27 = var0;
        var17 = (var17 + 4);
        ret = ((union { unsigned long long bits; double value; }){ .value = ((double)(*(float *)((var17 - 0x4))) - var16) }).bits;
        var0 = (var0 + ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value);
        var16 = ((var0 - var27) - ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value);
    } while ((var11 != var17));
    ret = ((union { unsigned long long bits; double value; }){ .value = var0 }).bits;
    return var0;
}
naive_sum_f64 pass 26 lines 
// glaurung: naive_sum_f64 @ 0x1140
double naive_sum_f64(const float * arg0, int32_t arg1) {
    int index;
    double total;
    long ret;
    long var5;
    long var6;
    long var7;
    double var8;
    ret = 0;
    if ((arg0 != 0)) {
        var5 = (unsigned long)((unsigned int)((arg1 - 1)));
        if (((unsigned long)(31) < (unsigned long)((unsigned long)((unsigned int)(var5))))) {
            return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
        }
        var6 = (long)((((long)arg0 + (var5 * 4)) + 4));
        var7 = (long)arg0;
        var8 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(0) }).value;
        do {
            var7 = (var7 + 4);
            ret = ((union { unsigned long long bits; double value; }){ .value = (var8 + (double)(*(float *)((var7 - 0x4)))) }).bits;
            var8 = ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
        } while ((var7 != var6));
    }
    return ((union { unsigned long long bits; double value; }){ .bits = (unsigned long long)(ret) }).value;
}
summation_disagrees pass 19 lines 
// glaurung: summation_disagrees @ 0x11e0
int32_t summation_disagrees(const float * arg0, int32_t arg1) {
    extern double kahan_sum_f64(int *, int);
    extern double naive_sum_f64(int *, int);
    double kahan;
    double naive;
    double local_20;
    long var0;
    int * var1;
    double var2;
    double var4;
    var0 = (unsigned long)((unsigned int)(arg1));
    var1 = (int *)arg0;
    var2 = naive_sum_f64((int *)(arg0), arg1);
    local_20 = var2;
    var4 = kahan_sum_f64(var1, (unsigned long)((unsigned int)(var0)));
    // x86-64 epilogue: tear down frame
    return ((((local_20 == var4) | ((local_20 != local_20) | (var4 != var4))) == 0) ? 1 : (((local_20 != local_20) | (var4 != var4)) & 255));
}

← 213 fixtures