Floating point division simplification yields bad results

[leventeBajczi]

The expression simplification of FpDivExprs sometimes yield bad results: https://github.com/ftsrg/theta/blob/b379e4ea0a77db26c040ad7798936fd1caeee7d3/subprojects/common/core/src/main/java/hu/bme/mit/theta/core/utils/ExprSimplifier.java#L1967-L1984

The underlying problem seems to be not with Theta, but with the BigFloat library. To see the wrong results:

val a = FpLitExpr.of(
        false,
        BvLitExpr.of(arrayOf(0, 0, 0, 0, 0, 0, 0, 0).map{it == 1}.toBooleanArray()),
        BvLitExpr.of(arrayOf(1,0,1,0,0,1,0,0,0,0,0,1,1,0,1,0,0,1,0,0,0,0,0).map{it == 1}.toBooleanArray()),
)
val b = FpLitExpr.of(
        false,
        BvLitExpr.of(arrayOf(0, 1, 1, 1, 1, 1, 1, 1).map{it == 1}.toBooleanArray()),
        BvLitExpr.of(arrayOf(1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,1).map{it == 1}.toBooleanArray()),
)

val expr = FpDivExpr.create(FpRoundingMode.RNE, a, b)

println("$a / $b [RNE] = ${ExprUtils.simplify(expr)}")

This will print

(#b0 #b00000000 #b10100100000110100100000) / (#b0 #b01111111 #b10011001100110011001101) [RNE] = (#b0 #b00000001 #b00000110100100000110100)

which is erroneous.

[leventeBajczi]

I tested this with native code calling MPFR functions (see code below), and it works surprisingly. This means that the java wrapper must be at fault.

#include <stdio.h>
#include <stdint.h>
#include <mpfr.h>

int main(int argc, char **argv) {
    mpfr_t a, b, result;
    uint32_t a_val, b_val;
    _Float32 res;

    a_val = 0x00520d20;
    b_val = 0x3fcccccd;

    mpfr_init2(a, 32);
    mpfr_init2(b, 32);
    mpfr_init2(result, 32);

    mpfr_set_flt(a, *((_Float32*)&a_val), MPFR_RNDN);
    mpfr_set_flt(b, *((_Float32*)&b_val), MPFR_RNDN);

    mpfr_div(result, a, b, MPFR_RNDN);

    res = mpfr_get_flt(result, MPFR_RNDN);
    printf("0x%08x / 0x%08x = 0x%08x\n", a_val, b_val, *((uint32_t*)&res));

    mpfr_clear(a);
    mpfr_clear(b);
    mpfr_clear(result);
    return 0;
}

Output: 0x00520d20 / 0x3fcccccd = 0x00334834. The correctness of this answer can be verified by the following native code:

#include <stdio.h>
#include <stdint.h>

int main(int argc, char **argv) {
    uint32_t a_val, b_val;
    _Float32 res;

    a_val = 0x00520d20;
    b_val = 0x3fcccccd;

    res = *((_Float32*)&a_val) / *((_Float32*)&b_val);
    printf("0x%08x / 0x%08x = 0x%08x\n", a_val, b_val, *((uint32_t*)&res));

    return 0;
}

By contrast, the output of the same program but from JVM via the wrapper results in 0x00834834.