1 /* Sum -- efficiently sum a list of floating-point numbers
3 Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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 3 of the License, or (at your
11 option) any later version.
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.
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.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 /* Reference: James Demmel and Yozo Hida, Fast and accurate floating-point
24 summation with application to computational geometry, Numerical Algorithms,
25 volume 37, number 1-4, pages 101--112, 2004. */
27 #define MPFR_NEED_LONGLONG_H
28 #include "mpfr-impl.h"
30 /* I would really like to use "mpfr_srcptr const []" but the norm is buggy:
31 it doesn't automaticaly cast a "mpfr_ptr []" to "mpfr_srcptr const []"
32 if necessary. So the choice are:
36 mpfr_s *const*const : ok
37 const mpfr_s *const* : no
38 const mpfr_s **const : no
39 const mpfr_s *const*const: no
40 VL: this is not a bug, but a feature. See the reason here:
41 http://c-faq.com/ansi/constmismatch.html
43 static void heap_sort (mpfr_srcptr
*const, unsigned long, mpfr_srcptr
*);
44 static void count_sort (mpfr_srcptr
*const, unsigned long, mpfr_srcptr
*,
45 mpfr_exp_t
, mpfr_uexp_t
);
47 /* Either sort the tab in perm and returns 0
48 Or returns 1 for +INF, -1 for -INF and 2 for NAN */
50 mpfr_sum_sort (mpfr_srcptr
*const tab
, unsigned long n
, mpfr_srcptr
*perm
)
60 for (i
= 0; i
< n
; i
++)
62 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tab
[i
])))
64 if (MPFR_IS_NAN (tab
[i
]))
65 return 2; /* Return NAN code */
66 else if (MPFR_IS_INF (tab
[i
]))
68 if (sign_inf
== 0) /* No previous INF */
69 sign_inf
= MPFR_SIGN (tab
[i
]);
70 else if (sign_inf
!= MPFR_SIGN (tab
[i
]))
71 return 2; /* Return NAN */
76 MPFR_ASSERTD (MPFR_IS_PURE_FP (tab
[i
]));
77 if (MPFR_GET_EXP (tab
[i
]) < min
)
78 min
= MPFR_GET_EXP(tab
[i
]);
79 if (MPFR_GET_EXP (tab
[i
]) > max
)
80 max
= MPFR_GET_EXP(tab
[i
]);
83 if (MPFR_UNLIKELY (sign_inf
!= 0))
86 exp_num
= max
- min
+ 1;
87 /* FIXME : better test */
88 if (exp_num
> n
* MPFR_INT_CEIL_LOG2 (n
))
89 heap_sort (tab
, n
, perm
);
91 count_sort (tab
, n
, perm
, min
, exp_num
);
95 #define GET_EXP1(x) (MPFR_IS_ZERO (x) ? min : MPFR_GET_EXP (x))
96 /* Performs a count sort of the entries */
98 count_sort (mpfr_srcptr
*const tab
, unsigned long n
,
99 mpfr_srcptr
*perm
, mpfr_exp_t min
, mpfr_uexp_t exp_num
)
101 unsigned long *account
;
102 unsigned long target_rank
, i
;
103 MPFR_TMP_DECL(marker
);
105 /* Reserve a place for potential 0 (with EXP min-1)
106 If there is no zero, we only lose one unused entry */
110 /* Performs a counting sort of the entries */
111 MPFR_TMP_MARK (marker
);
112 account
= (unsigned long *) MPFR_TMP_ALLOC (exp_num
* sizeof *account
);
113 for (i
= 0; i
< exp_num
; i
++)
115 for (i
= 0; i
< n
; i
++)
116 account
[GET_EXP1 (tab
[i
]) - min
]++;
117 for (i
= exp_num
- 1; i
>= 1; i
--)
118 account
[i
- 1] += account
[i
];
119 for (i
= 0; i
< n
; i
++)
121 target_rank
= --account
[GET_EXP1 (tab
[i
]) - min
];
122 perm
[target_rank
] = tab
[i
];
124 MPFR_TMP_FREE (marker
);
128 #define GET_EXP2(x) (MPFR_IS_ZERO (x) ? MPFR_EMIN_MIN : MPFR_GET_EXP (x))
130 /* Performs a heap sort of the entries */
132 heap_sort (mpfr_srcptr
*const tab
, unsigned long n
, mpfr_srcptr
*perm
)
134 unsigned long dernier_traite
;
135 unsigned long i
, pere
;
137 unsigned long fils_gauche
, fils_droit
, fils_indigne
;
138 /* Reminder of a heap structure :
139 node(i) has for left son node(2i +1) and right son node(2i)
140 and father(node(i)) = node((i - 1) / 2)
143 /* initialize the permutation to identity */
144 for (i
= 0; i
< n
; i
++)
147 /* insertion phase */
148 for (dernier_traite
= 1; dernier_traite
< n
; dernier_traite
++)
154 if (GET_EXP2 (perm
[pere
]) > GET_EXP2 (perm
[i
]))
157 perm
[pere
] = perm
[i
];
166 /* extraction phase */
167 for (dernier_traite
= n
- 1; dernier_traite
> 0; dernier_traite
--)
170 perm
[0] = perm
[dernier_traite
];
171 perm
[dernier_traite
] = tmp
;
176 fils_gauche
= 2 * i
+ 1;
177 fils_droit
= fils_gauche
+ 1;
178 if (fils_gauche
< dernier_traite
)
180 if (fils_droit
< dernier_traite
)
182 if (GET_EXP2(perm
[fils_droit
]) < GET_EXP2(perm
[fils_gauche
]))
183 fils_indigne
= fils_droit
;
185 fils_indigne
= fils_gauche
;
187 if (GET_EXP2 (perm
[i
]) > GET_EXP2 (perm
[fils_indigne
]))
190 perm
[i
] = perm
[fils_indigne
];
191 perm
[fils_indigne
] = tmp
;
197 else /* on a un fils gauche, pas de fils droit */
199 if (GET_EXP2 (perm
[i
]) > GET_EXP2 (perm
[fils_gauche
]))
202 perm
[i
] = perm
[fils_gauche
];
203 perm
[fils_gauche
] = tmp
;
208 else /* on n'a pas de fils */
215 /* Sum a list of float with order given by permutation perm,
216 * intermediate size set to F.
217 * Internal use function.
220 sum_once (mpfr_ptr ret
, mpfr_srcptr
*const tab
, unsigned long n
, mpfr_prec_t F
)
226 MPFR_ASSERTD (n
>= 2);
229 error_trap
= mpfr_set (sum
, tab
[0], MPFR_RNDN
);
230 for (i
= 1; i
< n
- 1; i
++)
232 MPFR_ASSERTD (!MPFR_IS_NAN (sum
) && !MPFR_IS_INF (sum
));
233 error_trap
|= mpfr_add (sum
, sum
, tab
[i
], MPFR_RNDN
);
235 error_trap
|= mpfr_add (ret
, sum
, tab
[n
- 1], MPFR_RNDN
);
240 /* Sum a list of floating-point numbers.
244 mpfr_sum (mpfr_ptr ret
, mpfr_ptr
*const tab_p
, unsigned long n
, mpfr_rnd_t rnd
)
248 mpfr_srcptr
*perm
, *const tab
= (mpfr_srcptr
*) tab_p
;
250 MPFR_ZIV_DECL (loop
);
251 MPFR_SAVE_EXPO_DECL (expo
);
252 MPFR_TMP_DECL (marker
);
254 if (MPFR_UNLIKELY (n
<= 1))
263 return mpfr_set (ret
, tab
[0], rnd
);
266 /* Sort and treat special cases */
267 MPFR_TMP_MARK (marker
);
268 perm
= (mpfr_srcptr
*) MPFR_TMP_ALLOC (n
* sizeof *perm
);
269 error_trap
= mpfr_sum_sort (tab
, n
, perm
);
270 /* Check if there was a NAN or a INF */
271 if (MPFR_UNLIKELY (error_trap
!= 0))
273 MPFR_TMP_FREE (marker
);
280 MPFR_SET_SIGN (ret
, error_trap
);
284 /* Initial precision */
285 prec
= MAX (MPFR_PREC (tab
[0]), MPFR_PREC (ret
));
286 k
= MPFR_INT_CEIL_LOG2 (n
) + 1;
288 mpfr_init2 (cur_sum
, prec
);
291 MPFR_SAVE_EXPO_MARK (expo
);
292 MPFR_ZIV_INIT (loop
, prec
);
295 error_trap
= sum_once (cur_sum
, perm
, n
, prec
+ k
);
296 if (MPFR_LIKELY (error_trap
== 0 ||
297 (!MPFR_IS_ZERO (cur_sum
) &&
298 mpfr_can_round (cur_sum
,
299 MPFR_GET_EXP (cur_sum
) - prec
+ 2,
300 MPFR_RNDN
, rnd
, MPFR_PREC (ret
)))))
302 MPFR_ZIV_NEXT (loop
, prec
);
303 mpfr_set_prec (cur_sum
, prec
);
305 MPFR_ZIV_FREE (loop
);
306 MPFR_TMP_FREE (marker
);
308 error_trap
|= mpfr_set (ret
, cur_sum
, rnd
);
309 mpfr_clear (cur_sum
);
311 MPFR_SAVE_EXPO_FREE (expo
);
312 error_trap
|= mpfr_check_range (ret
, 0, rnd
);
313 return error_trap
; /* It doesn't return the ternary value */