Added mpfr_sum function.

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

1 files changed, 231 insertions, 0 deletions

diff --git a/sum.c b/sum.c
new file mode 100644
index 000000000..9a5622167
--- /dev/null
+++ b/sum.c
@@ -0,0 +1,231 @@
+/* Sum -- efficiently sum a list of floating-point numbers
+
+Copyright 2004 Free Software Foundation, Inc.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the MPFR Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include "mpfr-impl.h"
+
+/* Performs a counting sort of the entries */
+
+static void heap_sort_exp_clean (mpfr_ptr const tab[], int n, 
+                                 mpfr_srcptr *perm);
+
+void mpfr_count_sort (mpfr_ptr const tab[], unsigned int n, 
+                 mpfr_srcptr *perm)
+{
+    mp_exp_t min, max;
+    unsigned int i;
+    unsigned int *account;
+    unsigned int exp_num;
+    unsigned int target_rank;
+    TMP_DECL(marker);
+    
+    TMP_MARK(marker);
+    min = max = MPFR_GET_EXP(tab[0]);
+
+    for (i = 1; i < n; i++)
+    {
+        if (MPFR_GET_EXP(tab[i]) < min)
+            min = MPFR_GET_EXP(tab[i]);
+        if (MPFR_GET_EXP(tab[i]) > max)
+            max = MPFR_GET_EXP(tab[i]);
+    }
+
+    exp_num = max - min + 1;
+    if (exp_num > 42 * __gmpfr_ceil_log2 (n)) /* FIXME : better test */
+    {
+        heap_sort_exp_clean (tab, n, perm);
+        return;
+    }
+    account = (unsigned int *) TMP_ALLOC(exp_num * sizeof(*account));
+    for (i = 0; i < exp_num; i++)
+        account[i] = 0;
+    for (i = 0; i < n; i++)
+        account[MPFR_GET_EXP(tab[i]) - min]++;
+    for (i = exp_num - 1; i >= 1; i--)
+        account[i - 1] += account[i];
+    for (i = 0; i < n; i++)
+    {
+        target_rank = --account[MPFR_GET_EXP(tab[i]) - min];
+        perm[target_rank] = tab[i];
+    }
+    
+    TMP_FREE(account);
+}
+
+/* Performs a heap sort of the entries */
+
+static void heap_sort_exp_clean (mpfr_ptr const tab[], int n, 
+                                 mpfr_srcptr *perm)
+{
+  unsigned int dernier_traite;
+  unsigned int i, pere;
+  mpfr_srcptr tmp;
+  unsigned int fils_gauche, fils_droit, fils_indigne;
+  /* Reminder of a heap structure :
+     node(i) has for left son node(2i +1) and right son node(2i)
+     and father(node(i)) = node((i - 1) / 2)
+  */
+  
+  /* initialize the permutation to identity */
+
+  for (i = 0; i < n; i++)
+    perm[i] = tab[i];
+
+  /* insertion phase */
+
+  for (dernier_traite = 1; dernier_traite < n; dernier_traite++)
+    {
+      i = dernier_traite;
+      while (i > 0)
+	{
+	  pere = (i - 1) / 2;
+	  if (MPFR_GET_EXP(perm[pere]) > MPFR_GET_EXP(perm[i]))
+	    {
+	      tmp = perm[pere];
+	      perm[pere] = perm[i];
+	      perm[i] = tmp;
+	      i = pere;
+	    }
+	  else
+	    break;
+	}
+    }
+
+  /* extraction phase */
+  
+  for (dernier_traite = n - 1; dernier_traite > 0; dernier_traite--)
+    {
+      tmp = perm[0];
+      perm[0] = perm[dernier_traite];
+      perm[dernier_traite] = tmp;
+
+      i = 0;
+      while (1)
+	{
+	  fils_gauche = 2 * i + 1;
+	  fils_droit = fils_gauche + 1;
+	  if (fils_gauche < dernier_traite)
+	    {
+	      if (fils_droit < dernier_traite)
+		{
+		  if (MPFR_GET_EXP(perm[fils_droit]) < MPFR_GET_EXP(perm[fils_gauche]))
+		    fils_indigne = fils_droit;
+		  else
+		    fils_indigne = fils_gauche;
+
+		  if (MPFR_GET_EXP(perm[i]) > MPFR_GET_EXP(perm[fils_indigne]))
+		    {
+		      tmp = perm[i];
+		      perm[i] = perm[fils_indigne];
+		      perm[fils_indigne] = tmp;
+		      i = fils_indigne;
+		    }
+		  else
+		    break;
+		}
+	      else /* on a un fils gauche, pas de fils droit */
+		{
+		  if (MPFR_GET_EXP(perm[i]) > MPFR_GET_EXP(perm[fils_gauche]))
+		    {
+		      tmp = perm[i];
+		      perm[i] = perm[fils_gauche];
+		      perm[fils_gauche] = tmp;
+		    }
+		  break;
+		}
+	    }
+	  else /* on n'a pas de fils */
+	    break;
+	}
+    }
+}
+
+
+/* Sum a list of float with order given by permutation perm,
+ * intermediate size set to F.
+ * Internal use function.
+ */
+
+static int mpfr_list_sum_once (mpfr_ptr ret, mpfr_srcptr const tab[], 
+                               unsigned int n, mp_prec_t F)
+{
+  unsigned int i;
+  mpfr_t sum;
+  int error_trap = 0;
+
+  mpfr_init2 (sum, F);
+  mpfr_set (sum, tab[0], GMP_RNDN);
+
+  for (i = 1; i < n - 1; i++)
+    {
+      error_trap |= mpfr_add (sum, sum, tab[i], GMP_RNDN);
+    }
+
+  error_trap |= mpfr_add (ret, sum, tab[n - 1], GMP_RNDN);
+  mpfr_clear (sum);
+  return error_trap;
+}
+
+/* Sum a list of floating-point numbers.
+ * FIXME : add reference to Demmel-Hida's paper.
+*/
+
+int mpfr_sum (mpfr_ptr ret, mpfr_ptr const tab[], unsigned int n, 
+              mp_rnd_t rnd)
+{
+  mp_prec_t initial_f, current_f;
+  unsigned int k;
+  mpfr_srcptr *perm;
+  unsigned int guard_digits;
+  unsigned int initial_guard_digits;
+  int error_trap;
+  mpfr_t cur_sum;
+  TMP_DECL(marker);
+    
+  TMP_MARK(marker);
+
+  perm = (mpfr_srcptr *) TMP_ALLOC(n * sizeof(mpfr_srcptr)); 
+
+  mpfr_count_sort (tab, n, perm);
+
+  initial_f = MPFR_PREC(tab[0]);
+  k = __gmpfr_ceil_log2 ((double) n) + 1;
+  mpfr_init2 (cur_sum, initial_f);
+  initial_guard_digits = k + 2;
+  guard_digits = initial_guard_digits;
+  do
+  {
+      current_f = initial_f + guard_digits;
+      mpfr_set_prec (cur_sum, current_f);
+      error_trap = mpfr_list_sum_once (cur_sum, perm, n, 
+                                       current_f + k);
+      guard_digits *= 2;
+  }
+  while ((error_trap != 0) &&
+          !(mpfr_can_round (cur_sum, MPFR_GET_EXP(cur_sum) - current_f + 2,
+                            GMP_RNDN, rnd, initial_f)));
+  error_trap |= mpfr_set (ret, cur_sum, rnd);
+  mpfr_clear (cur_sum);
+  TMP_FREE(marker);
+  return error_trap;
+}
+
+
+/* __END__ */