blob: 166fa7adbea354c4a4a890aba1b4828a1f4b1031 [file] [log] [blame]
Damien George04b91472014-05-03 23:27:38 +01001/*
2 * This file is part of the Micro Python 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
Damien George438c88d2014-02-22 19:25:23 +000027#include <string.h>
28#include <assert.h>
29
Damien George51dfcb42015-01-01 20:27:54 +000030#include "py/mpz.h"
Damien George438c88d2014-02-22 19:25:23 +000031
Damien George4c02e542015-10-01 17:57:36 +010032// this is only needed for mp_not_implemented, which should eventually be removed
33#include "py/runtime.h"
34
Damien George438c88d2014-02-22 19:25:23 +000035#if MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ
36
Damien George06201ff2014-03-01 19:50:50 +000037#define DIG_SIZE (MPZ_DIG_SIZE)
stijn0e557fa2014-10-30 14:39:22 +010038#define DIG_MASK ((MPZ_LONG_1 << DIG_SIZE) - 1)
39#define DIG_MSB (MPZ_LONG_1 << (DIG_SIZE - 1))
40#define DIG_BASE (MPZ_LONG_1 << DIG_SIZE)
Damien George438c88d2014-02-22 19:25:23 +000041
42/*
Damien George06201ff2014-03-01 19:50:50 +000043 mpz is an arbitrary precision integer type with a public API.
44
45 mpn functions act on non-negative integers represented by an array of generalised
46 digits (eg a word per digit). You also need to specify separately the length of the
47 array. There is no public API for mpn. Rather, the functions are used by mpz to
48 implement its features.
49
50 Integer values are stored little endian (first digit is first in memory).
51
52 Definition of normalise: ?
Damien George438c88d2014-02-22 19:25:23 +000053*/
54
55/* compares i with j
56 returns sign(i - j)
57 assumes i, j are normalised
58*/
Damien George42f3de92014-10-03 17:44:14 +000059STATIC int mpn_cmp(const mpz_dig_t *idig, mp_uint_t ilen, const mpz_dig_t *jdig, mp_uint_t jlen) {
Damien George438c88d2014-02-22 19:25:23 +000060 if (ilen < jlen) { return -1; }
61 if (ilen > jlen) { return 1; }
62
63 for (idig += ilen, jdig += ilen; ilen > 0; --ilen) {
Damien George9a21d2e2014-09-06 17:15:34 +010064 mpz_dbl_dig_signed_t cmp = (mpz_dbl_dig_t)*(--idig) - (mpz_dbl_dig_t)*(--jdig);
Damien George438c88d2014-02-22 19:25:23 +000065 if (cmp < 0) { return -1; }
66 if (cmp > 0) { return 1; }
67 }
68
69 return 0;
70}
71
Damien Georgec5ac2ac2014-02-26 16:56:30 +000072/* computes i = j << n
73 returns number of digits in i
Damien George06201ff2014-03-01 19:50:50 +000074 assumes enough memory in i; assumes normalised j; assumes n > 0
Damien Georgec5ac2ac2014-02-26 16:56:30 +000075 can have i, j pointing to same memory
76*/
Damien Georgeafb1cf72014-09-05 20:37:06 +010077STATIC mp_uint_t mpn_shl(mpz_dig_t *idig, mpz_dig_t *jdig, mp_uint_t jlen, mp_uint_t n) {
78 mp_uint_t n_whole = (n + DIG_SIZE - 1) / DIG_SIZE;
79 mp_uint_t n_part = n % DIG_SIZE;
Damien Georgebb4a43f2014-03-12 15:36:06 +000080 if (n_part == 0) {
81 n_part = DIG_SIZE;
82 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +000083
Damien George06201ff2014-03-01 19:50:50 +000084 // start from the high end of the digit arrays
85 idig += jlen + n_whole - 1;
86 jdig += jlen - 1;
Damien Georgec5ac2ac2014-02-26 16:56:30 +000087
Damien George06201ff2014-03-01 19:50:50 +000088 // shift the digits
89 mpz_dbl_dig_t d = 0;
Damien Georgeafb1cf72014-09-05 20:37:06 +010090 for (mp_uint_t i = jlen; i > 0; i--, idig--, jdig--) {
Damien George06201ff2014-03-01 19:50:50 +000091 d |= *jdig;
Damien George5d9b8162014-08-07 14:27:48 +000092 *idig = (d >> (DIG_SIZE - n_part)) & DIG_MASK;
Damien George06201ff2014-03-01 19:50:50 +000093 d <<= DIG_SIZE;
Damien Georgec5ac2ac2014-02-26 16:56:30 +000094 }
95
Damien George06201ff2014-03-01 19:50:50 +000096 // store remaining bits
Damien George5d9b8162014-08-07 14:27:48 +000097 *idig = (d >> (DIG_SIZE - n_part)) & DIG_MASK;
Damien George06201ff2014-03-01 19:50:50 +000098 idig -= n_whole - 1;
Damien Georgebb4a43f2014-03-12 15:36:06 +000099 memset(idig, 0, (n_whole - 1) * sizeof(mpz_dig_t));
Damien George06201ff2014-03-01 19:50:50 +0000100
101 // work out length of result
102 jlen += n_whole;
Damien Georgef66ee4d2015-04-22 23:16:49 +0100103 while (jlen != 0 && idig[jlen - 1] == 0) {
Damien George06201ff2014-03-01 19:50:50 +0000104 jlen--;
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000105 }
106
Damien George06201ff2014-03-01 19:50:50 +0000107 // return length of result
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000108 return jlen;
109}
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000110
Damien George438c88d2014-02-22 19:25:23 +0000111/* computes i = j >> n
112 returns number of digits in i
Damien George06201ff2014-03-01 19:50:50 +0000113 assumes enough memory in i; assumes normalised j; assumes n > 0
Damien George438c88d2014-02-22 19:25:23 +0000114 can have i, j pointing to same memory
115*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100116STATIC mp_uint_t mpn_shr(mpz_dig_t *idig, mpz_dig_t *jdig, mp_uint_t jlen, mp_uint_t n) {
117 mp_uint_t n_whole = n / DIG_SIZE;
118 mp_uint_t n_part = n % DIG_SIZE;
Damien George438c88d2014-02-22 19:25:23 +0000119
120 if (n_whole >= jlen) {
121 return 0;
122 }
123
124 jdig += n_whole;
125 jlen -= n_whole;
126
Damien Georgeafb1cf72014-09-05 20:37:06 +0100127 for (mp_uint_t i = jlen; i > 0; i--, idig++, jdig++) {
Damien George438c88d2014-02-22 19:25:23 +0000128 mpz_dbl_dig_t d = *jdig;
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000129 if (i > 1) {
Damien George9a21d2e2014-09-06 17:15:34 +0100130 d |= (mpz_dbl_dig_t)jdig[1] << DIG_SIZE;
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000131 }
Damien George438c88d2014-02-22 19:25:23 +0000132 d >>= n_part;
133 *idig = d & DIG_MASK;
134 }
135
136 if (idig[-1] == 0) {
Damien George06201ff2014-03-01 19:50:50 +0000137 jlen--;
Damien George438c88d2014-02-22 19:25:23 +0000138 }
139
140 return jlen;
141}
142
143/* computes i = j + k
144 returns number of digits in i
145 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
146 can have i, j, k pointing to same memory
147*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100148STATIC mp_uint_t mpn_add(mpz_dig_t *idig, const mpz_dig_t *jdig, mp_uint_t jlen, const mpz_dig_t *kdig, mp_uint_t klen) {
Damien George438c88d2014-02-22 19:25:23 +0000149 mpz_dig_t *oidig = idig;
150 mpz_dbl_dig_t carry = 0;
151
152 jlen -= klen;
153
154 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
Damien George9a21d2e2014-09-06 17:15:34 +0100155 carry += (mpz_dbl_dig_t)*jdig + (mpz_dbl_dig_t)*kdig;
Damien George438c88d2014-02-22 19:25:23 +0000156 *idig = carry & DIG_MASK;
157 carry >>= DIG_SIZE;
158 }
159
160 for (; jlen > 0; --jlen, ++idig, ++jdig) {
161 carry += *jdig;
162 *idig = carry & DIG_MASK;
163 carry >>= DIG_SIZE;
164 }
165
166 if (carry != 0) {
167 *idig++ = carry;
168 }
169
170 return idig - oidig;
171}
172
173/* computes i = j - k
174 returns number of digits in i
175 assumes enough memory in i; assumes normalised j, k; assumes j >= k
176 can have i, j, k pointing to same memory
177*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100178STATIC mp_uint_t mpn_sub(mpz_dig_t *idig, const mpz_dig_t *jdig, mp_uint_t jlen, const mpz_dig_t *kdig, mp_uint_t klen) {
Damien George438c88d2014-02-22 19:25:23 +0000179 mpz_dig_t *oidig = idig;
180 mpz_dbl_dig_signed_t borrow = 0;
181
182 jlen -= klen;
183
184 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
Damien George9a21d2e2014-09-06 17:15:34 +0100185 borrow += (mpz_dbl_dig_t)*jdig - (mpz_dbl_dig_t)*kdig;
Damien George438c88d2014-02-22 19:25:23 +0000186 *idig = borrow & DIG_MASK;
187 borrow >>= DIG_SIZE;
188 }
189
Damien Georgeaca14122014-02-24 21:32:52 +0000190 for (; jlen > 0; --jlen, ++idig, ++jdig) {
Damien George438c88d2014-02-22 19:25:23 +0000191 borrow += *jdig;
192 *idig = borrow & DIG_MASK;
193 borrow >>= DIG_SIZE;
194 }
195
196 for (--idig; idig >= oidig && *idig == 0; --idig) {
197 }
198
199 return idig + 1 - oidig;
200}
201
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200202/* computes i = j & k
203 returns number of digits in i
Damien Georgeff8dd3f2015-01-20 12:47:20 +0000204 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen (jlen argument not needed)
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200205 can have i, j, k pointing to same memory
206*/
Damien Georgeff8dd3f2015-01-20 12:47:20 +0000207STATIC mp_uint_t mpn_and(mpz_dig_t *idig, const mpz_dig_t *jdig, const mpz_dig_t *kdig, mp_uint_t klen) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200208 mpz_dig_t *oidig = idig;
209
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200210 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
211 *idig = *jdig & *kdig;
212 }
213
Damien George561e4252014-05-12 23:27:29 +0100214 // remove trailing zeros
Damien George51fab282014-05-13 22:58:00 +0100215 for (--idig; idig >= oidig && *idig == 0; --idig) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200216 }
217
Damien George51fab282014-05-13 22:58:00 +0100218 return idig + 1 - oidig;
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200219}
220
Damien Georgef55cf102014-05-29 15:01:49 +0100221/* computes i = j & -k = j & (~k + 1)
222 returns number of digits in i
223 assumes enough memory in i; assumes normalised j, k
224 can have i, j, k pointing to same memory
225*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100226STATIC mp_uint_t mpn_and_neg(mpz_dig_t *idig, const mpz_dig_t *jdig, mp_uint_t jlen, const mpz_dig_t *kdig, mp_uint_t klen) {
Damien Georgef55cf102014-05-29 15:01:49 +0100227 mpz_dig_t *oidig = idig;
228 mpz_dbl_dig_t carry = 1;
229
230 for (; jlen > 0 && klen > 0; --jlen, --klen, ++idig, ++jdig, ++kdig) {
231 carry += *kdig ^ DIG_MASK;
232 *idig = (*jdig & carry) & DIG_MASK;
233 carry >>= DIG_SIZE;
234 }
235
236 for (; jlen > 0; --jlen, ++idig, ++jdig) {
237 carry += DIG_MASK;
238 *idig = (*jdig & carry) & DIG_MASK;
239 carry >>= DIG_SIZE;
240 }
241
242 if (carry != 0) {
243 *idig = carry;
244 } else {
245 // remove trailing zeros
246 for (--idig; idig >= oidig && *idig == 0; --idig) {
247 }
248 }
249
250 return idig + 1 - oidig;
251}
252
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200253/* computes i = j | k
254 returns number of digits in i
255 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
256 can have i, j, k pointing to same memory
257*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100258STATIC mp_uint_t mpn_or(mpz_dig_t *idig, const mpz_dig_t *jdig, mp_uint_t jlen, const mpz_dig_t *kdig, mp_uint_t klen) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200259 mpz_dig_t *oidig = idig;
260
261 jlen -= klen;
262
263 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
264 *idig = *jdig | *kdig;
265 }
266
267 for (; jlen > 0; --jlen, ++idig, ++jdig) {
268 *idig = *jdig;
269 }
270
271 return idig - oidig;
272}
273
274/* computes i = j ^ k
275 returns number of digits in i
276 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
277 can have i, j, k pointing to same memory
278*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100279STATIC mp_uint_t mpn_xor(mpz_dig_t *idig, const mpz_dig_t *jdig, mp_uint_t jlen, const mpz_dig_t *kdig, mp_uint_t klen) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200280 mpz_dig_t *oidig = idig;
281
282 jlen -= klen;
283
284 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
285 *idig = *jdig ^ *kdig;
286 }
287
288 for (; jlen > 0; --jlen, ++idig, ++jdig) {
289 *idig = *jdig;
290 }
291
292 return idig - oidig;
293}
294
Damien George438c88d2014-02-22 19:25:23 +0000295/* computes i = i * d1 + d2
296 returns number of digits in i
297 assumes enough memory in i; assumes normalised i; assumes dmul != 0
298*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100299STATIC mp_uint_t mpn_mul_dig_add_dig(mpz_dig_t *idig, mp_uint_t ilen, mpz_dig_t dmul, mpz_dig_t dadd) {
Damien George438c88d2014-02-22 19:25:23 +0000300 mpz_dig_t *oidig = idig;
301 mpz_dbl_dig_t carry = dadd;
302
303 for (; ilen > 0; --ilen, ++idig) {
Damien George9a21d2e2014-09-06 17:15:34 +0100304 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
Damien George438c88d2014-02-22 19:25:23 +0000305 *idig = carry & DIG_MASK;
306 carry >>= DIG_SIZE;
307 }
308
309 if (carry != 0) {
310 *idig++ = carry;
311 }
312
313 return idig - oidig;
314}
315
316/* computes i = j * k
317 returns number of digits in i
318 assumes enough memory in i; assumes i is zeroed; assumes normalised j, k
319 can have j, k point to same memory
320*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100321STATIC mp_uint_t mpn_mul(mpz_dig_t *idig, mpz_dig_t *jdig, mp_uint_t jlen, mpz_dig_t *kdig, mp_uint_t klen) {
Damien George438c88d2014-02-22 19:25:23 +0000322 mpz_dig_t *oidig = idig;
Damien Georgeafb1cf72014-09-05 20:37:06 +0100323 mp_uint_t ilen = 0;
Damien George438c88d2014-02-22 19:25:23 +0000324
325 for (; klen > 0; --klen, ++idig, ++kdig) {
326 mpz_dig_t *id = idig;
327 mpz_dbl_dig_t carry = 0;
328
Damien Georgeafb1cf72014-09-05 20:37:06 +0100329 mp_uint_t jl = jlen;
Damien George438c88d2014-02-22 19:25:23 +0000330 for (mpz_dig_t *jd = jdig; jl > 0; --jl, ++jd, ++id) {
Damien George9a21d2e2014-09-06 17:15:34 +0100331 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
Damien George438c88d2014-02-22 19:25:23 +0000332 *id = carry & DIG_MASK;
333 carry >>= DIG_SIZE;
334 }
335
336 if (carry != 0) {
337 *id++ = carry;
338 }
339
340 ilen = id - oidig;
341 }
342
343 return ilen;
344}
345
346/* natural_div - quo * den + new_num = old_num (ie num is replaced with rem)
347 assumes den != 0
348 assumes num_dig has enough memory to be extended by 1 digit
349 assumes quo_dig has enough memory (as many digits as num)
350 assumes quo_dig is filled with zeros
351 modifies den_dig memory, but restors it to original state at end
352*/
353
Damien George40f3c022014-07-03 13:25:24 +0100354STATIC void mpn_div(mpz_dig_t *num_dig, mp_uint_t *num_len, mpz_dig_t *den_dig, mp_uint_t den_len, mpz_dig_t *quo_dig, mp_uint_t *quo_len) {
Damien George438c88d2014-02-22 19:25:23 +0000355 mpz_dig_t *orig_num_dig = num_dig;
356 mpz_dig_t *orig_quo_dig = quo_dig;
357 mpz_dig_t norm_shift = 0;
358 mpz_dbl_dig_t lead_den_digit;
359
360 // handle simple cases
361 {
Damien George42f3de92014-10-03 17:44:14 +0000362 int cmp = mpn_cmp(num_dig, *num_len, den_dig, den_len);
Damien George438c88d2014-02-22 19:25:23 +0000363 if (cmp == 0) {
364 *num_len = 0;
365 quo_dig[0] = 1;
366 *quo_len = 1;
367 return;
368 } else if (cmp < 0) {
369 // numerator remains the same
370 *quo_len = 0;
371 return;
372 }
373 }
374
375 // count number of leading zeros in leading digit of denominator
376 {
377 mpz_dig_t d = den_dig[den_len - 1];
Damien George9a21d2e2014-09-06 17:15:34 +0100378 while ((d & DIG_MSB) == 0) {
Damien George438c88d2014-02-22 19:25:23 +0000379 d <<= 1;
380 ++norm_shift;
381 }
382 }
383
384 // normalise denomenator (leading bit of leading digit is 1)
385 for (mpz_dig_t *den = den_dig, carry = 0; den < den_dig + den_len; ++den) {
386 mpz_dig_t d = *den;
387 *den = ((d << norm_shift) | carry) & DIG_MASK;
388 carry = d >> (DIG_SIZE - norm_shift);
389 }
390
391 // now need to shift numerator by same amount as denominator
392 // first, increase length of numerator in case we need more room to shift
393 num_dig[*num_len] = 0;
394 ++(*num_len);
395 for (mpz_dig_t *num = num_dig, carry = 0; num < num_dig + *num_len; ++num) {
396 mpz_dig_t n = *num;
397 *num = ((n << norm_shift) | carry) & DIG_MASK;
398 carry = n >> (DIG_SIZE - norm_shift);
399 }
400
401 // cache the leading digit of the denominator
402 lead_den_digit = den_dig[den_len - 1];
403
404 // point num_dig to last digit in numerator
405 num_dig += *num_len - 1;
406
407 // calculate number of digits in quotient
408 *quo_len = *num_len - den_len;
409
410 // point to last digit to store for quotient
411 quo_dig += *quo_len - 1;
412
413 // keep going while we have enough digits to divide
414 while (*num_len > den_len) {
Damien George9a21d2e2014-09-06 17:15:34 +0100415 mpz_dbl_dig_t quo = ((mpz_dbl_dig_t)*num_dig << DIG_SIZE) | num_dig[-1];
Damien George438c88d2014-02-22 19:25:23 +0000416
417 // get approximate quotient
418 quo /= lead_den_digit;
419
Damien George9a21d2e2014-09-06 17:15:34 +0100420 // Multiply quo by den and subtract from num to get remainder.
421 // We have different code here to handle different compile-time
422 // configurations of mpz:
423 //
424 // 1. DIG_SIZE is stricly less than half the number of bits
425 // available in mpz_dbl_dig_t. In this case we can use a
426 // slightly more optimal (in time and space) routine that
427 // uses the extra bits in mpz_dbl_dig_signed_t to store a
428 // sign bit.
429 //
430 // 2. DIG_SIZE is exactly half the number of bits available in
431 // mpz_dbl_dig_t. In this (common) case we need to be careful
432 // not to overflow the borrow variable. And the shifting of
433 // borrow needs some special logic (it's a shift right with
434 // round up).
435
436 if (DIG_SIZE < 8 * sizeof(mpz_dbl_dig_t) / 2) {
Damien George438c88d2014-02-22 19:25:23 +0000437 mpz_dbl_dig_signed_t borrow = 0;
438
439 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
Damien George9a21d2e2014-09-06 17:15:34 +0100440 borrow += (mpz_dbl_dig_t)*n - (mpz_dbl_dig_t)quo * (mpz_dbl_dig_t)*d; // will overflow if DIG_SIZE >= 8*sizeof(mpz_dbl_dig_t)/2
Damien George438c88d2014-02-22 19:25:23 +0000441 *n = borrow & DIG_MASK;
442 borrow >>= DIG_SIZE;
443 }
Damien George9a21d2e2014-09-06 17:15:34 +0100444 borrow += *num_dig; // will overflow if DIG_SIZE >= 8*sizeof(mpz_dbl_dig_t)/2
Damien George438c88d2014-02-22 19:25:23 +0000445 *num_dig = borrow & DIG_MASK;
446 borrow >>= DIG_SIZE;
447
448 // adjust quotient if it is too big
449 for (; borrow != 0; --quo) {
450 mpz_dbl_dig_t carry = 0;
451 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
Damien George9a21d2e2014-09-06 17:15:34 +0100452 carry += (mpz_dbl_dig_t)*n + (mpz_dbl_dig_t)*d;
Damien George438c88d2014-02-22 19:25:23 +0000453 *n = carry & DIG_MASK;
454 carry >>= DIG_SIZE;
455 }
456 carry += *num_dig;
457 *num_dig = carry & DIG_MASK;
458 carry >>= DIG_SIZE;
459
460 borrow += carry;
461 }
Damien George9a21d2e2014-09-06 17:15:34 +0100462 } else { // DIG_SIZE == 8 * sizeof(mpz_dbl_dig_t) / 2
463 mpz_dbl_dig_t borrow = 0;
464
465 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
466 mpz_dbl_dig_t x = (mpz_dbl_dig_t)quo * (mpz_dbl_dig_t)(*d);
467 if (x >= *n || *n - x <= borrow) {
468 borrow += (mpz_dbl_dig_t)x - (mpz_dbl_dig_t)*n;
469 *n = (-borrow) & DIG_MASK;
470 borrow = (borrow >> DIG_SIZE) + ((borrow & DIG_MASK) == 0 ? 0 : 1); // shift-right with round-up
471 } else {
472 *n = ((mpz_dbl_dig_t)*n - (mpz_dbl_dig_t)x - (mpz_dbl_dig_t)borrow) & DIG_MASK;
473 borrow = 0;
474 }
475 }
476 if (borrow >= *num_dig) {
477 borrow -= (mpz_dbl_dig_t)*num_dig;
478 *num_dig = (-borrow) & DIG_MASK;
479 borrow = (borrow >> DIG_SIZE) + ((borrow & DIG_MASK) == 0 ? 0 : 1); // shift-right with round-up
480 } else {
481 *num_dig = (*num_dig - borrow) & DIG_MASK;
482 borrow = 0;
483 }
484
485 // adjust quotient if it is too big
486 for (; borrow != 0; --quo) {
487 mpz_dbl_dig_t carry = 0;
488 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
489 carry += (mpz_dbl_dig_t)*n + (mpz_dbl_dig_t)*d;
490 *n = carry & DIG_MASK;
491 carry >>= DIG_SIZE;
492 }
493 carry += (mpz_dbl_dig_t)*num_dig;
494 *num_dig = carry & DIG_MASK;
495 carry >>= DIG_SIZE;
496
497 //assert(borrow >= carry); // enable this to check the logic
498 borrow -= carry;
499 }
Damien George438c88d2014-02-22 19:25:23 +0000500 }
501
502 // store this digit of the quotient
503 *quo_dig = quo & DIG_MASK;
504 --quo_dig;
505
506 // move down to next digit of numerator
507 --num_dig;
508 --(*num_len);
509 }
510
511 // unnormalise denomenator
512 for (mpz_dig_t *den = den_dig + den_len - 1, carry = 0; den >= den_dig; --den) {
513 mpz_dig_t d = *den;
514 *den = ((d >> norm_shift) | carry) & DIG_MASK;
515 carry = d << (DIG_SIZE - norm_shift);
516 }
517
518 // unnormalise numerator (remainder now)
519 for (mpz_dig_t *num = orig_num_dig + *num_len - 1, carry = 0; num >= orig_num_dig; --num) {
520 mpz_dig_t n = *num;
521 *num = ((n >> norm_shift) | carry) & DIG_MASK;
522 carry = n << (DIG_SIZE - norm_shift);
523 }
524
525 // strip trailing zeros
526
527 while (*quo_len > 0 && orig_quo_dig[*quo_len - 1] == 0) {
528 --(*quo_len);
529 }
530
531 while (*num_len > 0 && orig_num_dig[*num_len - 1] == 0) {
532 --(*num_len);
533 }
534}
535
Damien George06201ff2014-03-01 19:50:50 +0000536#define MIN_ALLOC (2)
Damien George438c88d2014-02-22 19:25:23 +0000537
Damien George1c70cbf2014-08-30 00:38:16 +0100538STATIC const uint8_t log_base2_floor[] = {
Damien George438c88d2014-02-22 19:25:23 +0000539 0,
540 0, 1, 1, 2,
541 2, 2, 2, 3,
542 3, 3, 3, 3,
543 3, 3, 3, 4,
544 4, 4, 4, 4,
545 4, 4, 4, 4,
546 4, 4, 4, 4,
547 4, 4, 4, 5
548};
549
Damien George438c88d2014-02-22 19:25:23 +0000550void mpz_init_zero(mpz_t *z) {
Damien George438c88d2014-02-22 19:25:23 +0000551 z->neg = 0;
Damien George06201ff2014-03-01 19:50:50 +0000552 z->fixed_dig = 0;
553 z->alloc = 0;
Damien George438c88d2014-02-22 19:25:23 +0000554 z->len = 0;
555 z->dig = NULL;
556}
557
Damien George40f3c022014-07-03 13:25:24 +0100558void mpz_init_from_int(mpz_t *z, mp_int_t val) {
Damien George438c88d2014-02-22 19:25:23 +0000559 mpz_init_zero(z);
560 mpz_set_from_int(z, val);
561}
562
Damien Georgeafb1cf72014-09-05 20:37:06 +0100563void mpz_init_fixed_from_int(mpz_t *z, mpz_dig_t *dig, mp_uint_t alloc, mp_int_t val) {
Damien George06201ff2014-03-01 19:50:50 +0000564 z->neg = 0;
565 z->fixed_dig = 1;
566 z->alloc = alloc;
567 z->len = 0;
568 z->dig = dig;
569 mpz_set_from_int(z, val);
570}
571
Damien George438c88d2014-02-22 19:25:23 +0000572void mpz_deinit(mpz_t *z) {
Damien George06201ff2014-03-01 19:50:50 +0000573 if (z != NULL && !z->fixed_dig) {
Damien George438c88d2014-02-22 19:25:23 +0000574 m_del(mpz_dig_t, z->dig, z->alloc);
575 }
576}
577
Damien George848dd0e2015-03-12 15:59:40 +0000578#if 0
579these functions are unused
580
Damien George438c88d2014-02-22 19:25:23 +0000581mpz_t *mpz_zero(void) {
582 mpz_t *z = m_new_obj(mpz_t);
583 mpz_init_zero(z);
584 return z;
585}
586
Damien George40f3c022014-07-03 13:25:24 +0100587mpz_t *mpz_from_int(mp_int_t val) {
Damien George438c88d2014-02-22 19:25:23 +0000588 mpz_t *z = mpz_zero();
589 mpz_set_from_int(z, val);
590 return z;
591}
592
Damien George95307432014-09-10 22:10:33 +0100593mpz_t *mpz_from_ll(long long val, bool is_signed) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000594 mpz_t *z = mpz_zero();
Damien George95307432014-09-10 22:10:33 +0100595 mpz_set_from_ll(z, val, is_signed);
Damien Georgebb4a43f2014-03-12 15:36:06 +0000596 return z;
597}
598
David Steinberg6e0b6d02015-01-02 12:39:22 +0000599#if MICROPY_PY_BUILTINS_FLOAT
600mpz_t *mpz_from_float(mp_float_t val) {
601 mpz_t *z = mpz_zero();
602 mpz_set_from_float(z, val);
603 return z;
604}
605#endif
606
Damien Georgeafb1cf72014-09-05 20:37:06 +0100607mpz_t *mpz_from_str(const char *str, mp_uint_t len, bool neg, mp_uint_t base) {
Damien George438c88d2014-02-22 19:25:23 +0000608 mpz_t *z = mpz_zero();
609 mpz_set_from_str(z, str, len, neg, base);
610 return z;
611}
Damien George848dd0e2015-03-12 15:59:40 +0000612#endif
Damien George438c88d2014-02-22 19:25:23 +0000613
Damien George848dd0e2015-03-12 15:59:40 +0000614STATIC void mpz_free(mpz_t *z) {
Damien George438c88d2014-02-22 19:25:23 +0000615 if (z != NULL) {
616 m_del(mpz_dig_t, z->dig, z->alloc);
617 m_del_obj(mpz_t, z);
618 }
619}
620
Damien Georgeafb1cf72014-09-05 20:37:06 +0100621STATIC void mpz_need_dig(mpz_t *z, mp_uint_t need) {
Damien George438c88d2014-02-22 19:25:23 +0000622 if (need < MIN_ALLOC) {
Damien George06201ff2014-03-01 19:50:50 +0000623 need = MIN_ALLOC;
Damien George438c88d2014-02-22 19:25:23 +0000624 }
625
Damien George06201ff2014-03-01 19:50:50 +0000626 if (z->dig == NULL || z->alloc < need) {
627 if (z->fixed_dig) {
628 // cannot reallocate fixed buffers
629 assert(0);
630 return;
631 }
632 z->dig = m_renew(mpz_dig_t, z->dig, z->alloc, need);
633 z->alloc = need;
Damien George438c88d2014-02-22 19:25:23 +0000634 }
635}
636
Damien George848dd0e2015-03-12 15:59:40 +0000637STATIC mpz_t *mpz_clone(const mpz_t *src) {
Damien George438c88d2014-02-22 19:25:23 +0000638 mpz_t *z = m_new_obj(mpz_t);
Damien George438c88d2014-02-22 19:25:23 +0000639 z->neg = src->neg;
Damien George06201ff2014-03-01 19:50:50 +0000640 z->fixed_dig = 0;
641 z->alloc = src->alloc;
Damien George438c88d2014-02-22 19:25:23 +0000642 z->len = src->len;
643 if (src->dig == NULL) {
644 z->dig = NULL;
645 } else {
646 z->dig = m_new(mpz_dig_t, z->alloc);
647 memcpy(z->dig, src->dig, src->alloc * sizeof(mpz_dig_t));
648 }
649 return z;
650}
651
Damien George06201ff2014-03-01 19:50:50 +0000652/* sets dest = src
653 can have dest, src the same
654*/
Damien George438c88d2014-02-22 19:25:23 +0000655void mpz_set(mpz_t *dest, const mpz_t *src) {
656 mpz_need_dig(dest, src->len);
657 dest->neg = src->neg;
658 dest->len = src->len;
659 memcpy(dest->dig, src->dig, src->len * sizeof(mpz_dig_t));
660}
661
Damien George40f3c022014-07-03 13:25:24 +0100662void mpz_set_from_int(mpz_t *z, mp_int_t val) {
Damien George58056b02015-01-09 20:58:58 +0000663 if (val == 0) {
664 z->len = 0;
665 return;
666 }
667
Damien George06201ff2014-03-01 19:50:50 +0000668 mpz_need_dig(z, MPZ_NUM_DIG_FOR_INT);
Damien George438c88d2014-02-22 19:25:23 +0000669
Damien George40f3c022014-07-03 13:25:24 +0100670 mp_uint_t uval;
Damien George438c88d2014-02-22 19:25:23 +0000671 if (val < 0) {
672 z->neg = 1;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000673 uval = -val;
Damien George438c88d2014-02-22 19:25:23 +0000674 } else {
675 z->neg = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000676 uval = val;
Damien George438c88d2014-02-22 19:25:23 +0000677 }
678
679 z->len = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000680 while (uval > 0) {
681 z->dig[z->len++] = uval & DIG_MASK;
682 uval >>= DIG_SIZE;
683 }
684}
685
Damien George95307432014-09-10 22:10:33 +0100686void mpz_set_from_ll(mpz_t *z, long long val, bool is_signed) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000687 mpz_need_dig(z, MPZ_NUM_DIG_FOR_LL);
688
689 unsigned long long uval;
Damien George95307432014-09-10 22:10:33 +0100690 if (is_signed && val < 0) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000691 z->neg = 1;
692 uval = -val;
693 } else {
694 z->neg = 0;
695 uval = val;
696 }
697
698 z->len = 0;
699 while (uval > 0) {
700 z->dig[z->len++] = uval & DIG_MASK;
701 uval >>= DIG_SIZE;
Damien George438c88d2014-02-22 19:25:23 +0000702 }
703}
704
David Steinberg6e0b6d02015-01-02 12:39:22 +0000705#if MICROPY_PY_BUILTINS_FLOAT
706void mpz_set_from_float(mpz_t *z, mp_float_t src) {
707#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
David Steinberg6e0b6d02015-01-02 12:39:22 +0000708typedef uint64_t mp_float_int_t;
David Steinbergc585ad12015-01-13 15:19:37 +0000709#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
David Steinberg6e0b6d02015-01-02 12:39:22 +0000710typedef uint32_t mp_float_int_t;
711#endif
712 union {
713 mp_float_t f;
David Steinbergc585ad12015-01-13 15:19:37 +0000714 struct { mp_float_int_t frc:MP_FLOAT_FRAC_BITS, exp:MP_FLOAT_EXP_BITS, sgn:1; } p;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000715 } u = {src};
716
717 z->neg = u.p.sgn;
718 if (u.p.exp == 0) {
719 // value == 0 || value < 1
Damien George58056b02015-01-09 20:58:58 +0000720 mpz_set_from_int(z, 0);
David Steinbergc585ad12015-01-13 15:19:37 +0000721 } else if (u.p.exp == ((1 << MP_FLOAT_EXP_BITS) - 1)) {
Damien George6fd4b362015-01-02 23:04:09 +0000722 // u.p.frc == 0 indicates inf, else NaN
723 // should be handled by caller
Damien George58056b02015-01-09 20:58:58 +0000724 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000725 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000726 const int adj_exp = (int)u.p.exp - MP_FLOAT_EXP_BIAS;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000727 if (adj_exp < 0) {
728 // value < 1 , truncates to 0
Damien George58056b02015-01-09 20:58:58 +0000729 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000730 } else if (adj_exp == 0) {
731 // 1 <= value < 2 , so truncates to 1
Damien George58056b02015-01-09 20:58:58 +0000732 mpz_set_from_int(z, 1);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000733 } else {
734 // 2 <= value
735 const int dig_cnt = (adj_exp + 1 + (DIG_SIZE - 1)) / DIG_SIZE;
736 const unsigned int rem = adj_exp % DIG_SIZE;
737 int dig_ind, shft;
David Steinbergc585ad12015-01-13 15:19:37 +0000738 mp_float_int_t frc = u.p.frc | ((mp_float_int_t)1 << MP_FLOAT_FRAC_BITS);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000739
David Steinbergc585ad12015-01-13 15:19:37 +0000740 if (adj_exp < MP_FLOAT_FRAC_BITS) {
David Steinberg6e0b6d02015-01-02 12:39:22 +0000741 shft = 0;
742 dig_ind = 0;
David Steinbergc585ad12015-01-13 15:19:37 +0000743 frc >>= MP_FLOAT_FRAC_BITS - adj_exp;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000744 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000745 shft = (rem - MP_FLOAT_FRAC_BITS) % DIG_SIZE;
746 dig_ind = (adj_exp - MP_FLOAT_FRAC_BITS) / DIG_SIZE;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000747 }
748 mpz_need_dig(z, dig_cnt);
749 z->len = dig_cnt;
750 if (dig_ind != 0) {
751 memset(z->dig, 0, dig_ind * sizeof(mpz_dig_t));
752 }
753 if (shft != 0) {
754 z->dig[dig_ind++] = (frc << shft) & DIG_MASK;
755 frc >>= DIG_SIZE - shft;
756 }
David Steinberg8d427b72015-01-13 15:20:32 +0000757#if DIG_SIZE < (MP_FLOAT_FRAC_BITS + 1)
David Steinberg6e0b6d02015-01-02 12:39:22 +0000758 while (dig_ind != dig_cnt) {
759 z->dig[dig_ind++] = frc & DIG_MASK;
760 frc >>= DIG_SIZE;
761 }
David Steinberg8d427b72015-01-13 15:20:32 +0000762#else
763 if (dig_ind != dig_cnt) {
764 z->dig[dig_ind] = frc;
765 }
766#endif
David Steinberg6e0b6d02015-01-02 12:39:22 +0000767 }
768 }
769}
David Steinberg6e0b6d02015-01-02 12:39:22 +0000770#endif
771
Damien George438c88d2014-02-22 19:25:23 +0000772// returns number of bytes from str that were processed
Damien Georgeafb1cf72014-09-05 20:37:06 +0100773mp_uint_t mpz_set_from_str(mpz_t *z, const char *str, mp_uint_t len, bool neg, mp_uint_t base) {
Damien George438c88d2014-02-22 19:25:23 +0000774 assert(base < 36);
775
776 const char *cur = str;
777 const char *top = str + len;
778
779 mpz_need_dig(z, len * 8 / DIG_SIZE + 1);
780
781 if (neg) {
782 z->neg = 1;
783 } else {
784 z->neg = 0;
785 }
786
787 z->len = 0;
788 for (; cur < top; ++cur) { // XXX UTF8 next char
Damien Georgeafb1cf72014-09-05 20:37:06 +0100789 //mp_uint_t v = char_to_numeric(cur#); // XXX UTF8 get char
790 mp_uint_t v = *cur;
Damien George438c88d2014-02-22 19:25:23 +0000791 if ('0' <= v && v <= '9') {
792 v -= '0';
793 } else if ('A' <= v && v <= 'Z') {
794 v -= 'A' - 10;
795 } else if ('a' <= v && v <= 'z') {
796 v -= 'a' - 10;
797 } else {
798 break;
799 }
800 if (v >= base) {
801 break;
802 }
803 z->len = mpn_mul_dig_add_dig(z->dig, z->len, base, v);
804 }
805
806 return cur - str;
807}
808
809bool mpz_is_zero(const mpz_t *z) {
810 return z->len == 0;
811}
812
Damien Georgea2e38382015-03-02 12:58:06 +0000813#if 0
814these functions are unused
815
Damien George438c88d2014-02-22 19:25:23 +0000816bool mpz_is_pos(const mpz_t *z) {
817 return z->len > 0 && z->neg == 0;
818}
819
820bool mpz_is_neg(const mpz_t *z) {
821 return z->len > 0 && z->neg != 0;
822}
823
824bool mpz_is_odd(const mpz_t *z) {
825 return z->len > 0 && (z->dig[0] & 1) != 0;
826}
827
828bool mpz_is_even(const mpz_t *z) {
829 return z->len == 0 || (z->dig[0] & 1) == 0;
830}
Damien Georgea2e38382015-03-02 12:58:06 +0000831#endif
Damien George438c88d2014-02-22 19:25:23 +0000832
Damien George42f3de92014-10-03 17:44:14 +0000833int mpz_cmp(const mpz_t *z1, const mpz_t *z2) {
Damien Georgea9dc9b82015-01-27 17:47:38 +0000834 // to catch comparison of -0 with +0
835 if (z1->len == 0 && z2->len == 0) {
836 return 0;
837 }
Damien George42f3de92014-10-03 17:44:14 +0000838 int cmp = (int)z2->neg - (int)z1->neg;
Damien George438c88d2014-02-22 19:25:23 +0000839 if (cmp != 0) {
840 return cmp;
841 }
842 cmp = mpn_cmp(z1->dig, z1->len, z2->dig, z2->len);
843 if (z1->neg != 0) {
844 cmp = -cmp;
845 }
846 return cmp;
847}
848
Damien George06201ff2014-03-01 19:50:50 +0000849#if 0
850// obsolete
851// compares mpz with an integer that fits within DIG_SIZE bits
Damien Georgeafb1cf72014-09-05 20:37:06 +0100852mp_int_t mpz_cmp_sml_int(const mpz_t *z, mp_int_t sml_int) {
853 mp_int_t cmp;
Damien George438c88d2014-02-22 19:25:23 +0000854 if (z->neg == 0) {
855 if (sml_int < 0) return 1;
856 if (sml_int == 0) {
857 if (z->len == 0) return 0;
858 return 1;
859 }
860 if (z->len == 0) return -1;
861 assert(sml_int < (1 << DIG_SIZE));
862 if (z->len != 1) return 1;
863 cmp = z->dig[0] - sml_int;
864 } else {
865 if (sml_int > 0) return -1;
866 if (sml_int == 0) {
867 if (z->len == 0) return 0;
868 return -1;
869 }
870 if (z->len == 0) return 1;
871 assert(sml_int > -(1 << DIG_SIZE));
872 if (z->len != 1) return -1;
873 cmp = -z->dig[0] - sml_int;
874 }
875 if (cmp < 0) return -1;
876 if (cmp > 0) return 1;
877 return 0;
878}
Damien George06201ff2014-03-01 19:50:50 +0000879#endif
Damien George438c88d2014-02-22 19:25:23 +0000880
Damien George438c88d2014-02-22 19:25:23 +0000881#if 0
882these functions are unused
883
884/* returns abs(z)
885*/
886mpz_t *mpz_abs(const mpz_t *z) {
887 mpz_t *z2 = mpz_clone(z);
888 z2->neg = 0;
889 return z2;
890}
891
892/* returns -z
893*/
894mpz_t *mpz_neg(const mpz_t *z) {
895 mpz_t *z2 = mpz_clone(z);
896 z2->neg = 1 - z2->neg;
897 return z2;
898}
899
900/* returns lhs + rhs
901 can have lhs, rhs the same
902*/
903mpz_t *mpz_add(const mpz_t *lhs, const mpz_t *rhs) {
904 mpz_t *z = mpz_zero();
905 mpz_add_inpl(z, lhs, rhs);
906 return z;
907}
908
909/* returns lhs - rhs
910 can have lhs, rhs the same
911*/
912mpz_t *mpz_sub(const mpz_t *lhs, const mpz_t *rhs) {
913 mpz_t *z = mpz_zero();
914 mpz_sub_inpl(z, lhs, rhs);
915 return z;
916}
917
918/* returns lhs * rhs
919 can have lhs, rhs the same
920*/
921mpz_t *mpz_mul(const mpz_t *lhs, const mpz_t *rhs) {
922 mpz_t *z = mpz_zero();
923 mpz_mul_inpl(z, lhs, rhs);
924 return z;
925}
926
927/* returns lhs ** rhs
928 can have lhs, rhs the same
929*/
930mpz_t *mpz_pow(const mpz_t *lhs, const mpz_t *rhs) {
931 mpz_t *z = mpz_zero();
932 mpz_pow_inpl(z, lhs, rhs);
933 return z;
934}
Damien Georgea2e38382015-03-02 12:58:06 +0000935
936/* computes new integers in quo and rem such that:
937 quo * rhs + rem = lhs
938 0 <= rem < rhs
939 can have lhs, rhs the same
940*/
941void mpz_divmod(const mpz_t *lhs, const mpz_t *rhs, mpz_t **quo, mpz_t **rem) {
942 *quo = mpz_zero();
943 *rem = mpz_zero();
944 mpz_divmod_inpl(*quo, *rem, lhs, rhs);
945}
Damien George438c88d2014-02-22 19:25:23 +0000946#endif
947
948/* computes dest = abs(z)
949 can have dest, z the same
950*/
951void mpz_abs_inpl(mpz_t *dest, const mpz_t *z) {
952 if (dest != z) {
953 mpz_set(dest, z);
954 }
955 dest->neg = 0;
956}
957
958/* computes dest = -z
959 can have dest, z the same
960*/
961void mpz_neg_inpl(mpz_t *dest, const mpz_t *z) {
962 if (dest != z) {
963 mpz_set(dest, z);
964 }
965 dest->neg = 1 - dest->neg;
966}
967
Damien George06201ff2014-03-01 19:50:50 +0000968/* computes dest = ~z (= -z - 1)
969 can have dest, z the same
970*/
971void mpz_not_inpl(mpz_t *dest, const mpz_t *z) {
972 if (dest != z) {
973 mpz_set(dest, z);
974 }
Damien Georgee0ac1942014-12-31 19:35:01 +0000975 if (dest->len == 0) {
976 mpz_need_dig(dest, 1);
977 dest->dig[0] = 1;
978 dest->len = 1;
979 dest->neg = 1;
980 } else if (dest->neg) {
Damien George06201ff2014-03-01 19:50:50 +0000981 dest->neg = 0;
982 mpz_dig_t k = 1;
983 dest->len = mpn_sub(dest->dig, dest->dig, dest->len, &k, 1);
984 } else {
Damien Georgee0ac1942014-12-31 19:35:01 +0000985 mpz_need_dig(dest, dest->len + 1);
Damien George06201ff2014-03-01 19:50:50 +0000986 mpz_dig_t k = 1;
987 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &k, 1);
988 dest->neg = 1;
989 }
990}
991
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000992/* computes dest = lhs << rhs
993 can have dest, lhs the same
994*/
Damien George2f4e8512015-10-01 18:01:37 +0100995void mpz_shl_inpl(mpz_t *dest, const mpz_t *lhs, mp_uint_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +0000996 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000997 mpz_set(dest, lhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000998 } else {
Damien George06201ff2014-03-01 19:50:50 +0000999 mpz_need_dig(dest, lhs->len + (rhs + DIG_SIZE - 1) / DIG_SIZE);
1000 dest->len = mpn_shl(dest->dig, lhs->dig, lhs->len, rhs);
1001 dest->neg = lhs->neg;
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001002 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001003}
1004
1005/* computes dest = lhs >> rhs
1006 can have dest, lhs the same
1007*/
Damien George2f4e8512015-10-01 18:01:37 +01001008void mpz_shr_inpl(mpz_t *dest, const mpz_t *lhs, mp_uint_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +00001009 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001010 mpz_set(dest, lhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001011 } else {
Damien George06201ff2014-03-01 19:50:50 +00001012 mpz_need_dig(dest, lhs->len);
1013 dest->len = mpn_shr(dest->dig, lhs->dig, lhs->len, rhs);
1014 dest->neg = lhs->neg;
1015 if (dest->neg) {
1016 // arithmetic shift right, rounding to negative infinity
Damien Georgeafb1cf72014-09-05 20:37:06 +01001017 mp_uint_t n_whole = rhs / DIG_SIZE;
1018 mp_uint_t n_part = rhs % DIG_SIZE;
Damien George06201ff2014-03-01 19:50:50 +00001019 mpz_dig_t round_up = 0;
Damien Georgeafb1cf72014-09-05 20:37:06 +01001020 for (mp_uint_t i = 0; i < lhs->len && i < n_whole; i++) {
Damien George06201ff2014-03-01 19:50:50 +00001021 if (lhs->dig[i] != 0) {
1022 round_up = 1;
1023 break;
1024 }
1025 }
1026 if (n_whole < lhs->len && (lhs->dig[n_whole] & ((1 << n_part) - 1)) != 0) {
1027 round_up = 1;
1028 }
1029 if (round_up) {
Damien Georgee0ac1942014-12-31 19:35:01 +00001030 if (dest->len == 0) {
1031 // dest == 0, so need to add 1 by hand (answer will be -1)
1032 dest->dig[0] = 1;
1033 dest->len = 1;
1034 } else {
1035 // dest > 0, so can use mpn_add to add 1
1036 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &round_up, 1);
1037 }
Damien George06201ff2014-03-01 19:50:50 +00001038 }
1039 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001040 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001041}
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001042
Damien George438c88d2014-02-22 19:25:23 +00001043/* computes dest = lhs + rhs
1044 can have dest, lhs, rhs the same
1045*/
1046void mpz_add_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1047 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1048 const mpz_t *temp = lhs;
1049 lhs = rhs;
1050 rhs = temp;
1051 }
1052
1053 if (lhs->neg == rhs->neg) {
1054 mpz_need_dig(dest, lhs->len + 1);
1055 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1056 } else {
1057 mpz_need_dig(dest, lhs->len);
1058 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1059 }
1060
1061 dest->neg = lhs->neg;
1062}
1063
1064/* computes dest = lhs - rhs
1065 can have dest, lhs, rhs the same
1066*/
1067void mpz_sub_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1068 bool neg = false;
1069
1070 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1071 const mpz_t *temp = lhs;
1072 lhs = rhs;
1073 rhs = temp;
1074 neg = true;
1075 }
1076
1077 if (lhs->neg != rhs->neg) {
1078 mpz_need_dig(dest, lhs->len + 1);
1079 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1080 } else {
1081 mpz_need_dig(dest, lhs->len);
1082 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1083 }
1084
1085 if (neg) {
1086 dest->neg = 1 - lhs->neg;
1087 } else {
1088 dest->neg = lhs->neg;
1089 }
1090}
1091
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001092/* computes dest = lhs & rhs
1093 can have dest, lhs, rhs the same
1094*/
1095void mpz_and_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001096 if (lhs->neg == rhs->neg) {
Damien Georgef55cf102014-05-29 15:01:49 +01001097 if (lhs->neg == 0) {
1098 // make sure lhs has the most digits
1099 if (lhs->len < rhs->len) {
1100 const mpz_t *temp = lhs;
1101 lhs = rhs;
1102 rhs = temp;
1103 }
1104 // do the and'ing
1105 mpz_need_dig(dest, rhs->len);
Damien Georgeff8dd3f2015-01-20 12:47:20 +00001106 dest->len = mpn_and(dest->dig, lhs->dig, rhs->dig, rhs->len);
Damien Georgef55cf102014-05-29 15:01:49 +01001107 dest->neg = 0;
1108 } else {
1109 // TODO both args are negative
Damien George4c02e542015-10-01 17:57:36 +01001110 mp_not_implemented("bignum and with negative args");
Damien Georgef55cf102014-05-29 15:01:49 +01001111 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001112 } else {
Damien Georgef55cf102014-05-29 15:01:49 +01001113 // args have different sign
1114 // make sure lhs is the positive arg
1115 if (rhs->neg == 0) {
1116 const mpz_t *temp = lhs;
1117 lhs = rhs;
1118 rhs = temp;
1119 }
1120 mpz_need_dig(dest, lhs->len + 1);
1121 dest->len = mpn_and_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1122 assert(dest->len <= dest->alloc);
1123 dest->neg = 0;
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001124 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001125}
1126
1127/* computes dest = lhs | rhs
1128 can have dest, lhs, rhs the same
1129*/
1130void mpz_or_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1131 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1132 const mpz_t *temp = lhs;
1133 lhs = rhs;
1134 rhs = temp;
1135 }
1136
1137 if (lhs->neg == rhs->neg) {
1138 mpz_need_dig(dest, lhs->len);
1139 dest->len = mpn_or(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1140 } else {
1141 mpz_need_dig(dest, lhs->len);
1142 // TODO
Damien George4c02e542015-10-01 17:57:36 +01001143 mp_not_implemented("bignum or with negative args");
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001144// dest->len = mpn_or_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1145 }
1146
1147 dest->neg = lhs->neg;
1148}
1149
1150/* computes dest = lhs ^ rhs
1151 can have dest, lhs, rhs the same
1152*/
1153void mpz_xor_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1154 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1155 const mpz_t *temp = lhs;
1156 lhs = rhs;
1157 rhs = temp;
1158 }
1159
1160 if (lhs->neg == rhs->neg) {
1161 mpz_need_dig(dest, lhs->len);
1162 dest->len = mpn_xor(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1163 } else {
1164 mpz_need_dig(dest, lhs->len);
1165 // TODO
Damien George4c02e542015-10-01 17:57:36 +01001166 mp_not_implemented("bignum xor with negative args");
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001167// dest->len = mpn_xor_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1168 }
1169
1170 dest->neg = 0;
1171}
1172
Damien George438c88d2014-02-22 19:25:23 +00001173/* computes dest = lhs * rhs
1174 can have dest, lhs, rhs the same
1175*/
Damien George4dea9222015-04-09 15:29:54 +00001176void mpz_mul_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Damien George438c88d2014-02-22 19:25:23 +00001177 if (lhs->len == 0 || rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001178 mpz_set_from_int(dest, 0);
1179 return;
Damien George438c88d2014-02-22 19:25:23 +00001180 }
1181
1182 mpz_t *temp = NULL;
1183 if (lhs == dest) {
1184 lhs = temp = mpz_clone(lhs);
1185 if (rhs == dest) {
1186 rhs = lhs;
1187 }
1188 } else if (rhs == dest) {
1189 rhs = temp = mpz_clone(rhs);
1190 }
1191
1192 mpz_need_dig(dest, lhs->len + rhs->len); // min mem l+r-1, max mem l+r
1193 memset(dest->dig, 0, dest->alloc * sizeof(mpz_dig_t));
1194 dest->len = mpn_mul(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1195
1196 if (lhs->neg == rhs->neg) {
1197 dest->neg = 0;
1198 } else {
1199 dest->neg = 1;
1200 }
1201
1202 mpz_free(temp);
1203}
1204
1205/* computes dest = lhs ** rhs
1206 can have dest, lhs, rhs the same
1207*/
1208void mpz_pow_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1209 if (lhs->len == 0 || rhs->neg != 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001210 mpz_set_from_int(dest, 0);
1211 return;
Damien George438c88d2014-02-22 19:25:23 +00001212 }
1213
1214 if (rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001215 mpz_set_from_int(dest, 1);
1216 return;
Damien George438c88d2014-02-22 19:25:23 +00001217 }
1218
1219 mpz_t *x = mpz_clone(lhs);
1220 mpz_t *n = mpz_clone(rhs);
1221
1222 mpz_set_from_int(dest, 1);
1223
1224 while (n->len > 0) {
Damien Georgea2e38382015-03-02 12:58:06 +00001225 if ((n->dig[0] & 1) != 0) {
Damien George438c88d2014-02-22 19:25:23 +00001226 mpz_mul_inpl(dest, dest, x);
1227 }
Damien George438c88d2014-02-22 19:25:23 +00001228 n->len = mpn_shr(n->dig, n->dig, n->len, 1);
Damien George5bf565e2014-04-04 00:16:32 +01001229 if (n->len == 0) {
1230 break;
1231 }
1232 mpz_mul_inpl(x, x, x);
Damien George438c88d2014-02-22 19:25:23 +00001233 }
1234
1235 mpz_free(x);
1236 mpz_free(n);
1237}
1238
Damien Georgea2e38382015-03-02 12:58:06 +00001239#if 0
1240these functions are unused
1241
Damien George438c88d2014-02-22 19:25:23 +00001242/* computes gcd(z1, z2)
1243 based on Knuth's modified gcd algorithm (I think?)
1244 gcd(z1, z2) >= 0
1245 gcd(0, 0) = 0
1246 gcd(z, 0) = abs(z)
1247*/
1248mpz_t *mpz_gcd(const mpz_t *z1, const mpz_t *z2) {
1249 if (z1->len == 0) {
1250 mpz_t *a = mpz_clone(z2);
1251 a->neg = 0;
1252 return a;
1253 } else if (z2->len == 0) {
1254 mpz_t *a = mpz_clone(z1);
1255 a->neg = 0;
1256 return a;
1257 }
1258
1259 mpz_t *a = mpz_clone(z1);
1260 mpz_t *b = mpz_clone(z2);
1261 mpz_t c; mpz_init_zero(&c);
1262 a->neg = 0;
1263 b->neg = 0;
1264
1265 for (;;) {
1266 if (mpz_cmp(a, b) < 0) {
1267 if (a->len == 0) {
1268 mpz_free(a);
1269 mpz_deinit(&c);
1270 return b;
1271 }
1272 mpz_t *t = a; a = b; b = t;
1273 }
1274 if (!(b->len >= 2 || (b->len == 1 && b->dig[0] > 1))) { // compute b > 0; could be mpz_cmp_small_int(b, 1) > 0
1275 break;
1276 }
1277 mpz_set(&c, b);
1278 do {
1279 mpz_add_inpl(&c, &c, &c);
1280 } while (mpz_cmp(&c, a) <= 0);
1281 c.len = mpn_shr(c.dig, c.dig, c.len, 1);
1282 mpz_sub_inpl(a, a, &c);
1283 }
1284
1285 mpz_deinit(&c);
1286
1287 if (b->len == 1 && b->dig[0] == 1) { // compute b == 1; could be mpz_cmp_small_int(b, 1) == 0
1288 mpz_free(a);
1289 return b;
1290 } else {
1291 mpz_free(b);
1292 return a;
1293 }
1294}
1295
1296/* computes lcm(z1, z2)
1297 = abs(z1) / gcd(z1, z2) * abs(z2)
1298 lcm(z1, z1) >= 0
1299 lcm(0, 0) = 0
1300 lcm(z, 0) = 0
1301*/
Damien George5d9b8162014-08-07 14:27:48 +00001302mpz_t *mpz_lcm(const mpz_t *z1, const mpz_t *z2) {
1303 if (z1->len == 0 || z2->len == 0) {
1304 return mpz_zero();
stijn01d6be42014-05-05 12:18:27 +02001305 }
Damien George438c88d2014-02-22 19:25:23 +00001306
1307 mpz_t *gcd = mpz_gcd(z1, z2);
1308 mpz_t *quo = mpz_zero();
1309 mpz_t *rem = mpz_zero();
1310 mpz_divmod_inpl(quo, rem, z1, gcd);
1311 mpz_mul_inpl(rem, quo, z2);
1312 mpz_free(gcd);
1313 mpz_free(quo);
1314 rem->neg = 0;
1315 return rem;
1316}
Damien Georgea2e38382015-03-02 12:58:06 +00001317#endif
Damien George438c88d2014-02-22 19:25:23 +00001318
1319/* computes new integers in quo and rem such that:
1320 quo * rhs + rem = lhs
1321 0 <= rem < rhs
1322 can have lhs, rhs the same
1323*/
1324void mpz_divmod_inpl(mpz_t *dest_quo, mpz_t *dest_rem, const mpz_t *lhs, const mpz_t *rhs) {
1325 if (rhs->len == 0) {
1326 mpz_set_from_int(dest_quo, 0);
1327 mpz_set_from_int(dest_rem, 0);
1328 return;
1329 }
1330
1331 mpz_need_dig(dest_quo, lhs->len + 1); // +1 necessary?
1332 memset(dest_quo->dig, 0, (lhs->len + 1) * sizeof(mpz_dig_t));
1333 dest_quo->len = 0;
1334 mpz_need_dig(dest_rem, lhs->len + 1); // +1 necessary?
1335 mpz_set(dest_rem, lhs);
1336 //rhs->dig[rhs->len] = 0;
1337 mpn_div(dest_rem->dig, &dest_rem->len, rhs->dig, rhs->len, dest_quo->dig, &dest_quo->len);
1338
1339 if (lhs->neg != rhs->neg) {
1340 dest_quo->neg = 1;
1341 }
1342}
1343
1344#if 0
1345these functions are unused
1346
1347/* computes floor(lhs / rhs)
1348 can have lhs, rhs the same
1349*/
1350mpz_t *mpz_div(const mpz_t *lhs, const mpz_t *rhs) {
1351 mpz_t *quo = mpz_zero();
1352 mpz_t rem; mpz_init_zero(&rem);
1353 mpz_divmod_inpl(quo, &rem, lhs, rhs);
1354 mpz_deinit(&rem);
1355 return quo;
1356}
1357
1358/* computes lhs % rhs ( >= 0)
1359 can have lhs, rhs the same
1360*/
1361mpz_t *mpz_mod(const mpz_t *lhs, const mpz_t *rhs) {
1362 mpz_t quo; mpz_init_zero(&quo);
1363 mpz_t *rem = mpz_zero();
1364 mpz_divmod_inpl(&quo, rem, lhs, rhs);
1365 mpz_deinit(&quo);
1366 return rem;
1367}
1368#endif
1369
Damien Georgeffe911d2014-07-24 14:21:37 +01001370// must return actual int value if it fits in mp_int_t
1371mp_int_t mpz_hash(const mpz_t *z) {
1372 mp_int_t val = 0;
1373 mpz_dig_t *d = z->dig + z->len;
1374
Damien George58056b02015-01-09 20:58:58 +00001375 while (d-- > z->dig) {
Damien Georgeffe911d2014-07-24 14:21:37 +01001376 val = (val << DIG_SIZE) | *d;
1377 }
1378
1379 if (z->neg != 0) {
1380 val = -val;
1381 }
1382
1383 return val;
1384}
1385
Damien George40f3c022014-07-03 13:25:24 +01001386bool mpz_as_int_checked(const mpz_t *i, mp_int_t *value) {
Damien George963a5a32015-01-16 17:47:07 +00001387 mp_uint_t val = 0;
Damien George8270e382014-04-03 11:00:54 +00001388 mpz_dig_t *d = i->dig + i->len;
1389
Damien George58056b02015-01-09 20:58:58 +00001390 while (d-- > i->dig) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001391 if (val > (~(WORD_MSBIT_HIGH) >> DIG_SIZE)) {
1392 // will overflow
Damien George8270e382014-04-03 11:00:54 +00001393 return false;
1394 }
Damien Georgec9aa58e2014-07-31 13:41:43 +00001395 val = (val << DIG_SIZE) | *d;
Damien George8270e382014-04-03 11:00:54 +00001396 }
1397
1398 if (i->neg != 0) {
1399 val = -val;
1400 }
1401
1402 *value = val;
1403 return true;
1404}
1405
Damien Georgec9aa58e2014-07-31 13:41:43 +00001406bool mpz_as_uint_checked(const mpz_t *i, mp_uint_t *value) {
1407 if (i->neg != 0) {
1408 // can't represent signed values
1409 return false;
1410 }
1411
1412 mp_uint_t val = 0;
1413 mpz_dig_t *d = i->dig + i->len;
1414
Damien George58056b02015-01-09 20:58:58 +00001415 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001416 if (val > (~(WORD_MSBIT_HIGH) >> (DIG_SIZE - 1))) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001417 // will overflow
1418 return false;
1419 }
1420 val = (val << DIG_SIZE) | *d;
1421 }
1422
1423 *value = val;
1424 return true;
1425}
1426
Damien George271d18e2015-04-25 23:16:39 +01001427// writes at most len bytes to buf (so buf should be zeroed before calling)
1428void mpz_as_bytes(const mpz_t *z, bool big_endian, mp_uint_t len, byte *buf) {
1429 byte *b = buf;
1430 if (big_endian) {
1431 b += len;
1432 }
1433 mpz_dig_t *zdig = z->dig;
1434 int bits = 0;
1435 mpz_dbl_dig_t d = 0;
1436 mpz_dbl_dig_t carry = 1;
1437 for (mp_uint_t zlen = z->len; zlen > 0; --zlen) {
1438 bits += DIG_SIZE;
1439 d = (d << DIG_SIZE) | *zdig++;
1440 for (; bits >= 8; bits -= 8, d >>= 8) {
1441 mpz_dig_t val = d;
1442 if (z->neg) {
Damien George94729072015-04-25 23:51:14 +01001443 val = (~val & 0xff) + carry;
1444 carry = val >> 8;
Damien George271d18e2015-04-25 23:16:39 +01001445 }
1446 if (big_endian) {
1447 *--b = val;
1448 if (b == buf) {
1449 return;
1450 }
1451 } else {
1452 *b++ = val;
1453 if (b == buf + len) {
1454 return;
1455 }
1456 }
1457 }
1458 }
1459}
1460
Damien Georgefb510b32014-06-01 13:32:54 +01001461#if MICROPY_PY_BUILTINS_FLOAT
Damien George52608102014-03-08 15:04:54 +00001462mp_float_t mpz_as_float(const mpz_t *i) {
1463 mp_float_t val = 0;
Damien George438c88d2014-02-22 19:25:23 +00001464 mpz_dig_t *d = i->dig + i->len;
1465
Damien George58056b02015-01-09 20:58:58 +00001466 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001467 val = val * DIG_BASE + *d;
Damien George438c88d2014-02-22 19:25:23 +00001468 }
1469
1470 if (i->neg != 0) {
1471 val = -val;
1472 }
1473
1474 return val;
1475}
Damien George52608102014-03-08 15:04:54 +00001476#endif
Damien George438c88d2014-02-22 19:25:23 +00001477
Damien Georgeafb1cf72014-09-05 20:37:06 +01001478mp_uint_t mpz_as_str_size(const mpz_t *i, mp_uint_t base, const char *prefix, char comma) {
Damien George438c88d2014-02-22 19:25:23 +00001479 if (base < 2 || base > 32) {
1480 return 0;
1481 }
1482
Damien Georgeafb1cf72014-09-05 20:37:06 +01001483 mp_uint_t num_digits = i->len * DIG_SIZE / log_base2_floor[base] + 1;
1484 mp_uint_t num_commas = comma ? num_digits / 3: 0;
1485 mp_uint_t prefix_len = prefix ? strlen(prefix) : 0;
Dave Hylandsc4029e52014-04-07 11:19:51 -07001486
1487 return num_digits + num_commas + prefix_len + 2; // +1 for sign, +1 for null byte
1488}
1489
Damien Georgeafb1cf72014-09-05 20:37:06 +01001490#if 0
1491this function is unused
1492char *mpz_as_str(const mpz_t *i, mp_uint_t base) {
1493 char *s = m_new(char, mpz_as_str_size(i, base, NULL, '\0'));
1494 mpz_as_str_inpl(i, base, NULL, 'a', '\0', s);
Damien George438c88d2014-02-22 19:25:23 +00001495 return s;
1496}
Damien Georgeafb1cf72014-09-05 20:37:06 +01001497#endif
Damien George438c88d2014-02-22 19:25:23 +00001498
1499// assumes enough space as calculated by mpz_as_str_size
1500// returns length of string, not including null byte
Damien Georgeafb1cf72014-09-05 20:37:06 +01001501mp_uint_t mpz_as_str_inpl(const mpz_t *i, mp_uint_t base, const char *prefix, char base_char, char comma, char *str) {
Damien George438c88d2014-02-22 19:25:23 +00001502 if (str == NULL || base < 2 || base > 32) {
1503 str[0] = 0;
1504 return 0;
1505 }
1506
Damien Georgeafb1cf72014-09-05 20:37:06 +01001507 mp_uint_t ilen = i->len;
Damien George438c88d2014-02-22 19:25:23 +00001508
Dave Hylandsc4029e52014-04-07 11:19:51 -07001509 char *s = str;
Damien George438c88d2014-02-22 19:25:23 +00001510 if (ilen == 0) {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001511 if (prefix) {
1512 while (*prefix)
1513 *s++ = *prefix++;
1514 }
1515 *s++ = '0';
1516 *s = '\0';
1517 return s - str;
Damien George438c88d2014-02-22 19:25:23 +00001518 }
1519
Damien Georgeeec91052014-04-08 23:11:00 +01001520 // make a copy of mpz digits, so we can do the div/mod calculation
Damien George438c88d2014-02-22 19:25:23 +00001521 mpz_dig_t *dig = m_new(mpz_dig_t, ilen);
1522 memcpy(dig, i->dig, ilen * sizeof(mpz_dig_t));
1523
1524 // convert
Dave Hylandsc4029e52014-04-07 11:19:51 -07001525 char *last_comma = str;
Damien George438c88d2014-02-22 19:25:23 +00001526 bool done;
1527 do {
1528 mpz_dig_t *d = dig + ilen;
1529 mpz_dbl_dig_t a = 0;
1530
1531 // compute next remainder
1532 while (--d >= dig) {
1533 a = (a << DIG_SIZE) | *d;
1534 *d = a / base;
1535 a %= base;
1536 }
1537
1538 // convert to character
1539 a += '0';
1540 if (a > '9') {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001541 a += base_char - '9' - 1;
Damien George438c88d2014-02-22 19:25:23 +00001542 }
1543 *s++ = a;
1544
1545 // check if number is zero
1546 done = true;
1547 for (d = dig; d < dig + ilen; ++d) {
1548 if (*d != 0) {
1549 done = false;
1550 break;
1551 }
1552 }
Dave Hylandsc4029e52014-04-07 11:19:51 -07001553 if (comma && (s - last_comma) == 3) {
1554 *s++ = comma;
1555 last_comma = s;
1556 }
1557 }
1558 while (!done);
Damien George438c88d2014-02-22 19:25:23 +00001559
Damien Georgeeec91052014-04-08 23:11:00 +01001560 // free the copy of the digits array
1561 m_del(mpz_dig_t, dig, ilen);
1562
Dave Hylandsc4029e52014-04-07 11:19:51 -07001563 if (prefix) {
1564 const char *p = &prefix[strlen(prefix)];
1565 while (p > prefix) {
1566 *s++ = *--p;
1567 }
1568 }
Damien George438c88d2014-02-22 19:25:23 +00001569 if (i->neg != 0) {
1570 *s++ = '-';
1571 }
1572
1573 // reverse string
1574 for (char *u = str, *v = s - 1; u < v; ++u, --v) {
1575 char temp = *u;
1576 *u = *v;
1577 *v = temp;
1578 }
1579
Dave Hylandsc4029e52014-04-07 11:19:51 -07001580 *s = '\0'; // null termination
Damien George438c88d2014-02-22 19:25:23 +00001581
1582 return s - str;
1583}
1584
1585#endif // MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ