blob: b3f8b15b60c5b7f510f903be5aea7aedba437241 [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
Paul Sokolovskyb3be4712015-11-22 22:03:18 +0200288 // remove trailing zeros
289 for (--idig; idig >= oidig && *idig == 0; --idig) {
290 }
291
292 return idig + 1 - oidig;
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200293}
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;
Damien George2adf7ec2016-01-08 17:56:58 +0000714 #if MP_ENDIANNESS_LITTLE
David Steinbergc585ad12015-01-13 15:19:37 +0000715 struct { mp_float_int_t frc:MP_FLOAT_FRAC_BITS, exp:MP_FLOAT_EXP_BITS, sgn:1; } p;
Damien George2adf7ec2016-01-08 17:56:58 +0000716 #else
717 struct { mp_float_int_t sgn:1, exp:MP_FLOAT_EXP_BITS, frc:MP_FLOAT_FRAC_BITS; } p;
718 #endif
David Steinberg6e0b6d02015-01-02 12:39:22 +0000719 } u = {src};
720
721 z->neg = u.p.sgn;
722 if (u.p.exp == 0) {
723 // value == 0 || value < 1
Damien George58056b02015-01-09 20:58:58 +0000724 mpz_set_from_int(z, 0);
David Steinbergc585ad12015-01-13 15:19:37 +0000725 } else if (u.p.exp == ((1 << MP_FLOAT_EXP_BITS) - 1)) {
Damien George6fd4b362015-01-02 23:04:09 +0000726 // u.p.frc == 0 indicates inf, else NaN
727 // should be handled by caller
Damien George58056b02015-01-09 20:58:58 +0000728 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000729 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000730 const int adj_exp = (int)u.p.exp - MP_FLOAT_EXP_BIAS;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000731 if (adj_exp < 0) {
732 // value < 1 , truncates to 0
Damien George58056b02015-01-09 20:58:58 +0000733 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000734 } else if (adj_exp == 0) {
735 // 1 <= value < 2 , so truncates to 1
Damien George58056b02015-01-09 20:58:58 +0000736 mpz_set_from_int(z, 1);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000737 } else {
738 // 2 <= value
739 const int dig_cnt = (adj_exp + 1 + (DIG_SIZE - 1)) / DIG_SIZE;
740 const unsigned int rem = adj_exp % DIG_SIZE;
741 int dig_ind, shft;
David Steinbergc585ad12015-01-13 15:19:37 +0000742 mp_float_int_t frc = u.p.frc | ((mp_float_int_t)1 << MP_FLOAT_FRAC_BITS);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000743
David Steinbergc585ad12015-01-13 15:19:37 +0000744 if (adj_exp < MP_FLOAT_FRAC_BITS) {
David Steinberg6e0b6d02015-01-02 12:39:22 +0000745 shft = 0;
746 dig_ind = 0;
David Steinbergc585ad12015-01-13 15:19:37 +0000747 frc >>= MP_FLOAT_FRAC_BITS - adj_exp;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000748 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000749 shft = (rem - MP_FLOAT_FRAC_BITS) % DIG_SIZE;
750 dig_ind = (adj_exp - MP_FLOAT_FRAC_BITS) / DIG_SIZE;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000751 }
752 mpz_need_dig(z, dig_cnt);
753 z->len = dig_cnt;
754 if (dig_ind != 0) {
755 memset(z->dig, 0, dig_ind * sizeof(mpz_dig_t));
756 }
757 if (shft != 0) {
758 z->dig[dig_ind++] = (frc << shft) & DIG_MASK;
759 frc >>= DIG_SIZE - shft;
760 }
David Steinberg8d427b72015-01-13 15:20:32 +0000761#if DIG_SIZE < (MP_FLOAT_FRAC_BITS + 1)
David Steinberg6e0b6d02015-01-02 12:39:22 +0000762 while (dig_ind != dig_cnt) {
763 z->dig[dig_ind++] = frc & DIG_MASK;
764 frc >>= DIG_SIZE;
765 }
David Steinberg8d427b72015-01-13 15:20:32 +0000766#else
767 if (dig_ind != dig_cnt) {
768 z->dig[dig_ind] = frc;
769 }
770#endif
David Steinberg6e0b6d02015-01-02 12:39:22 +0000771 }
772 }
773}
David Steinberg6e0b6d02015-01-02 12:39:22 +0000774#endif
775
Damien George438c88d2014-02-22 19:25:23 +0000776// returns number of bytes from str that were processed
Damien Georgeafb1cf72014-09-05 20:37:06 +0100777mp_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 +0000778 assert(base < 36);
779
780 const char *cur = str;
781 const char *top = str + len;
782
783 mpz_need_dig(z, len * 8 / DIG_SIZE + 1);
784
785 if (neg) {
786 z->neg = 1;
787 } else {
788 z->neg = 0;
789 }
790
791 z->len = 0;
792 for (; cur < top; ++cur) { // XXX UTF8 next char
Damien Georgeafb1cf72014-09-05 20:37:06 +0100793 //mp_uint_t v = char_to_numeric(cur#); // XXX UTF8 get char
794 mp_uint_t v = *cur;
Damien George438c88d2014-02-22 19:25:23 +0000795 if ('0' <= v && v <= '9') {
796 v -= '0';
797 } else if ('A' <= v && v <= 'Z') {
798 v -= 'A' - 10;
799 } else if ('a' <= v && v <= 'z') {
800 v -= 'a' - 10;
801 } else {
802 break;
803 }
804 if (v >= base) {
805 break;
806 }
807 z->len = mpn_mul_dig_add_dig(z->dig, z->len, base, v);
808 }
809
810 return cur - str;
811}
812
813bool mpz_is_zero(const mpz_t *z) {
814 return z->len == 0;
815}
816
Damien Georgea2e38382015-03-02 12:58:06 +0000817#if 0
818these functions are unused
819
Damien George438c88d2014-02-22 19:25:23 +0000820bool mpz_is_pos(const mpz_t *z) {
821 return z->len > 0 && z->neg == 0;
822}
823
824bool mpz_is_neg(const mpz_t *z) {
825 return z->len > 0 && z->neg != 0;
826}
827
828bool mpz_is_odd(const mpz_t *z) {
829 return z->len > 0 && (z->dig[0] & 1) != 0;
830}
831
832bool mpz_is_even(const mpz_t *z) {
833 return z->len == 0 || (z->dig[0] & 1) == 0;
834}
Damien Georgea2e38382015-03-02 12:58:06 +0000835#endif
Damien George438c88d2014-02-22 19:25:23 +0000836
Damien George42f3de92014-10-03 17:44:14 +0000837int mpz_cmp(const mpz_t *z1, const mpz_t *z2) {
Damien Georgea9dc9b82015-01-27 17:47:38 +0000838 // to catch comparison of -0 with +0
839 if (z1->len == 0 && z2->len == 0) {
840 return 0;
841 }
Damien George42f3de92014-10-03 17:44:14 +0000842 int cmp = (int)z2->neg - (int)z1->neg;
Damien George438c88d2014-02-22 19:25:23 +0000843 if (cmp != 0) {
844 return cmp;
845 }
846 cmp = mpn_cmp(z1->dig, z1->len, z2->dig, z2->len);
847 if (z1->neg != 0) {
848 cmp = -cmp;
849 }
850 return cmp;
851}
852
Damien George06201ff2014-03-01 19:50:50 +0000853#if 0
854// obsolete
855// compares mpz with an integer that fits within DIG_SIZE bits
Damien Georgeafb1cf72014-09-05 20:37:06 +0100856mp_int_t mpz_cmp_sml_int(const mpz_t *z, mp_int_t sml_int) {
857 mp_int_t cmp;
Damien George438c88d2014-02-22 19:25:23 +0000858 if (z->neg == 0) {
859 if (sml_int < 0) return 1;
860 if (sml_int == 0) {
861 if (z->len == 0) return 0;
862 return 1;
863 }
864 if (z->len == 0) return -1;
865 assert(sml_int < (1 << DIG_SIZE));
866 if (z->len != 1) return 1;
867 cmp = z->dig[0] - sml_int;
868 } else {
869 if (sml_int > 0) return -1;
870 if (sml_int == 0) {
871 if (z->len == 0) return 0;
872 return -1;
873 }
874 if (z->len == 0) return 1;
875 assert(sml_int > -(1 << DIG_SIZE));
876 if (z->len != 1) return -1;
877 cmp = -z->dig[0] - sml_int;
878 }
879 if (cmp < 0) return -1;
880 if (cmp > 0) return 1;
881 return 0;
882}
Damien George06201ff2014-03-01 19:50:50 +0000883#endif
Damien George438c88d2014-02-22 19:25:23 +0000884
Damien George438c88d2014-02-22 19:25:23 +0000885#if 0
886these functions are unused
887
888/* returns abs(z)
889*/
890mpz_t *mpz_abs(const mpz_t *z) {
891 mpz_t *z2 = mpz_clone(z);
892 z2->neg = 0;
893 return z2;
894}
895
896/* returns -z
897*/
898mpz_t *mpz_neg(const mpz_t *z) {
899 mpz_t *z2 = mpz_clone(z);
900 z2->neg = 1 - z2->neg;
901 return z2;
902}
903
904/* returns lhs + rhs
905 can have lhs, rhs the same
906*/
907mpz_t *mpz_add(const mpz_t *lhs, const mpz_t *rhs) {
908 mpz_t *z = mpz_zero();
909 mpz_add_inpl(z, lhs, rhs);
910 return z;
911}
912
913/* returns lhs - rhs
914 can have lhs, rhs the same
915*/
916mpz_t *mpz_sub(const mpz_t *lhs, const mpz_t *rhs) {
917 mpz_t *z = mpz_zero();
918 mpz_sub_inpl(z, lhs, rhs);
919 return z;
920}
921
922/* returns lhs * rhs
923 can have lhs, rhs the same
924*/
925mpz_t *mpz_mul(const mpz_t *lhs, const mpz_t *rhs) {
926 mpz_t *z = mpz_zero();
927 mpz_mul_inpl(z, lhs, rhs);
928 return z;
929}
930
931/* returns lhs ** rhs
932 can have lhs, rhs the same
933*/
934mpz_t *mpz_pow(const mpz_t *lhs, const mpz_t *rhs) {
935 mpz_t *z = mpz_zero();
936 mpz_pow_inpl(z, lhs, rhs);
937 return z;
938}
Damien Georgea2e38382015-03-02 12:58:06 +0000939
940/* computes new integers in quo and rem such that:
941 quo * rhs + rem = lhs
942 0 <= rem < rhs
943 can have lhs, rhs the same
944*/
945void mpz_divmod(const mpz_t *lhs, const mpz_t *rhs, mpz_t **quo, mpz_t **rem) {
946 *quo = mpz_zero();
947 *rem = mpz_zero();
948 mpz_divmod_inpl(*quo, *rem, lhs, rhs);
949}
Damien George438c88d2014-02-22 19:25:23 +0000950#endif
951
952/* computes dest = abs(z)
953 can have dest, z the same
954*/
955void mpz_abs_inpl(mpz_t *dest, const mpz_t *z) {
956 if (dest != z) {
957 mpz_set(dest, z);
958 }
959 dest->neg = 0;
960}
961
962/* computes dest = -z
963 can have dest, z the same
964*/
965void mpz_neg_inpl(mpz_t *dest, const mpz_t *z) {
966 if (dest != z) {
967 mpz_set(dest, z);
968 }
969 dest->neg = 1 - dest->neg;
970}
971
Damien George06201ff2014-03-01 19:50:50 +0000972/* computes dest = ~z (= -z - 1)
973 can have dest, z the same
974*/
975void mpz_not_inpl(mpz_t *dest, const mpz_t *z) {
976 if (dest != z) {
977 mpz_set(dest, z);
978 }
Damien Georgee0ac1942014-12-31 19:35:01 +0000979 if (dest->len == 0) {
980 mpz_need_dig(dest, 1);
981 dest->dig[0] = 1;
982 dest->len = 1;
983 dest->neg = 1;
984 } else if (dest->neg) {
Damien George06201ff2014-03-01 19:50:50 +0000985 dest->neg = 0;
986 mpz_dig_t k = 1;
987 dest->len = mpn_sub(dest->dig, dest->dig, dest->len, &k, 1);
988 } else {
Damien Georgee0ac1942014-12-31 19:35:01 +0000989 mpz_need_dig(dest, dest->len + 1);
Damien George06201ff2014-03-01 19:50:50 +0000990 mpz_dig_t k = 1;
991 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &k, 1);
992 dest->neg = 1;
993 }
994}
995
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000996/* computes dest = lhs << rhs
997 can have dest, lhs the same
998*/
Damien George2f4e8512015-10-01 18:01:37 +0100999void mpz_shl_inpl(mpz_t *dest, const mpz_t *lhs, mp_uint_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +00001000 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001001 mpz_set(dest, lhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001002 } else {
Damien George06201ff2014-03-01 19:50:50 +00001003 mpz_need_dig(dest, lhs->len + (rhs + DIG_SIZE - 1) / DIG_SIZE);
1004 dest->len = mpn_shl(dest->dig, lhs->dig, lhs->len, rhs);
1005 dest->neg = lhs->neg;
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001006 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001007}
1008
1009/* computes dest = lhs >> rhs
1010 can have dest, lhs the same
1011*/
Damien George2f4e8512015-10-01 18:01:37 +01001012void mpz_shr_inpl(mpz_t *dest, const mpz_t *lhs, mp_uint_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +00001013 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001014 mpz_set(dest, lhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001015 } else {
Damien George06201ff2014-03-01 19:50:50 +00001016 mpz_need_dig(dest, lhs->len);
1017 dest->len = mpn_shr(dest->dig, lhs->dig, lhs->len, rhs);
1018 dest->neg = lhs->neg;
1019 if (dest->neg) {
1020 // arithmetic shift right, rounding to negative infinity
Damien Georgeafb1cf72014-09-05 20:37:06 +01001021 mp_uint_t n_whole = rhs / DIG_SIZE;
1022 mp_uint_t n_part = rhs % DIG_SIZE;
Damien George06201ff2014-03-01 19:50:50 +00001023 mpz_dig_t round_up = 0;
Damien Georgeafb1cf72014-09-05 20:37:06 +01001024 for (mp_uint_t i = 0; i < lhs->len && i < n_whole; i++) {
Damien George06201ff2014-03-01 19:50:50 +00001025 if (lhs->dig[i] != 0) {
1026 round_up = 1;
1027 break;
1028 }
1029 }
1030 if (n_whole < lhs->len && (lhs->dig[n_whole] & ((1 << n_part) - 1)) != 0) {
1031 round_up = 1;
1032 }
1033 if (round_up) {
Damien Georgee0ac1942014-12-31 19:35:01 +00001034 if (dest->len == 0) {
1035 // dest == 0, so need to add 1 by hand (answer will be -1)
1036 dest->dig[0] = 1;
1037 dest->len = 1;
1038 } else {
1039 // dest > 0, so can use mpn_add to add 1
1040 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &round_up, 1);
1041 }
Damien George06201ff2014-03-01 19:50:50 +00001042 }
1043 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001044 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001045}
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001046
Damien George438c88d2014-02-22 19:25:23 +00001047/* computes dest = lhs + rhs
1048 can have dest, lhs, rhs the same
1049*/
1050void mpz_add_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1051 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1052 const mpz_t *temp = lhs;
1053 lhs = rhs;
1054 rhs = temp;
1055 }
1056
1057 if (lhs->neg == rhs->neg) {
1058 mpz_need_dig(dest, lhs->len + 1);
1059 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1060 } else {
1061 mpz_need_dig(dest, lhs->len);
1062 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1063 }
1064
1065 dest->neg = lhs->neg;
1066}
1067
1068/* computes dest = lhs - rhs
1069 can have dest, lhs, rhs the same
1070*/
1071void mpz_sub_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1072 bool neg = false;
1073
1074 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1075 const mpz_t *temp = lhs;
1076 lhs = rhs;
1077 rhs = temp;
1078 neg = true;
1079 }
1080
1081 if (lhs->neg != rhs->neg) {
1082 mpz_need_dig(dest, lhs->len + 1);
1083 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1084 } else {
1085 mpz_need_dig(dest, lhs->len);
1086 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1087 }
1088
1089 if (neg) {
1090 dest->neg = 1 - lhs->neg;
1091 } else {
1092 dest->neg = lhs->neg;
1093 }
1094}
1095
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001096/* computes dest = lhs & rhs
1097 can have dest, lhs, rhs the same
1098*/
1099void mpz_and_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001100 if (lhs->neg == rhs->neg) {
Damien Georgef55cf102014-05-29 15:01:49 +01001101 if (lhs->neg == 0) {
1102 // make sure lhs has the most digits
1103 if (lhs->len < rhs->len) {
1104 const mpz_t *temp = lhs;
1105 lhs = rhs;
1106 rhs = temp;
1107 }
1108 // do the and'ing
1109 mpz_need_dig(dest, rhs->len);
Damien Georgeff8dd3f2015-01-20 12:47:20 +00001110 dest->len = mpn_and(dest->dig, lhs->dig, rhs->dig, rhs->len);
Damien Georgef55cf102014-05-29 15:01:49 +01001111 dest->neg = 0;
1112 } else {
1113 // TODO both args are negative
Damien George4c02e542015-10-01 17:57:36 +01001114 mp_not_implemented("bignum and with negative args");
Damien Georgef55cf102014-05-29 15:01:49 +01001115 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001116 } else {
Damien Georgef55cf102014-05-29 15:01:49 +01001117 // args have different sign
1118 // make sure lhs is the positive arg
1119 if (rhs->neg == 0) {
1120 const mpz_t *temp = lhs;
1121 lhs = rhs;
1122 rhs = temp;
1123 }
1124 mpz_need_dig(dest, lhs->len + 1);
1125 dest->len = mpn_and_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1126 assert(dest->len <= dest->alloc);
1127 dest->neg = 0;
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001128 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001129}
1130
1131/* computes dest = lhs | rhs
1132 can have dest, lhs, rhs the same
1133*/
1134void mpz_or_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1135 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1136 const mpz_t *temp = lhs;
1137 lhs = rhs;
1138 rhs = temp;
1139 }
1140
1141 if (lhs->neg == rhs->neg) {
1142 mpz_need_dig(dest, lhs->len);
1143 dest->len = mpn_or(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1144 } else {
1145 mpz_need_dig(dest, lhs->len);
1146 // TODO
Damien George4c02e542015-10-01 17:57:36 +01001147 mp_not_implemented("bignum or with negative args");
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001148// dest->len = mpn_or_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1149 }
1150
1151 dest->neg = lhs->neg;
1152}
1153
1154/* computes dest = lhs ^ rhs
1155 can have dest, lhs, rhs the same
1156*/
1157void mpz_xor_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1158 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1159 const mpz_t *temp = lhs;
1160 lhs = rhs;
1161 rhs = temp;
1162 }
1163
1164 if (lhs->neg == rhs->neg) {
1165 mpz_need_dig(dest, lhs->len);
1166 dest->len = mpn_xor(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1167 } else {
1168 mpz_need_dig(dest, lhs->len);
1169 // TODO
Damien George4c02e542015-10-01 17:57:36 +01001170 mp_not_implemented("bignum xor with negative args");
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001171// dest->len = mpn_xor_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1172 }
1173
1174 dest->neg = 0;
1175}
1176
Damien George438c88d2014-02-22 19:25:23 +00001177/* computes dest = lhs * rhs
1178 can have dest, lhs, rhs the same
1179*/
Damien George4dea9222015-04-09 15:29:54 +00001180void mpz_mul_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Damien George438c88d2014-02-22 19:25:23 +00001181 if (lhs->len == 0 || rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001182 mpz_set_from_int(dest, 0);
1183 return;
Damien George438c88d2014-02-22 19:25:23 +00001184 }
1185
1186 mpz_t *temp = NULL;
1187 if (lhs == dest) {
1188 lhs = temp = mpz_clone(lhs);
1189 if (rhs == dest) {
1190 rhs = lhs;
1191 }
1192 } else if (rhs == dest) {
1193 rhs = temp = mpz_clone(rhs);
1194 }
1195
1196 mpz_need_dig(dest, lhs->len + rhs->len); // min mem l+r-1, max mem l+r
1197 memset(dest->dig, 0, dest->alloc * sizeof(mpz_dig_t));
1198 dest->len = mpn_mul(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1199
1200 if (lhs->neg == rhs->neg) {
1201 dest->neg = 0;
1202 } else {
1203 dest->neg = 1;
1204 }
1205
1206 mpz_free(temp);
1207}
1208
1209/* computes dest = lhs ** rhs
1210 can have dest, lhs, rhs the same
1211*/
1212void mpz_pow_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1213 if (lhs->len == 0 || rhs->neg != 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001214 mpz_set_from_int(dest, 0);
1215 return;
Damien George438c88d2014-02-22 19:25:23 +00001216 }
1217
1218 if (rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001219 mpz_set_from_int(dest, 1);
1220 return;
Damien George438c88d2014-02-22 19:25:23 +00001221 }
1222
1223 mpz_t *x = mpz_clone(lhs);
1224 mpz_t *n = mpz_clone(rhs);
1225
1226 mpz_set_from_int(dest, 1);
1227
1228 while (n->len > 0) {
Damien Georgea2e38382015-03-02 12:58:06 +00001229 if ((n->dig[0] & 1) != 0) {
Damien George438c88d2014-02-22 19:25:23 +00001230 mpz_mul_inpl(dest, dest, x);
1231 }
Damien George438c88d2014-02-22 19:25:23 +00001232 n->len = mpn_shr(n->dig, n->dig, n->len, 1);
Damien George5bf565e2014-04-04 00:16:32 +01001233 if (n->len == 0) {
1234 break;
1235 }
1236 mpz_mul_inpl(x, x, x);
Damien George438c88d2014-02-22 19:25:23 +00001237 }
1238
1239 mpz_free(x);
1240 mpz_free(n);
1241}
1242
Damien Georgea2e38382015-03-02 12:58:06 +00001243#if 0
1244these functions are unused
1245
Damien George438c88d2014-02-22 19:25:23 +00001246/* computes gcd(z1, z2)
1247 based on Knuth's modified gcd algorithm (I think?)
1248 gcd(z1, z2) >= 0
1249 gcd(0, 0) = 0
1250 gcd(z, 0) = abs(z)
1251*/
1252mpz_t *mpz_gcd(const mpz_t *z1, const mpz_t *z2) {
1253 if (z1->len == 0) {
1254 mpz_t *a = mpz_clone(z2);
1255 a->neg = 0;
1256 return a;
1257 } else if (z2->len == 0) {
1258 mpz_t *a = mpz_clone(z1);
1259 a->neg = 0;
1260 return a;
1261 }
1262
1263 mpz_t *a = mpz_clone(z1);
1264 mpz_t *b = mpz_clone(z2);
1265 mpz_t c; mpz_init_zero(&c);
1266 a->neg = 0;
1267 b->neg = 0;
1268
1269 for (;;) {
1270 if (mpz_cmp(a, b) < 0) {
1271 if (a->len == 0) {
1272 mpz_free(a);
1273 mpz_deinit(&c);
1274 return b;
1275 }
1276 mpz_t *t = a; a = b; b = t;
1277 }
1278 if (!(b->len >= 2 || (b->len == 1 && b->dig[0] > 1))) { // compute b > 0; could be mpz_cmp_small_int(b, 1) > 0
1279 break;
1280 }
1281 mpz_set(&c, b);
1282 do {
1283 mpz_add_inpl(&c, &c, &c);
1284 } while (mpz_cmp(&c, a) <= 0);
1285 c.len = mpn_shr(c.dig, c.dig, c.len, 1);
1286 mpz_sub_inpl(a, a, &c);
1287 }
1288
1289 mpz_deinit(&c);
1290
1291 if (b->len == 1 && b->dig[0] == 1) { // compute b == 1; could be mpz_cmp_small_int(b, 1) == 0
1292 mpz_free(a);
1293 return b;
1294 } else {
1295 mpz_free(b);
1296 return a;
1297 }
1298}
1299
1300/* computes lcm(z1, z2)
1301 = abs(z1) / gcd(z1, z2) * abs(z2)
1302 lcm(z1, z1) >= 0
1303 lcm(0, 0) = 0
1304 lcm(z, 0) = 0
1305*/
Damien George5d9b8162014-08-07 14:27:48 +00001306mpz_t *mpz_lcm(const mpz_t *z1, const mpz_t *z2) {
1307 if (z1->len == 0 || z2->len == 0) {
1308 return mpz_zero();
stijn01d6be42014-05-05 12:18:27 +02001309 }
Damien George438c88d2014-02-22 19:25:23 +00001310
1311 mpz_t *gcd = mpz_gcd(z1, z2);
1312 mpz_t *quo = mpz_zero();
1313 mpz_t *rem = mpz_zero();
1314 mpz_divmod_inpl(quo, rem, z1, gcd);
1315 mpz_mul_inpl(rem, quo, z2);
1316 mpz_free(gcd);
1317 mpz_free(quo);
1318 rem->neg = 0;
1319 return rem;
1320}
Damien Georgea2e38382015-03-02 12:58:06 +00001321#endif
Damien George438c88d2014-02-22 19:25:23 +00001322
1323/* computes new integers in quo and rem such that:
1324 quo * rhs + rem = lhs
1325 0 <= rem < rhs
1326 can have lhs, rhs the same
1327*/
1328void mpz_divmod_inpl(mpz_t *dest_quo, mpz_t *dest_rem, const mpz_t *lhs, const mpz_t *rhs) {
1329 if (rhs->len == 0) {
1330 mpz_set_from_int(dest_quo, 0);
1331 mpz_set_from_int(dest_rem, 0);
1332 return;
1333 }
1334
1335 mpz_need_dig(dest_quo, lhs->len + 1); // +1 necessary?
1336 memset(dest_quo->dig, 0, (lhs->len + 1) * sizeof(mpz_dig_t));
1337 dest_quo->len = 0;
1338 mpz_need_dig(dest_rem, lhs->len + 1); // +1 necessary?
1339 mpz_set(dest_rem, lhs);
1340 //rhs->dig[rhs->len] = 0;
1341 mpn_div(dest_rem->dig, &dest_rem->len, rhs->dig, rhs->len, dest_quo->dig, &dest_quo->len);
1342
1343 if (lhs->neg != rhs->neg) {
1344 dest_quo->neg = 1;
1345 }
1346}
1347
1348#if 0
1349these functions are unused
1350
1351/* computes floor(lhs / rhs)
1352 can have lhs, rhs the same
1353*/
1354mpz_t *mpz_div(const mpz_t *lhs, const mpz_t *rhs) {
1355 mpz_t *quo = mpz_zero();
1356 mpz_t rem; mpz_init_zero(&rem);
1357 mpz_divmod_inpl(quo, &rem, lhs, rhs);
1358 mpz_deinit(&rem);
1359 return quo;
1360}
1361
1362/* computes lhs % rhs ( >= 0)
1363 can have lhs, rhs the same
1364*/
1365mpz_t *mpz_mod(const mpz_t *lhs, const mpz_t *rhs) {
1366 mpz_t quo; mpz_init_zero(&quo);
1367 mpz_t *rem = mpz_zero();
1368 mpz_divmod_inpl(&quo, rem, lhs, rhs);
1369 mpz_deinit(&quo);
1370 return rem;
1371}
1372#endif
1373
Damien Georgeffe911d2014-07-24 14:21:37 +01001374// must return actual int value if it fits in mp_int_t
1375mp_int_t mpz_hash(const mpz_t *z) {
1376 mp_int_t val = 0;
1377 mpz_dig_t *d = z->dig + z->len;
1378
Damien George58056b02015-01-09 20:58:58 +00001379 while (d-- > z->dig) {
Damien Georgeffe911d2014-07-24 14:21:37 +01001380 val = (val << DIG_SIZE) | *d;
1381 }
1382
1383 if (z->neg != 0) {
1384 val = -val;
1385 }
1386
1387 return val;
1388}
1389
Damien George40f3c022014-07-03 13:25:24 +01001390bool mpz_as_int_checked(const mpz_t *i, mp_int_t *value) {
Damien George963a5a32015-01-16 17:47:07 +00001391 mp_uint_t val = 0;
Damien George8270e382014-04-03 11:00:54 +00001392 mpz_dig_t *d = i->dig + i->len;
1393
Damien George58056b02015-01-09 20:58:58 +00001394 while (d-- > i->dig) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001395 if (val > (~(WORD_MSBIT_HIGH) >> DIG_SIZE)) {
1396 // will overflow
Damien George8270e382014-04-03 11:00:54 +00001397 return false;
1398 }
Damien Georgec9aa58e2014-07-31 13:41:43 +00001399 val = (val << DIG_SIZE) | *d;
Damien George8270e382014-04-03 11:00:54 +00001400 }
1401
1402 if (i->neg != 0) {
1403 val = -val;
1404 }
1405
1406 *value = val;
1407 return true;
1408}
1409
Damien Georgec9aa58e2014-07-31 13:41:43 +00001410bool mpz_as_uint_checked(const mpz_t *i, mp_uint_t *value) {
1411 if (i->neg != 0) {
1412 // can't represent signed values
1413 return false;
1414 }
1415
1416 mp_uint_t val = 0;
1417 mpz_dig_t *d = i->dig + i->len;
1418
Damien George58056b02015-01-09 20:58:58 +00001419 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001420 if (val > (~(WORD_MSBIT_HIGH) >> (DIG_SIZE - 1))) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001421 // will overflow
1422 return false;
1423 }
1424 val = (val << DIG_SIZE) | *d;
1425 }
1426
1427 *value = val;
1428 return true;
1429}
1430
Damien George271d18e2015-04-25 23:16:39 +01001431// writes at most len bytes to buf (so buf should be zeroed before calling)
1432void mpz_as_bytes(const mpz_t *z, bool big_endian, mp_uint_t len, byte *buf) {
1433 byte *b = buf;
1434 if (big_endian) {
1435 b += len;
1436 }
1437 mpz_dig_t *zdig = z->dig;
1438 int bits = 0;
1439 mpz_dbl_dig_t d = 0;
1440 mpz_dbl_dig_t carry = 1;
1441 for (mp_uint_t zlen = z->len; zlen > 0; --zlen) {
1442 bits += DIG_SIZE;
1443 d = (d << DIG_SIZE) | *zdig++;
1444 for (; bits >= 8; bits -= 8, d >>= 8) {
1445 mpz_dig_t val = d;
1446 if (z->neg) {
Damien George94729072015-04-25 23:51:14 +01001447 val = (~val & 0xff) + carry;
1448 carry = val >> 8;
Damien George271d18e2015-04-25 23:16:39 +01001449 }
1450 if (big_endian) {
1451 *--b = val;
1452 if (b == buf) {
1453 return;
1454 }
1455 } else {
1456 *b++ = val;
1457 if (b == buf + len) {
1458 return;
1459 }
1460 }
1461 }
1462 }
1463}
1464
Damien Georgefb510b32014-06-01 13:32:54 +01001465#if MICROPY_PY_BUILTINS_FLOAT
Damien George52608102014-03-08 15:04:54 +00001466mp_float_t mpz_as_float(const mpz_t *i) {
1467 mp_float_t val = 0;
Damien George438c88d2014-02-22 19:25:23 +00001468 mpz_dig_t *d = i->dig + i->len;
1469
Damien George58056b02015-01-09 20:58:58 +00001470 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001471 val = val * DIG_BASE + *d;
Damien George438c88d2014-02-22 19:25:23 +00001472 }
1473
1474 if (i->neg != 0) {
1475 val = -val;
1476 }
1477
1478 return val;
1479}
Damien George52608102014-03-08 15:04:54 +00001480#endif
Damien George438c88d2014-02-22 19:25:23 +00001481
Damien Georgeafb1cf72014-09-05 20:37:06 +01001482mp_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 +00001483 if (base < 2 || base > 32) {
1484 return 0;
1485 }
1486
Damien Georgeafb1cf72014-09-05 20:37:06 +01001487 mp_uint_t num_digits = i->len * DIG_SIZE / log_base2_floor[base] + 1;
1488 mp_uint_t num_commas = comma ? num_digits / 3: 0;
1489 mp_uint_t prefix_len = prefix ? strlen(prefix) : 0;
Dave Hylandsc4029e52014-04-07 11:19:51 -07001490
1491 return num_digits + num_commas + prefix_len + 2; // +1 for sign, +1 for null byte
1492}
1493
Damien Georgeafb1cf72014-09-05 20:37:06 +01001494#if 0
1495this function is unused
1496char *mpz_as_str(const mpz_t *i, mp_uint_t base) {
1497 char *s = m_new(char, mpz_as_str_size(i, base, NULL, '\0'));
1498 mpz_as_str_inpl(i, base, NULL, 'a', '\0', s);
Damien George438c88d2014-02-22 19:25:23 +00001499 return s;
1500}
Damien Georgeafb1cf72014-09-05 20:37:06 +01001501#endif
Damien George438c88d2014-02-22 19:25:23 +00001502
1503// assumes enough space as calculated by mpz_as_str_size
1504// returns length of string, not including null byte
Damien Georgeafb1cf72014-09-05 20:37:06 +01001505mp_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 +00001506 if (str == NULL || base < 2 || base > 32) {
1507 str[0] = 0;
1508 return 0;
1509 }
1510
Damien Georgeafb1cf72014-09-05 20:37:06 +01001511 mp_uint_t ilen = i->len;
Damien George438c88d2014-02-22 19:25:23 +00001512
Dave Hylandsc4029e52014-04-07 11:19:51 -07001513 char *s = str;
Damien George438c88d2014-02-22 19:25:23 +00001514 if (ilen == 0) {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001515 if (prefix) {
1516 while (*prefix)
1517 *s++ = *prefix++;
1518 }
1519 *s++ = '0';
1520 *s = '\0';
1521 return s - str;
Damien George438c88d2014-02-22 19:25:23 +00001522 }
1523
Damien Georgeeec91052014-04-08 23:11:00 +01001524 // make a copy of mpz digits, so we can do the div/mod calculation
Damien George438c88d2014-02-22 19:25:23 +00001525 mpz_dig_t *dig = m_new(mpz_dig_t, ilen);
1526 memcpy(dig, i->dig, ilen * sizeof(mpz_dig_t));
1527
1528 // convert
Dave Hylandsc4029e52014-04-07 11:19:51 -07001529 char *last_comma = str;
Damien George438c88d2014-02-22 19:25:23 +00001530 bool done;
1531 do {
1532 mpz_dig_t *d = dig + ilen;
1533 mpz_dbl_dig_t a = 0;
1534
1535 // compute next remainder
1536 while (--d >= dig) {
1537 a = (a << DIG_SIZE) | *d;
1538 *d = a / base;
1539 a %= base;
1540 }
1541
1542 // convert to character
1543 a += '0';
1544 if (a > '9') {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001545 a += base_char - '9' - 1;
Damien George438c88d2014-02-22 19:25:23 +00001546 }
1547 *s++ = a;
1548
1549 // check if number is zero
1550 done = true;
1551 for (d = dig; d < dig + ilen; ++d) {
1552 if (*d != 0) {
1553 done = false;
1554 break;
1555 }
1556 }
Dave Hylandsc4029e52014-04-07 11:19:51 -07001557 if (comma && (s - last_comma) == 3) {
1558 *s++ = comma;
1559 last_comma = s;
1560 }
1561 }
1562 while (!done);
Damien George438c88d2014-02-22 19:25:23 +00001563
Damien Georgeeec91052014-04-08 23:11:00 +01001564 // free the copy of the digits array
1565 m_del(mpz_dig_t, dig, ilen);
1566
Dave Hylandsc4029e52014-04-07 11:19:51 -07001567 if (prefix) {
1568 const char *p = &prefix[strlen(prefix)];
1569 while (p > prefix) {
1570 *s++ = *--p;
1571 }
1572 }
Damien George438c88d2014-02-22 19:25:23 +00001573 if (i->neg != 0) {
1574 *s++ = '-';
1575 }
1576
1577 // reverse string
1578 for (char *u = str, *v = s - 1; u < v; ++u, --v) {
1579 char temp = *u;
1580 *u = *v;
1581 *v = temp;
1582 }
1583
Dave Hylandsc4029e52014-04-07 11:19:51 -07001584 *s = '\0'; // null termination
Damien George438c88d2014-02-22 19:25:23 +00001585
1586 return s - str;
1587}
1588
1589#endif // MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ