From 6a6c4807451bb89890f0a7d6a8902f00daedfd8e Mon Sep 17 00:00:00 2001 From: Michael Tautschnig Date: Thu, 11 Jul 2024 10:38:40 +0000 Subject: [PATCH] Fix multiplication and division of complex numbers Multiplication and division of complex numbers are not just pointwise applications of those operations. Fixes: #8375 --- regression/cbmc/complex2/main.c | 15 ++++ regression/cbmc/complex2/test.desc | 10 +++ src/goto-programs/remove_complex.cpp | 56 ++++++++++-- src/solvers/flattening/boolbv_floatbv_op.cpp | 90 ++++++++++++-------- 4 files changed, 131 insertions(+), 40 deletions(-) create mode 100644 regression/cbmc/complex2/main.c create mode 100644 regression/cbmc/complex2/test.desc diff --git a/regression/cbmc/complex2/main.c b/regression/cbmc/complex2/main.c new file mode 100644 index 00000000000..d7d128319ff --- /dev/null +++ b/regression/cbmc/complex2/main.c @@ -0,0 +1,15 @@ +#include + +int main() +{ + char choice; + float re = choice ? 1.3f : 2.1f; // a non-constant well-behaved float + float complex z1 = I + re; + float complex z2 = z1 * z1; + float complex expected = 2 * I * re + re * re - 1; // (a+i)^2 = 2ai + a^2 - 1 + float complex actual = + re * re + I; // (a1 + b1*i)*(a2 + b2*i) = (a1*a2 + b1*b2*i) + __CPROVER_assert(z2 == expected, "right"); + __CPROVER_assert(z2 == actual, "wrong"); + return 0; +} diff --git a/regression/cbmc/complex2/test.desc b/regression/cbmc/complex2/test.desc new file mode 100644 index 00000000000..557a34ce2a2 --- /dev/null +++ b/regression/cbmc/complex2/test.desc @@ -0,0 +1,10 @@ +CORE no-new-smt +main.c + +^\[main.assertion.1\] line 12 right: SUCCESS$ +^\[main.assertion.2\] line 13 wrong: FAILURE$ +^VERIFICATION FAILED$ +^EXIT=10$ +^SIGNAL=0$ +-- +^warning: ignoring diff --git a/src/goto-programs/remove_complex.cpp b/src/goto-programs/remove_complex.cpp index c983692cfb1..45d470e2e68 100644 --- a/src/goto-programs/remove_complex.cpp +++ b/src/goto-programs/remove_complex.cpp @@ -127,13 +127,10 @@ static void remove_complex(exprt &expr) if(expr.type().id()==ID_complex) { - if(expr.id()==ID_plus || expr.id()==ID_minus || - expr.id()==ID_mult || expr.id()==ID_div) + if(expr.id() == ID_plus || expr.id() == ID_minus) { - // FIXME plus and mult are defined as n-ary operations - // rather than binary. This code assumes that they - // can only have exactly 2 operands, and it is not clear - // that it is safe to do so in this context + // plus and mult are n-ary expressions, but front-ends currently ensure + // that we only see them as binary ones PRECONDITION(expr.operands().size() == 2); // do component-wise: // x+y -> complex(x.r+y.r,x.i+y.i) @@ -153,6 +150,53 @@ static void remove_complex(exprt &expr) expr=struct_expr; } + else if(expr.id() == ID_mult) + { + // plus and mult are n-ary expressions, but front-ends currently ensure + // that we only see them as binary ones + PRECONDITION(expr.operands().size() == 2); + exprt lhs_real = complex_member(to_binary_expr(expr).op0(), ID_real); + exprt lhs_imag = complex_member(to_binary_expr(expr).op0(), ID_imag); + exprt rhs_real = complex_member(to_binary_expr(expr).op1(), ID_real); + exprt rhs_imag = complex_member(to_binary_expr(expr).op1(), ID_imag); + + struct_exprt struct_expr{ + {minus_exprt{ + mult_exprt{lhs_real, rhs_real}, mult_exprt{lhs_imag, rhs_imag}}, + plus_exprt{ + mult_exprt{lhs_imag, rhs_real}, mult_exprt{lhs_real, rhs_imag}}}, + expr.type()}; + + struct_expr.op0().add_source_location() = expr.source_location(); + struct_expr.op1().add_source_location() = expr.source_location(); + + expr = struct_expr; + } + else if(expr.id() == ID_div) + { + exprt lhs_real = complex_member(to_binary_expr(expr).op0(), ID_real); + exprt lhs_imag = complex_member(to_binary_expr(expr).op0(), ID_imag); + exprt rhs_real = complex_member(to_binary_expr(expr).op1(), ID_real); + exprt rhs_imag = complex_member(to_binary_expr(expr).op1(), ID_imag); + + plus_exprt numerator_real{ + mult_exprt{lhs_real, rhs_real}, mult_exprt{lhs_imag, rhs_imag}}; + minus_exprt numerator_imag{ + mult_exprt{lhs_imag, rhs_real}, mult_exprt{lhs_real, rhs_imag}}; + + plus_exprt denominator{ + mult_exprt{rhs_real, rhs_real}, mult_exprt{rhs_imag, rhs_imag}}; + + struct_exprt struct_expr{ + {div_exprt{numerator_real, denominator}, + div_exprt{numerator_imag, denominator}}, + expr.type()}; + + struct_expr.op0().add_source_location() = expr.source_location(); + struct_expr.op1().add_source_location() = expr.source_location(); + + expr = struct_expr; + } else if(expr.id()==ID_unary_minus) { auto const &unary_minus_expr = to_unary_minus_expr(expr); diff --git a/src/solvers/flattening/boolbv_floatbv_op.cpp b/src/solvers/flattening/boolbv_floatbv_op.cpp index e42479573f5..7a18e0ca5f4 100644 --- a/src/solvers/flattening/boolbv_floatbv_op.cpp +++ b/src/solvers/flattening/boolbv_floatbv_op.cpp @@ -131,44 +131,66 @@ bvt boolbvt::convert_floatbv_op(const ieee_float_op_exprt &expr) sub_width > 0 && width % sub_width == 0, "width of a complex subtype must be positive and evenly divide the " "width of the complex expression"); + DATA_INVARIANT( + sub_width * 2 == width, "a complex type consists of exactly two parts"); + + bvt lhs_real{lhs_as_bv.begin(), lhs_as_bv.begin() + sub_width}; + bvt rhs_real{rhs_as_bv.begin(), rhs_as_bv.begin() + sub_width}; - std::size_t size=width/sub_width; - bvt result_bv; - result_bv.resize(width); + bvt lhs_imag{lhs_as_bv.begin() + sub_width, lhs_as_bv.end()}; + bvt rhs_imag{rhs_as_bv.begin() + sub_width, rhs_as_bv.end()}; - for(std::size_t i=0; i