contrib/mpfr/sum.c

   1 /* Sum -- efficiently sum a list of floating-point numbers
   2
   3 Copyright 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 2.1 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* I would really like to use "mpfr_srcptr const []" but the norm is buggy:
  27    it doesn't automaticaly cast a "mpfr_ptr []" to "mpfr_srcptr const []"
  28    if necessary. So the choice are:
  29      mpfr_s **                : ok
  30      mpfr_s *const*           : ok
  31      mpfr_s **const           : ok
  32      mpfr_s *const*const      : ok
  33      const mpfr_s *const*     : no
  34      const mpfr_s **const     : no
  35      const mpfr_s *const*const: no
  36    VL: this is not a bug, but a feature. See the reason here:
  37      http://c-faq.com/ansi/constmismatch.html
  38 */
  39 static void heap_sort (mpfr_srcptr *const, unsigned long, mpfr_srcptr *);
  40 static void count_sort (mpfr_srcptr *const, unsigned long, mpfr_srcptr *,
  41                         mp_exp_t, mpfr_uexp_t);
  42
  43 /* Either sort the tab in perm and returns 0
  44    Or returns 1 for +INF, -1 for -INF and 2 for NAN */
  45 int
  46 mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm)
  47 {
  48   mp_exp_t min, max;
  49   mpfr_uexp_t exp_num;
  50   unsigned long i;
  51   int sign_inf;
  52
  53   sign_inf = 0;
  54   min = MPFR_EMIN_MAX;
  55   max = MPFR_EMAX_MIN;
  56   for (i = 0; i < n; i++)
  57     {
  58       if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tab[i])))
  59         {
  60           if (MPFR_IS_NAN (tab[i]))
  61             return 2; /* Return NAN code */
  62           else if (MPFR_IS_INF (tab[i]))
  63             {
  64               if (sign_inf == 0) /* No previous INF */
  65                 sign_inf = MPFR_SIGN (tab[i]);
  66               else if (sign_inf != MPFR_SIGN (tab[i]))
  67                 return 2; /* Return NAN */
  68             }
  69         }
  70       else
  71         {
  72           if (MPFR_GET_EXP (tab[i]) < min)
  73             min = MPFR_GET_EXP(tab[i]);
  74           if (MPFR_GET_EXP (tab[i]) > max)
  75             max = MPFR_GET_EXP(tab[i]);
  76         }
  77     }
  78   if (MPFR_UNLIKELY (sign_inf != 0))
  79     return sign_inf;
  80
  81   exp_num = max - min + 1;
  82   /* FIXME : better test */
  83   if (exp_num > n * MPFR_INT_CEIL_LOG2 (n))
  84     heap_sort (tab, n, perm);
  85   else
  86     count_sort (tab, n, perm, min, exp_num);
  87   return 0;
  88 }
  89
  90 #define GET_EXP1(x) (MPFR_IS_ZERO (x) ? min : MPFR_GET_EXP (x))
  91 /* Performs a count sort of the entries */
  92 static void
  93 count_sort (mpfr_srcptr *const tab, unsigned long n,
  94             mpfr_srcptr *perm, mp_exp_t min, mpfr_uexp_t exp_num)
  95 {
  96   unsigned long *account;
  97   unsigned long target_rank, i;
  98   MPFR_TMP_DECL(marker);
  99
 100   /* Reserve a place for potential 0 (with EXP min-1)
 101      If there is no zero, we only lose one unused entry */
 102   min--;
 103   exp_num++;
 104
 105   /* Performs a counting sort of the entries */
 106   MPFR_TMP_MARK (marker);
 107   account = (unsigned long *) MPFR_TMP_ALLOC (exp_num * sizeof * account);
 108   for (i = 0; i < exp_num; i++)
 109     account[i] = 0;
 110   for (i = 0; i < n; i++)
 111     account[GET_EXP1 (tab[i]) - min]++;
 112   for (i = exp_num - 1; i >= 1; i--)
 113     account[i - 1] += account[i];
 114   for (i = 0; i < n; i++)
 115     {
 116       target_rank = --account[GET_EXP1 (tab[i]) - min];
 117       perm[target_rank] = tab[i];
 118     }
 119   MPFR_TMP_FREE (marker);
 120 }
 121
 122
 123 #define GET_EXP2(x) (MPFR_IS_ZERO (x) ? MPFR_EMIN_MIN : MPFR_GET_EXP (x))
 124
 125 /* Performs a heap sort of the entries */
 126 static void
 127 heap_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm)
 128 {
 129   unsigned long dernier_traite;
 130   unsigned long i, pere;
 131   mpfr_srcptr tmp;
 132   unsigned long fils_gauche, fils_droit, fils_indigne;
 133   /* Reminder of a heap structure :
 134      node(i) has for left son node(2i +1) and right son node(2i)
 135      and father(node(i)) = node((i - 1) / 2)
 136   */
 137
 138   /* initialize the permutation to identity */
 139   for (i = 0; i < n; i++)
 140     perm[i] = tab[i];
 141
 142   /* insertion phase */
 143   for (dernier_traite = 1; dernier_traite < n; dernier_traite++)
 144     {
 145       i = dernier_traite;
 146       while (i > 0)
 147         {
 148           pere = (i - 1) / 2;
 149           if (GET_EXP2 (perm[pere]) > GET_EXP2 (perm[i]))
 150             {
 151               tmp = perm[pere];
 152               perm[pere] = perm[i];
 153               perm[i] = tmp;
 154               i = pere;
 155             }
 156           else
 157             break;
 158         }
 159     }
 160
 161   /* extraction phase */
 162   for (dernier_traite = n - 1; dernier_traite > 0; dernier_traite--)
 163     {
 164       tmp = perm[0];
 165       perm[0] = perm[dernier_traite];
 166       perm[dernier_traite] = tmp;
 167
 168       i = 0;
 169       while (1)
 170         {
 171           fils_gauche = 2 * i + 1;
 172           fils_droit = fils_gauche + 1;
 173           if (fils_gauche < dernier_traite)
 174             {
 175               if (fils_droit < dernier_traite)
 176                 {
 177                   if (GET_EXP2(perm[fils_droit]) < GET_EXP2(perm[fils_gauche]))
 178                     fils_indigne = fils_droit;
 179                   else
 180                     fils_indigne = fils_gauche;
 181
 182                   if (GET_EXP2 (perm[i]) > GET_EXP2 (perm[fils_indigne]))
 183                     {
 184                       tmp = perm[i];
 185                       perm[i] = perm[fils_indigne];
 186                       perm[fils_indigne] = tmp;
 187                       i = fils_indigne;
 188                     }
 189                   else
 190                     break;
 191                 }
 192               else /* on a un fils gauche, pas de fils droit */
 193                 {
 194                   if (GET_EXP2 (perm[i]) > GET_EXP2 (perm[fils_gauche]))
 195                     {
 196                       tmp = perm[i];
 197                       perm[i] = perm[fils_gauche];
 198                       perm[fils_gauche] = tmp;
 199                     }
 200                   break;
 201                 }
 202             }
 203           else /* on n'a pas de fils */
 204             break;
 205         }
 206     }
 207 }
 208
 209
 210 /* Sum a list of float with order given by permutation perm,
 211  * intermediate size set to F.
 212  * Internal use function.
 213  */
 214 static int
 215 sum_once (mpfr_ptr ret, mpfr_srcptr *const tab, unsigned long n, mp_prec_t F)
 216 {
 217   mpfr_t sum;
 218   unsigned long i;
 219   int error_trap;
 220
 221   MPFR_ASSERTD (n >= 2);
 222
 223   mpfr_init2 (sum, F);
 224   error_trap = mpfr_set (sum, tab[0], GMP_RNDN);
 225   for (i = 1; i < n - 1; i++)
 226     {
 227       MPFR_ASSERTD (!MPFR_IS_NAN (sum) && !MPFR_IS_INF (sum));
 228       error_trap |= mpfr_add (sum, sum, tab[i], GMP_RNDN);
 229     }
 230   error_trap |= mpfr_add (ret, sum, tab[n - 1], GMP_RNDN);
 231   mpfr_clear (sum);
 232   return error_trap;
 233 }
 234
 235 /* Sum a list of floating-point numbers.
 236  * FIXME : add reference to Demmel-Hida's paper.
 237  */
 238
 239 int
 240 mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mp_rnd_t rnd)
 241 {
 242   mpfr_t cur_sum;
 243   mp_prec_t prec;
 244   mpfr_srcptr *perm, *const tab = (mpfr_srcptr *) tab_p;
 245   int k, error_trap;
 246   MPFR_ZIV_DECL (loop);
 247   MPFR_SAVE_EXPO_DECL (expo);
 248   MPFR_TMP_DECL (marker);
 249
 250   if (MPFR_UNLIKELY (n <= 1))
 251     {
 252       if (n < 1)
 253         {
 254           MPFR_SET_ZERO (ret);
 255           MPFR_SET_POS (ret);
 256           return 0;
 257         }
 258       else
 259         return mpfr_set (ret, tab[0], rnd);
 260     }
 261
 262   /* Sort and treat special cases */
 263   MPFR_TMP_MARK (marker);
 264   perm = (mpfr_srcptr *) MPFR_TMP_ALLOC (n * sizeof *perm);
 265   error_trap = mpfr_sum_sort (tab, n, perm);
 266   /* Check if there was a NAN or a INF */
 267   if (MPFR_UNLIKELY (error_trap != 0))
 268     {
 269       MPFR_TMP_FREE (marker);
 270       if (error_trap == 2)
 271         {
 272           MPFR_SET_NAN (ret);
 273           MPFR_RET_NAN;
 274         }
 275       MPFR_SET_INF (ret);
 276       MPFR_SET_SIGN (ret, error_trap);
 277       MPFR_RET (0);
 278     }
 279
 280   /* Initial precision */
 281   prec = MAX (MPFR_PREC (tab[0]), MPFR_PREC (ret));
 282   k = MPFR_INT_CEIL_LOG2 (n) + 1;
 283   prec += k + 2;
 284   mpfr_init2 (cur_sum, prec);
 285
 286   /* Ziv Loop */
 287   MPFR_SAVE_EXPO_MARK (expo);
 288   MPFR_ZIV_INIT (loop, prec);
 289   for (;;)
 290     {
 291       error_trap = sum_once (cur_sum, perm, n, prec + k);
 292       if (MPFR_LIKELY (error_trap == 0 ||
 293                        (!MPFR_IS_ZERO (cur_sum) &&
 294                         mpfr_can_round (cur_sum,
 295                                         MPFR_GET_EXP (cur_sum) - prec + 2,
 296                                         GMP_RNDN, rnd, MPFR_PREC (ret)))))
 297         break;
 298       MPFR_ZIV_NEXT (loop, prec);
 299       mpfr_set_prec (cur_sum, prec);
 300     }
 301   MPFR_ZIV_FREE (loop);
 302   MPFR_TMP_FREE (marker);
 303
 304   error_trap |= mpfr_set (ret, cur_sum, rnd);
 305   mpfr_clear (cur_sum);
 306
 307   MPFR_SAVE_EXPO_FREE (expo);
 308   error_trap |= mpfr_check_range (ret, 0, rnd);
 309   return error_trap; /* It doesn't return the ternary value */
 310 }
 311
 312 /* __END__ */