Python code coverage for Modules/_decimal/libmpdec/crt.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 <assert.h>
32	n/a	#include "numbertheory.h"
33	n/a	#include "umodarith.h"
34	n/a	#include "crt.h"
35	n/a
36	n/a
37	n/a	/* Bignum: Chinese Remainder Theorem, extends the maximum transform length. */
38	n/a
39	n/a
40	n/a	/* Multiply P1P2 by v, store result in w. */
41	n/a	static inline void
42	n/a	_crt_mulP1P2_3(mpd_uint_t w[3], mpd_uint_t v)
43	n/a	{
44	n/a	mpd_uint_t hi1, hi2, lo;
45	n/a
46	n/a	_mpd_mul_words(&hi1, &lo, LH_P1P2, v);
47	n/a	w[0] = lo;
48	n/a
49	n/a	_mpd_mul_words(&hi2, &lo, UH_P1P2, v);
50	n/a	lo = hi1 + lo;
51	n/a	if (lo < hi1) hi2++;
52	n/a
53	n/a	w[1] = lo;
54	n/a	w[2] = hi2;
55	n/a	}
56	n/a
57	n/a	/* Add 3 words from v to w. The result is known to fit in w. */
58	n/a	static inline void
59	n/a	_crt_add3(mpd_uint_t w[3], mpd_uint_t v[3])
60	n/a	{
61	n/a	mpd_uint_t carry;
62	n/a	mpd_uint_t s;
63	n/a
64	n/a	s = w[0] + v[0];
65	n/a	carry = (s < w[0]);
66	n/a	w[0] = s;
67	n/a
68	n/a	s = w[1] + (v[1] + carry);
69	n/a	carry = (s < w[1]);
70	n/a	w[1] = s;
71	n/a
72	n/a	w[2] = w[2] + (v[2] + carry);
73	n/a	}
74	n/a
75	n/a	/* Divide 3 words in u by v, store result in w, return remainder. */
76	n/a	static inline mpd_uint_t
77	n/a	_crt_div3(mpd_uint_t w, const mpd_uint_t u, mpd_uint_t v)
78	n/a	{
79	n/a	mpd_uint_t r1 = u[2];
80	n/a	mpd_uint_t r2;
81	n/a
82	n/a	if (r1 < v) {
83	n/a	w[2] = 0;
84	n/a	}
85	n/a	else {
86	n/a	_mpd_div_word(&w[2], &r1, u[2], v); /* GCOV_NOT_REACHED */
87	n/a	}
88	n/a
89	n/a	_mpd_div_words(&w[1], &r2, r1, u[1], v);
90	n/a	_mpd_div_words(&w[0], &r1, r2, u[0], v);
91	n/a
92	n/a	return r1;
93	n/a	}
94	n/a
95	n/a
96	n/a	/*
97	n/a	* Chinese Remainder Theorem:
98	n/a	* Algorithm from Joerg Arndt, "Matters Computational",
99	n/a	* Chapter 37.4.1 [http://www.jjj.de/fxt/]
100	n/a	*
101	n/a	* See also Knuth, TAOCP, Volume 2, 4.3.2, exercise 7.
102	n/a	*/
103	n/a
104	n/a	/*
105	n/a	* CRT with carry: x1, x2, x3 contain numbers modulo p1, p2, p3. For each
106	n/a	* triple of members of the arrays, find the unique z modulo p1p2p3, with
107	n/a	* zmax = p1p2p3 - 1.
108	n/a	*
109	n/a	* In each iteration of the loop, split z into result[i] = z % MPD_RADIX
110	n/a	* and carry = z / MPD_RADIX. Let N be the size of carry[] and cmax the
111	n/a	* maximum carry.
112	n/a	*
113	n/a	* Limits for the 32-bit build:
114	n/a	*
115	n/a	* N = 2**96
116	n/a	* cmax = 7711435591312380274
117	n/a	*
118	n/a	* Limits for the 64 bit build:
119	n/a	*
120	n/a	* N = 2**192
121	n/a	* cmax = 627710135393475385904124401220046371710
122	n/a	*
123	n/a	* The following statements hold for both versions:
124	n/a	*
125	n/a	* 1) cmax + zmax < N, so the addition does not overflow.
126	n/a	*
127	n/a	* 2) (cmax + zmax) / MPD_RADIX == cmax.
128	n/a	*
129	n/a	* 3) If c <= cmax, then c_next = (c + zmax) / MPD_RADIX <= cmax.
130	n/a	*/
131	n/a	void
132	n/a	crt3(mpd_uint_t x1, mpd_uint_t x2, mpd_uint_t *x3, mpd_size_t rsize)
133	n/a	{
134	n/a	mpd_uint_t p1 = mpd_moduli[P1];
135	n/a	mpd_uint_t umod;
136	n/a	#ifdef PPRO
137	n/a	double dmod;
138	n/a	uint32_t dinvmod[3];
139	n/a	#endif
140	n/a	mpd_uint_t a1, a2, a3;
141	n/a	mpd_uint_t s;
142	n/a	mpd_uint_t z[3], t[3];
143	n/a	mpd_uint_t carry[3] = {0,0,0};
144	n/a	mpd_uint_t hi, lo;
145	n/a	mpd_size_t i;
146	n/a
147	n/a	for (i = 0; i < rsize; i++) {
148	n/a
149	n/a	a1 = x1[i];
150	n/a	a2 = x2[i];
151	n/a	a3 = x3[i];
152	n/a
153	n/a	SETMODULUS(P2);
154	n/a	s = ext_submod(a2, a1, umod);
155	n/a	s = MULMOD(s, INV_P1_MOD_P2);
156	n/a
157	n/a	_mpd_mul_words(&hi, &lo, s, p1);
158	n/a	lo = lo + a1;
159	n/a	if (lo < a1) hi++;
160	n/a
161	n/a	SETMODULUS(P3);
162	n/a	s = dw_submod(a3, hi, lo, umod);
163	n/a	s = MULMOD(s, INV_P1P2_MOD_P3);
164	n/a
165	n/a	z[0] = lo;
166	n/a	z[1] = hi;
167	n/a	z[2] = 0;
168	n/a
169	n/a	_crt_mulP1P2_3(t, s);
170	n/a	_crt_add3(z, t);
171	n/a	_crt_add3(carry, z);
172	n/a
173	n/a	x1[i] = _crt_div3(carry, carry, MPD_RADIX);
174	n/a	}
175	n/a
176	n/a	assert(carry[0] == 0 && carry[1] == 0 && carry[2] == 0);
177	n/a	}
178	n/a
179	n/a