[src/dot.c] implement a suggestion from Trevor Spiteri

https://lists.gforge.inria.fr/pipermail/mpc-discuss/2018-April/001454.html

1 files changed, 18 insertions, 8 deletions

diff --git a/src/dot.c b/src/dot.c
index 67f97b3..7c61018 100644
--- a/src/dot.c
+++ b/src/dot.c
@@ -41,11 +41,19 @@ mpc_dot (mpc_ptr res, const mpc_ptr *x, const mpc_ptr *y,
      and in z[n+i] that of -Im(x[i])*Im(y[i]) */
   for (i = 0; i < n; i++)
     {
-      mpfr_init2 (z[i], mpfr_get_prec (mpc_realref (x[i]))
-                  + mpfr_get_prec (mpc_realref (y[i])));
+      mpfr_prec_t prec_x_re = mpfr_get_prec (mpc_realref (x[i]));
+      mpfr_prec_t prec_x_im = mpfr_get_prec (mpc_imagref (x[i]));
+      mpfr_prec_t prec_y_re = mpfr_get_prec (mpc_realref (y[i]));
+      mpfr_prec_t prec_y_im = mpfr_get_prec (mpc_imagref (y[i]));
+      mpfr_prec_t prec_y_max = MPC_MAX (prec_y_re, prec_y_im);
+      /* we allocate z[i] with prec_x_re + prec_y_max bits
+         so that the second loop below does not reallocate */
+      mpfr_init2 (z[i], prec_x_re + prec_y_max);
+      mpfr_set_prec (z[i], prec_x_re + prec_y_re);
       mpfr_mul (z[i], mpc_realref (x[i]), mpc_realref (y[i]), MPFR_RNDZ);
-      mpfr_init2 (z[n+i], mpfr_get_prec (mpc_imagref (x[i]))
-                  + mpfr_get_prec (mpc_imagref (y[i])));
+      /* idem for z[n+i]: we allocate with prec_x_im + prec_y_max bits */
+      mpfr_init2 (z[n+i], prec_x_im + prec_y_max);
+      mpfr_set_prec (z[i], prec_x_im + prec_y_im);
       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);
     }
@@ -54,11 +62,13 @@ mpc_dot (mpc_ptr res, const mpc_ptr *x, const mpc_ptr *y,
      and in z[n+i] that of Im(x[i])*Re(y[i]) */
   for (i = 0; i < n; i++)
     {
-      mpfr_set_prec (z[i], mpfr_get_prec (mpc_realref (x[i]))
-                     + mpfr_get_prec (mpc_imagref (y[i])));
+      mpfr_prec_t prec_x_re = mpfr_get_prec (mpc_realref (x[i]));
+      mpfr_prec_t prec_x_im = mpfr_get_prec (mpc_imagref (x[i]));
+      mpfr_prec_t prec_y_re = mpfr_get_prec (mpc_realref (y[i]));
+      mpfr_prec_t prec_y_im = mpfr_get_prec (mpc_imagref (y[i]));
+      mpfr_set_prec (z[i], prec_x_re + prec_y_im);
       mpfr_mul (z[i], mpc_realref (x[i]), mpc_imagref (y[i]), MPFR_RNDZ);
-      mpfr_set_prec (z[n+i], mpfr_get_prec (mpc_imagref (x[i]))
-                  + mpfr_get_prec (mpc_realref (y[i])));
+      mpfr_set_prec (z[n+i], prec_x_im + prec_y_re);
       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));