blob: e0475d60a1d52b65425dbdfaff171e4712589332 [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
Damien George20653732015-10-01 22:35:06 +0100242 // remove trailing zeros
243 for (--idig; idig >= oidig && *idig == 0; --idig) {
Damien Georgef55cf102014-05-29 15:01:49 +0100244 }
245
246 return idig + 1 - oidig;
247}
248
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200249/* computes i = j | k
250 returns number of digits in i
251 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
252 can have i, j, k pointing to same memory
253*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100254STATIC 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 +0200255 mpz_dig_t *oidig = idig;
256
257 jlen -= klen;
258
259 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
260 *idig = *jdig | *kdig;
261 }
262
263 for (; jlen > 0; --jlen, ++idig, ++jdig) {
264 *idig = *jdig;
265 }
266
267 return idig - oidig;
268}
269
270/* computes i = j ^ k
271 returns number of digits in i
272 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
273 can have i, j, k pointing to same memory
274*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100275STATIC 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 +0200276 mpz_dig_t *oidig = idig;
277
278 jlen -= klen;
279
280 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
281 *idig = *jdig ^ *kdig;
282 }
283
284 for (; jlen > 0; --jlen, ++idig, ++jdig) {
285 *idig = *jdig;
286 }
287
288 return idig - oidig;
289}
290
Damien George438c88d2014-02-22 19:25:23 +0000291/* computes i = i * d1 + d2
292 returns number of digits in i
293 assumes enough memory in i; assumes normalised i; assumes dmul != 0
294*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100295STATIC 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 +0000296 mpz_dig_t *oidig = idig;
297 mpz_dbl_dig_t carry = dadd;
298
299 for (; ilen > 0; --ilen, ++idig) {
Damien George9a21d2e2014-09-06 17:15:34 +0100300 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 +0000301 *idig = carry & DIG_MASK;
302 carry >>= DIG_SIZE;
303 }
304
305 if (carry != 0) {
306 *idig++ = carry;
307 }
308
309 return idig - oidig;
310}
311
312/* computes i = j * k
313 returns number of digits in i
314 assumes enough memory in i; assumes i is zeroed; assumes normalised j, k
315 can have j, k point to same memory
316*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100317STATIC 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 +0000318 mpz_dig_t *oidig = idig;
Damien Georgeafb1cf72014-09-05 20:37:06 +0100319 mp_uint_t ilen = 0;
Damien George438c88d2014-02-22 19:25:23 +0000320
321 for (; klen > 0; --klen, ++idig, ++kdig) {
322 mpz_dig_t *id = idig;
323 mpz_dbl_dig_t carry = 0;
324
Damien Georgeafb1cf72014-09-05 20:37:06 +0100325 mp_uint_t jl = jlen;
Damien George438c88d2014-02-22 19:25:23 +0000326 for (mpz_dig_t *jd = jdig; jl > 0; --jl, ++jd, ++id) {
Damien George9a21d2e2014-09-06 17:15:34 +0100327 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 +0000328 *id = carry & DIG_MASK;
329 carry >>= DIG_SIZE;
330 }
331
332 if (carry != 0) {
333 *id++ = carry;
334 }
335
336 ilen = id - oidig;
337 }
338
339 return ilen;
340}
341
342/* natural_div - quo * den + new_num = old_num (ie num is replaced with rem)
343 assumes den != 0
344 assumes num_dig has enough memory to be extended by 1 digit
345 assumes quo_dig has enough memory (as many digits as num)
346 assumes quo_dig is filled with zeros
347 modifies den_dig memory, but restors it to original state at end
348*/
349
Damien George40f3c022014-07-03 13:25:24 +0100350STATIC 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 +0000351 mpz_dig_t *orig_num_dig = num_dig;
352 mpz_dig_t *orig_quo_dig = quo_dig;
353 mpz_dig_t norm_shift = 0;
354 mpz_dbl_dig_t lead_den_digit;
355
356 // handle simple cases
357 {
Damien George42f3de92014-10-03 17:44:14 +0000358 int cmp = mpn_cmp(num_dig, *num_len, den_dig, den_len);
Damien George438c88d2014-02-22 19:25:23 +0000359 if (cmp == 0) {
360 *num_len = 0;
361 quo_dig[0] = 1;
362 *quo_len = 1;
363 return;
364 } else if (cmp < 0) {
365 // numerator remains the same
366 *quo_len = 0;
367 return;
368 }
369 }
370
371 // count number of leading zeros in leading digit of denominator
372 {
373 mpz_dig_t d = den_dig[den_len - 1];
Damien George9a21d2e2014-09-06 17:15:34 +0100374 while ((d & DIG_MSB) == 0) {
Damien George438c88d2014-02-22 19:25:23 +0000375 d <<= 1;
376 ++norm_shift;
377 }
378 }
379
380 // normalise denomenator (leading bit of leading digit is 1)
381 for (mpz_dig_t *den = den_dig, carry = 0; den < den_dig + den_len; ++den) {
382 mpz_dig_t d = *den;
383 *den = ((d << norm_shift) | carry) & DIG_MASK;
384 carry = d >> (DIG_SIZE - norm_shift);
385 }
386
387 // now need to shift numerator by same amount as denominator
388 // first, increase length of numerator in case we need more room to shift
389 num_dig[*num_len] = 0;
390 ++(*num_len);
391 for (mpz_dig_t *num = num_dig, carry = 0; num < num_dig + *num_len; ++num) {
392 mpz_dig_t n = *num;
393 *num = ((n << norm_shift) | carry) & DIG_MASK;
394 carry = n >> (DIG_SIZE - norm_shift);
395 }
396
397 // cache the leading digit of the denominator
398 lead_den_digit = den_dig[den_len - 1];
399
400 // point num_dig to last digit in numerator
401 num_dig += *num_len - 1;
402
403 // calculate number of digits in quotient
404 *quo_len = *num_len - den_len;
405
406 // point to last digit to store for quotient
407 quo_dig += *quo_len - 1;
408
409 // keep going while we have enough digits to divide
410 while (*num_len > den_len) {
Damien George9a21d2e2014-09-06 17:15:34 +0100411 mpz_dbl_dig_t quo = ((mpz_dbl_dig_t)*num_dig << DIG_SIZE) | num_dig[-1];
Damien George438c88d2014-02-22 19:25:23 +0000412
413 // get approximate quotient
414 quo /= lead_den_digit;
415
Damien George9a21d2e2014-09-06 17:15:34 +0100416 // Multiply quo by den and subtract from num to get remainder.
417 // We have different code here to handle different compile-time
418 // configurations of mpz:
419 //
420 // 1. DIG_SIZE is stricly less than half the number of bits
421 // available in mpz_dbl_dig_t. In this case we can use a
422 // slightly more optimal (in time and space) routine that
423 // uses the extra bits in mpz_dbl_dig_signed_t to store a
424 // sign bit.
425 //
426 // 2. DIG_SIZE is exactly half the number of bits available in
427 // mpz_dbl_dig_t. In this (common) case we need to be careful
428 // not to overflow the borrow variable. And the shifting of
429 // borrow needs some special logic (it's a shift right with
430 // round up).
431
432 if (DIG_SIZE < 8 * sizeof(mpz_dbl_dig_t) / 2) {
Damien George438c88d2014-02-22 19:25:23 +0000433 mpz_dbl_dig_signed_t borrow = 0;
434
435 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
Damien George9a21d2e2014-09-06 17:15:34 +0100436 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 +0000437 *n = borrow & DIG_MASK;
438 borrow >>= DIG_SIZE;
439 }
Damien George9a21d2e2014-09-06 17:15:34 +0100440 borrow += *num_dig; // will overflow if DIG_SIZE >= 8*sizeof(mpz_dbl_dig_t)/2
Damien George438c88d2014-02-22 19:25:23 +0000441 *num_dig = borrow & DIG_MASK;
442 borrow >>= DIG_SIZE;
443
444 // adjust quotient if it is too big
445 for (; borrow != 0; --quo) {
446 mpz_dbl_dig_t carry = 0;
447 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
Damien George9a21d2e2014-09-06 17:15:34 +0100448 carry += (mpz_dbl_dig_t)*n + (mpz_dbl_dig_t)*d;
Damien George438c88d2014-02-22 19:25:23 +0000449 *n = carry & DIG_MASK;
450 carry >>= DIG_SIZE;
451 }
452 carry += *num_dig;
453 *num_dig = carry & DIG_MASK;
454 carry >>= DIG_SIZE;
455
456 borrow += carry;
457 }
Damien George9a21d2e2014-09-06 17:15:34 +0100458 } else { // DIG_SIZE == 8 * sizeof(mpz_dbl_dig_t) / 2
459 mpz_dbl_dig_t borrow = 0;
460
461 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
462 mpz_dbl_dig_t x = (mpz_dbl_dig_t)quo * (mpz_dbl_dig_t)(*d);
463 if (x >= *n || *n - x <= borrow) {
464 borrow += (mpz_dbl_dig_t)x - (mpz_dbl_dig_t)*n;
465 *n = (-borrow) & DIG_MASK;
466 borrow = (borrow >> DIG_SIZE) + ((borrow & DIG_MASK) == 0 ? 0 : 1); // shift-right with round-up
467 } else {
468 *n = ((mpz_dbl_dig_t)*n - (mpz_dbl_dig_t)x - (mpz_dbl_dig_t)borrow) & DIG_MASK;
469 borrow = 0;
470 }
471 }
472 if (borrow >= *num_dig) {
473 borrow -= (mpz_dbl_dig_t)*num_dig;
474 *num_dig = (-borrow) & DIG_MASK;
475 borrow = (borrow >> DIG_SIZE) + ((borrow & DIG_MASK) == 0 ? 0 : 1); // shift-right with round-up
476 } else {
477 *num_dig = (*num_dig - borrow) & DIG_MASK;
478 borrow = 0;
479 }
480
481 // adjust quotient if it is too big
482 for (; borrow != 0; --quo) {
483 mpz_dbl_dig_t carry = 0;
484 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
485 carry += (mpz_dbl_dig_t)*n + (mpz_dbl_dig_t)*d;
486 *n = carry & DIG_MASK;
487 carry >>= DIG_SIZE;
488 }
489 carry += (mpz_dbl_dig_t)*num_dig;
490 *num_dig = carry & DIG_MASK;
491 carry >>= DIG_SIZE;
492
493 //assert(borrow >= carry); // enable this to check the logic
494 borrow -= carry;
495 }
Damien George438c88d2014-02-22 19:25:23 +0000496 }
497
498 // store this digit of the quotient
499 *quo_dig = quo & DIG_MASK;
500 --quo_dig;
501
502 // move down to next digit of numerator
503 --num_dig;
504 --(*num_len);
505 }
506
507 // unnormalise denomenator
508 for (mpz_dig_t *den = den_dig + den_len - 1, carry = 0; den >= den_dig; --den) {
509 mpz_dig_t d = *den;
510 *den = ((d >> norm_shift) | carry) & DIG_MASK;
511 carry = d << (DIG_SIZE - norm_shift);
512 }
513
514 // unnormalise numerator (remainder now)
515 for (mpz_dig_t *num = orig_num_dig + *num_len - 1, carry = 0; num >= orig_num_dig; --num) {
516 mpz_dig_t n = *num;
517 *num = ((n >> norm_shift) | carry) & DIG_MASK;
518 carry = n << (DIG_SIZE - norm_shift);
519 }
520
521 // strip trailing zeros
522
523 while (*quo_len > 0 && orig_quo_dig[*quo_len - 1] == 0) {
524 --(*quo_len);
525 }
526
527 while (*num_len > 0 && orig_num_dig[*num_len - 1] == 0) {
528 --(*num_len);
529 }
530}
531
Damien George06201ff2014-03-01 19:50:50 +0000532#define MIN_ALLOC (2)
Damien George438c88d2014-02-22 19:25:23 +0000533
Damien George1c70cbf2014-08-30 00:38:16 +0100534STATIC const uint8_t log_base2_floor[] = {
Damien George438c88d2014-02-22 19:25:23 +0000535 0,
536 0, 1, 1, 2,
537 2, 2, 2, 3,
538 3, 3, 3, 3,
539 3, 3, 3, 4,
540 4, 4, 4, 4,
541 4, 4, 4, 4,
542 4, 4, 4, 4,
543 4, 4, 4, 5
544};
545
Damien George438c88d2014-02-22 19:25:23 +0000546void mpz_init_zero(mpz_t *z) {
Damien George438c88d2014-02-22 19:25:23 +0000547 z->neg = 0;
Damien George06201ff2014-03-01 19:50:50 +0000548 z->fixed_dig = 0;
549 z->alloc = 0;
Damien George438c88d2014-02-22 19:25:23 +0000550 z->len = 0;
551 z->dig = NULL;
552}
553
Damien George40f3c022014-07-03 13:25:24 +0100554void mpz_init_from_int(mpz_t *z, mp_int_t val) {
Damien George438c88d2014-02-22 19:25:23 +0000555 mpz_init_zero(z);
556 mpz_set_from_int(z, val);
557}
558
Damien Georgeafb1cf72014-09-05 20:37:06 +0100559void 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 +0000560 z->neg = 0;
561 z->fixed_dig = 1;
562 z->alloc = alloc;
563 z->len = 0;
564 z->dig = dig;
565 mpz_set_from_int(z, val);
566}
567
Damien George438c88d2014-02-22 19:25:23 +0000568void mpz_deinit(mpz_t *z) {
Damien George06201ff2014-03-01 19:50:50 +0000569 if (z != NULL && !z->fixed_dig) {
Damien George438c88d2014-02-22 19:25:23 +0000570 m_del(mpz_dig_t, z->dig, z->alloc);
571 }
572}
573
Damien George848dd0e2015-03-12 15:59:40 +0000574#if 0
575these functions are unused
576
Damien George438c88d2014-02-22 19:25:23 +0000577mpz_t *mpz_zero(void) {
578 mpz_t *z = m_new_obj(mpz_t);
579 mpz_init_zero(z);
580 return z;
581}
582
Damien George40f3c022014-07-03 13:25:24 +0100583mpz_t *mpz_from_int(mp_int_t val) {
Damien George438c88d2014-02-22 19:25:23 +0000584 mpz_t *z = mpz_zero();
585 mpz_set_from_int(z, val);
586 return z;
587}
588
Damien George95307432014-09-10 22:10:33 +0100589mpz_t *mpz_from_ll(long long val, bool is_signed) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000590 mpz_t *z = mpz_zero();
Damien George95307432014-09-10 22:10:33 +0100591 mpz_set_from_ll(z, val, is_signed);
Damien Georgebb4a43f2014-03-12 15:36:06 +0000592 return z;
593}
594
David Steinberg6e0b6d02015-01-02 12:39:22 +0000595#if MICROPY_PY_BUILTINS_FLOAT
596mpz_t *mpz_from_float(mp_float_t val) {
597 mpz_t *z = mpz_zero();
598 mpz_set_from_float(z, val);
599 return z;
600}
601#endif
602
Damien Georgeafb1cf72014-09-05 20:37:06 +0100603mpz_t *mpz_from_str(const char *str, mp_uint_t len, bool neg, mp_uint_t base) {
Damien George438c88d2014-02-22 19:25:23 +0000604 mpz_t *z = mpz_zero();
605 mpz_set_from_str(z, str, len, neg, base);
606 return z;
607}
Damien George848dd0e2015-03-12 15:59:40 +0000608#endif
Damien George438c88d2014-02-22 19:25:23 +0000609
Damien George848dd0e2015-03-12 15:59:40 +0000610STATIC void mpz_free(mpz_t *z) {
Damien George438c88d2014-02-22 19:25:23 +0000611 if (z != NULL) {
612 m_del(mpz_dig_t, z->dig, z->alloc);
613 m_del_obj(mpz_t, z);
614 }
615}
616
Damien Georgeafb1cf72014-09-05 20:37:06 +0100617STATIC void mpz_need_dig(mpz_t *z, mp_uint_t need) {
Damien George438c88d2014-02-22 19:25:23 +0000618 if (need < MIN_ALLOC) {
Damien George06201ff2014-03-01 19:50:50 +0000619 need = MIN_ALLOC;
Damien George438c88d2014-02-22 19:25:23 +0000620 }
621
Damien George06201ff2014-03-01 19:50:50 +0000622 if (z->dig == NULL || z->alloc < need) {
623 if (z->fixed_dig) {
624 // cannot reallocate fixed buffers
625 assert(0);
626 return;
627 }
628 z->dig = m_renew(mpz_dig_t, z->dig, z->alloc, need);
629 z->alloc = need;
Damien George438c88d2014-02-22 19:25:23 +0000630 }
631}
632
Damien George848dd0e2015-03-12 15:59:40 +0000633STATIC mpz_t *mpz_clone(const mpz_t *src) {
Damien George438c88d2014-02-22 19:25:23 +0000634 mpz_t *z = m_new_obj(mpz_t);
Damien George438c88d2014-02-22 19:25:23 +0000635 z->neg = src->neg;
Damien George06201ff2014-03-01 19:50:50 +0000636 z->fixed_dig = 0;
637 z->alloc = src->alloc;
Damien George438c88d2014-02-22 19:25:23 +0000638 z->len = src->len;
639 if (src->dig == NULL) {
640 z->dig = NULL;
641 } else {
642 z->dig = m_new(mpz_dig_t, z->alloc);
643 memcpy(z->dig, src->dig, src->alloc * sizeof(mpz_dig_t));
644 }
645 return z;
646}
647
Damien George06201ff2014-03-01 19:50:50 +0000648/* sets dest = src
649 can have dest, src the same
650*/
Damien George438c88d2014-02-22 19:25:23 +0000651void mpz_set(mpz_t *dest, const mpz_t *src) {
652 mpz_need_dig(dest, src->len);
653 dest->neg = src->neg;
654 dest->len = src->len;
655 memcpy(dest->dig, src->dig, src->len * sizeof(mpz_dig_t));
656}
657
Damien George40f3c022014-07-03 13:25:24 +0100658void mpz_set_from_int(mpz_t *z, mp_int_t val) {
Damien George58056b02015-01-09 20:58:58 +0000659 if (val == 0) {
660 z->len = 0;
661 return;
662 }
663
Damien George06201ff2014-03-01 19:50:50 +0000664 mpz_need_dig(z, MPZ_NUM_DIG_FOR_INT);
Damien George438c88d2014-02-22 19:25:23 +0000665
Damien George40f3c022014-07-03 13:25:24 +0100666 mp_uint_t uval;
Damien George438c88d2014-02-22 19:25:23 +0000667 if (val < 0) {
668 z->neg = 1;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000669 uval = -val;
Damien George438c88d2014-02-22 19:25:23 +0000670 } else {
671 z->neg = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000672 uval = val;
Damien George438c88d2014-02-22 19:25:23 +0000673 }
674
675 z->len = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000676 while (uval > 0) {
677 z->dig[z->len++] = uval & DIG_MASK;
678 uval >>= DIG_SIZE;
679 }
680}
681
Damien George95307432014-09-10 22:10:33 +0100682void mpz_set_from_ll(mpz_t *z, long long val, bool is_signed) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000683 mpz_need_dig(z, MPZ_NUM_DIG_FOR_LL);
684
685 unsigned long long uval;
Damien George95307432014-09-10 22:10:33 +0100686 if (is_signed && val < 0) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000687 z->neg = 1;
688 uval = -val;
689 } else {
690 z->neg = 0;
691 uval = val;
692 }
693
694 z->len = 0;
695 while (uval > 0) {
696 z->dig[z->len++] = uval & DIG_MASK;
697 uval >>= DIG_SIZE;
Damien George438c88d2014-02-22 19:25:23 +0000698 }
699}
700
David Steinberg6e0b6d02015-01-02 12:39:22 +0000701#if MICROPY_PY_BUILTINS_FLOAT
702void mpz_set_from_float(mpz_t *z, mp_float_t src) {
703#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
David Steinberg6e0b6d02015-01-02 12:39:22 +0000704typedef uint64_t mp_float_int_t;
David Steinbergc585ad12015-01-13 15:19:37 +0000705#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
David Steinberg6e0b6d02015-01-02 12:39:22 +0000706typedef uint32_t mp_float_int_t;
707#endif
708 union {
709 mp_float_t f;
David Steinbergc585ad12015-01-13 15:19:37 +0000710 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 +0000711 } u = {src};
712
713 z->neg = u.p.sgn;
714 if (u.p.exp == 0) {
715 // value == 0 || value < 1
Damien George58056b02015-01-09 20:58:58 +0000716 mpz_set_from_int(z, 0);
David Steinbergc585ad12015-01-13 15:19:37 +0000717 } else if (u.p.exp == ((1 << MP_FLOAT_EXP_BITS) - 1)) {
Damien George6fd4b362015-01-02 23:04:09 +0000718 // u.p.frc == 0 indicates inf, else NaN
719 // should be handled by caller
Damien George58056b02015-01-09 20:58:58 +0000720 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000721 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000722 const int adj_exp = (int)u.p.exp - MP_FLOAT_EXP_BIAS;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000723 if (adj_exp < 0) {
724 // value < 1 , truncates to 0
Damien George58056b02015-01-09 20:58:58 +0000725 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000726 } else if (adj_exp == 0) {
727 // 1 <= value < 2 , so truncates to 1
Damien George58056b02015-01-09 20:58:58 +0000728 mpz_set_from_int(z, 1);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000729 } else {
730 // 2 <= value
731 const int dig_cnt = (adj_exp + 1 + (DIG_SIZE - 1)) / DIG_SIZE;
732 const unsigned int rem = adj_exp % DIG_SIZE;
733 int dig_ind, shft;
David Steinbergc585ad12015-01-13 15:19:37 +0000734 mp_float_int_t frc = u.p.frc | ((mp_float_int_t)1 << MP_FLOAT_FRAC_BITS);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000735
David Steinbergc585ad12015-01-13 15:19:37 +0000736 if (adj_exp < MP_FLOAT_FRAC_BITS) {
David Steinberg6e0b6d02015-01-02 12:39:22 +0000737 shft = 0;
738 dig_ind = 0;
David Steinbergc585ad12015-01-13 15:19:37 +0000739 frc >>= MP_FLOAT_FRAC_BITS - adj_exp;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000740 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000741 shft = (rem - MP_FLOAT_FRAC_BITS) % DIG_SIZE;
742 dig_ind = (adj_exp - MP_FLOAT_FRAC_BITS) / DIG_SIZE;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000743 }
744 mpz_need_dig(z, dig_cnt);
745 z->len = dig_cnt;
746 if (dig_ind != 0) {
747 memset(z->dig, 0, dig_ind * sizeof(mpz_dig_t));
748 }
749 if (shft != 0) {
750 z->dig[dig_ind++] = (frc << shft) & DIG_MASK;
751 frc >>= DIG_SIZE - shft;
752 }
David Steinberg8d427b72015-01-13 15:20:32 +0000753#if DIG_SIZE < (MP_FLOAT_FRAC_BITS + 1)
David Steinberg6e0b6d02015-01-02 12:39:22 +0000754 while (dig_ind != dig_cnt) {
755 z->dig[dig_ind++] = frc & DIG_MASK;
756 frc >>= DIG_SIZE;
757 }
David Steinberg8d427b72015-01-13 15:20:32 +0000758#else
759 if (dig_ind != dig_cnt) {
760 z->dig[dig_ind] = frc;
761 }
762#endif
David Steinberg6e0b6d02015-01-02 12:39:22 +0000763 }
764 }
765}
David Steinberg6e0b6d02015-01-02 12:39:22 +0000766#endif
767
Damien George438c88d2014-02-22 19:25:23 +0000768// returns number of bytes from str that were processed
Damien Georgeafb1cf72014-09-05 20:37:06 +0100769mp_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 +0000770 assert(base < 36);
771
772 const char *cur = str;
773 const char *top = str + len;
774
775 mpz_need_dig(z, len * 8 / DIG_SIZE + 1);
776
777 if (neg) {
778 z->neg = 1;
779 } else {
780 z->neg = 0;
781 }
782
783 z->len = 0;
784 for (; cur < top; ++cur) { // XXX UTF8 next char
Damien Georgeafb1cf72014-09-05 20:37:06 +0100785 //mp_uint_t v = char_to_numeric(cur#); // XXX UTF8 get char
786 mp_uint_t v = *cur;
Damien George438c88d2014-02-22 19:25:23 +0000787 if ('0' <= v && v <= '9') {
788 v -= '0';
789 } else if ('A' <= v && v <= 'Z') {
790 v -= 'A' - 10;
791 } else if ('a' <= v && v <= 'z') {
792 v -= 'a' - 10;
793 } else {
794 break;
795 }
796 if (v >= base) {
797 break;
798 }
799 z->len = mpn_mul_dig_add_dig(z->dig, z->len, base, v);
800 }
801
802 return cur - str;
803}
804
805bool mpz_is_zero(const mpz_t *z) {
806 return z->len == 0;
807}
808
Damien Georgea2e38382015-03-02 12:58:06 +0000809#if 0
810these functions are unused
811
Damien George438c88d2014-02-22 19:25:23 +0000812bool mpz_is_pos(const mpz_t *z) {
813 return z->len > 0 && z->neg == 0;
814}
815
816bool mpz_is_neg(const mpz_t *z) {
817 return z->len > 0 && z->neg != 0;
818}
819
820bool mpz_is_odd(const mpz_t *z) {
821 return z->len > 0 && (z->dig[0] & 1) != 0;
822}
823
824bool mpz_is_even(const mpz_t *z) {
825 return z->len == 0 || (z->dig[0] & 1) == 0;
826}
Damien Georgea2e38382015-03-02 12:58:06 +0000827#endif
Damien George438c88d2014-02-22 19:25:23 +0000828
Damien George42f3de92014-10-03 17:44:14 +0000829int mpz_cmp(const mpz_t *z1, const mpz_t *z2) {
Damien Georgea9dc9b82015-01-27 17:47:38 +0000830 // to catch comparison of -0 with +0
831 if (z1->len == 0 && z2->len == 0) {
832 return 0;
833 }
Damien George42f3de92014-10-03 17:44:14 +0000834 int cmp = (int)z2->neg - (int)z1->neg;
Damien George438c88d2014-02-22 19:25:23 +0000835 if (cmp != 0) {
836 return cmp;
837 }
838 cmp = mpn_cmp(z1->dig, z1->len, z2->dig, z2->len);
839 if (z1->neg != 0) {
840 cmp = -cmp;
841 }
842 return cmp;
843}
844
Damien George06201ff2014-03-01 19:50:50 +0000845#if 0
846// obsolete
847// compares mpz with an integer that fits within DIG_SIZE bits
Damien Georgeafb1cf72014-09-05 20:37:06 +0100848mp_int_t mpz_cmp_sml_int(const mpz_t *z, mp_int_t sml_int) {
849 mp_int_t cmp;
Damien George438c88d2014-02-22 19:25:23 +0000850 if (z->neg == 0) {
851 if (sml_int < 0) return 1;
852 if (sml_int == 0) {
853 if (z->len == 0) return 0;
854 return 1;
855 }
856 if (z->len == 0) return -1;
857 assert(sml_int < (1 << DIG_SIZE));
858 if (z->len != 1) return 1;
859 cmp = z->dig[0] - sml_int;
860 } else {
861 if (sml_int > 0) return -1;
862 if (sml_int == 0) {
863 if (z->len == 0) return 0;
864 return -1;
865 }
866 if (z->len == 0) return 1;
867 assert(sml_int > -(1 << DIG_SIZE));
868 if (z->len != 1) return -1;
869 cmp = -z->dig[0] - sml_int;
870 }
871 if (cmp < 0) return -1;
872 if (cmp > 0) return 1;
873 return 0;
874}
Damien George06201ff2014-03-01 19:50:50 +0000875#endif
Damien George438c88d2014-02-22 19:25:23 +0000876
Damien George438c88d2014-02-22 19:25:23 +0000877#if 0
878these functions are unused
879
880/* returns abs(z)
881*/
882mpz_t *mpz_abs(const mpz_t *z) {
883 mpz_t *z2 = mpz_clone(z);
884 z2->neg = 0;
885 return z2;
886}
887
888/* returns -z
889*/
890mpz_t *mpz_neg(const mpz_t *z) {
891 mpz_t *z2 = mpz_clone(z);
892 z2->neg = 1 - z2->neg;
893 return z2;
894}
895
896/* returns lhs + rhs
897 can have lhs, rhs the same
898*/
899mpz_t *mpz_add(const mpz_t *lhs, const mpz_t *rhs) {
900 mpz_t *z = mpz_zero();
901 mpz_add_inpl(z, lhs, rhs);
902 return z;
903}
904
905/* returns lhs - rhs
906 can have lhs, rhs the same
907*/
908mpz_t *mpz_sub(const mpz_t *lhs, const mpz_t *rhs) {
909 mpz_t *z = mpz_zero();
910 mpz_sub_inpl(z, lhs, rhs);
911 return z;
912}
913
914/* returns lhs * rhs
915 can have lhs, rhs the same
916*/
917mpz_t *mpz_mul(const mpz_t *lhs, const mpz_t *rhs) {
918 mpz_t *z = mpz_zero();
919 mpz_mul_inpl(z, lhs, rhs);
920 return z;
921}
922
923/* returns lhs ** rhs
924 can have lhs, rhs the same
925*/
926mpz_t *mpz_pow(const mpz_t *lhs, const mpz_t *rhs) {
927 mpz_t *z = mpz_zero();
928 mpz_pow_inpl(z, lhs, rhs);
929 return z;
930}
Damien Georgea2e38382015-03-02 12:58:06 +0000931
932/* computes new integers in quo and rem such that:
933 quo * rhs + rem = lhs
934 0 <= rem < rhs
935 can have lhs, rhs the same
936*/
937void mpz_divmod(const mpz_t *lhs, const mpz_t *rhs, mpz_t **quo, mpz_t **rem) {
938 *quo = mpz_zero();
939 *rem = mpz_zero();
940 mpz_divmod_inpl(*quo, *rem, lhs, rhs);
941}
Damien George438c88d2014-02-22 19:25:23 +0000942#endif
943
944/* computes dest = abs(z)
945 can have dest, z the same
946*/
947void mpz_abs_inpl(mpz_t *dest, const mpz_t *z) {
948 if (dest != z) {
949 mpz_set(dest, z);
950 }
951 dest->neg = 0;
952}
953
954/* computes dest = -z
955 can have dest, z the same
956*/
957void mpz_neg_inpl(mpz_t *dest, const mpz_t *z) {
958 if (dest != z) {
959 mpz_set(dest, z);
960 }
961 dest->neg = 1 - dest->neg;
962}
963
Damien George06201ff2014-03-01 19:50:50 +0000964/* computes dest = ~z (= -z - 1)
965 can have dest, z the same
966*/
967void mpz_not_inpl(mpz_t *dest, const mpz_t *z) {
968 if (dest != z) {
969 mpz_set(dest, z);
970 }
Damien Georgee0ac1942014-12-31 19:35:01 +0000971 if (dest->len == 0) {
972 mpz_need_dig(dest, 1);
973 dest->dig[0] = 1;
974 dest->len = 1;
975 dest->neg = 1;
976 } else if (dest->neg) {
Damien George06201ff2014-03-01 19:50:50 +0000977 dest->neg = 0;
978 mpz_dig_t k = 1;
979 dest->len = mpn_sub(dest->dig, dest->dig, dest->len, &k, 1);
980 } else {
Damien Georgee0ac1942014-12-31 19:35:01 +0000981 mpz_need_dig(dest, dest->len + 1);
Damien George06201ff2014-03-01 19:50:50 +0000982 mpz_dig_t k = 1;
983 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &k, 1);
984 dest->neg = 1;
985 }
986}
987
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000988/* computes dest = lhs << rhs
989 can have dest, lhs the same
990*/
Damien George2f4e8512015-10-01 18:01:37 +0100991void mpz_shl_inpl(mpz_t *dest, const mpz_t *lhs, mp_uint_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +0000992 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000993 mpz_set(dest, lhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000994 } else {
Damien George06201ff2014-03-01 19:50:50 +0000995 mpz_need_dig(dest, lhs->len + (rhs + DIG_SIZE - 1) / DIG_SIZE);
996 dest->len = mpn_shl(dest->dig, lhs->dig, lhs->len, rhs);
997 dest->neg = lhs->neg;
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000998 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000999}
1000
1001/* computes dest = lhs >> rhs
1002 can have dest, lhs the same
1003*/
Damien George2f4e8512015-10-01 18:01:37 +01001004void mpz_shr_inpl(mpz_t *dest, const mpz_t *lhs, mp_uint_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +00001005 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001006 mpz_set(dest, lhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001007 } else {
Damien George06201ff2014-03-01 19:50:50 +00001008 mpz_need_dig(dest, lhs->len);
1009 dest->len = mpn_shr(dest->dig, lhs->dig, lhs->len, rhs);
1010 dest->neg = lhs->neg;
1011 if (dest->neg) {
1012 // arithmetic shift right, rounding to negative infinity
Damien Georgeafb1cf72014-09-05 20:37:06 +01001013 mp_uint_t n_whole = rhs / DIG_SIZE;
1014 mp_uint_t n_part = rhs % DIG_SIZE;
Damien George06201ff2014-03-01 19:50:50 +00001015 mpz_dig_t round_up = 0;
Damien Georgeafb1cf72014-09-05 20:37:06 +01001016 for (mp_uint_t i = 0; i < lhs->len && i < n_whole; i++) {
Damien George06201ff2014-03-01 19:50:50 +00001017 if (lhs->dig[i] != 0) {
1018 round_up = 1;
1019 break;
1020 }
1021 }
1022 if (n_whole < lhs->len && (lhs->dig[n_whole] & ((1 << n_part) - 1)) != 0) {
1023 round_up = 1;
1024 }
1025 if (round_up) {
Damien Georgee0ac1942014-12-31 19:35:01 +00001026 if (dest->len == 0) {
1027 // dest == 0, so need to add 1 by hand (answer will be -1)
1028 dest->dig[0] = 1;
1029 dest->len = 1;
1030 } else {
1031 // dest > 0, so can use mpn_add to add 1
1032 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &round_up, 1);
1033 }
Damien George06201ff2014-03-01 19:50:50 +00001034 }
1035 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001036 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001037}
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001038
Damien George438c88d2014-02-22 19:25:23 +00001039/* computes dest = lhs + rhs
1040 can have dest, lhs, rhs the same
1041*/
1042void mpz_add_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1043 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1044 const mpz_t *temp = lhs;
1045 lhs = rhs;
1046 rhs = temp;
1047 }
1048
1049 if (lhs->neg == rhs->neg) {
1050 mpz_need_dig(dest, lhs->len + 1);
1051 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1052 } else {
1053 mpz_need_dig(dest, lhs->len);
1054 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1055 }
1056
1057 dest->neg = lhs->neg;
1058}
1059
1060/* computes dest = lhs - rhs
1061 can have dest, lhs, rhs the same
1062*/
1063void mpz_sub_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1064 bool neg = false;
1065
1066 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1067 const mpz_t *temp = lhs;
1068 lhs = rhs;
1069 rhs = temp;
1070 neg = true;
1071 }
1072
1073 if (lhs->neg != rhs->neg) {
1074 mpz_need_dig(dest, lhs->len + 1);
1075 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1076 } else {
1077 mpz_need_dig(dest, lhs->len);
1078 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1079 }
1080
1081 if (neg) {
1082 dest->neg = 1 - lhs->neg;
1083 } else {
1084 dest->neg = lhs->neg;
1085 }
1086}
1087
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001088/* computes dest = lhs & rhs
1089 can have dest, lhs, rhs the same
1090*/
1091void mpz_and_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001092 if (lhs->neg == rhs->neg) {
Damien Georgef55cf102014-05-29 15:01:49 +01001093 if (lhs->neg == 0) {
1094 // make sure lhs has the most digits
1095 if (lhs->len < rhs->len) {
1096 const mpz_t *temp = lhs;
1097 lhs = rhs;
1098 rhs = temp;
1099 }
1100 // do the and'ing
1101 mpz_need_dig(dest, rhs->len);
Damien Georgeff8dd3f2015-01-20 12:47:20 +00001102 dest->len = mpn_and(dest->dig, lhs->dig, rhs->dig, rhs->len);
Damien Georgef55cf102014-05-29 15:01:49 +01001103 dest->neg = 0;
1104 } else {
1105 // TODO both args are negative
Damien George4c02e542015-10-01 17:57:36 +01001106 mp_not_implemented("bignum and with negative args");
Damien Georgef55cf102014-05-29 15:01:49 +01001107 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001108 } else {
Damien Georgef55cf102014-05-29 15:01:49 +01001109 // args have different sign
1110 // make sure lhs is the positive arg
1111 if (rhs->neg == 0) {
1112 const mpz_t *temp = lhs;
1113 lhs = rhs;
1114 rhs = temp;
1115 }
1116 mpz_need_dig(dest, lhs->len + 1);
1117 dest->len = mpn_and_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1118 assert(dest->len <= dest->alloc);
1119 dest->neg = 0;
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001120 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001121}
1122
1123/* computes dest = lhs | rhs
1124 can have dest, lhs, rhs the same
1125*/
1126void mpz_or_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1127 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1128 const mpz_t *temp = lhs;
1129 lhs = rhs;
1130 rhs = temp;
1131 }
1132
1133 if (lhs->neg == rhs->neg) {
1134 mpz_need_dig(dest, lhs->len);
1135 dest->len = mpn_or(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1136 } else {
1137 mpz_need_dig(dest, lhs->len);
1138 // TODO
Damien George4c02e542015-10-01 17:57:36 +01001139 mp_not_implemented("bignum or with negative args");
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001140// dest->len = mpn_or_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1141 }
1142
1143 dest->neg = lhs->neg;
1144}
1145
1146/* computes dest = lhs ^ rhs
1147 can have dest, lhs, rhs the same
1148*/
1149void mpz_xor_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1150 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1151 const mpz_t *temp = lhs;
1152 lhs = rhs;
1153 rhs = temp;
1154 }
1155
1156 if (lhs->neg == rhs->neg) {
1157 mpz_need_dig(dest, lhs->len);
1158 dest->len = mpn_xor(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1159 } else {
1160 mpz_need_dig(dest, lhs->len);
1161 // TODO
Damien George4c02e542015-10-01 17:57:36 +01001162 mp_not_implemented("bignum xor with negative args");
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001163// dest->len = mpn_xor_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1164 }
1165
1166 dest->neg = 0;
1167}
1168
Damien George438c88d2014-02-22 19:25:23 +00001169/* computes dest = lhs * rhs
1170 can have dest, lhs, rhs the same
1171*/
Damien George4dea9222015-04-09 15:29:54 +00001172void mpz_mul_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Damien George438c88d2014-02-22 19:25:23 +00001173 if (lhs->len == 0 || rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001174 mpz_set_from_int(dest, 0);
1175 return;
Damien George438c88d2014-02-22 19:25:23 +00001176 }
1177
1178 mpz_t *temp = NULL;
1179 if (lhs == dest) {
1180 lhs = temp = mpz_clone(lhs);
1181 if (rhs == dest) {
1182 rhs = lhs;
1183 }
1184 } else if (rhs == dest) {
1185 rhs = temp = mpz_clone(rhs);
1186 }
1187
1188 mpz_need_dig(dest, lhs->len + rhs->len); // min mem l+r-1, max mem l+r
1189 memset(dest->dig, 0, dest->alloc * sizeof(mpz_dig_t));
1190 dest->len = mpn_mul(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1191
1192 if (lhs->neg == rhs->neg) {
1193 dest->neg = 0;
1194 } else {
1195 dest->neg = 1;
1196 }
1197
1198 mpz_free(temp);
1199}
1200
1201/* computes dest = lhs ** rhs
1202 can have dest, lhs, rhs the same
1203*/
1204void mpz_pow_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1205 if (lhs->len == 0 || rhs->neg != 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001206 mpz_set_from_int(dest, 0);
1207 return;
Damien George438c88d2014-02-22 19:25:23 +00001208 }
1209
1210 if (rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001211 mpz_set_from_int(dest, 1);
1212 return;
Damien George438c88d2014-02-22 19:25:23 +00001213 }
1214
1215 mpz_t *x = mpz_clone(lhs);
1216 mpz_t *n = mpz_clone(rhs);
1217
1218 mpz_set_from_int(dest, 1);
1219
1220 while (n->len > 0) {
Damien Georgea2e38382015-03-02 12:58:06 +00001221 if ((n->dig[0] & 1) != 0) {
Damien George438c88d2014-02-22 19:25:23 +00001222 mpz_mul_inpl(dest, dest, x);
1223 }
Damien George438c88d2014-02-22 19:25:23 +00001224 n->len = mpn_shr(n->dig, n->dig, n->len, 1);
Damien George5bf565e2014-04-04 00:16:32 +01001225 if (n->len == 0) {
1226 break;
1227 }
1228 mpz_mul_inpl(x, x, x);
Damien George438c88d2014-02-22 19:25:23 +00001229 }
1230
1231 mpz_free(x);
1232 mpz_free(n);
1233}
1234
Damien Georgea2e38382015-03-02 12:58:06 +00001235#if 0
1236these functions are unused
1237
Damien George438c88d2014-02-22 19:25:23 +00001238/* computes gcd(z1, z2)
1239 based on Knuth's modified gcd algorithm (I think?)
1240 gcd(z1, z2) >= 0
1241 gcd(0, 0) = 0
1242 gcd(z, 0) = abs(z)
1243*/
1244mpz_t *mpz_gcd(const mpz_t *z1, const mpz_t *z2) {
1245 if (z1->len == 0) {
1246 mpz_t *a = mpz_clone(z2);
1247 a->neg = 0;
1248 return a;
1249 } else if (z2->len == 0) {
1250 mpz_t *a = mpz_clone(z1);
1251 a->neg = 0;
1252 return a;
1253 }
1254
1255 mpz_t *a = mpz_clone(z1);
1256 mpz_t *b = mpz_clone(z2);
1257 mpz_t c; mpz_init_zero(&c);
1258 a->neg = 0;
1259 b->neg = 0;
1260
1261 for (;;) {
1262 if (mpz_cmp(a, b) < 0) {
1263 if (a->len == 0) {
1264 mpz_free(a);
1265 mpz_deinit(&c);
1266 return b;
1267 }
1268 mpz_t *t = a; a = b; b = t;
1269 }
1270 if (!(b->len >= 2 || (b->len == 1 && b->dig[0] > 1))) { // compute b > 0; could be mpz_cmp_small_int(b, 1) > 0
1271 break;
1272 }
1273 mpz_set(&c, b);
1274 do {
1275 mpz_add_inpl(&c, &c, &c);
1276 } while (mpz_cmp(&c, a) <= 0);
1277 c.len = mpn_shr(c.dig, c.dig, c.len, 1);
1278 mpz_sub_inpl(a, a, &c);
1279 }
1280
1281 mpz_deinit(&c);
1282
1283 if (b->len == 1 && b->dig[0] == 1) { // compute b == 1; could be mpz_cmp_small_int(b, 1) == 0
1284 mpz_free(a);
1285 return b;
1286 } else {
1287 mpz_free(b);
1288 return a;
1289 }
1290}
1291
1292/* computes lcm(z1, z2)
1293 = abs(z1) / gcd(z1, z2) * abs(z2)
1294 lcm(z1, z1) >= 0
1295 lcm(0, 0) = 0
1296 lcm(z, 0) = 0
1297*/
Damien George5d9b8162014-08-07 14:27:48 +00001298mpz_t *mpz_lcm(const mpz_t *z1, const mpz_t *z2) {
1299 if (z1->len == 0 || z2->len == 0) {
1300 return mpz_zero();
stijn01d6be42014-05-05 12:18:27 +02001301 }
Damien George438c88d2014-02-22 19:25:23 +00001302
1303 mpz_t *gcd = mpz_gcd(z1, z2);
1304 mpz_t *quo = mpz_zero();
1305 mpz_t *rem = mpz_zero();
1306 mpz_divmod_inpl(quo, rem, z1, gcd);
1307 mpz_mul_inpl(rem, quo, z2);
1308 mpz_free(gcd);
1309 mpz_free(quo);
1310 rem->neg = 0;
1311 return rem;
1312}
Damien Georgea2e38382015-03-02 12:58:06 +00001313#endif
Damien George438c88d2014-02-22 19:25:23 +00001314
1315/* computes new integers in quo and rem such that:
1316 quo * rhs + rem = lhs
1317 0 <= rem < rhs
1318 can have lhs, rhs the same
1319*/
1320void mpz_divmod_inpl(mpz_t *dest_quo, mpz_t *dest_rem, const mpz_t *lhs, const mpz_t *rhs) {
1321 if (rhs->len == 0) {
1322 mpz_set_from_int(dest_quo, 0);
1323 mpz_set_from_int(dest_rem, 0);
1324 return;
1325 }
1326
1327 mpz_need_dig(dest_quo, lhs->len + 1); // +1 necessary?
1328 memset(dest_quo->dig, 0, (lhs->len + 1) * sizeof(mpz_dig_t));
1329 dest_quo->len = 0;
1330 mpz_need_dig(dest_rem, lhs->len + 1); // +1 necessary?
1331 mpz_set(dest_rem, lhs);
1332 //rhs->dig[rhs->len] = 0;
1333 mpn_div(dest_rem->dig, &dest_rem->len, rhs->dig, rhs->len, dest_quo->dig, &dest_quo->len);
1334
1335 if (lhs->neg != rhs->neg) {
1336 dest_quo->neg = 1;
1337 }
1338}
1339
1340#if 0
1341these functions are unused
1342
1343/* computes floor(lhs / rhs)
1344 can have lhs, rhs the same
1345*/
1346mpz_t *mpz_div(const mpz_t *lhs, const mpz_t *rhs) {
1347 mpz_t *quo = mpz_zero();
1348 mpz_t rem; mpz_init_zero(&rem);
1349 mpz_divmod_inpl(quo, &rem, lhs, rhs);
1350 mpz_deinit(&rem);
1351 return quo;
1352}
1353
1354/* computes lhs % rhs ( >= 0)
1355 can have lhs, rhs the same
1356*/
1357mpz_t *mpz_mod(const mpz_t *lhs, const mpz_t *rhs) {
1358 mpz_t quo; mpz_init_zero(&quo);
1359 mpz_t *rem = mpz_zero();
1360 mpz_divmod_inpl(&quo, rem, lhs, rhs);
1361 mpz_deinit(&quo);
1362 return rem;
1363}
1364#endif
1365
Damien Georgeffe911d2014-07-24 14:21:37 +01001366// must return actual int value if it fits in mp_int_t
1367mp_int_t mpz_hash(const mpz_t *z) {
1368 mp_int_t val = 0;
1369 mpz_dig_t *d = z->dig + z->len;
1370
Damien George58056b02015-01-09 20:58:58 +00001371 while (d-- > z->dig) {
Damien Georgeffe911d2014-07-24 14:21:37 +01001372 val = (val << DIG_SIZE) | *d;
1373 }
1374
1375 if (z->neg != 0) {
1376 val = -val;
1377 }
1378
1379 return val;
1380}
1381
Damien George40f3c022014-07-03 13:25:24 +01001382bool mpz_as_int_checked(const mpz_t *i, mp_int_t *value) {
Damien George963a5a32015-01-16 17:47:07 +00001383 mp_uint_t val = 0;
Damien George8270e382014-04-03 11:00:54 +00001384 mpz_dig_t *d = i->dig + i->len;
1385
Damien George58056b02015-01-09 20:58:58 +00001386 while (d-- > i->dig) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001387 if (val > (~(WORD_MSBIT_HIGH) >> DIG_SIZE)) {
1388 // will overflow
Damien George8270e382014-04-03 11:00:54 +00001389 return false;
1390 }
Damien Georgec9aa58e2014-07-31 13:41:43 +00001391 val = (val << DIG_SIZE) | *d;
Damien George8270e382014-04-03 11:00:54 +00001392 }
1393
1394 if (i->neg != 0) {
1395 val = -val;
1396 }
1397
1398 *value = val;
1399 return true;
1400}
1401
Damien Georgec9aa58e2014-07-31 13:41:43 +00001402bool mpz_as_uint_checked(const mpz_t *i, mp_uint_t *value) {
1403 if (i->neg != 0) {
1404 // can't represent signed values
1405 return false;
1406 }
1407
1408 mp_uint_t val = 0;
1409 mpz_dig_t *d = i->dig + i->len;
1410
Damien George58056b02015-01-09 20:58:58 +00001411 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001412 if (val > (~(WORD_MSBIT_HIGH) >> (DIG_SIZE - 1))) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001413 // will overflow
1414 return false;
1415 }
1416 val = (val << DIG_SIZE) | *d;
1417 }
1418
1419 *value = val;
1420 return true;
1421}
1422
Damien George271d18e2015-04-25 23:16:39 +01001423// writes at most len bytes to buf (so buf should be zeroed before calling)
1424void mpz_as_bytes(const mpz_t *z, bool big_endian, mp_uint_t len, byte *buf) {
1425 byte *b = buf;
1426 if (big_endian) {
1427 b += len;
1428 }
1429 mpz_dig_t *zdig = z->dig;
1430 int bits = 0;
1431 mpz_dbl_dig_t d = 0;
1432 mpz_dbl_dig_t carry = 1;
1433 for (mp_uint_t zlen = z->len; zlen > 0; --zlen) {
1434 bits += DIG_SIZE;
1435 d = (d << DIG_SIZE) | *zdig++;
1436 for (; bits >= 8; bits -= 8, d >>= 8) {
1437 mpz_dig_t val = d;
1438 if (z->neg) {
Damien George94729072015-04-25 23:51:14 +01001439 val = (~val & 0xff) + carry;
1440 carry = val >> 8;
Damien George271d18e2015-04-25 23:16:39 +01001441 }
1442 if (big_endian) {
1443 *--b = val;
1444 if (b == buf) {
1445 return;
1446 }
1447 } else {
1448 *b++ = val;
1449 if (b == buf + len) {
1450 return;
1451 }
1452 }
1453 }
1454 }
1455}
1456
Damien Georgefb510b32014-06-01 13:32:54 +01001457#if MICROPY_PY_BUILTINS_FLOAT
Damien George52608102014-03-08 15:04:54 +00001458mp_float_t mpz_as_float(const mpz_t *i) {
1459 mp_float_t val = 0;
Damien George438c88d2014-02-22 19:25:23 +00001460 mpz_dig_t *d = i->dig + i->len;
1461
Damien George58056b02015-01-09 20:58:58 +00001462 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001463 val = val * DIG_BASE + *d;
Damien George438c88d2014-02-22 19:25:23 +00001464 }
1465
1466 if (i->neg != 0) {
1467 val = -val;
1468 }
1469
1470 return val;
1471}
Damien George52608102014-03-08 15:04:54 +00001472#endif
Damien George438c88d2014-02-22 19:25:23 +00001473
Damien Georgeafb1cf72014-09-05 20:37:06 +01001474mp_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 +00001475 if (base < 2 || base > 32) {
1476 return 0;
1477 }
1478
Damien Georgeafb1cf72014-09-05 20:37:06 +01001479 mp_uint_t num_digits = i->len * DIG_SIZE / log_base2_floor[base] + 1;
1480 mp_uint_t num_commas = comma ? num_digits / 3: 0;
1481 mp_uint_t prefix_len = prefix ? strlen(prefix) : 0;
Dave Hylandsc4029e52014-04-07 11:19:51 -07001482
1483 return num_digits + num_commas + prefix_len + 2; // +1 for sign, +1 for null byte
1484}
1485
Damien Georgeafb1cf72014-09-05 20:37:06 +01001486#if 0
1487this function is unused
1488char *mpz_as_str(const mpz_t *i, mp_uint_t base) {
1489 char *s = m_new(char, mpz_as_str_size(i, base, NULL, '\0'));
1490 mpz_as_str_inpl(i, base, NULL, 'a', '\0', s);
Damien George438c88d2014-02-22 19:25:23 +00001491 return s;
1492}
Damien Georgeafb1cf72014-09-05 20:37:06 +01001493#endif
Damien George438c88d2014-02-22 19:25:23 +00001494
1495// assumes enough space as calculated by mpz_as_str_size
1496// returns length of string, not including null byte
Damien Georgeafb1cf72014-09-05 20:37:06 +01001497mp_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 +00001498 if (str == NULL || base < 2 || base > 32) {
1499 str[0] = 0;
1500 return 0;
1501 }
1502
Damien Georgeafb1cf72014-09-05 20:37:06 +01001503 mp_uint_t ilen = i->len;
Damien George438c88d2014-02-22 19:25:23 +00001504
Dave Hylandsc4029e52014-04-07 11:19:51 -07001505 char *s = str;
Damien George438c88d2014-02-22 19:25:23 +00001506 if (ilen == 0) {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001507 if (prefix) {
1508 while (*prefix)
1509 *s++ = *prefix++;
1510 }
1511 *s++ = '0';
1512 *s = '\0';
1513 return s - str;
Damien George438c88d2014-02-22 19:25:23 +00001514 }
1515
Damien Georgeeec91052014-04-08 23:11:00 +01001516 // make a copy of mpz digits, so we can do the div/mod calculation
Damien George438c88d2014-02-22 19:25:23 +00001517 mpz_dig_t *dig = m_new(mpz_dig_t, ilen);
1518 memcpy(dig, i->dig, ilen * sizeof(mpz_dig_t));
1519
1520 // convert
Dave Hylandsc4029e52014-04-07 11:19:51 -07001521 char *last_comma = str;
Damien George438c88d2014-02-22 19:25:23 +00001522 bool done;
1523 do {
1524 mpz_dig_t *d = dig + ilen;
1525 mpz_dbl_dig_t a = 0;
1526
1527 // compute next remainder
1528 while (--d >= dig) {
1529 a = (a << DIG_SIZE) | *d;
1530 *d = a / base;
1531 a %= base;
1532 }
1533
1534 // convert to character
1535 a += '0';
1536 if (a > '9') {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001537 a += base_char - '9' - 1;
Damien George438c88d2014-02-22 19:25:23 +00001538 }
1539 *s++ = a;
1540
1541 // check if number is zero
1542 done = true;
1543 for (d = dig; d < dig + ilen; ++d) {
1544 if (*d != 0) {
1545 done = false;
1546 break;
1547 }
1548 }
Dave Hylandsc4029e52014-04-07 11:19:51 -07001549 if (comma && (s - last_comma) == 3) {
1550 *s++ = comma;
1551 last_comma = s;
1552 }
1553 }
1554 while (!done);
Damien George438c88d2014-02-22 19:25:23 +00001555
Damien Georgeeec91052014-04-08 23:11:00 +01001556 // free the copy of the digits array
1557 m_del(mpz_dig_t, dig, ilen);
1558
Dave Hylandsc4029e52014-04-07 11:19:51 -07001559 if (prefix) {
1560 const char *p = &prefix[strlen(prefix)];
1561 while (p > prefix) {
1562 *s++ = *--p;
1563 }
1564 }
Damien George438c88d2014-02-22 19:25:23 +00001565 if (i->neg != 0) {
1566 *s++ = '-';
1567 }
1568
1569 // reverse string
1570 for (char *u = str, *v = s - 1; u < v; ++u, --v) {
1571 char temp = *u;
1572 *u = *v;
1573 *v = temp;
1574 }
1575
Dave Hylandsc4029e52014-04-07 11:19:51 -07001576 *s = '\0'; // null termination
Damien George438c88d2014-02-22 19:25:23 +00001577
1578 return s - str;
1579}
1580
1581#endif // MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ