Python code coverage for Modules/_decimal/libmpdec/basearith.c

#	count	content
1	n/a	/*
2	n/a	* Copyright (c) 2008-2016 Stefan Krah. All rights reserved.
3	n/a	*
4	n/a	* Redistribution and use in source and binary forms, with or without
5	n/a	* modification, are permitted provided that the following conditions
6	n/a	* are met:
7	n/a	*
8	n/a	* 1. Redistributions of source code must retain the above copyright
9	n/a	* notice, this list of conditions and the following disclaimer.
10	n/a	*
11	n/a	* 2. Redistributions in binary form must reproduce the above copyright
12	n/a	* notice, this list of conditions and the following disclaimer in the
13	n/a	* documentation and/or other materials provided with the distribution.
14	n/a	*
15	n/a	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
16	n/a	* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17	n/a	* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18	n/a	* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19	n/a	* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20	n/a	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21	n/a	* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22	n/a	* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23	n/a	* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24	n/a	* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25	n/a	* SUCH DAMAGE.
26	n/a	*/
27	n/a
28	n/a
29	n/a	#include "mpdecimal.h"
30	n/a	#include <stdio.h>
31	n/a	#include <stdlib.h>
32	n/a	#include <string.h>
33	n/a	#include <assert.h>
34	n/a	#include "constants.h"
35	n/a	#include "typearith.h"
36	n/a	#include "basearith.h"
37	n/a
38	n/a
39	n/a	/*********************************************************************/
40	n/a	/* Calculations in base MPD_RADIX */
41	n/a	/*********************************************************************/
42	n/a
43	n/a
44	n/a	/*
45	n/a	* Knuth, TAOCP, Volume 2, 4.3.1:
46	n/a	* w := sum of u (len m) and v (len n)
47	n/a	* n > 0 and m >= n
48	n/a	* The calling function has to handle a possible final carry.
49	n/a	*/
50	n/a	mpd_uint_t
51	n/a	_mpd_baseadd(mpd_uint_t w, const mpd_uint_t u, const mpd_uint_t *v,
52	n/a	mpd_size_t m, mpd_size_t n)
53	n/a	{
54	n/a	mpd_uint_t s;
55	n/a	mpd_uint_t carry = 0;
56	n/a	mpd_size_t i;
57	n/a
58	n/a	assert(n > 0 && m >= n);
59	n/a
60	n/a	/* add n members of u and v */
61	n/a	for (i = 0; i < n; i++) {
62	n/a	s = u[i] + (v[i] + carry);
63	n/a	carry = (s < u[i]) \| (s >= MPD_RADIX);
64	n/a	w[i] = carry ? s-MPD_RADIX : s;
65	n/a	}
66	n/a	/* if there is a carry, propagate it */
67	n/a	for (; carry && i < m; i++) {
68	n/a	s = u[i] + carry;
69	n/a	carry = (s == MPD_RADIX);
70	n/a	w[i] = carry ? 0 : s;
71	n/a	}
72	n/a	/* copy the rest of u */
73	n/a	for (; i < m; i++) {
74	n/a	w[i] = u[i];
75	n/a	}
76	n/a
77	n/a	return carry;
78	n/a	}
79	n/a
80	n/a	/*
81	n/a	* Add the contents of u to w. Carries are propagated further. The caller
82	n/a	* has to make sure that w is big enough.
83	n/a	*/
84	n/a	void
85	n/a	_mpd_baseaddto(mpd_uint_t w, const mpd_uint_t u, mpd_size_t n)
86	n/a	{
87	n/a	mpd_uint_t s;
88	n/a	mpd_uint_t carry = 0;
89	n/a	mpd_size_t i;
90	n/a
91	n/a	if (n == 0) return;
92	n/a
93	n/a	/* add n members of u to w */
94	n/a	for (i = 0; i < n; i++) {
95	n/a	s = w[i] + (u[i] + carry);
96	n/a	carry = (s < w[i]) \| (s >= MPD_RADIX);
97	n/a	w[i] = carry ? s-MPD_RADIX : s;
98	n/a	}
99	n/a	/* if there is a carry, propagate it */
100	n/a	for (; carry; i++) {
101	n/a	s = w[i] + carry;
102	n/a	carry = (s == MPD_RADIX);
103	n/a	w[i] = carry ? 0 : s;
104	n/a	}
105	n/a	}
106	n/a
107	n/a	/*
108	n/a	* Add v to w (len m). The calling function has to handle a possible
109	n/a	* final carry. Assumption: m > 0.
110	n/a	*/
111	n/a	mpd_uint_t
112	n/a	_mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v)
113	n/a	{
114	n/a	mpd_uint_t s;
115	n/a	mpd_uint_t carry;
116	n/a	mpd_size_t i;
117	n/a
118	n/a	assert(m > 0);
119	n/a
120	n/a	/* add v to w */
121	n/a	s = w[0] + v;
122	n/a	carry = (s < v) \| (s >= MPD_RADIX);
123	n/a	w[0] = carry ? s-MPD_RADIX : s;
124	n/a
125	n/a	/* if there is a carry, propagate it */
126	n/a	for (i = 1; carry && i < m; i++) {
127	n/a	s = w[i] + carry;
128	n/a	carry = (s == MPD_RADIX);
129	n/a	w[i] = carry ? 0 : s;
130	n/a	}
131	n/a
132	n/a	return carry;
133	n/a	}
134	n/a
135	n/a	/* Increment u. The calling function has to handle a possible carry. */
136	n/a	mpd_uint_t
137	n/a	_mpd_baseincr(mpd_uint_t *u, mpd_size_t n)
138	n/a	{
139	n/a	mpd_uint_t s;
140	n/a	mpd_uint_t carry = 1;
141	n/a	mpd_size_t i;
142	n/a
143	n/a	assert(n > 0);
144	n/a
145	n/a	/* if there is a carry, propagate it */
146	n/a	for (i = 0; carry && i < n; i++) {
147	n/a	s = u[i] + carry;
148	n/a	carry = (s == MPD_RADIX);
149	n/a	u[i] = carry ? 0 : s;
150	n/a	}
151	n/a
152	n/a	return carry;
153	n/a	}
154	n/a
155	n/a	/*
156	n/a	* Knuth, TAOCP, Volume 2, 4.3.1:
157	n/a	* w := difference of u (len m) and v (len n).
158	n/a	* number in u >= number in v;
159	n/a	*/
160	n/a	void
161	n/a	_mpd_basesub(mpd_uint_t w, const mpd_uint_t u, const mpd_uint_t *v,
162	n/a	mpd_size_t m, mpd_size_t n)
163	n/a	{
164	n/a	mpd_uint_t d;
165	n/a	mpd_uint_t borrow = 0;
166	n/a	mpd_size_t i;
167	n/a
168	n/a	assert(m > 0 && n > 0);
169	n/a
170	n/a	/* subtract n members of v from u */
171	n/a	for (i = 0; i < n; i++) {
172	n/a	d = u[i] - (v[i] + borrow);
173	n/a	borrow = (u[i] < d);
174	n/a	w[i] = borrow ? d + MPD_RADIX : d;
175	n/a	}
176	n/a	/* if there is a borrow, propagate it */
177	n/a	for (; borrow && i < m; i++) {
178	n/a	d = u[i] - borrow;
179	n/a	borrow = (u[i] == 0);
180	n/a	w[i] = borrow ? MPD_RADIX-1 : d;
181	n/a	}
182	n/a	/* copy the rest of u */
183	n/a	for (; i < m; i++) {
184	n/a	w[i] = u[i];
185	n/a	}
186	n/a	}
187	n/a
188	n/a	/*
189	n/a	* Subtract the contents of u from w. w is larger than u. Borrows are
190	n/a	* propagated further, but eventually w can absorb the final borrow.
191	n/a	*/
192	n/a	void
193	n/a	_mpd_basesubfrom(mpd_uint_t w, const mpd_uint_t u, mpd_size_t n)
194	n/a	{
195	n/a	mpd_uint_t d;
196	n/a	mpd_uint_t borrow = 0;
197	n/a	mpd_size_t i;
198	n/a
199	n/a	if (n == 0) return;
200	n/a
201	n/a	/* subtract n members of u from w */
202	n/a	for (i = 0; i < n; i++) {
203	n/a	d = w[i] - (u[i] + borrow);
204	n/a	borrow = (w[i] < d);
205	n/a	w[i] = borrow ? d + MPD_RADIX : d;
206	n/a	}
207	n/a	/* if there is a borrow, propagate it */
208	n/a	for (; borrow; i++) {
209	n/a	d = w[i] - borrow;
210	n/a	borrow = (w[i] == 0);
211	n/a	w[i] = borrow ? MPD_RADIX-1 : d;
212	n/a	}
213	n/a	}
214	n/a
215	n/a	/* w := product of u (len n) and v (single word) */
216	n/a	void
217	n/a	_mpd_shortmul(mpd_uint_t w, const mpd_uint_t u, mpd_size_t n, mpd_uint_t v)
218	n/a	{
219	n/a	mpd_uint_t hi, lo;
220	n/a	mpd_uint_t carry = 0;
221	n/a	mpd_size_t i;
222	n/a
223	n/a	assert(n > 0);
224	n/a
225	n/a	for (i=0; i < n; i++) {
226	n/a
227	n/a	_mpd_mul_words(&hi, &lo, u[i], v);
228	n/a	lo = carry + lo;
229	n/a	if (lo < carry) hi++;
230	n/a
231	n/a	_mpd_div_words_r(&carry, &w[i], hi, lo);
232	n/a	}
233	n/a	w[i] = carry;
234	n/a	}
235	n/a
236	n/a	/*
237	n/a	* Knuth, TAOCP, Volume 2, 4.3.1:
238	n/a	* w := product of u (len m) and v (len n)
239	n/a	* w must be initialized to zero
240	n/a	*/
241	n/a	void
242	n/a	_mpd_basemul(mpd_uint_t w, const mpd_uint_t u, const mpd_uint_t *v,
243	n/a	mpd_size_t m, mpd_size_t n)
244	n/a	{
245	n/a	mpd_uint_t hi, lo;
246	n/a	mpd_uint_t carry;
247	n/a	mpd_size_t i, j;
248	n/a
249	n/a	assert(m > 0 && n > 0);
250	n/a
251	n/a	for (j=0; j < n; j++) {
252	n/a	carry = 0;
253	n/a	for (i=0; i < m; i++) {
254	n/a
255	n/a	_mpd_mul_words(&hi, &lo, u[i], v[j]);
256	n/a	lo = w[i+j] + lo;
257	n/a	if (lo < w[i+j]) hi++;
258	n/a	lo = carry + lo;
259	n/a	if (lo < carry) hi++;
260	n/a
261	n/a	_mpd_div_words_r(&carry, &w[i+j], hi, lo);
262	n/a	}
263	n/a	w[j+m] = carry;
264	n/a	}
265	n/a	}
266	n/a
267	n/a	/*
268	n/a	* Knuth, TAOCP Volume 2, 4.3.1, exercise 16:
269	n/a	* w := quotient of u (len n) divided by a single word v
270	n/a	*/
271	n/a	mpd_uint_t
272	n/a	_mpd_shortdiv(mpd_uint_t w, const mpd_uint_t u, mpd_size_t n, mpd_uint_t v)
273	n/a	{
274	n/a	mpd_uint_t hi, lo;
275	n/a	mpd_uint_t rem = 0;
276	n/a	mpd_size_t i;
277	n/a
278	n/a	assert(n > 0);
279	n/a
280	n/a	for (i=n-1; i != MPD_SIZE_MAX; i--) {
281	n/a
282	n/a	_mpd_mul_words(&hi, &lo, rem, MPD_RADIX);
283	n/a	lo = u[i] + lo;
284	n/a	if (lo < u[i]) hi++;
285	n/a
286	n/a	_mpd_div_words(&w[i], &rem, hi, lo, v);
287	n/a	}
288	n/a
289	n/a	return rem;
290	n/a	}
291	n/a
292	n/a	/*
293	n/a	* Knuth, TAOCP Volume 2, 4.3.1:
294	n/a	* q, r := quotient and remainder of uconst (len nplusm)
295	n/a	* divided by vconst (len n)
296	n/a	* nplusm >= n
297	n/a	*
298	n/a	* If r is not NULL, r will contain the remainder. If r is NULL, the
299	n/a	* return value indicates if there is a remainder: 1 for true, 0 for
300	n/a	* false. A return value of -1 indicates an error.
301	n/a	*/
302	n/a	int
303	n/a	_mpd_basedivmod(mpd_uint_t q, mpd_uint_t r,
304	n/a	const mpd_uint_t uconst, const mpd_uint_t vconst,
305	n/a	mpd_size_t nplusm, mpd_size_t n)
306	n/a	{
307	n/a	mpd_uint_t ustatic[MPD_MINALLOC_MAX];
308	n/a	mpd_uint_t vstatic[MPD_MINALLOC_MAX];
309	n/a	mpd_uint_t *u = ustatic;
310	n/a	mpd_uint_t *v = vstatic;
311	n/a	mpd_uint_t d, qhat, rhat, w2[2];
312	n/a	mpd_uint_t hi, lo, x;
313	n/a	mpd_uint_t carry;
314	n/a	mpd_size_t i, j, m;
315	n/a	int retval = 0;
316	n/a
317	n/a	assert(n > 1 && nplusm >= n);
318	n/a	m = sub_size_t(nplusm, n);
319	n/a
320	n/a	/* D1: normalize */
321	n/a	d = MPD_RADIX / (vconst[n-1] + 1);
322	n/a
323	n/a	if (nplusm >= MPD_MINALLOC_MAX) {
324	n/a	if ((u = mpd_alloc(nplusm+1, sizeof *u)) == NULL) {
325	n/a	return -1;
326	n/a	}
327	n/a	}
328	n/a	if (n >= MPD_MINALLOC_MAX) {
329	n/a	if ((v = mpd_alloc(n+1, sizeof *v)) == NULL) {
330	n/a	mpd_free(u);
331	n/a	return -1;
332	n/a	}
333	n/a	}
334	n/a
335	n/a	_mpd_shortmul(u, uconst, nplusm, d);
336	n/a	_mpd_shortmul(v, vconst, n, d);
337	n/a
338	n/a	/* D2: loop */
339	n/a	for (j=m; j != MPD_SIZE_MAX; j--) {
340	n/a
341	n/a	/* D3: calculate qhat and rhat */
342	n/a	rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]);
343	n/a	qhat = w2[1] * MPD_RADIX + w2[0];
344	n/a
345	n/a	while (1) {
346	n/a	if (qhat < MPD_RADIX) {
347	n/a	_mpd_singlemul(w2, qhat, v[n-2]);
348	n/a	if (w2[1] <= rhat) {
349	n/a	if (w2[1] != rhat \|\| w2[0] <= u[j+n-2]) {
350	n/a	break;
351	n/a	}
352	n/a	}
353	n/a	}
354	n/a	qhat -= 1;
355	n/a	rhat += v[n-1];
356	n/a	if (rhat < v[n-1] \|\| rhat >= MPD_RADIX) {
357	n/a	break;
358	n/a	}
359	n/a	}
360	n/a	/* D4: multiply and subtract */
361	n/a	carry = 0;
362	n/a	for (i=0; i <= n; i++) {
363	n/a
364	n/a	_mpd_mul_words(&hi, &lo, qhat, v[i]);
365	n/a
366	n/a	lo = carry + lo;
367	n/a	if (lo < carry) hi++;
368	n/a
369	n/a	_mpd_div_words_r(&hi, &lo, hi, lo);
370	n/a
371	n/a	x = u[i+j] - lo;
372	n/a	carry = (u[i+j] < x);
373	n/a	u[i+j] = carry ? x+MPD_RADIX : x;
374	n/a	carry += hi;
375	n/a	}
376	n/a	q[j] = qhat;
377	n/a	/* D5: test remainder */
378	n/a	if (carry) {
379	n/a	q[j] -= 1;
380	n/a	/* D6: add back */
381	n/a	(void)_mpd_baseadd(u+j, u+j, v, n+1, n);
382	n/a	}
383	n/a	}
384	n/a
385	n/a	/* D8: unnormalize */
386	n/a	if (r != NULL) {
387	n/a	_mpd_shortdiv(r, u, n, d);
388	n/a	/* we are not interested in the return value here */
389	n/a	retval = 0;
390	n/a	}
391	n/a	else {
392	n/a	retval = !_mpd_isallzero(u, n);
393	n/a	}
394	n/a
395	n/a
396	n/a	if (u != ustatic) mpd_free(u);
397	n/a	if (v != vstatic) mpd_free(v);
398	n/a	return retval;
399	n/a	}
400	n/a
401	n/a	/*
402	n/a	* Left shift of src by 'shift' digits; src may equal dest.
403	n/a	*
404	n/a	* dest := area of n mpd_uint_t with space for srcdigits+shift digits.
405	n/a	* src := coefficient with length m.
406	n/a	*
407	n/a	* The case splits in the function are non-obvious. The following
408	n/a	* equations might help:
409	n/a	*
410	n/a	* Let msdigits denote the number of digits in the most significant
411	n/a	* word of src. Then 1 <= msdigits <= rdigits.
412	n/a	*
413	n/a	* 1) shift = q * rdigits + r
414	n/a	* 2) srcdigits = qsrc * rdigits + msdigits
415	n/a	* 3) destdigits = shift + srcdigits
416	n/a	* = q * rdigits + r + qsrc * rdigits + msdigits
417	n/a	* = q * rdigits + (qsrc * rdigits + (r + msdigits))
418	n/a	*
419	n/a	* The result has q zero words, followed by the coefficient that
420	n/a	* is left-shifted by r. The case r == 0 is trivial. For r > 0, it
421	n/a	* is important to keep in mind that we always read m source words,
422	n/a	* but write m+1 destination words if r + msdigits > rdigits, m words
423	n/a	* otherwise.
424	n/a	*/
425	n/a	void
426	n/a	_mpd_baseshiftl(mpd_uint_t dest, mpd_uint_t src, mpd_size_t n, mpd_size_t m,
427	n/a	mpd_size_t shift)
428	n/a	{
429	n/a	#if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
430	n/a	/* spurious uninitialized warnings */
431	n/a	mpd_uint_t l=l, lprev=lprev, h=h;
432	n/a	#else
433	n/a	mpd_uint_t l, lprev, h;
434	n/a	#endif
435	n/a	mpd_uint_t q, r;
436	n/a	mpd_uint_t ph;
437	n/a
438	n/a	assert(m > 0 && n >= m);
439	n/a
440	n/a	_mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
441	n/a
442	n/a	if (r != 0) {
443	n/a
444	n/a	ph = mpd_pow10[r];
445	n/a
446	n/a	--m; --n;
447	n/a	_mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r);
448	n/a	if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */
449	n/a	dest[n--] = h;
450	n/a	}
451	n/a	/* write m-1 shifted words */
452	n/a	for (; m != MPD_SIZE_MAX; m--,n--) {
453	n/a	_mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r);
454	n/a	dest[n] = ph * lprev + h;
455	n/a	lprev = l;
456	n/a	}
457	n/a	/* write least significant word */
458	n/a	dest[q] = ph * lprev;
459	n/a	}
460	n/a	else {
461	n/a	while (--m != MPD_SIZE_MAX) {
462	n/a	dest[m+q] = src[m];
463	n/a	}
464	n/a	}
465	n/a
466	n/a	mpd_uint_zero(dest, q);
467	n/a	}
468	n/a
469	n/a	/*
470	n/a	* Right shift of src by 'shift' digits; src may equal dest.
471	n/a	* Assumption: srcdigits-shift > 0.
472	n/a	*
473	n/a	* dest := area with space for srcdigits-shift digits.
474	n/a	* src := coefficient with length 'slen'.
475	n/a	*
476	n/a	* The case splits in the function rely on the following equations:
477	n/a	*
478	n/a	* Let msdigits denote the number of digits in the most significant
479	n/a	* word of src. Then 1 <= msdigits <= rdigits.
480	n/a	*
481	n/a	* 1) shift = q * rdigits + r
482	n/a	* 2) srcdigits = qsrc * rdigits + msdigits
483	n/a	* 3) destdigits = srcdigits - shift
484	n/a	* = qsrc * rdigits + msdigits - (q * rdigits + r)
485	n/a	* = (qsrc - q) * rdigits + msdigits - r
486	n/a	*
487	n/a	* Since destdigits > 0 and 1 <= msdigits <= rdigits:
488	n/a	*
489	n/a	* 4) qsrc >= q
490	n/a	* 5) qsrc == q ==> msdigits > r
491	n/a	*
492	n/a	* The result has slen-q words if msdigits > r, slen-q-1 words otherwise.
493	n/a	*/
494	n/a	mpd_uint_t
495	n/a	_mpd_baseshiftr(mpd_uint_t dest, mpd_uint_t src, mpd_size_t slen,
496	n/a	mpd_size_t shift)
497	n/a	{
498	n/a	#if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
499	n/a	/* spurious uninitialized warnings */
500	n/a	mpd_uint_t l=l, h=h, hprev=hprev; /* low, high, previous high */
501	n/a	#else
502	n/a	mpd_uint_t l, h, hprev; /* low, high, previous high */
503	n/a	#endif
504	n/a	mpd_uint_t rnd, rest; /* rounding digit, rest */
505	n/a	mpd_uint_t q, r;
506	n/a	mpd_size_t i, j;
507	n/a	mpd_uint_t ph;
508	n/a
509	n/a	assert(slen > 0);
510	n/a
511	n/a	_mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
512	n/a
513	n/a	rnd = rest = 0;
514	n/a	if (r != 0) {
515	n/a
516	n/a	ph = mpd_pow10[MPD_RDIGITS-r];
517	n/a
518	n/a	_mpd_divmod_pow10(&hprev, &rest, src[q], r);
519	n/a	_mpd_divmod_pow10(&rnd, &rest, rest, r-1);
520	n/a
521	n/a	if (rest == 0 && q > 0) {
522	n/a	rest = !_mpd_isallzero(src, q);
523	n/a	}
524	n/a	/* write slen-q-1 words */
525	n/a	for (j=0,i=q+1; i<slen; i++,j++) {
526	n/a	_mpd_divmod_pow10(&h, &l, src[i], r);
527	n/a	dest[j] = ph * l + hprev;
528	n/a	hprev = h;
529	n/a	}
530	n/a	/* write most significant word */
531	n/a	if (hprev != 0) { /* always the case if slen==q-1 */
532	n/a	dest[j] = hprev;
533	n/a	}
534	n/a	}
535	n/a	else {
536	n/a	if (q > 0) {
537	n/a	_mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1);
538	n/a	/* is there any non-zero digit below rnd? */
539	n/a	if (rest == 0) rest = !_mpd_isallzero(src, q-1);
540	n/a	}
541	n/a	for (j = 0; j < slen-q; j++) {
542	n/a	dest[j] = src[q+j];
543	n/a	}
544	n/a	}
545	n/a
546	n/a	/* 0-4 ==> rnd+rest < 0.5 */
547	n/a	/* 5 ==> rnd+rest == 0.5 */
548	n/a	/* 6-9 ==> rnd+rest > 0.5 */
549	n/a	return (rnd == 0 \|\| rnd == 5) ? rnd + !!rest : rnd;
550	n/a	}
551	n/a
552	n/a
553	n/a	/*********************************************************************/
554	n/a	/* Calculations in base b */
555	n/a	/*********************************************************************/
556	n/a
557	n/a	/*
558	n/a	* Add v to w (len m). The calling function has to handle a possible
559	n/a	* final carry. Assumption: m > 0.
560	n/a	*/
561	n/a	mpd_uint_t
562	n/a	_mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, mpd_uint_t b)
563	n/a	{
564	n/a	mpd_uint_t s;
565	n/a	mpd_uint_t carry;
566	n/a	mpd_size_t i;
567	n/a
568	n/a	assert(m > 0);
569	n/a
570	n/a	/* add v to w */
571	n/a	s = w[0] + v;
572	n/a	carry = (s < v) \| (s >= b);
573	n/a	w[0] = carry ? s-b : s;
574	n/a
575	n/a	/* if there is a carry, propagate it */
576	n/a	for (i = 1; carry && i < m; i++) {
577	n/a	s = w[i] + carry;
578	n/a	carry = (s == b);
579	n/a	w[i] = carry ? 0 : s;
580	n/a	}
581	n/a
582	n/a	return carry;
583	n/a	}
584	n/a
585	n/a	/* w := product of u (len n) and v (single word). Return carry. */
586	n/a	mpd_uint_t
587	n/a	_mpd_shortmul_c(mpd_uint_t w, const mpd_uint_t u, mpd_size_t n, mpd_uint_t v)
588	n/a	{
589	n/a	mpd_uint_t hi, lo;
590	n/a	mpd_uint_t carry = 0;
591	n/a	mpd_size_t i;
592	n/a
593	n/a	assert(n > 0);
594	n/a
595	n/a	for (i=0; i < n; i++) {
596	n/a
597	n/a	_mpd_mul_words(&hi, &lo, u[i], v);
598	n/a	lo = carry + lo;
599	n/a	if (lo < carry) hi++;
600	n/a
601	n/a	_mpd_div_words_r(&carry, &w[i], hi, lo);
602	n/a	}
603	n/a
604	n/a	return carry;
605	n/a	}
606	n/a
607	n/a	/* w := product of u (len n) and v (single word) */
608	n/a	mpd_uint_t
609	n/a	_mpd_shortmul_b(mpd_uint_t w, const mpd_uint_t u, mpd_size_t n,
610	n/a	mpd_uint_t v, mpd_uint_t b)
611	n/a	{
612	n/a	mpd_uint_t hi, lo;
613	n/a	mpd_uint_t carry = 0;
614	n/a	mpd_size_t i;
615	n/a
616	n/a	assert(n > 0);
617	n/a
618	n/a	for (i=0; i < n; i++) {
619	n/a
620	n/a	_mpd_mul_words(&hi, &lo, u[i], v);
621	n/a	lo = carry + lo;
622	n/a	if (lo < carry) hi++;
623	n/a
624	n/a	_mpd_div_words(&carry, &w[i], hi, lo, b);
625	n/a	}
626	n/a
627	n/a	return carry;
628	n/a	}
629	n/a
630	n/a	/*
631	n/a	* Knuth, TAOCP Volume 2, 4.3.1, exercise 16:
632	n/a	* w := quotient of u (len n) divided by a single word v
633	n/a	*/
634	n/a	mpd_uint_t
635	n/a	_mpd_shortdiv_b(mpd_uint_t w, const mpd_uint_t u, mpd_size_t n,
636	n/a	mpd_uint_t v, mpd_uint_t b)
637	n/a	{
638	n/a	mpd_uint_t hi, lo;
639	n/a	mpd_uint_t rem = 0;
640	n/a	mpd_size_t i;
641	n/a
642	n/a	assert(n > 0);
643	n/a
644	n/a	for (i=n-1; i != MPD_SIZE_MAX; i--) {
645	n/a
646	n/a	_mpd_mul_words(&hi, &lo, rem, b);
647	n/a	lo = u[i] + lo;
648	n/a	if (lo < u[i]) hi++;
649	n/a
650	n/a	_mpd_div_words(&w[i], &rem, hi, lo, v);
651	n/a	}
652	n/a
653	n/a	return rem;
654	n/a	}
655	n/a
656	n/a
657	n/a