Branch data Line data Source code
1 : : /*
2 : : * This file is part of the MicroPython project, http://micropython.org/
3 : : *
4 : : * The MIT License (MIT)
5 : : *
6 : : * Copyright (c) 2013, 2014 Damien P. George
7 : : *
8 : : * Permission is hereby granted, free of charge, to any person obtaining a copy
9 : : * of this software and associated documentation files (the "Software"), to deal
10 : : * in the Software without restriction, including without limitation the rights
11 : : * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 : : * copies of the Software, and to permit persons to whom the Software is
13 : : * furnished to do so, subject to the following conditions:
14 : : *
15 : : * The above copyright notice and this permission notice shall be included in
16 : : * all copies or substantial portions of the Software.
17 : : *
18 : : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 : : * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 : : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 : : * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 : : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 : : * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24 : : * THE SOFTWARE.
25 : : */
26 : :
27 : : #include <string.h>
28 : : #include <assert.h>
29 : :
30 : : #include "py/mpz.h"
31 : :
32 : : #if MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ
33 : :
34 : : #define DIG_SIZE (MPZ_DIG_SIZE)
35 : : #define DIG_MASK ((MPZ_LONG_1 << DIG_SIZE) - 1)
36 : : #define DIG_MSB (MPZ_LONG_1 << (DIG_SIZE - 1))
37 : : #define DIG_BASE (MPZ_LONG_1 << DIG_SIZE)
38 : :
39 : : /*
40 : : mpz is an arbitrary precision integer type with a public API.
41 : :
42 : : mpn functions act on non-negative integers represented by an array of generalised
43 : : digits (eg a word per digit). You also need to specify separately the length of the
44 : : array. There is no public API for mpn. Rather, the functions are used by mpz to
45 : : implement its features.
46 : :
47 : : Integer values are stored little endian (first digit is first in memory).
48 : :
49 : : Definition of normalise: ?
50 : : */
51 : :
52 : 914 : static size_t mpn_remove_trailing_zeros(mpz_dig_t *oidig, mpz_dig_t *idig) {
53 [ + + + + ]: 2208 : for (--idig; idig >= oidig && *idig == 0; --idig) {
54 : 1294 : }
55 : 914 : return idig + 1 - oidig;
56 : : }
57 : :
58 : : /* compares i with j
59 : : returns sign(i - j)
60 : : assumes i, j are normalised
61 : : */
62 : 44234 : static int mpn_cmp(const mpz_dig_t *idig, size_t ilen, const mpz_dig_t *jdig, size_t jlen) {
63 [ + + ]: 44234 : if (ilen < jlen) {
64 : : return -1;
65 : : }
66 [ + + ]: 39658 : if (ilen > jlen) {
67 : : return 1;
68 : : }
69 : :
70 [ + + ]: 13618 : for (idig += ilen, jdig += ilen; ilen > 0; --ilen) {
71 : 12986 : mpz_dbl_dig_signed_t cmp = (mpz_dbl_dig_t)*(--idig) - (mpz_dbl_dig_t)*(--jdig);
72 [ + + ]: 12986 : if (cmp < 0) {
73 : : return -1;
74 : : }
75 [ + + ]: 10500 : if (cmp > 0) {
76 : : return 1;
77 : : }
78 : : }
79 : :
80 : : return 0;
81 : : }
82 : :
83 : : /* computes i = j << n
84 : : returns number of digits in i
85 : : assumes enough memory in i; assumes normalised j; assumes n > 0
86 : : can have i, j pointing to same memory
87 : : */
88 : 1258 : static size_t mpn_shl(mpz_dig_t *idig, mpz_dig_t *jdig, size_t jlen, mp_uint_t n) {
89 : 1258 : mp_uint_t n_whole = (n + DIG_SIZE - 1) / DIG_SIZE;
90 : 1258 : mp_uint_t n_part = n % DIG_SIZE;
91 [ + + ]: 1258 : if (n_part == 0) {
92 : 126 : n_part = DIG_SIZE;
93 : : }
94 : :
95 : : // start from the high end of the digit arrays
96 : 1258 : idig += jlen + n_whole - 1;
97 : 1258 : jdig += jlen - 1;
98 : :
99 : : // shift the digits
100 : 1258 : mpz_dbl_dig_t d = 0;
101 [ + + ]: 4836 : for (size_t i = jlen; i > 0; i--, idig--, jdig--) {
102 : 3578 : d |= *jdig;
103 : 3578 : *idig = (d >> (DIG_SIZE - n_part)) & DIG_MASK;
104 : 3578 : d <<= DIG_SIZE;
105 : : }
106 : :
107 : : // store remaining bits
108 : 1258 : *idig = (d >> (DIG_SIZE - n_part)) & DIG_MASK;
109 : 1258 : idig -= n_whole - 1;
110 : 1258 : memset(idig, 0, (n_whole - 1) * sizeof(mpz_dig_t));
111 : :
112 : : // work out length of result
113 : 1258 : jlen += n_whole;
114 [ + - + + ]: 2230 : while (jlen != 0 && idig[jlen - 1] == 0) {
115 : 972 : jlen--;
116 : : }
117 : :
118 : : // return length of result
119 : 1258 : return jlen;
120 : : }
121 : :
122 : : /* computes i = j >> n
123 : : returns number of digits in i
124 : : assumes enough memory in i; assumes normalised j; assumes n > 0
125 : : can have i, j pointing to same memory
126 : : */
127 : 16866 : static size_t mpn_shr(mpz_dig_t *idig, mpz_dig_t *jdig, size_t jlen, mp_uint_t n) {
128 : 16866 : mp_uint_t n_whole = n / DIG_SIZE;
129 : 16866 : mp_uint_t n_part = n % DIG_SIZE;
130 : :
131 [ + + ]: 16866 : if (n_whole >= jlen) {
132 : : return 0;
133 : : }
134 : :
135 : 16858 : jdig += n_whole;
136 : 16858 : jlen -= n_whole;
137 : :
138 [ + + ]: 435676 : for (size_t i = jlen; i > 0; i--, idig++, jdig++) {
139 : 418818 : mpz_dbl_dig_t d = *jdig;
140 [ + + ]: 418818 : if (i > 1) {
141 : 401960 : d |= (mpz_dbl_dig_t)jdig[1] << DIG_SIZE;
142 : : }
143 : 418818 : d >>= n_part;
144 : 418818 : *idig = d & DIG_MASK;
145 : : }
146 : :
147 [ + + ]: 16858 : if (idig[-1] == 0) {
148 : 1512 : jlen--;
149 : : }
150 : :
151 : : return jlen;
152 : : }
153 : :
154 : : /* computes i = j + k
155 : : returns number of digits in i
156 : : assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
157 : : can have i, j, k pointing to same memory
158 : : */
159 : 378 : static size_t mpn_add(mpz_dig_t *idig, const mpz_dig_t *jdig, size_t jlen, const mpz_dig_t *kdig, size_t klen) {
160 : 378 : mpz_dig_t *oidig = idig;
161 : 378 : mpz_dbl_dig_t carry = 0;
162 : :
163 : 378 : jlen -= klen;
164 : :
165 [ + + ]: 844 : for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
166 : 466 : carry += (mpz_dbl_dig_t)*jdig + (mpz_dbl_dig_t)*kdig;
167 : 466 : *idig = carry & DIG_MASK;
168 : 466 : carry >>= DIG_SIZE;
169 : : }
170 : :
171 [ + + ]: 1788 : for (; jlen > 0; --jlen, ++idig, ++jdig) {
172 : 1410 : carry += *jdig;
173 : 1410 : *idig = carry & DIG_MASK;
174 : 1410 : carry >>= DIG_SIZE;
175 : : }
176 : :
177 [ + + ]: 378 : if (carry != 0) {
178 : 6 : *idig++ = carry;
179 : : }
180 : :
181 : 378 : return idig - oidig;
182 : : }
183 : :
184 : : /* computes i = j - k
185 : : returns number of digits in i
186 : : assumes enough memory in i; assumes normalised j, k; assumes j >= k
187 : : can have i, j, k pointing to same memory
188 : : */
189 : 396 : static size_t mpn_sub(mpz_dig_t *idig, const mpz_dig_t *jdig, size_t jlen, const mpz_dig_t *kdig, size_t klen) {
190 : 396 : mpz_dig_t *oidig = idig;
191 : 396 : mpz_dbl_dig_signed_t borrow = 0;
192 : :
193 : 396 : jlen -= klen;
194 : :
195 [ + + ]: 2030 : for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
196 : 1634 : borrow += (mpz_dbl_dig_t)*jdig - (mpz_dbl_dig_t)*kdig;
197 : 1634 : *idig = borrow & DIG_MASK;
198 : 1634 : borrow >>= DIG_SIZE;
199 : : }
200 : :
201 [ + + ]: 1610 : for (; jlen > 0; --jlen, ++idig, ++jdig) {
202 : 1214 : borrow += *jdig;
203 : 1214 : *idig = borrow & DIG_MASK;
204 : 1214 : borrow >>= DIG_SIZE;
205 : : }
206 : :
207 : 396 : return mpn_remove_trailing_zeros(oidig, idig);
208 : : }
209 : :
210 : : #if MICROPY_OPT_MPZ_BITWISE
211 : :
212 : : /* computes i = j & k
213 : : returns number of digits in i
214 : : assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen (jlen argument not needed)
215 : : can have i, j, k pointing to same memory
216 : : */
217 : 92 : static size_t mpn_and(mpz_dig_t *idig, const mpz_dig_t *jdig, const mpz_dig_t *kdig, size_t klen) {
218 : 92 : mpz_dig_t *oidig = idig;
219 : :
220 [ + + ]: 468 : for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
221 : 376 : *idig = *jdig & *kdig;
222 : : }
223 : :
224 : 92 : return mpn_remove_trailing_zeros(oidig, idig);
225 : : }
226 : :
227 : : #endif
228 : :
229 : : /* i = -((-j) & (-k)) = ~((~j + 1) & (~k + 1)) + 1
230 : : i = (j & (-k)) = (j & (~k + 1)) = ( j & (~k + 1))
231 : : i = ((-j) & k) = ((~j + 1) & k) = ((~j + 1) & k )
232 : : computes general form:
233 : : i = (im ^ (((j ^ jm) + jc) & ((k ^ km) + kc))) + ic where Xm = Xc == 0 ? 0 : DIG_MASK
234 : : returns number of digits in i
235 : : assumes enough memory in i; assumes normalised j, k; assumes length j >= length k
236 : : can have i, j, k pointing to same memory
237 : : */
238 : 140 : static size_t mpn_and_neg(mpz_dig_t *idig, const mpz_dig_t *jdig, size_t jlen, const mpz_dig_t *kdig, size_t klen,
239 : : mpz_dbl_dig_t carryi, mpz_dbl_dig_t carryj, mpz_dbl_dig_t carryk) {
240 : 140 : mpz_dig_t *oidig = idig;
241 [ + + ]: 140 : mpz_dig_t imask = (0 == carryi) ? 0 : DIG_MASK;
242 [ + + ]: 140 : mpz_dig_t jmask = (0 == carryj) ? 0 : DIG_MASK;
243 [ + + ]: 140 : mpz_dig_t kmask = (0 == carryk) ? 0 : DIG_MASK;
244 : :
245 [ + + ]: 1446 : for (; jlen > 0; ++idig, ++jdig) {
246 : 1306 : carryj += *jdig ^ jmask;
247 [ + + ]: 1306 : carryk += (--klen <= --jlen) ? (*kdig++ ^ kmask) : kmask;
248 : 1306 : carryi += ((carryj & carryk) ^ imask) & DIG_MASK;
249 : 1306 : *idig = carryi & DIG_MASK;
250 : 1306 : carryk >>= DIG_SIZE;
251 : 1306 : carryj >>= DIG_SIZE;
252 : 1306 : carryi >>= DIG_SIZE;
253 : : }
254 : :
255 [ + + ]: 140 : if (0 != carryi) {
256 : 4 : *idig++ = carryi;
257 : : }
258 : :
259 : 140 : return mpn_remove_trailing_zeros(oidig, idig);
260 : : }
261 : :
262 : : #if MICROPY_OPT_MPZ_BITWISE
263 : :
264 : : /* computes i = j | k
265 : : returns number of digits in i
266 : : assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
267 : : can have i, j, k pointing to same memory
268 : : */
269 : 54 : static size_t mpn_or(mpz_dig_t *idig, const mpz_dig_t *jdig, size_t jlen, const mpz_dig_t *kdig, size_t klen) {
270 : 54 : mpz_dig_t *oidig = idig;
271 : :
272 : 54 : jlen -= klen;
273 : :
274 [ + + ]: 426 : for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
275 : 372 : *idig = *jdig | *kdig;
276 : : }
277 : :
278 [ + + ]: 200 : for (; jlen > 0; --jlen, ++idig, ++jdig) {
279 : 146 : *idig = *jdig;
280 : : }
281 : :
282 : 54 : return idig - oidig;
283 : : }
284 : :
285 : : #endif
286 : :
287 : : /* i = -((-j) | (-k)) = ~((~j + 1) | (~k + 1)) + 1
288 : : i = -(j | (-k)) = -(j | (~k + 1)) = ~( j | (~k + 1)) + 1
289 : : i = -((-j) | k) = -((~j + 1) | k) = ~((~j + 1) | k ) + 1
290 : : computes general form:
291 : : i = ~(((j ^ jm) + jc) | ((k ^ km) + kc)) + 1 where Xm = Xc == 0 ? 0 : DIG_MASK
292 : : returns number of digits in i
293 : : assumes enough memory in i; assumes normalised j, k; assumes length j >= length k
294 : : can have i, j, k pointing to same memory
295 : : */
296 : :
297 : : #if MICROPY_OPT_MPZ_BITWISE
298 : :
299 : 106 : static size_t mpn_or_neg(mpz_dig_t *idig, const mpz_dig_t *jdig, size_t jlen, const mpz_dig_t *kdig, size_t klen,
300 : : mpz_dbl_dig_t carryj, mpz_dbl_dig_t carryk) {
301 : 106 : mpz_dig_t *oidig = idig;
302 : 106 : mpz_dbl_dig_t carryi = 1;
303 [ + + ]: 106 : mpz_dig_t jmask = (0 == carryj) ? 0 : DIG_MASK;
304 [ + + ]: 106 : mpz_dig_t kmask = (0 == carryk) ? 0 : DIG_MASK;
305 : :
306 [ + + ]: 1282 : for (; jlen > 0; ++idig, ++jdig) {
307 : 1176 : carryj += *jdig ^ jmask;
308 [ + + ]: 1176 : carryk += (--klen <= --jlen) ? (*kdig++ ^ kmask) : kmask;
309 : 1176 : carryi += ((carryj | carryk) ^ DIG_MASK) & DIG_MASK;
310 : 1176 : *idig = carryi & DIG_MASK;
311 : 1176 : carryk >>= DIG_SIZE;
312 : 1176 : carryj >>= DIG_SIZE;
313 : 1176 : carryi >>= DIG_SIZE;
314 : : }
315 : :
316 : : // At least one of j,k must be negative so the above for-loop runs at least
317 : : // once. For carryi to be non-zero here it must be equal to 1 at the end of
318 : : // each iteration of the loop. So the accumulation of carryi must overflow
319 : : // each time, ie carryi += 0xff..ff. So carryj|carryk must be 0 in the
320 : : // DIG_MASK bits on each iteration. But considering all cases of signs of
321 : : // j,k one sees that this is not possible.
322 [ - + ]: 106 : assert(carryi == 0);
323 : :
324 : 106 : return mpn_remove_trailing_zeros(oidig, idig);
325 : : }
326 : :
327 : : #else
328 : :
329 : : static size_t mpn_or_neg(mpz_dig_t *idig, const mpz_dig_t *jdig, size_t jlen, const mpz_dig_t *kdig, size_t klen,
330 : : mpz_dbl_dig_t carryi, mpz_dbl_dig_t carryj, mpz_dbl_dig_t carryk) {
331 : : mpz_dig_t *oidig = idig;
332 : : mpz_dig_t imask = (0 == carryi) ? 0 : DIG_MASK;
333 : : mpz_dig_t jmask = (0 == carryj) ? 0 : DIG_MASK;
334 : : mpz_dig_t kmask = (0 == carryk) ? 0 : DIG_MASK;
335 : :
336 : : for (; jlen > 0; ++idig, ++jdig) {
337 : : carryj += *jdig ^ jmask;
338 : : carryk += (--klen <= --jlen) ? (*kdig++ ^ kmask) : kmask;
339 : : carryi += ((carryj | carryk) ^ imask) & DIG_MASK;
340 : : *idig = carryi & DIG_MASK;
341 : : carryk >>= DIG_SIZE;
342 : : carryj >>= DIG_SIZE;
343 : : carryi >>= DIG_SIZE;
344 : : }
345 : :
346 : : // See comment in above mpn_or_neg for why carryi must be 0.
347 : : assert(carryi == 0);
348 : :
349 : : return mpn_remove_trailing_zeros(oidig, idig);
350 : : }
351 : :
352 : : #endif
353 : :
354 : : #if MICROPY_OPT_MPZ_BITWISE
355 : :
356 : : /* computes i = j ^ k
357 : : returns number of digits in i
358 : : assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
359 : : can have i, j, k pointing to same memory
360 : : */
361 : 42 : static size_t mpn_xor(mpz_dig_t *idig, const mpz_dig_t *jdig, size_t jlen, const mpz_dig_t *kdig, size_t klen) {
362 : 42 : mpz_dig_t *oidig = idig;
363 : :
364 : 42 : jlen -= klen;
365 : :
366 [ + + ]: 376 : for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
367 : 334 : *idig = *jdig ^ *kdig;
368 : : }
369 : :
370 [ + + ]: 136 : for (; jlen > 0; --jlen, ++idig, ++jdig) {
371 : 94 : *idig = *jdig;
372 : : }
373 : :
374 : 42 : return mpn_remove_trailing_zeros(oidig, idig);
375 : : }
376 : :
377 : : #endif
378 : :
379 : : /* i = (-j) ^ (-k) = ~(j - 1) ^ ~(k - 1) = (j - 1) ^ (k - 1)
380 : : i = -(j ^ (-k)) = -(j ^ ~(k - 1)) = ~(j ^ ~(k - 1)) + 1 = (j ^ (k - 1)) + 1
381 : : i = -((-j) ^ k) = -(~(j - 1) ^ k) = ~(~(j - 1) ^ k) + 1 = ((j - 1) ^ k) + 1
382 : : computes general form:
383 : : i = ((j - 1 + jc) ^ (k - 1 + kc)) + ic
384 : : returns number of digits in i
385 : : assumes enough memory in i; assumes normalised j, k; assumes length j >= length k
386 : : can have i, j, k pointing to same memory
387 : : */
388 : 114 : static size_t mpn_xor_neg(mpz_dig_t *idig, const mpz_dig_t *jdig, size_t jlen, const mpz_dig_t *kdig, size_t klen,
389 : : mpz_dbl_dig_t carryi, mpz_dbl_dig_t carryj, mpz_dbl_dig_t carryk) {
390 : 114 : mpz_dig_t *oidig = idig;
391 : :
392 [ + + ]: 1326 : for (; jlen > 0; ++idig, ++jdig) {
393 : 1212 : carryj += *jdig + DIG_MASK;
394 [ + + ]: 1212 : carryk += (--klen <= --jlen) ? (*kdig++ + DIG_MASK) : DIG_MASK;
395 : 1212 : carryi += (carryj ^ carryk) & DIG_MASK;
396 : 1212 : *idig = carryi & DIG_MASK;
397 : 1212 : carryk >>= DIG_SIZE;
398 : 1212 : carryj >>= DIG_SIZE;
399 : 1212 : carryi >>= DIG_SIZE;
400 : : }
401 : :
402 [ + + ]: 114 : if (0 != carryi) {
403 : 4 : *idig++ = carryi;
404 : : }
405 : :
406 : 114 : return mpn_remove_trailing_zeros(oidig, idig);
407 : : }
408 : :
409 : : /* computes i = i * d1 + d2
410 : : returns number of digits in i
411 : : assumes enough memory in i; assumes normalised i; assumes dmul != 0
412 : : */
413 : 47746 : static size_t mpn_mul_dig_add_dig(mpz_dig_t *idig, size_t ilen, mpz_dig_t dmul, mpz_dig_t dadd) {
414 : 47746 : mpz_dig_t *oidig = idig;
415 : 47746 : mpz_dbl_dig_t carry = dadd;
416 : :
417 [ + + ]: 415941 : for (; ilen > 0; --ilen, ++idig) {
418 : 368195 : carry += (mpz_dbl_dig_t)*idig * (mpz_dbl_dig_t)dmul; // will never overflow so long as DIG_SIZE <= 8*sizeof(mpz_dbl_dig_t)/2
419 : 368195 : *idig = carry & DIG_MASK;
420 : 368195 : carry >>= DIG_SIZE;
421 : : }
422 : :
423 [ + + ]: 47746 : if (carry != 0) {
424 : 10522 : *idig++ = carry;
425 : : }
426 : :
427 : 47746 : return idig - oidig;
428 : : }
429 : :
430 : : /* computes i = j * k
431 : : returns number of digits in i
432 : : assumes enough memory in i; assumes i is zeroed; assumes normalised j, k
433 : : can have j, k point to same memory
434 : : */
435 : 25550 : static size_t mpn_mul(mpz_dig_t *idig, mpz_dig_t *jdig, size_t jlen, mpz_dig_t *kdig, size_t klen) {
436 : 25550 : mpz_dig_t *oidig = idig;
437 : 25550 : size_t ilen = 0;
438 : :
439 [ + + ]: 1204548 : for (; klen > 0; --klen, ++idig, ++kdig) {
440 : : mpz_dig_t *id = idig;
441 : : mpz_dbl_dig_t carry = 0;
442 : :
443 : : size_t jl = jlen;
444 [ + + ]: 81484372 : for (mpz_dig_t *jd = jdig; jl > 0; --jl, ++jd, ++id) {
445 : 80305374 : carry += (mpz_dbl_dig_t)*id + (mpz_dbl_dig_t)*jd * (mpz_dbl_dig_t)*kdig; // will never overflow so long as DIG_SIZE <= 8*sizeof(mpz_dbl_dig_t)/2
446 : 80305374 : *id = carry & DIG_MASK;
447 : 80305374 : carry >>= DIG_SIZE;
448 : : }
449 : :
450 [ + + ]: 1178998 : if (carry != 0) {
451 : 1164595 : *id++ = carry;
452 : : }
453 : :
454 : 1178998 : ilen = id - oidig;
455 : : }
456 : :
457 : 25550 : return ilen;
458 : : }
459 : :
460 : : /* natural_div - quo * den + new_num = old_num (ie num is replaced with rem)
461 : : assumes den != 0
462 : : assumes num_dig has enough memory to be extended by 1 digit
463 : : assumes quo_dig has enough memory (as many digits as num)
464 : : assumes quo_dig is filled with zeros
465 : : */
466 : 30330 : static void mpn_div(mpz_dig_t *num_dig, size_t *num_len, const mpz_dig_t *den_dig, size_t den_len, mpz_dig_t *quo_dig, size_t *quo_len) {
467 : 30330 : mpz_dig_t *orig_num_dig = num_dig;
468 : 30330 : mpz_dig_t *orig_quo_dig = quo_dig;
469 : 30330 : mpz_dig_t norm_shift = 0;
470 : 30330 : mpz_dbl_dig_t lead_den_digit;
471 : :
472 : : // handle simple cases
473 : : {
474 : 30330 : int cmp = mpn_cmp(num_dig, *num_len, den_dig, den_len);
475 [ + + ]: 30330 : if (cmp == 0) {
476 : 80 : *num_len = 0;
477 : 80 : quo_dig[0] = 1;
478 : 80 : *quo_len = 1;
479 : 80 : return;
480 [ + + ]: 30250 : } else if (cmp < 0) {
481 : : // numerator remains the same
482 : 4116 : *quo_len = 0;
483 : 4116 : return;
484 : : }
485 : : }
486 : :
487 : : // We need to normalise the denominator (leading bit of leading digit is 1)
488 : : // so that the division routine works. Since the denominator memory is
489 : : // read-only we do the normalisation on the fly, each time a digit of the
490 : : // denominator is needed. We need to know is how many bits to shift by.
491 : :
492 : : // count number of leading zeros in leading digit of denominator
493 : : {
494 : 26134 : mpz_dig_t d = den_dig[den_len - 1];
495 [ + + ]: 167314 : while ((d & DIG_MSB) == 0) {
496 : 141180 : d <<= 1;
497 : 141180 : ++norm_shift;
498 : : }
499 : : }
500 : :
501 : : // now need to shift numerator by same amount as denominator
502 : : // first, increase length of numerator in case we need more room to shift
503 : 26134 : num_dig[*num_len] = 0;
504 : 26134 : ++(*num_len);
505 [ + + ]: 2436876 : for (mpz_dig_t *num = num_dig, carry = 0; num < num_dig + *num_len; ++num) {
506 : 2410742 : mpz_dig_t n = *num;
507 : 2410742 : *num = ((n << norm_shift) | carry) & DIG_MASK;
508 : 2410742 : carry = (mpz_dbl_dig_t)n >> (DIG_SIZE - norm_shift);
509 : : }
510 : :
511 : : // cache the leading digit of the denominator
512 : 26134 : lead_den_digit = (mpz_dbl_dig_t)den_dig[den_len - 1] << norm_shift;
513 [ + + ]: 26134 : if (den_len >= 2) {
514 : 17198 : lead_den_digit |= (mpz_dbl_dig_t)den_dig[den_len - 2] >> (DIG_SIZE - norm_shift);
515 : : }
516 : :
517 : : // point num_dig to last digit in numerator
518 : 26134 : num_dig += *num_len - 1;
519 : :
520 : : // calculate number of digits in quotient
521 : 26134 : *quo_len = *num_len - den_len;
522 : :
523 : : // point to last digit to store for quotient
524 : 26134 : quo_dig += *quo_len - 1;
525 : :
526 : : // keep going while we have enough digits to divide
527 [ + + ]: 1262926 : while (*num_len > den_len) {
528 : 1236792 : mpz_dbl_dig_t quo = ((mpz_dbl_dig_t)*num_dig << DIG_SIZE) | num_dig[-1];
529 : :
530 : : // get approximate quotient
531 : 1236792 : quo /= lead_den_digit;
532 : :
533 : : // Multiply quo by den and subtract from num to get remainder.
534 : : // Must be careful with overflow of the borrow variable. Both
535 : : // borrow and low_digs are signed values and need signed right-shift,
536 : : // but x is unsigned and may take a full-range value.
537 : 1236792 : const mpz_dig_t *d = den_dig;
538 : 1236792 : mpz_dbl_dig_t d_norm = 0;
539 : 1236792 : mpz_dbl_dig_signed_t borrow = 0;
540 [ + + ]: 82587344 : for (mpz_dig_t *n = num_dig - den_len; n < num_dig; ++n, ++d) {
541 : : // Get the next digit in (den).
542 : 81350552 : d_norm = ((mpz_dbl_dig_t)*d << norm_shift) | (d_norm >> DIG_SIZE);
543 : : // Multiply the next digit in (quo * den).
544 : 81350552 : mpz_dbl_dig_t x = (mpz_dbl_dig_t)quo * (d_norm & DIG_MASK);
545 : : // Compute the low DIG_MASK bits of the next digit in (num - quo * den)
546 : 81350552 : mpz_dbl_dig_signed_t low_digs = (borrow & DIG_MASK) + *n - (x & DIG_MASK);
547 : : // Store the digit result for (num).
548 : 81350552 : *n = low_digs & DIG_MASK;
549 : : // Compute the borrow, shifted right before summing to avoid overflow.
550 : 81350552 : borrow = (borrow >> DIG_SIZE) - (x >> DIG_SIZE) + (low_digs >> DIG_SIZE);
551 : : }
552 : :
553 : : // At this point we have either:
554 : : //
555 : : // 1. quo was the correct value and the most-sig-digit of num is exactly
556 : : // cancelled by borrow (borrow + *num_dig == 0). In this case there is
557 : : // nothing more to do.
558 : : //
559 : : // 2. quo was too large, we subtracted too many den from num, and the
560 : : // most-sig-digit of num is less than needed (borrow + *num_dig < 0).
561 : : // In this case we must reduce quo and add back den to num until the
562 : : // carry from this operation cancels out the borrow.
563 : : //
564 : 1236792 : borrow += *num_dig;
565 [ + + ]: 1318470 : for (; borrow != 0; --quo) {
566 : : d = den_dig;
567 : : d_norm = 0;
568 : : mpz_dbl_dig_t carry = 0;
569 [ + + ]: 5709942 : for (mpz_dig_t *n = num_dig - den_len; n < num_dig; ++n, ++d) {
570 : 5628264 : d_norm = ((mpz_dbl_dig_t)*d << norm_shift) | (d_norm >> DIG_SIZE);
571 : 5628264 : carry += (mpz_dbl_dig_t)*n + (d_norm & DIG_MASK);
572 : 5628264 : *n = carry & DIG_MASK;
573 : 5628264 : carry >>= DIG_SIZE;
574 : : }
575 : 81678 : borrow += carry;
576 : : }
577 : :
578 : : // store this digit of the quotient
579 : 1236792 : *quo_dig = quo & DIG_MASK;
580 : 1236792 : --quo_dig;
581 : :
582 : : // move down to next digit of numerator
583 : 1236792 : --num_dig;
584 : 1236792 : --(*num_len);
585 : : }
586 : :
587 : : // unnormalise numerator (remainder now)
588 [ + + ]: 1200084 : for (mpz_dig_t *num = orig_num_dig + *num_len - 1, carry = 0; num >= orig_num_dig; --num) {
589 : 1173950 : mpz_dig_t n = *num;
590 : 1173950 : *num = ((n >> norm_shift) | carry) & DIG_MASK;
591 : 1173950 : carry = (mpz_dbl_dig_t)n << (DIG_SIZE - norm_shift);
592 : : }
593 : :
594 : : // strip trailing zeros
595 : :
596 [ + - + + ]: 46472 : while (*quo_len > 0 && orig_quo_dig[*quo_len - 1] == 0) {
597 : 20338 : --(*quo_len);
598 : : }
599 : :
600 [ + + + + ]: 28061 : while (*num_len > 0 && orig_num_dig[*num_len - 1] == 0) {
601 : 1927 : --(*num_len);
602 : : }
603 : : }
604 : :
605 : : #define MIN_ALLOC (2)
606 : :
607 : 35188 : void mpz_init_zero(mpz_t *z) {
608 : 35188 : z->neg = 0;
609 : 35188 : z->fixed_dig = 0;
610 : 35188 : z->alloc = 0;
611 : 35188 : z->len = 0;
612 : 35188 : z->dig = NULL;
613 : 35188 : }
614 : :
615 : 242 : void mpz_init_from_int(mpz_t *z, mp_int_t val) {
616 : 242 : mpz_init_zero(z);
617 : 242 : mpz_set_from_int(z, val);
618 : 242 : }
619 : :
620 : 25444 : void mpz_init_fixed_from_int(mpz_t *z, mpz_dig_t *dig, size_t alloc, mp_int_t val) {
621 : 25444 : z->neg = 0;
622 : 25444 : z->fixed_dig = 1;
623 : 25444 : z->alloc = alloc;
624 : 25444 : z->len = 0;
625 : 25444 : z->dig = dig;
626 : 25444 : mpz_set_from_int(z, val);
627 : 25444 : }
628 : :
629 : 9566 : void mpz_deinit(mpz_t *z) {
630 [ + - + - ]: 9566 : if (z != NULL && !z->fixed_dig) {
631 : 9566 : m_del(mpz_dig_t, z->dig, z->alloc);
632 : : }
633 : 9566 : }
634 : :
635 : : #if 0
636 : : these functions are unused
637 : :
638 : : mpz_t *mpz_zero(void) {
639 : : mpz_t *z = m_new_obj(mpz_t);
640 : : mpz_init_zero(z);
641 : : return z;
642 : : }
643 : :
644 : : mpz_t *mpz_from_int(mp_int_t val) {
645 : : mpz_t *z = mpz_zero();
646 : : mpz_set_from_int(z, val);
647 : : return z;
648 : : }
649 : :
650 : : mpz_t *mpz_from_ll(long long val, bool is_signed) {
651 : : mpz_t *z = mpz_zero();
652 : : mpz_set_from_ll(z, val, is_signed);
653 : : return z;
654 : : }
655 : :
656 : : #if MICROPY_PY_BUILTINS_FLOAT
657 : : mpz_t *mpz_from_float(mp_float_t val) {
658 : : mpz_t *z = mpz_zero();
659 : : mpz_set_from_float(z, val);
660 : : return z;
661 : : }
662 : : #endif
663 : :
664 : : mpz_t *mpz_from_str(const char *str, size_t len, bool neg, unsigned int base) {
665 : : mpz_t *z = mpz_zero();
666 : : mpz_set_from_str(z, str, len, neg, base);
667 : : return z;
668 : : }
669 : : #endif
670 : :
671 : 26362 : static void mpz_free(mpz_t *z) {
672 [ + + ]: 26362 : if (z != NULL) {
673 : 21068 : m_del(mpz_dig_t, z->dig, z->alloc);
674 : 21068 : m_del_obj(mpz_t, z);
675 : : }
676 : 26362 : }
677 : :
678 : 152183 : static void mpz_need_dig(mpz_t *z, size_t need) {
679 [ + + ]: 152183 : if (need < MIN_ALLOC) {
680 : 11660 : need = MIN_ALLOC;
681 : : }
682 : :
683 [ + + + + ]: 152183 : if (z->dig == NULL || z->alloc < need) {
684 : : // if z has fixed digit buffer there's not much we can do as the caller will
685 : : // be expecting a buffer with at least "need" bytes (but it shouldn't happen)
686 [ - + ]: 36222 : assert(!z->fixed_dig);
687 : 36222 : z->dig = m_renew(mpz_dig_t, z->dig, z->alloc, need);
688 : 36222 : z->alloc = need;
689 : : }
690 : 152183 : }
691 : :
692 : 21068 : static mpz_t *mpz_clone(const mpz_t *src) {
693 [ - + ]: 21068 : assert(src->alloc != 0);
694 : 21068 : mpz_t *z = m_new_obj(mpz_t);
695 : 21068 : z->neg = src->neg;
696 : 21068 : z->fixed_dig = 0;
697 : 21068 : z->alloc = src->alloc;
698 : 21068 : z->len = src->len;
699 : 21068 : z->dig = m_new(mpz_dig_t, z->alloc);
700 : 21068 : memcpy(z->dig, src->dig, src->alloc * sizeof(mpz_dig_t));
701 : 21068 : return z;
702 : : }
703 : :
704 : : /* sets dest = src
705 : : can have dest, src the same
706 : : */
707 : 30964 : void mpz_set(mpz_t *dest, const mpz_t *src) {
708 : 30964 : mpz_need_dig(dest, src->len);
709 : 30964 : dest->neg = src->neg;
710 : 30964 : dest->len = src->len;
711 : 30964 : memcpy(dest->dig, src->dig, src->len * sizeof(mpz_dig_t));
712 : 30964 : }
713 : :
714 : 30014 : void mpz_set_from_int(mpz_t *z, mp_int_t val) {
715 [ + + ]: 30014 : if (val == 0) {
716 : 4871 : z->neg = 0;
717 : 4871 : z->len = 0;
718 : 4871 : return;
719 : : }
720 : :
721 : 25143 : mpz_need_dig(z, MPZ_NUM_DIG_FOR_INT);
722 : :
723 : 25143 : mp_uint_t uval;
724 [ + + ]: 25143 : if (val < 0) {
725 : 932 : z->neg = 1;
726 : 932 : uval = -val;
727 : : } else {
728 : 24211 : z->neg = 0;
729 : 24211 : uval = val;
730 : : }
731 : :
732 : 25143 : z->len = 0;
733 [ + + ]: 50662 : while (uval > 0) {
734 : 25519 : z->dig[z->len++] = uval & DIG_MASK;
735 : 25519 : uval >>= DIG_SIZE;
736 : : }
737 : : }
738 : :
739 : 1617 : void mpz_set_from_ll(mpz_t *z, long long val, bool is_signed) {
740 : 1617 : mpz_need_dig(z, MPZ_NUM_DIG_FOR_LL);
741 : :
742 : 1617 : unsigned long long uval;
743 [ + + ]: 1617 : if (is_signed && val < 0) {
744 : 66 : z->neg = 1;
745 : 66 : uval = -(unsigned long long)val;
746 : : } else {
747 : 1551 : z->neg = 0;
748 : 1551 : uval = val;
749 : : }
750 : :
751 : 1617 : z->len = 0;
752 [ + + ]: 3973 : while (uval > 0) {
753 : 2356 : z->dig[z->len++] = uval & DIG_MASK;
754 : 2356 : uval >>= DIG_SIZE;
755 : : }
756 : 1617 : }
757 : :
758 : : #if MICROPY_PY_BUILTINS_FLOAT
759 : 4466 : void mpz_set_from_float(mpz_t *z, mp_float_t src) {
760 : 4466 : mp_float_union_t u = {src};
761 : 4466 : z->neg = u.p.sgn;
762 [ + + ]: 4466 : if (u.p.exp == 0) {
763 : : // value == 0 || value < 1
764 : 2 : mpz_set_from_int(z, 0);
765 [ + + ]: 4464 : } else if (u.p.exp == ((1 << MP_FLOAT_EXP_BITS) - 1)) {
766 : : // u.p.frc == 0 indicates inf, else NaN
767 : : // should be handled by caller
768 : 2 : mpz_set_from_int(z, 0);
769 : : } else {
770 : 4462 : const int adj_exp = (int)u.p.exp - MP_FLOAT_EXP_BIAS;
771 [ + + ]: 4462 : if (adj_exp < 0) {
772 : : // value < 1 , truncates to 0
773 : 2 : mpz_set_from_int(z, 0);
774 [ + + ]: 4460 : } else if (adj_exp == 0) {
775 : : // 1 <= value < 2 , so truncates to 1
776 : 2 : mpz_set_from_int(z, 1);
777 : : } else {
778 : : // 2 <= value
779 : 4458 : const int dig_cnt = (adj_exp + 1 + (DIG_SIZE - 1)) / DIG_SIZE;
780 : 4458 : const unsigned int rem = adj_exp % DIG_SIZE;
781 : 4458 : int dig_ind, shft;
782 : 4458 : mp_float_uint_t frc = u.p.frc | ((mp_float_uint_t)1 << MP_FLOAT_FRAC_BITS);
783 : :
784 [ + + ]: 4458 : if (adj_exp < MP_FLOAT_FRAC_BITS) {
785 : 2 : shft = 0;
786 : 2 : dig_ind = 0;
787 : 2 : frc >>= MP_FLOAT_FRAC_BITS - adj_exp;
788 : : } else {
789 : 4456 : shft = (rem - MP_FLOAT_FRAC_BITS) % DIG_SIZE;
790 : 4456 : dig_ind = (adj_exp - MP_FLOAT_FRAC_BITS) / DIG_SIZE;
791 : : }
792 : 4458 : mpz_need_dig(z, dig_cnt);
793 : 4458 : z->len = dig_cnt;
794 [ + + ]: 4458 : if (dig_ind != 0) {
795 : 4376 : memset(z->dig, 0, dig_ind * sizeof(mpz_dig_t));
796 : : }
797 [ + + ]: 4458 : if (shft != 0) {
798 : 4168 : z->dig[dig_ind++] = (frc << shft) & DIG_MASK;
799 : 4168 : frc >>= DIG_SIZE - shft;
800 : : }
801 : : #if DIG_SIZE < (MP_FLOAT_FRAC_BITS + 1)
802 [ + + ]: 19240 : while (dig_ind != dig_cnt) {
803 : 14782 : z->dig[dig_ind++] = frc & DIG_MASK;
804 : 14782 : frc >>= DIG_SIZE;
805 : : }
806 : : #else
807 : : if (dig_ind != dig_cnt) {
808 : : z->dig[dig_ind] = frc;
809 : : }
810 : : #endif
811 : : }
812 : : }
813 : 4466 : }
814 : : #endif
815 : :
816 : : // returns number of bytes from str that were processed
817 : 957 : size_t mpz_set_from_str(mpz_t *z, const char *str, size_t len, bool neg, unsigned int base) {
818 [ - + ]: 957 : assert(base <= 36);
819 : :
820 : 957 : const char *cur = str;
821 : 957 : const char *top = str + len;
822 : :
823 : 957 : mpz_need_dig(z, len * 8 / DIG_SIZE + 1);
824 : :
825 [ + + ]: 957 : if (neg) {
826 : 154 : z->neg = 1;
827 : : } else {
828 : 803 : z->neg = 0;
829 : : }
830 : :
831 : 957 : z->len = 0;
832 [ + + ]: 48703 : for (; cur < top; ++cur) { // XXX UTF8 next char
833 : : // mp_uint_t v = char_to_numeric(cur#); // XXX UTF8 get char
834 : 47754 : mp_uint_t v = *cur;
835 [ + + ]: 47754 : if ('0' <= v && v <= '9') {
836 : : v -= '0';
837 [ + + ]: 1594 : } else if ('A' <= v && v <= 'Z') {
838 : 352 : v -= 'A' - 10;
839 [ + + ]: 1242 : } else if ('a' <= v && v <= 'z') {
840 : 1238 : v -= 'a' - 10;
841 : : } else {
842 : : break;
843 : : }
844 [ + + ]: 47750 : if (v >= base) {
845 : : break;
846 : : }
847 : 47746 : z->len = mpn_mul_dig_add_dig(z->dig, z->len, base, v);
848 : : }
849 : :
850 : 957 : return cur - str;
851 : : }
852 : :
853 : 24 : void mpz_set_from_bytes(mpz_t *z, bool big_endian, size_t len, const byte *buf) {
854 : 24 : int delta = 1;
855 [ + + ]: 24 : if (big_endian) {
856 : 8 : buf += len - 1;
857 : 8 : delta = -1;
858 : : }
859 : :
860 : 24 : mpz_need_dig(z, (len * 8 + DIG_SIZE - 1) / DIG_SIZE);
861 : :
862 : 24 : mpz_dig_t d = 0;
863 : 24 : int num_bits = 0;
864 : 24 : z->neg = 0;
865 : 24 : z->len = 0;
866 [ + + ]: 256 : while (len) {
867 [ + + ]: 696 : while (len && num_bits < DIG_SIZE) {
868 : 464 : d |= *buf << num_bits;
869 : 464 : num_bits += 8;
870 : 464 : buf += delta;
871 : 464 : len--;
872 : : }
873 : 232 : z->dig[z->len++] = d & DIG_MASK;
874 : : // Need this #if because it's C undefined behavior to do: uint32_t >> 32
875 : : #if DIG_SIZE != 8 && DIG_SIZE != 16 && DIG_SIZE != 32
876 : : d >>= DIG_SIZE;
877 : : #else
878 : 232 : d = 0;
879 : : #endif
880 : 232 : num_bits -= DIG_SIZE;
881 : : }
882 : :
883 : 24 : z->len = mpn_remove_trailing_zeros(z->dig, z->dig + z->len);
884 : 24 : }
885 : :
886 : : #if 0
887 : : these functions are unused
888 : :
889 : : bool mpz_is_pos(const mpz_t *z) {
890 : : return z->len > 0 && z->neg == 0;
891 : : }
892 : :
893 : : bool mpz_is_odd(const mpz_t *z) {
894 : : return z->len > 0 && (z->dig[0] & 1) != 0;
895 : : }
896 : :
897 : : bool mpz_is_even(const mpz_t *z) {
898 : : return z->len == 0 || (z->dig[0] & 1) == 0;
899 : : }
900 : : #endif
901 : :
902 : 16234 : int mpz_cmp(const mpz_t *z1, const mpz_t *z2) {
903 : 16234 : int cmp = (int)z2->neg - (int)z1->neg;
904 [ + + ]: 16234 : if (cmp != 0) {
905 : : return cmp;
906 : : }
907 : 13346 : cmp = mpn_cmp(z1->dig, z1->len, z2->dig, z2->len);
908 [ + + ]: 13346 : if (z1->neg != 0) {
909 : 2968 : cmp = -cmp;
910 : : }
911 : : return cmp;
912 : : }
913 : :
914 : : #if 0
915 : : // obsolete
916 : : // compares mpz with an integer that fits within DIG_SIZE bits
917 : : mp_int_t mpz_cmp_sml_int(const mpz_t *z, mp_int_t sml_int) {
918 : : mp_int_t cmp;
919 : : if (z->neg == 0) {
920 : : if (sml_int < 0) {
921 : : return 1;
922 : : }
923 : : if (sml_int == 0) {
924 : : if (z->len == 0) {
925 : : return 0;
926 : : }
927 : : return 1;
928 : : }
929 : : if (z->len == 0) {
930 : : return -1;
931 : : }
932 : : assert(sml_int < (1 << DIG_SIZE));
933 : : if (z->len != 1) {
934 : : return 1;
935 : : }
936 : : cmp = z->dig[0] - sml_int;
937 : : } else {
938 : : if (sml_int > 0) {
939 : : return -1;
940 : : }
941 : : if (sml_int == 0) {
942 : : if (z->len == 0) {
943 : : return 0;
944 : : }
945 : : return -1;
946 : : }
947 : : if (z->len == 0) {
948 : : return 1;
949 : : }
950 : : assert(sml_int > -(1 << DIG_SIZE));
951 : : if (z->len != 1) {
952 : : return -1;
953 : : }
954 : : cmp = -z->dig[0] - sml_int;
955 : : }
956 : : if (cmp < 0) {
957 : : return -1;
958 : : }
959 : : if (cmp > 0) {
960 : : return 1;
961 : : }
962 : : return 0;
963 : : }
964 : : #endif
965 : :
966 : : #if 0
967 : : these functions are unused
968 : :
969 : : /* returns abs(z)
970 : : */
971 : : mpz_t *mpz_abs(const mpz_t *z) {
972 : : // TODO: handle case of z->alloc=0
973 : : mpz_t *z2 = mpz_clone(z);
974 : : z2->neg = 0;
975 : : return z2;
976 : : }
977 : :
978 : : /* returns -z
979 : : */
980 : : mpz_t *mpz_neg(const mpz_t *z) {
981 : : // TODO: handle case of z->alloc=0
982 : : mpz_t *z2 = mpz_clone(z);
983 : : z2->neg = 1 - z2->neg;
984 : : return z2;
985 : : }
986 : :
987 : : /* returns lhs + rhs
988 : : can have lhs, rhs the same
989 : : */
990 : : mpz_t *mpz_add(const mpz_t *lhs, const mpz_t *rhs) {
991 : : mpz_t *z = mpz_zero();
992 : : mpz_add_inpl(z, lhs, rhs);
993 : : return z;
994 : : }
995 : :
996 : : /* returns lhs - rhs
997 : : can have lhs, rhs the same
998 : : */
999 : : mpz_t *mpz_sub(const mpz_t *lhs, const mpz_t *rhs) {
1000 : : mpz_t *z = mpz_zero();
1001 : : mpz_sub_inpl(z, lhs, rhs);
1002 : : return z;
1003 : : }
1004 : :
1005 : : /* returns lhs * rhs
1006 : : can have lhs, rhs the same
1007 : : */
1008 : : mpz_t *mpz_mul(const mpz_t *lhs, const mpz_t *rhs) {
1009 : : mpz_t *z = mpz_zero();
1010 : : mpz_mul_inpl(z, lhs, rhs);
1011 : : return z;
1012 : : }
1013 : :
1014 : : /* returns lhs ** rhs
1015 : : can have lhs, rhs the same
1016 : : */
1017 : : mpz_t *mpz_pow(const mpz_t *lhs, const mpz_t *rhs) {
1018 : : mpz_t *z = mpz_zero();
1019 : : mpz_pow_inpl(z, lhs, rhs);
1020 : : return z;
1021 : : }
1022 : :
1023 : : /* computes new integers in quo and rem such that:
1024 : : quo * rhs + rem = lhs
1025 : : 0 <= rem < rhs
1026 : : can have lhs, rhs the same
1027 : : */
1028 : : void mpz_divmod(const mpz_t *lhs, const mpz_t *rhs, mpz_t **quo, mpz_t **rem) {
1029 : : *quo = mpz_zero();
1030 : : *rem = mpz_zero();
1031 : : mpz_divmod_inpl(*quo, *rem, lhs, rhs);
1032 : : }
1033 : : #endif
1034 : :
1035 : : /* computes dest = abs(z)
1036 : : can have dest, z the same
1037 : : */
1038 : 6 : void mpz_abs_inpl(mpz_t *dest, const mpz_t *z) {
1039 [ + - ]: 6 : if (dest != z) {
1040 : 6 : mpz_set(dest, z);
1041 : : }
1042 : 6 : dest->neg = 0;
1043 : 6 : }
1044 : :
1045 : : /* computes dest = -z
1046 : : can have dest, z the same
1047 : : */
1048 : 534 : void mpz_neg_inpl(mpz_t *dest, const mpz_t *z) {
1049 [ + - ]: 534 : if (dest != z) {
1050 : 534 : mpz_set(dest, z);
1051 : : }
1052 [ + - ]: 534 : if (dest->len) {
1053 : 534 : dest->neg = 1 - dest->neg;
1054 : : }
1055 : 534 : }
1056 : :
1057 : : /* computes dest = ~z (= -z - 1)
1058 : : can have dest, z the same
1059 : : */
1060 : 14 : void mpz_not_inpl(mpz_t *dest, const mpz_t *z) {
1061 [ + + ]: 14 : if (dest != z) {
1062 : 12 : mpz_set(dest, z);
1063 : : }
1064 [ + + ]: 14 : if (dest->len == 0) {
1065 : 2 : mpz_need_dig(dest, 1);
1066 : 2 : dest->dig[0] = 1;
1067 : 2 : dest->len = 1;
1068 : 2 : dest->neg = 1;
1069 [ + + ]: 12 : } else if (dest->neg) {
1070 : 4 : dest->neg = 0;
1071 : 4 : mpz_dig_t k = 1;
1072 : 4 : dest->len = mpn_sub(dest->dig, dest->dig, dest->len, &k, 1);
1073 : : } else {
1074 : 8 : mpz_need_dig(dest, dest->len + 1);
1075 : 8 : mpz_dig_t k = 1;
1076 : 8 : dest->len = mpn_add(dest->dig, dest->dig, dest->len, &k, 1);
1077 : 8 : dest->neg = 1;
1078 : : }
1079 : 14 : }
1080 : :
1081 : : /* computes dest = lhs << rhs
1082 : : can have dest, lhs the same
1083 : : */
1084 : 1300 : void mpz_shl_inpl(mpz_t *dest, const mpz_t *lhs, mp_uint_t rhs) {
1085 [ + - + + ]: 1300 : if (lhs->len == 0 || rhs == 0) {
1086 : 42 : mpz_set(dest, lhs);
1087 : : } else {
1088 : 1258 : mpz_need_dig(dest, lhs->len + (rhs + DIG_SIZE - 1) / DIG_SIZE);
1089 : 1258 : dest->len = mpn_shl(dest->dig, lhs->dig, lhs->len, rhs);
1090 : 1258 : dest->neg = lhs->neg;
1091 : : }
1092 : 1300 : }
1093 : :
1094 : : /* computes dest = lhs >> rhs
1095 : : can have dest, lhs the same
1096 : : */
1097 : 476 : void mpz_shr_inpl(mpz_t *dest, const mpz_t *lhs, mp_uint_t rhs) {
1098 [ + - + + ]: 476 : if (lhs->len == 0 || rhs == 0) {
1099 : 40 : mpz_set(dest, lhs);
1100 : : } else {
1101 : 436 : mpz_need_dig(dest, lhs->len);
1102 : 436 : dest->len = mpn_shr(dest->dig, lhs->dig, lhs->len, rhs);
1103 : 436 : dest->neg = lhs->neg;
1104 [ + + ]: 436 : if (dest->neg) {
1105 : : // arithmetic shift right, rounding to negative infinity
1106 : 288 : mp_uint_t n_whole = rhs / DIG_SIZE;
1107 : 288 : mp_uint_t n_part = rhs % DIG_SIZE;
1108 : 288 : mpz_dig_t round_up = 0;
1109 [ + - + + ]: 304 : for (size_t i = 0; i < lhs->len && i < n_whole; i++) {
1110 [ + + ]: 24 : if (lhs->dig[i] != 0) {
1111 : 8 : round_up = 1;
1112 : 8 : break;
1113 : : }
1114 : : }
1115 [ + + + + ]: 288 : if (n_whole < lhs->len && (lhs->dig[n_whole] & ((1 << n_part) - 1)) != 0) {
1116 : 200 : round_up = 1;
1117 : : }
1118 [ + + ]: 288 : if (round_up) {
1119 [ + + ]: 208 : if (dest->len == 0) {
1120 : : // dest == 0, so need to add 1 by hand (answer will be -1)
1121 : 4 : dest->dig[0] = 1;
1122 : 4 : dest->len = 1;
1123 : : } else {
1124 : : // dest > 0, so can use mpn_add to add 1
1125 : 204 : dest->len = mpn_add(dest->dig, dest->dig, dest->len, &round_up, 1);
1126 : : }
1127 : : }
1128 : : }
1129 : : }
1130 : 476 : }
1131 : :
1132 : : /* computes dest = lhs + rhs
1133 : : can have dest, lhs, rhs the same
1134 : : */
1135 : 412 : void mpz_add_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1136 [ + + ]: 412 : if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1137 : 236 : const mpz_t *temp = lhs;
1138 : 236 : lhs = rhs;
1139 : 236 : rhs = temp;
1140 : : }
1141 : :
1142 [ + + ]: 412 : if (lhs->neg == rhs->neg) {
1143 : 154 : mpz_need_dig(dest, lhs->len + 1);
1144 : 154 : dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1145 : : } else {
1146 : 258 : mpz_need_dig(dest, lhs->len);
1147 : 258 : dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1148 : : }
1149 : :
1150 : 412 : dest->neg = lhs->neg & !!dest->len;
1151 : 412 : }
1152 : :
1153 : : /* computes dest = lhs - rhs
1154 : : can have dest, lhs, rhs the same
1155 : : */
1156 : 146 : void mpz_sub_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1157 : 146 : bool neg = false;
1158 : :
1159 [ + + ]: 146 : if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1160 : 16 : const mpz_t *temp = lhs;
1161 : 16 : lhs = rhs;
1162 : 16 : rhs = temp;
1163 : 16 : neg = true;
1164 : : }
1165 : :
1166 [ + + ]: 146 : if (lhs->neg != rhs->neg) {
1167 : 12 : mpz_need_dig(dest, lhs->len + 1);
1168 : 12 : dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1169 : : } else {
1170 : 134 : mpz_need_dig(dest, lhs->len);
1171 : 134 : dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1172 : : }
1173 : :
1174 [ + + ]: 146 : if (dest->len == 0) {
1175 : 40 : dest->neg = 0;
1176 [ + + ]: 106 : } else if (neg) {
1177 : 16 : dest->neg = 1 - lhs->neg;
1178 : : } else {
1179 : 90 : dest->neg = lhs->neg;
1180 : : }
1181 : 146 : }
1182 : :
1183 : : /* computes dest = lhs & rhs
1184 : : can have dest, lhs, rhs the same
1185 : : */
1186 : 232 : void mpz_and_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1187 : : // make sure lhs has the most digits
1188 [ + + ]: 232 : if (lhs->len < rhs->len) {
1189 : 48 : const mpz_t *temp = lhs;
1190 : 48 : lhs = rhs;
1191 : 48 : rhs = temp;
1192 : : }
1193 : :
1194 : : #if MICROPY_OPT_MPZ_BITWISE
1195 : :
1196 [ + + + + ]: 232 : if ((0 == lhs->neg) && (0 == rhs->neg)) {
1197 : 92 : mpz_need_dig(dest, lhs->len);
1198 : 92 : dest->len = mpn_and(dest->dig, lhs->dig, rhs->dig, rhs->len);
1199 : 92 : dest->neg = 0;
1200 : : } else {
1201 : 140 : mpz_need_dig(dest, lhs->len + 1);
1202 : 280 : dest->len = mpn_and_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len,
1203 : 140 : lhs->neg == rhs->neg, 0 != lhs->neg, 0 != rhs->neg);
1204 : 140 : dest->neg = lhs->neg & rhs->neg;
1205 : : }
1206 : :
1207 : : #else
1208 : :
1209 : : mpz_need_dig(dest, lhs->len + (lhs->neg || rhs->neg));
1210 : : dest->len = mpn_and_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len,
1211 : : (lhs->neg == rhs->neg) ? lhs->neg : 0, lhs->neg, rhs->neg);
1212 : : dest->neg = lhs->neg & rhs->neg;
1213 : :
1214 : : #endif
1215 : 232 : }
1216 : :
1217 : : /* computes dest = lhs | rhs
1218 : : can have dest, lhs, rhs the same
1219 : : */
1220 : 160 : void mpz_or_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1221 : : // make sure lhs has the most digits
1222 [ + + ]: 160 : if (lhs->len < rhs->len) {
1223 : 34 : const mpz_t *temp = lhs;
1224 : 34 : lhs = rhs;
1225 : 34 : rhs = temp;
1226 : : }
1227 : :
1228 : : #if MICROPY_OPT_MPZ_BITWISE
1229 : :
1230 [ + + + + ]: 160 : if ((0 == lhs->neg) && (0 == rhs->neg)) {
1231 : 54 : mpz_need_dig(dest, lhs->len);
1232 : 54 : dest->len = mpn_or(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1233 : 54 : dest->neg = 0;
1234 : : } else {
1235 : 106 : mpz_need_dig(dest, lhs->len + 1);
1236 : 212 : dest->len = mpn_or_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len,
1237 : 106 : 0 != lhs->neg, 0 != rhs->neg);
1238 : 106 : dest->neg = 1;
1239 : : }
1240 : :
1241 : : #else
1242 : :
1243 : : mpz_need_dig(dest, lhs->len + (lhs->neg || rhs->neg));
1244 : : dest->len = mpn_or_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len,
1245 : : (lhs->neg || rhs->neg), lhs->neg, rhs->neg);
1246 : : dest->neg = lhs->neg | rhs->neg;
1247 : :
1248 : : #endif
1249 : 160 : }
1250 : :
1251 : : /* computes dest = lhs ^ rhs
1252 : : can have dest, lhs, rhs the same
1253 : : */
1254 : 156 : void mpz_xor_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1255 : : // make sure lhs has the most digits
1256 [ + + ]: 156 : if (lhs->len < rhs->len) {
1257 : 32 : const mpz_t *temp = lhs;
1258 : 32 : lhs = rhs;
1259 : 32 : rhs = temp;
1260 : : }
1261 : :
1262 : : #if MICROPY_OPT_MPZ_BITWISE
1263 : :
1264 [ + + ]: 156 : if (lhs->neg == rhs->neg) {
1265 : 80 : mpz_need_dig(dest, lhs->len);
1266 [ + + ]: 80 : if (lhs->neg == 0) {
1267 : 42 : dest->len = mpn_xor(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1268 : : } else {
1269 : 38 : dest->len = mpn_xor_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len, 0, 0, 0);
1270 : : }
1271 : 80 : dest->neg = 0;
1272 : : } else {
1273 : 76 : mpz_need_dig(dest, lhs->len + 1);
1274 : 152 : dest->len = mpn_xor_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len, 1,
1275 : 76 : 0 == lhs->neg, 0 == rhs->neg);
1276 : 76 : dest->neg = 1;
1277 : : }
1278 : :
1279 : : #else
1280 : :
1281 : : mpz_need_dig(dest, lhs->len + (lhs->neg || rhs->neg));
1282 : : dest->len = mpn_xor_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len,
1283 : : (lhs->neg != rhs->neg), 0 == lhs->neg, 0 == rhs->neg);
1284 : : dest->neg = lhs->neg ^ rhs->neg;
1285 : :
1286 : : #endif
1287 : 156 : }
1288 : :
1289 : : /* computes dest = lhs * rhs
1290 : : can have dest, lhs, rhs the same
1291 : : */
1292 : 29430 : void mpz_mul_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1293 [ + + + + ]: 29430 : if (lhs->len == 0 || rhs->len == 0) {
1294 : 3880 : mpz_set_from_int(dest, 0);
1295 : 3880 : return;
1296 : : }
1297 : :
1298 : 25550 : mpz_t *temp = NULL;
1299 [ + + ]: 25550 : if (lhs == dest) {
1300 : 20254 : lhs = temp = mpz_clone(lhs);
1301 [ + + ]: 20254 : if (rhs == dest) {
1302 : 13448 : rhs = lhs;
1303 : : }
1304 [ + + ]: 5296 : } else if (rhs == dest) {
1305 : 2 : rhs = temp = mpz_clone(rhs);
1306 : : }
1307 : :
1308 : 25550 : mpz_need_dig(dest, lhs->len + rhs->len); // min mem l+r-1, max mem l+r
1309 : 25550 : memset(dest->dig, 0, dest->alloc * sizeof(mpz_dig_t));
1310 : 25550 : dest->len = mpn_mul(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1311 : :
1312 [ + + ]: 25550 : if (lhs->neg == rhs->neg) {
1313 : 25094 : dest->neg = 0;
1314 : : } else {
1315 : 456 : dest->neg = 1;
1316 : : }
1317 : :
1318 : 25550 : mpz_free(temp);
1319 : : }
1320 : :
1321 : : /* computes dest = lhs ** rhs
1322 : : can have dest, lhs, rhs the same
1323 : : */
1324 : 382 : void mpz_pow_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1325 [ + + - + ]: 382 : if (lhs->len == 0 || rhs->neg != 0) {
1326 : 4 : mpz_set_from_int(dest, 0);
1327 : 4 : return;
1328 : : }
1329 : :
1330 [ + + ]: 378 : if (rhs->len == 0) {
1331 : 4 : mpz_set_from_int(dest, 1);
1332 : 4 : return;
1333 : : }
1334 : :
1335 : 374 : mpz_t *x = mpz_clone(lhs);
1336 : 374 : mpz_t *n = mpz_clone(rhs);
1337 : :
1338 : 374 : mpz_set_from_int(dest, 1);
1339 : :
1340 [ + - ]: 2326 : while (n->len > 0) {
1341 [ + + ]: 2326 : if ((n->dig[0] & 1) != 0) {
1342 : 1214 : mpz_mul_inpl(dest, dest, x);
1343 : : }
1344 : 2326 : n->len = mpn_shr(n->dig, n->dig, n->len, 1);
1345 [ + + ]: 2326 : if (n->len == 0) {
1346 : : break;
1347 : : }
1348 : 1952 : mpz_mul_inpl(x, x, x);
1349 : : }
1350 : :
1351 : 374 : mpz_free(x);
1352 : 374 : mpz_free(n);
1353 : : }
1354 : :
1355 : : /* computes dest = (lhs ** rhs) % mod
1356 : : can have dest, lhs, rhs the same; mod can't be the same as dest
1357 : : */
1358 : 48 : void mpz_pow3_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs, const mpz_t *mod) {
1359 [ + + + - : 48 : if (lhs->len == 0 || rhs->neg != 0 || (mod->len == 1 && mod->dig[0] == 1)) {
+ + + + ]
1360 : 12 : mpz_set_from_int(dest, 0);
1361 : 28 : return;
1362 : : }
1363 : :
1364 : 36 : mpz_set_from_int(dest, 1);
1365 : :
1366 [ + + ]: 36 : if (rhs->len == 0) {
1367 : : return;
1368 : : }
1369 : :
1370 : 32 : mpz_t *x = mpz_clone(lhs);
1371 : 32 : mpz_t *n = mpz_clone(rhs);
1372 : 32 : mpz_t quo;
1373 : 32 : mpz_init_zero(&quo);
1374 : :
1375 [ + - ]: 14104 : while (n->len > 0) {
1376 [ + + ]: 14104 : if ((n->dig[0] & 1) != 0) {
1377 : 6804 : mpz_mul_inpl(dest, dest, x);
1378 : 6804 : mpz_divmod_inpl(&quo, dest, dest, mod);
1379 : : }
1380 : 14104 : n->len = mpn_shr(n->dig, n->dig, n->len, 1);
1381 [ + + ]: 14104 : if (n->len == 0) {
1382 : : break;
1383 : : }
1384 : 14072 : mpz_mul_inpl(x, x, x);
1385 : 14072 : mpz_divmod_inpl(&quo, x, x, mod);
1386 : : }
1387 : :
1388 : 32 : mpz_deinit(&quo);
1389 : 32 : mpz_free(x);
1390 : 32 : mpz_free(n);
1391 : : }
1392 : :
1393 : : #if 0
1394 : : these functions are unused
1395 : :
1396 : : /* computes gcd(z1, z2)
1397 : : based on Knuth's modified gcd algorithm (I think?)
1398 : : gcd(z1, z2) >= 0
1399 : : gcd(0, 0) = 0
1400 : : gcd(z, 0) = abs(z)
1401 : : */
1402 : : mpz_t *mpz_gcd(const mpz_t *z1, const mpz_t *z2) {
1403 : : if (z1->len == 0) {
1404 : : // TODO: handle case of z2->alloc=0
1405 : : mpz_t *a = mpz_clone(z2);
1406 : : a->neg = 0;
1407 : : return a;
1408 : : } else if (z2->len == 0) {
1409 : : mpz_t *a = mpz_clone(z1);
1410 : : a->neg = 0;
1411 : : return a;
1412 : : }
1413 : :
1414 : : mpz_t *a = mpz_clone(z1);
1415 : : mpz_t *b = mpz_clone(z2);
1416 : : mpz_t c;
1417 : : mpz_init_zero(&c);
1418 : : a->neg = 0;
1419 : : b->neg = 0;
1420 : :
1421 : : for (;;) {
1422 : : if (mpz_cmp(a, b) < 0) {
1423 : : if (a->len == 0) {
1424 : : mpz_free(a);
1425 : : mpz_deinit(&c);
1426 : : return b;
1427 : : }
1428 : : mpz_t *t = a;
1429 : : a = b;
1430 : : b = t;
1431 : : }
1432 : : if (!(b->len >= 2 || (b->len == 1 && b->dig[0] > 1))) { // compute b > 0; could be mpz_cmp_small_int(b, 1) > 0
1433 : : break;
1434 : : }
1435 : : mpz_set(&c, b);
1436 : : do {
1437 : : mpz_add_inpl(&c, &c, &c);
1438 : : } while (mpz_cmp(&c, a) <= 0);
1439 : : c.len = mpn_shr(c.dig, c.dig, c.len, 1);
1440 : : mpz_sub_inpl(a, a, &c);
1441 : : }
1442 : :
1443 : : mpz_deinit(&c);
1444 : :
1445 : : if (b->len == 1 && b->dig[0] == 1) { // compute b == 1; could be mpz_cmp_small_int(b, 1) == 0
1446 : : mpz_free(a);
1447 : : return b;
1448 : : } else {
1449 : : mpz_free(b);
1450 : : return a;
1451 : : }
1452 : : }
1453 : :
1454 : : /* computes lcm(z1, z2)
1455 : : = abs(z1) / gcd(z1, z2) * abs(z2)
1456 : : lcm(z1, z1) >= 0
1457 : : lcm(0, 0) = 0
1458 : : lcm(z, 0) = 0
1459 : : */
1460 : : mpz_t *mpz_lcm(const mpz_t *z1, const mpz_t *z2) {
1461 : : if (z1->len == 0 || z2->len == 0) {
1462 : : return mpz_zero();
1463 : : }
1464 : :
1465 : : mpz_t *gcd = mpz_gcd(z1, z2);
1466 : : mpz_t *quo = mpz_zero();
1467 : : mpz_t *rem = mpz_zero();
1468 : : mpz_divmod_inpl(quo, rem, z1, gcd);
1469 : : mpz_mul_inpl(rem, quo, z2);
1470 : : mpz_free(gcd);
1471 : : mpz_free(quo);
1472 : : rem->neg = 0;
1473 : : return rem;
1474 : : }
1475 : : #endif
1476 : :
1477 : : /* computes new integers in quo and rem such that:
1478 : : quo * rhs + rem = lhs
1479 : : 0 <= rem < rhs
1480 : : can have lhs, rhs the same
1481 : : assumes rhs != 0 (undefined behaviour if it is)
1482 : : */
1483 : 30330 : void mpz_divmod_inpl(mpz_t *dest_quo, mpz_t *dest_rem, const mpz_t *lhs, const mpz_t *rhs) {
1484 [ - + ]: 30330 : assert(!mpz_is_zero(rhs));
1485 : :
1486 : 30330 : mpz_need_dig(dest_quo, lhs->len + 1); // +1 necessary?
1487 : 30330 : memset(dest_quo->dig, 0, (lhs->len + 1) * sizeof(mpz_dig_t));
1488 : 30330 : dest_quo->neg = 0;
1489 : 30330 : dest_quo->len = 0;
1490 : 30330 : mpz_need_dig(dest_rem, lhs->len + 1); // +1 necessary?
1491 : 30330 : mpz_set(dest_rem, lhs);
1492 : 30330 : mpn_div(dest_rem->dig, &dest_rem->len, rhs->dig, rhs->len, dest_quo->dig, &dest_quo->len);
1493 : 30330 : dest_rem->neg &= !!dest_rem->len;
1494 : :
1495 : : // check signs and do Python style modulo
1496 [ + + ]: 30330 : if (lhs->neg != rhs->neg) {
1497 : 264 : dest_quo->neg = !!dest_quo->len;
1498 [ + + ]: 264 : if (!mpz_is_zero(dest_rem)) {
1499 : 144 : mpz_t mpzone;
1500 : 144 : mpz_init_from_int(&mpzone, -1);
1501 : 144 : mpz_add_inpl(dest_quo, dest_quo, &mpzone);
1502 : 144 : mpz_add_inpl(dest_rem, dest_rem, rhs);
1503 : : }
1504 : : }
1505 : 30330 : }
1506 : :
1507 : : #if 0
1508 : : these functions are unused
1509 : :
1510 : : /* computes floor(lhs / rhs)
1511 : : can have lhs, rhs the same
1512 : : */
1513 : : mpz_t *mpz_div(const mpz_t *lhs, const mpz_t *rhs) {
1514 : : mpz_t *quo = mpz_zero();
1515 : : mpz_t rem;
1516 : : mpz_init_zero(&rem);
1517 : : mpz_divmod_inpl(quo, &rem, lhs, rhs);
1518 : : mpz_deinit(&rem);
1519 : : return quo;
1520 : : }
1521 : :
1522 : : /* computes lhs % rhs ( >= 0)
1523 : : can have lhs, rhs the same
1524 : : */
1525 : : mpz_t *mpz_mod(const mpz_t *lhs, const mpz_t *rhs) {
1526 : : mpz_t quo;
1527 : : mpz_init_zero(&quo);
1528 : : mpz_t *rem = mpz_zero();
1529 : : mpz_divmod_inpl(&quo, rem, lhs, rhs);
1530 : : mpz_deinit(&quo);
1531 : : return rem;
1532 : : }
1533 : : #endif
1534 : :
1535 : : // must return actual int value if it fits in mp_int_t
1536 : 16 : mp_int_t mpz_hash(const mpz_t *z) {
1537 : 16 : mp_uint_t val = 0;
1538 : 16 : mpz_dig_t *d = z->dig + z->len;
1539 : :
1540 [ + + ]: 136 : while (d-- > z->dig) {
1541 : 120 : val = (val << DIG_SIZE) | *d;
1542 : : }
1543 : :
1544 [ + + ]: 16 : if (z->neg != 0) {
1545 : 4 : val = -val;
1546 : : }
1547 : :
1548 : 16 : return val;
1549 : : }
1550 : :
1551 : 17806 : bool mpz_as_int_checked(const mpz_t *i, mp_int_t *value) {
1552 : 17806 : mp_uint_t val = 0;
1553 : 17806 : mpz_dig_t *d = i->dig + i->len;
1554 : :
1555 [ + + ]: 86010 : while (d-- > i->dig) {
1556 [ + + ]: 85198 : if (val > (~(MP_OBJ_WORD_MSBIT_HIGH) >> DIG_SIZE)) {
1557 : : // will overflow
1558 : : return false;
1559 : : }
1560 : 68204 : val = (val << DIG_SIZE) | *d;
1561 : : }
1562 : :
1563 [ + + ]: 812 : if (i->neg != 0) {
1564 : 74 : val = -val;
1565 : : }
1566 : :
1567 : 812 : *value = val;
1568 : 812 : return true;
1569 : : }
1570 : :
1571 : 22 : bool mpz_as_uint_checked(const mpz_t *i, mp_uint_t *value) {
1572 [ + + ]: 22 : if (i->neg != 0) {
1573 : : // can't represent signed values
1574 : : return false;
1575 : : }
1576 : :
1577 : 18 : mp_uint_t val = 0;
1578 : 18 : mpz_dig_t *d = i->dig + i->len;
1579 : :
1580 [ + + ]: 38 : while (d-- > i->dig) {
1581 [ + + ]: 22 : if (val > (~(MP_OBJ_WORD_MSBIT_HIGH) >> (DIG_SIZE - 1))) {
1582 : : // will overflow
1583 : : return false;
1584 : : }
1585 : 20 : val = (val << DIG_SIZE) | *d;
1586 : : }
1587 : :
1588 : 16 : *value = val;
1589 : 16 : return true;
1590 : : }
1591 : :
1592 : 132 : bool mpz_as_bytes(const mpz_t *z, bool big_endian, bool as_signed, size_t len, byte *buf) {
1593 : 132 : byte *b = buf;
1594 [ + + ]: 132 : if (big_endian) {
1595 : 52 : b += len;
1596 : : }
1597 : 132 : mpz_dig_t *zdig = z->dig;
1598 : 132 : int bits = 0;
1599 : 132 : mpz_dbl_dig_t d = 0;
1600 : 132 : mpz_dbl_dig_t carry = 1;
1601 : 132 : size_t olen = len; // bytes in output buffer
1602 : 132 : bool ok = true;
1603 [ + + ]: 896 : for (size_t zlen = z->len; zlen > 0; --zlen) {
1604 : 764 : bits += DIG_SIZE;
1605 : 764 : d = (d << DIG_SIZE) | *zdig++;
1606 [ + + ]: 2292 : for (; bits >= 8; bits -= 8, d >>= 8) {
1607 : 1528 : mpz_dig_t val = d;
1608 [ + + ]: 1528 : if (z->neg) {
1609 : 456 : val = (~val & 0xff) + carry;
1610 : 456 : carry = val >> 8;
1611 : : }
1612 : :
1613 [ + + ]: 1528 : if (!olen) {
1614 : : // Buffer is full, only OK if all remaining bytes are zeroes
1615 [ + + + + ]: 132 : ok = ok && ((byte)val == 0);
1616 : 132 : continue;
1617 : : }
1618 : :
1619 [ + + ]: 1396 : if (big_endian) {
1620 : 720 : *--b = val;
1621 : : } else {
1622 : 676 : *b++ = val;
1623 : : }
1624 : 1396 : olen--;
1625 : : }
1626 : : }
1627 : :
1628 [ + + + - ]: 132 : if (as_signed && olen == 0 && len > 0) {
1629 : : // If output exhausted then ensure there was enough space for the sign bit
1630 [ + + ]: 16 : byte most_sig = big_endian ? buf[0] : buf[len - 1];
1631 [ + - + + ]: 20 : ok = ok && (bool)(most_sig & 0x80) == (bool)z->neg;
1632 : : } else {
1633 : : // fill remainder of buf with zero/sign extension of the integer
1634 [ + + + + ]: 280 : memset(big_endian ? buf : b, z->neg ? 0xff : 0x00, olen);
1635 : : }
1636 : :
1637 : 132 : return ok;
1638 : : }
1639 : :
1640 : : #if MICROPY_PY_BUILTINS_FLOAT
1641 : 21648 : mp_float_t mpz_as_float(const mpz_t *i) {
1642 : 21648 : mp_float_t val = 0;
1643 : 21648 : mpz_dig_t *d = i->dig + i->len;
1644 : :
1645 [ + + ]: 301504 : while (d-- > i->dig) {
1646 : 279856 : val = val * DIG_BASE + *d;
1647 : : }
1648 : :
1649 [ + + ]: 21648 : if (i->neg != 0) {
1650 : 20 : val = -val;
1651 : : }
1652 : :
1653 : 21648 : return val;
1654 : : }
1655 : : #endif
1656 : :
1657 : : #if 0
1658 : : this function is unused
1659 : : char *mpz_as_str(const mpz_t *i, unsigned int base) {
1660 : : char *s = m_new(char, mp_int_format_size(mpz_max_num_bits(i), base, NULL, '\0'));
1661 : : mpz_as_str_inpl(i, base, NULL, 'a', '\0', s);
1662 : : return s;
1663 : : }
1664 : : #endif
1665 : :
1666 : : // assumes enough space in str as calculated by mp_int_format_size
1667 : : // base must be between 2 and 32 inclusive
1668 : : // returns length of string, not including null byte
1669 : 22670 : size_t mpz_as_str_inpl(const mpz_t *i, unsigned int base, const char *prefix, char base_char, char comma, char *str) {
1670 [ - + ]: 22670 : assert(str != NULL);
1671 [ - + ]: 22670 : assert(2 <= base && base <= 32);
1672 : :
1673 : 22670 : size_t ilen = i->len;
1674 : :
1675 : 22670 : char *s = str;
1676 [ + + ]: 22670 : if (ilen == 0) {
1677 [ + + ]: 44 : if (prefix) {
1678 [ + + ]: 12 : while (*prefix) {
1679 : 8 : *s++ = *prefix++;
1680 : : }
1681 : : }
1682 : 44 : *s++ = '0';
1683 : 44 : *s = '\0';
1684 : 44 : return s - str;
1685 : : }
1686 : :
1687 : : // make a copy of mpz digits, so we can do the div/mod calculation
1688 : 22626 : mpz_dig_t *dig = m_new(mpz_dig_t, ilen);
1689 : 22626 : memcpy(dig, i->dig, ilen * sizeof(mpz_dig_t));
1690 : :
1691 : : // convert
1692 : 22626 : char *last_comma = str;
1693 : 2752311 : bool done;
1694 : 2752311 : do {
1695 : 2752311 : mpz_dig_t *d = dig + ilen;
1696 : 2752311 : mpz_dbl_dig_t a = 0;
1697 : :
1698 : : // compute next remainder
1699 [ + + ]: 99331959 : while (--d >= dig) {
1700 : 96579648 : a = (a << DIG_SIZE) | *d;
1701 : 96579648 : *d = a / base;
1702 : 96579648 : a %= base;
1703 : : }
1704 : :
1705 : : // convert to character
1706 : 2752311 : a += '0';
1707 [ + + ]: 2752311 : if (a > '9') {
1708 : 620 : a += base_char - '9' - 1;
1709 : : }
1710 : 2752311 : *s++ = a;
1711 : :
1712 : : // check if number is zero
1713 : 2752311 : done = true;
1714 [ + + ]: 46854463 : for (d = dig; d < dig + ilen; ++d) {
1715 [ + + ]: 46831837 : if (*d != 0) {
1716 : : done = false;
1717 : : break;
1718 : : }
1719 : : }
1720 [ + + + + ]: 2752311 : if (!done && comma && (s - last_comma) == 3) {
1721 : 232 : *s++ = comma;
1722 : 232 : last_comma = s;
1723 : : }
1724 : : }
1725 [ + + ]: 2752311 : while (!done);
1726 : :
1727 : : // free the copy of the digits array
1728 : 22626 : m_del(mpz_dig_t, dig, ilen);
1729 : :
1730 [ + + ]: 22626 : if (prefix) {
1731 : 4496 : const char *p = &prefix[strlen(prefix)];
1732 [ + + ]: 13344 : while (p > prefix) {
1733 : 8848 : *s++ = *--p;
1734 : : }
1735 : : }
1736 [ + + ]: 22626 : if (i->neg != 0) {
1737 : 1818 : *s++ = '-';
1738 : : }
1739 : :
1740 : : // reverse string
1741 [ + + ]: 1397907 : for (char *u = str, *v = s - 1; u < v; ++u, --v) {
1742 : 1375281 : char temp = *u;
1743 : 1375281 : *u = *v;
1744 : 1375281 : *v = temp;
1745 : : }
1746 : :
1747 : 22626 : *s = '\0'; // null termination
1748 : :
1749 : 22626 : return s - str;
1750 : : }
1751 : :
1752 : : #endif // MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ
|
