Floating point division simplification yields bad results

[leventeBajczi]

The expression simplification of FpDivExprs sometimes yield bad results:

theta/subprojects/common/core/src/main/java/hu/bme/mit/theta/core/utils/ExprSimplifier.java

Lines 1967 to 1984 in b379e4e

	private static Expr<FpType> simplifyFpDiv(final FpDivExpr expr, final Valuation val) {
	final Expr<FpType> leftOp = simplify(expr.getLeftOp(), val);
	final Expr<FpType> rightOp = simplify(expr.getRightOp(), val);
	if (leftOp instanceof FpLitExpr && rightOp instanceof FpLitExpr) {
	final FpLitExpr leftLit = (FpLitExpr) leftOp;
	final FpLitExpr rightLit = (FpLitExpr) rightOp;
	return leftLit.div(expr.getRoundingMode(), rightLit);
	}
	if (leftOp instanceof RefExpr && rightOp instanceof RefExpr) {
	if (leftOp.equals(rightOp)) {
	return FpUtils.bigFloatToFpLitExpr(new BigFloat(1.0f, FpUtils.getMathContext(expr.getType(), expr.getRoundingMode())), expr.getType());
	}
	}
	return expr.with(leftOp, rightOp);
	}

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.