[sum.c] mpfr_sum should consider the maximal precision of its inputs/output

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@8699 280ebfd0-de03-0410-8827-d642c229c3f4

3 files changed, 15 insertions, 7 deletions

diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h
index e71dd877f..63b12e278 100644
--- a/src/mpfr-impl.h
+++ b/src/mpfr-impl.h
@@ -1958,7 +1958,8 @@ __MPFR_DECLSPEC int mpfr_round_raw_4 _MPFR_PROTO ((mp_limb_t *,
 __MPFR_DECLSPEC int mpfr_check _MPFR_PROTO ((mpfr_srcptr));
 
 __MPFR_DECLSPEC int mpfr_sum_sort _MPFR_PROTO ((mpfr_srcptr *const,
-                                                unsigned long, mpfr_srcptr *));
+                                                unsigned long, mpfr_srcptr *,
+                                                mpfr_prec_t *));
 
 __MPFR_DECLSPEC int mpfr_get_cputime _MPFR_PROTO ((void));
 
diff --git a/src/sum.c b/src/sum.c
index 86d21243d..6d345a385 100644
--- a/src/sum.c
+++ b/src/sum.c
@@ -45,9 +45,13 @@ static void count_sort (mpfr_srcptr *const, unsigned long, mpfr_srcptr *,
                         mpfr_exp_t, mpfr_uexp_t);
 
 /* Either sort the tab in perm and returns 0
-   Or returns 1 for +INF, -1 for -INF and 2 for NAN */
+   Or returns 1 for +INF, -1 for -INF and 2 for NAN.
+   Also set *maxprec to the maximal precision of tab[0..n-1] and of the
+   initial value of *maxprec.
+*/
 int
-mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm)
+mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm,
+               mpfr_prec_t *maxprec)
 {
   mpfr_exp_t min, max;
   mpfr_uexp_t exp_num;
@@ -79,6 +83,8 @@ mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm)
           if (MPFR_GET_EXP (tab[i]) > max)
             max = MPFR_GET_EXP(tab[i]);
         }
+      if (MPFR_PREC (tab[i]) > *maxprec)
+        *maxprec = MPFR_PREC (tab[i]);
     }
   if (MPFR_UNLIKELY (sign_inf != 0))
     return sign_inf;
@@ -266,7 +272,8 @@ mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mpfr_rnd_t rnd)
   /* Sort and treat special cases */
   MPFR_TMP_MARK (marker);
   perm = (mpfr_srcptr *) MPFR_TMP_ALLOC (n * sizeof *perm);
-  error_trap = mpfr_sum_sort (tab, n, perm);
+  prec = MPFR_PREC (ret);
+  error_trap = mpfr_sum_sort (tab, n, perm, &prec);
   /* Check if there was a NAN or a INF */
   if (MPFR_UNLIKELY (error_trap != 0))
     {
@@ -281,8 +288,7 @@ mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mpfr_rnd_t rnd)
       MPFR_RET (0);
     }
 
-  /* Initial precision */
-  prec = MAX (MPFR_PREC (tab[0]), MPFR_PREC (ret));
+  /* Initial precision is max(prec(ret),prec(tab[0]),...,prec(tab[n-1])) */
   k = MPFR_INT_CEIL_LOG2 (n) + 1;
   prec += k + 2;
   mpfr_init2 (cur_sum, prec);
diff --git a/tests/tsum.c b/tests/tsum.c
index ae235540e..b03eb1fb7 100644
--- a/tests/tsum.c
+++ b/tests/tsum.c
@@ -136,6 +136,7 @@ test_sort (mpfr_prec_t f, unsigned long n)
   mpfr_ptr *tabtmp;
   mpfr_srcptr *perm;
   unsigned long i;
+  mpfr_prec_t prec = MPFR_PREC_MIN;
 
   /* Init stuff */
   tab = (mpfr_t *) (*__gmp_allocate_func) (n * sizeof (mpfr_t));
@@ -150,7 +151,7 @@ test_sort (mpfr_prec_t f, unsigned long n)
       tabtmp[i] = tab[i];
     }
 
-  mpfr_sum_sort ((mpfr_srcptr *)tabtmp, n, perm);
+  mpfr_sum_sort ((mpfr_srcptr *)tabtmp, n, perm, &prec);
 
   if (check_is_sorted (n, perm) == 0)
     {