fixed bug reported by Trevor Spiteri

2 files changed, 23 insertions, 1 deletions

diff --git a/src/dot.c b/src/dot.c
index 6f655c0..1e270af 100644
--- a/src/dot.c
+++ b/src/dot.c
@@ -30,6 +30,7 @@ mpc_dot (mpc_ptr res, const mpc_ptr *x, const mpc_ptr *y,
   mpfr_ptr *t;
   mpfr_t *z;
   unsigned long i;
+  mpfr_t re_res;
 
   z = (mpfr_t *) malloc (2 * n * sizeof (mpfr_t));
   MPC_ASSERT(z != NULL);
@@ -57,7 +58,10 @@ mpc_dot (mpc_ptr res, const mpc_ptr *x, const mpc_ptr *y,
       mpfr_mul (z[n+i], mpc_imagref (x[i]), mpc_imagref (y[i]), MPFR_RNDZ);
       mpfr_neg (z[n+i], z[n+i], MPFR_RNDZ);
     }
-  inex_re = mpfr_sum (mpc_realref (res), t, 2 * n, MPC_RND_RE (rnd));
+  /* copy the real part in a temporary variable, since it might be in the
+     input array */
+  mpfr_init2 (re_res, mpfr_get_prec (mpc_realref (res)));
+  inex_re = mpfr_sum (re_res, t, 2 * n, MPC_RND_RE (rnd));
   /* we then store in z[i] the value of Re(x[i])*Im(y[i])
      and in z[n+i] that of Im(x[i])*Re(y[i]) */
   for (i = 0; i < n; i++)
@@ -72,6 +76,8 @@ mpc_dot (mpc_ptr res, const mpc_ptr *x, const mpc_ptr *y,
       mpfr_mul (z[n+i], mpc_imagref (x[i]), mpc_realref (y[i]), MPFR_RNDZ);
     }
   inex_im = mpfr_sum (mpc_imagref (res), t, 2 * n, MPC_RND_IM (rnd));
+  mpfr_swap (mpc_realref (res), re_res);
+  mpfr_clear (re_res);
   for (i = 0; i < 2 * n; i++)
     mpfr_clear (z[i]);
   free (t);
diff --git a/tests/tdot.c b/tests/tdot.c
index c6d1f39..ce31482 100644
--- a/tests/tdot.c
+++ b/tests/tdot.c
@@ -68,11 +68,27 @@ check_special (void)
   mpc_clear (res);
 }
 
+/* bug reported by Trevor Spiteri */
+static void
+bug20200717 (void)
+{
+  mpc_t a;
+  mpc_ptr p[1];
+  mpc_init2 (a, 53);
+  mpc_set_ui_ui (a, 1, 2, MPC_RNDNN);
+  p[0] = a;
+  mpc_dot (a, p, p, 1, MPC_RNDNN);
+  MPC_ASSERT (mpfr_cmp_si (mpc_realref (a), -3) == 0);
+  MPC_ASSERT (mpfr_cmp_ui (mpc_imagref (a), 4) == 0);
+  mpc_clear (a);
+}
+
 int
 main (void)
 {
   test_start ();
 
+  bug20200717 ();
   check_special ();
 
   test_end ();