From 5d069df45240102db3fd1b6bafbcffc79353f347 Mon Sep 17 00:00:00 2001
From: Michael Tautschnig <tautschn@amazon.com>
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           |  8 ++
 src/solvers/flattening/boolbv_floatbv_op.cpp | 90 ++++++++++++--------
 3 files changed, 79 insertions(+), 34 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 <complex.h>
+
+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..9efefbc7362
--- /dev/null
+++ b/regression/cbmc/complex2/test.desc
@@ -0,0 +1,8 @@
+CORE
+main.c
+
+^EXIT=0$
+^SIGNAL=0$
+^VERIFICATION SUCCESSFUL$
+--
+^warning: ignoring
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<size; i++)
+      bvt result_real, result_imag;
+
+      if(expr.id() == ID_floatbv_plus || expr.id() == ID_floatbv_minus)
+      {
+        result_real = float_utils.add_sub(
+          lhs_real, rhs_real, expr.id() == ID_floatbv_minus);
+        result_imag = float_utils.add_sub(
+          lhs_imag, rhs_imag, expr.id() == ID_floatbv_minus);
+      }
+      else if(expr.id() == ID_floatbv_mult)
+      {
+        // Could be optimised to just three multiplications with more additions
+        // instead, but then we'd have to worry about the impact of possible
+        // overflows. So we use the naive approach for now:
+        result_real = float_utils.add_sub(
+          float_utils.mul(lhs_real, rhs_real),
+          float_utils.mul(lhs_imag, rhs_imag),
+          true);
+        result_imag = float_utils.add_sub(
+          float_utils.mul(lhs_real, rhs_imag),
+          float_utils.mul(lhs_imag, rhs_real),
+          false);
+      }
+      else if(expr.id() == ID_floatbv_div)
       {
-        bvt lhs_sub_bv, rhs_sub_bv, sub_result_bv;
-
-        lhs_sub_bv.assign(
-          lhs_as_bv.begin() + i * sub_width,
-          lhs_as_bv.begin() + (i + 1) * sub_width);
-        rhs_sub_bv.assign(
-          rhs_as_bv.begin() + i * sub_width,
-          rhs_as_bv.begin() + (i + 1) * sub_width);
-
-        if(expr.id()==ID_floatbv_plus)
-          sub_result_bv = float_utils.add_sub(lhs_sub_bv, rhs_sub_bv, false);
-        else if(expr.id()==ID_floatbv_minus)
-          sub_result_bv = float_utils.add_sub(lhs_sub_bv, rhs_sub_bv, true);
-        else if(expr.id()==ID_floatbv_mult)
-          sub_result_bv = float_utils.mul(lhs_sub_bv, rhs_sub_bv);
-        else if(expr.id()==ID_floatbv_div)
-          sub_result_bv = float_utils.div(lhs_sub_bv, rhs_sub_bv);
-        else
-          UNREACHABLE;
-
-        INVARIANT(
-          sub_result_bv.size() == sub_width,
-          "we constructed a new complex of the right size");
-        INVARIANT(
-          i * sub_width + sub_width - 1 < result_bv.size(),
-          "the sub-bitvector fits into the result bitvector");
-        std::copy(
-          sub_result_bv.begin(),
-          sub_result_bv.end(),
-          result_bv.begin() + i * sub_width);
+        bvt numerator_real = float_utils.add_sub(
+          float_utils.mul(lhs_real, rhs_real),
+          float_utils.mul(lhs_imag, rhs_imag),
+          false);
+        bvt numerator_imag = float_utils.add_sub(
+          float_utils.mul(lhs_imag, rhs_real),
+          float_utils.mul(lhs_real, rhs_imag),
+          true);
+
+        bvt denominator = float_utils.add_sub(
+          float_utils.mul(rhs_real, rhs_real),
+          float_utils.mul(rhs_imag, rhs_imag),
+          false);
+
+        result_real = float_utils.div(numerator_real, denominator);
+        result_imag = float_utils.div(numerator_imag, denominator);
       }
+      else
+        UNREACHABLE;
+
+      bvt result_bv = std::move(result_real);
+      result_bv.reserve(width);
+      result_bv.insert(
+        result_bv.end(),
+        std::make_move_iterator(result_imag.begin()),
+        std::make_move_iterator(result_imag.end()));
 
       return result_bv;
     }