| /* Sum -- efficiently sum a list of floating-point numbers |
| |
| Copyright 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc. |
| Contributed by the Arenaire and Cacao projects, INRIA. |
| |
| This file is part of the GNU MPFR Library. |
| |
| The GNU 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 GNU 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 GNU MPFR Library; see the file COPYING.LIB. If not, write to |
| the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, |
| MA 02110-1301, USA. */ |
| |
| #define MPFR_NEED_LONGLONG_H |
| #include "mpfr-impl.h" |
| |
| /* I would really like to use "mpfr_srcptr const []" but the norm is buggy: |
| it doesn't automaticaly cast a "mpfr_ptr []" to "mpfr_srcptr const []" |
| if necessary. So the choice are: |
| mpfr_s ** : ok |
| mpfr_s *const* : ok |
| mpfr_s **const : ok |
| mpfr_s *const*const : ok |
| const mpfr_s *const* : no |
| const mpfr_s **const : no |
| const mpfr_s *const*const: no |
| VL: this is not a bug, but a feature. See the reason here: |
| http://c-faq.com/ansi/constmismatch.html |
| */ |
| static void heap_sort (mpfr_srcptr *const, unsigned long, mpfr_srcptr *); |
| static void count_sort (mpfr_srcptr *const, unsigned long, mpfr_srcptr *, |
| mp_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 */ |
| int |
| mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm) |
| { |
| mp_exp_t min, max; |
| mpfr_uexp_t exp_num; |
| unsigned long i; |
| int sign_inf; |
| |
| sign_inf = 0; |
| min = MPFR_EMIN_MAX; |
| max = MPFR_EMAX_MIN; |
| for (i = 0; i < n; i++) |
| { |
| if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tab[i]))) |
| { |
| if (MPFR_IS_NAN (tab[i])) |
| return 2; /* Return NAN code */ |
| else if (MPFR_IS_INF (tab[i])) |
| { |
| if (sign_inf == 0) /* No previous INF */ |
| sign_inf = MPFR_SIGN (tab[i]); |
| else if (sign_inf != MPFR_SIGN (tab[i])) |
| return 2; /* Return NAN */ |
| } |
| } |
| else |
| { |
| 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]); |
| } |
| } |
| if (MPFR_UNLIKELY (sign_inf != 0)) |
| return sign_inf; |
| |
| exp_num = max - min + 1; |
| /* FIXME : better test */ |
| if (exp_num > n * MPFR_INT_CEIL_LOG2 (n)) |
| heap_sort (tab, n, perm); |
| else |
| count_sort (tab, n, perm, min, exp_num); |
| return 0; |
| } |
| |
| #define GET_EXP1(x) (MPFR_IS_ZERO (x) ? min : MPFR_GET_EXP (x)) |
| /* Performs a count sort of the entries */ |
| static void |
| count_sort (mpfr_srcptr *const tab, unsigned long n, |
| mpfr_srcptr *perm, mp_exp_t min, mpfr_uexp_t exp_num) |
| { |
| unsigned long *account; |
| unsigned long target_rank, i; |
| MPFR_TMP_DECL(marker); |
| |
| /* Reserve a place for potential 0 (with EXP min-1) |
| If there is no zero, we only lose one unused entry */ |
| min--; |
| exp_num++; |
| |
| /* Performs a counting sort of the entries */ |
| MPFR_TMP_MARK (marker); |
| account = (unsigned long *) MPFR_TMP_ALLOC (exp_num * sizeof * account); |
| for (i = 0; i < exp_num; i++) |
| account[i] = 0; |
| for (i = 0; i < n; i++) |
| account[GET_EXP1 (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[GET_EXP1 (tab[i]) - min]; |
| perm[target_rank] = tab[i]; |
| } |
| MPFR_TMP_FREE (marker); |
| } |
| |
| |
| #define GET_EXP2(x) (MPFR_IS_ZERO (x) ? MPFR_EMIN_MIN : MPFR_GET_EXP (x)) |
| |
| /* Performs a heap sort of the entries */ |
| static void |
| heap_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm) |
| { |
| unsigned long dernier_traite; |
| unsigned long i, pere; |
| mpfr_srcptr tmp; |
| unsigned long 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 (GET_EXP2 (perm[pere]) > GET_EXP2 (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 (GET_EXP2(perm[fils_droit]) < GET_EXP2(perm[fils_gauche])) |
| fils_indigne = fils_droit; |
| else |
| fils_indigne = fils_gauche; |
| |
| if (GET_EXP2 (perm[i]) > GET_EXP2 (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 (GET_EXP2 (perm[i]) > GET_EXP2 (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 |
| sum_once (mpfr_ptr ret, mpfr_srcptr *const tab, unsigned long n, mp_prec_t F) |
| { |
| mpfr_t sum; |
| unsigned long i; |
| int error_trap; |
| |
| MPFR_ASSERTD (n >= 2); |
| |
| mpfr_init2 (sum, F); |
| error_trap = mpfr_set (sum, tab[0], GMP_RNDN); |
| for (i = 1; i < n - 1; i++) |
| { |
| MPFR_ASSERTD (!MPFR_IS_NAN (sum) && !MPFR_IS_INF (sum)); |
| 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_p, unsigned long n, mp_rnd_t rnd) |
| { |
| mpfr_t cur_sum; |
| mp_prec_t prec; |
| mpfr_srcptr *perm, *const tab = (mpfr_srcptr *) tab_p; |
| int k, error_trap; |
| MPFR_ZIV_DECL (loop); |
| MPFR_SAVE_EXPO_DECL (expo); |
| MPFR_TMP_DECL (marker); |
| |
| if (MPFR_UNLIKELY (n <= 1)) |
| { |
| if (n < 1) |
| { |
| MPFR_SET_ZERO (ret); |
| MPFR_SET_POS (ret); |
| return 0; |
| } |
| else |
| return mpfr_set (ret, tab[0], 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); |
| /* Check if there was a NAN or a INF */ |
| if (MPFR_UNLIKELY (error_trap != 0)) |
| { |
| MPFR_TMP_FREE (marker); |
| if (error_trap == 2) |
| { |
| MPFR_SET_NAN (ret); |
| MPFR_RET_NAN; |
| } |
| MPFR_SET_INF (ret); |
| MPFR_SET_SIGN (ret, error_trap); |
| MPFR_RET (0); |
| } |
| |
| /* Initial precision */ |
| prec = MAX (MPFR_PREC (tab[0]), MPFR_PREC (ret)); |
| k = MPFR_INT_CEIL_LOG2 (n) + 1; |
| prec += k + 2; |
| mpfr_init2 (cur_sum, prec); |
| |
| /* Ziv Loop */ |
| MPFR_SAVE_EXPO_MARK (expo); |
| MPFR_ZIV_INIT (loop, prec); |
| for (;;) |
| { |
| error_trap = sum_once (cur_sum, perm, n, prec + k); |
| if (MPFR_LIKELY (error_trap == 0 || |
| (!MPFR_IS_ZERO (cur_sum) && |
| mpfr_can_round (cur_sum, |
| MPFR_GET_EXP (cur_sum) - prec + 2, |
| GMP_RNDN, rnd, MPFR_PREC (ret))))) |
| break; |
| MPFR_ZIV_NEXT (loop, prec); |
| mpfr_set_prec (cur_sum, prec); |
| } |
| MPFR_ZIV_FREE (loop); |
| MPFR_TMP_FREE (marker); |
| |
| error_trap |= mpfr_set (ret, cur_sum, rnd); |
| mpfr_clear (cur_sum); |
| |
| MPFR_SAVE_EXPO_FREE (expo); |
| error_trap |= mpfr_check_range (ret, 0, rnd); |
| return error_trap; /* It doesn't return the ternary value */ |
| } |
| |
| /* __END__ */ |