blob: a056a6e8ad620624840ea8d0fd2cc47666c20587 [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
32#if MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ
33
Damien George06201ff2014-03-01 19:50:50 +000034#define DIG_SIZE (MPZ_DIG_SIZE)
stijn0e557fa2014-10-30 14:39:22 +010035#define DIG_MASK ((MPZ_LONG_1 << DIG_SIZE) - 1)
36#define DIG_MSB (MPZ_LONG_1 << (DIG_SIZE - 1))
37#define DIG_BASE (MPZ_LONG_1 << DIG_SIZE)
Damien George438c88d2014-02-22 19:25:23 +000038
39/*
Damien George06201ff2014-03-01 19:50:50 +000040 mpz is an arbitrary precision integer type with a public API.
41
42 mpn functions act on non-negative integers represented by an array of generalised
43 digits (eg a word per digit). You also need to specify separately the length of the
44 array. There is no public API for mpn. Rather, the functions are used by mpz to
45 implement its features.
46
47 Integer values are stored little endian (first digit is first in memory).
48
49 Definition of normalise: ?
Damien George438c88d2014-02-22 19:25:23 +000050*/
51
52/* compares i with j
53 returns sign(i - j)
54 assumes i, j are normalised
55*/
Damien George42f3de92014-10-03 17:44:14 +000056STATIC 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 +000057 if (ilen < jlen) { return -1; }
58 if (ilen > jlen) { return 1; }
59
60 for (idig += ilen, jdig += ilen; ilen > 0; --ilen) {
Damien George9a21d2e2014-09-06 17:15:34 +010061 mpz_dbl_dig_signed_t cmp = (mpz_dbl_dig_t)*(--idig) - (mpz_dbl_dig_t)*(--jdig);
Damien George438c88d2014-02-22 19:25:23 +000062 if (cmp < 0) { return -1; }
63 if (cmp > 0) { return 1; }
64 }
65
66 return 0;
67}
68
Damien Georgec5ac2ac2014-02-26 16:56:30 +000069/* computes i = j << n
70 returns number of digits in i
Damien George06201ff2014-03-01 19:50:50 +000071 assumes enough memory in i; assumes normalised j; assumes n > 0
Damien Georgec5ac2ac2014-02-26 16:56:30 +000072 can have i, j pointing to same memory
73*/
Damien Georgeafb1cf72014-09-05 20:37:06 +010074STATIC mp_uint_t mpn_shl(mpz_dig_t *idig, mpz_dig_t *jdig, mp_uint_t jlen, mp_uint_t n) {
75 mp_uint_t n_whole = (n + DIG_SIZE - 1) / DIG_SIZE;
76 mp_uint_t n_part = n % DIG_SIZE;
Damien Georgebb4a43f2014-03-12 15:36:06 +000077 if (n_part == 0) {
78 n_part = DIG_SIZE;
79 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +000080
Damien George06201ff2014-03-01 19:50:50 +000081 // start from the high end of the digit arrays
82 idig += jlen + n_whole - 1;
83 jdig += jlen - 1;
Damien Georgec5ac2ac2014-02-26 16:56:30 +000084
Damien George06201ff2014-03-01 19:50:50 +000085 // shift the digits
86 mpz_dbl_dig_t d = 0;
Damien Georgeafb1cf72014-09-05 20:37:06 +010087 for (mp_uint_t i = jlen; i > 0; i--, idig--, jdig--) {
Damien George06201ff2014-03-01 19:50:50 +000088 d |= *jdig;
Damien George5d9b8162014-08-07 14:27:48 +000089 *idig = (d >> (DIG_SIZE - n_part)) & DIG_MASK;
Damien George06201ff2014-03-01 19:50:50 +000090 d <<= DIG_SIZE;
Damien Georgec5ac2ac2014-02-26 16:56:30 +000091 }
92
Damien George06201ff2014-03-01 19:50:50 +000093 // store remaining bits
Damien George5d9b8162014-08-07 14:27:48 +000094 *idig = (d >> (DIG_SIZE - n_part)) & DIG_MASK;
Damien George06201ff2014-03-01 19:50:50 +000095 idig -= n_whole - 1;
Damien Georgebb4a43f2014-03-12 15:36:06 +000096 memset(idig, 0, (n_whole - 1) * sizeof(mpz_dig_t));
Damien George06201ff2014-03-01 19:50:50 +000097
98 // work out length of result
99 jlen += n_whole;
100 if (idig[jlen - 1] == 0) {
101 jlen--;
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000102 }
103
Damien George06201ff2014-03-01 19:50:50 +0000104 // return length of result
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000105 return jlen;
106}
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000107
Damien George438c88d2014-02-22 19:25:23 +0000108/* computes i = j >> n
109 returns number of digits in i
Damien George06201ff2014-03-01 19:50:50 +0000110 assumes enough memory in i; assumes normalised j; assumes n > 0
Damien George438c88d2014-02-22 19:25:23 +0000111 can have i, j pointing to same memory
112*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100113STATIC mp_uint_t mpn_shr(mpz_dig_t *idig, mpz_dig_t *jdig, mp_uint_t jlen, mp_uint_t n) {
114 mp_uint_t n_whole = n / DIG_SIZE;
115 mp_uint_t n_part = n % DIG_SIZE;
Damien George438c88d2014-02-22 19:25:23 +0000116
117 if (n_whole >= jlen) {
118 return 0;
119 }
120
121 jdig += n_whole;
122 jlen -= n_whole;
123
Damien Georgeafb1cf72014-09-05 20:37:06 +0100124 for (mp_uint_t i = jlen; i > 0; i--, idig++, jdig++) {
Damien George438c88d2014-02-22 19:25:23 +0000125 mpz_dbl_dig_t d = *jdig;
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000126 if (i > 1) {
Damien George9a21d2e2014-09-06 17:15:34 +0100127 d |= (mpz_dbl_dig_t)jdig[1] << DIG_SIZE;
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000128 }
Damien George438c88d2014-02-22 19:25:23 +0000129 d >>= n_part;
130 *idig = d & DIG_MASK;
131 }
132
133 if (idig[-1] == 0) {
Damien George06201ff2014-03-01 19:50:50 +0000134 jlen--;
Damien George438c88d2014-02-22 19:25:23 +0000135 }
136
137 return jlen;
138}
139
140/* computes i = j + k
141 returns number of digits in i
142 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
143 can have i, j, k pointing to same memory
144*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100145STATIC 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 +0000146 mpz_dig_t *oidig = idig;
147 mpz_dbl_dig_t carry = 0;
148
149 jlen -= klen;
150
151 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
Damien George9a21d2e2014-09-06 17:15:34 +0100152 carry += (mpz_dbl_dig_t)*jdig + (mpz_dbl_dig_t)*kdig;
Damien George438c88d2014-02-22 19:25:23 +0000153 *idig = carry & DIG_MASK;
154 carry >>= DIG_SIZE;
155 }
156
157 for (; jlen > 0; --jlen, ++idig, ++jdig) {
158 carry += *jdig;
159 *idig = carry & DIG_MASK;
160 carry >>= DIG_SIZE;
161 }
162
163 if (carry != 0) {
164 *idig++ = carry;
165 }
166
167 return idig - oidig;
168}
169
170/* computes i = j - k
171 returns number of digits in i
172 assumes enough memory in i; assumes normalised j, k; assumes j >= k
173 can have i, j, k pointing to same memory
174*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100175STATIC 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 +0000176 mpz_dig_t *oidig = idig;
177 mpz_dbl_dig_signed_t borrow = 0;
178
179 jlen -= klen;
180
181 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
Damien George9a21d2e2014-09-06 17:15:34 +0100182 borrow += (mpz_dbl_dig_t)*jdig - (mpz_dbl_dig_t)*kdig;
Damien George438c88d2014-02-22 19:25:23 +0000183 *idig = borrow & DIG_MASK;
184 borrow >>= DIG_SIZE;
185 }
186
Damien Georgeaca14122014-02-24 21:32:52 +0000187 for (; jlen > 0; --jlen, ++idig, ++jdig) {
Damien George438c88d2014-02-22 19:25:23 +0000188 borrow += *jdig;
189 *idig = borrow & DIG_MASK;
190 borrow >>= DIG_SIZE;
191 }
192
193 for (--idig; idig >= oidig && *idig == 0; --idig) {
194 }
195
196 return idig + 1 - oidig;
197}
198
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200199/* computes i = j & k
200 returns number of digits in i
Damien Georgeff8dd3f2015-01-20 12:47:20 +0000201 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen (jlen argument not needed)
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200202 can have i, j, k pointing to same memory
203*/
Damien Georgeff8dd3f2015-01-20 12:47:20 +0000204STATIC 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 +0200205 mpz_dig_t *oidig = idig;
206
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200207 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
208 *idig = *jdig & *kdig;
209 }
210
Damien George561e4252014-05-12 23:27:29 +0100211 // remove trailing zeros
Damien George51fab282014-05-13 22:58:00 +0100212 for (--idig; idig >= oidig && *idig == 0; --idig) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200213 }
214
Damien George51fab282014-05-13 22:58:00 +0100215 return idig + 1 - oidig;
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200216}
217
Damien Georgef55cf102014-05-29 15:01:49 +0100218/* computes i = j & -k = j & (~k + 1)
219 returns number of digits in i
220 assumes enough memory in i; assumes normalised j, k
221 can have i, j, k pointing to same memory
222*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100223STATIC 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 +0100224 mpz_dig_t *oidig = idig;
225 mpz_dbl_dig_t carry = 1;
226
227 for (; jlen > 0 && klen > 0; --jlen, --klen, ++idig, ++jdig, ++kdig) {
228 carry += *kdig ^ DIG_MASK;
229 *idig = (*jdig & carry) & DIG_MASK;
230 carry >>= DIG_SIZE;
231 }
232
233 for (; jlen > 0; --jlen, ++idig, ++jdig) {
234 carry += DIG_MASK;
235 *idig = (*jdig & carry) & DIG_MASK;
236 carry >>= DIG_SIZE;
237 }
238
239 if (carry != 0) {
240 *idig = carry;
241 } else {
242 // remove trailing zeros
243 for (--idig; idig >= oidig && *idig == 0; --idig) {
244 }
245 }
246
247 return idig + 1 - oidig;
248}
249
Paul Sokolovsky57207b82014-03-23 01:52:36 +0200250/* computes i = j | k
251 returns number of digits in i
252 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
253 can have i, j, k pointing to same memory
254*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100255STATIC 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 +0200256 mpz_dig_t *oidig = idig;
257
258 jlen -= klen;
259
260 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
261 *idig = *jdig | *kdig;
262 }
263
264 for (; jlen > 0; --jlen, ++idig, ++jdig) {
265 *idig = *jdig;
266 }
267
268 return idig - oidig;
269}
270
271/* computes i = j ^ k
272 returns number of digits in i
273 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
274 can have i, j, k pointing to same memory
275*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100276STATIC 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 +0200277 mpz_dig_t *oidig = idig;
278
279 jlen -= klen;
280
281 for (; klen > 0; --klen, ++idig, ++jdig, ++kdig) {
282 *idig = *jdig ^ *kdig;
283 }
284
285 for (; jlen > 0; --jlen, ++idig, ++jdig) {
286 *idig = *jdig;
287 }
288
289 return idig - oidig;
290}
291
Damien George438c88d2014-02-22 19:25:23 +0000292/* computes i = i * d1 + d2
293 returns number of digits in i
294 assumes enough memory in i; assumes normalised i; assumes dmul != 0
295*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100296STATIC 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 +0000297 mpz_dig_t *oidig = idig;
298 mpz_dbl_dig_t carry = dadd;
299
300 for (; ilen > 0; --ilen, ++idig) {
Damien George9a21d2e2014-09-06 17:15:34 +0100301 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 +0000302 *idig = carry & DIG_MASK;
303 carry >>= DIG_SIZE;
304 }
305
306 if (carry != 0) {
307 *idig++ = carry;
308 }
309
310 return idig - oidig;
311}
312
313/* computes i = j * k
314 returns number of digits in i
315 assumes enough memory in i; assumes i is zeroed; assumes normalised j, k
316 can have j, k point to same memory
317*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100318STATIC 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 +0000319 mpz_dig_t *oidig = idig;
Damien Georgeafb1cf72014-09-05 20:37:06 +0100320 mp_uint_t ilen = 0;
Damien George438c88d2014-02-22 19:25:23 +0000321
322 for (; klen > 0; --klen, ++idig, ++kdig) {
323 mpz_dig_t *id = idig;
324 mpz_dbl_dig_t carry = 0;
325
Damien Georgeafb1cf72014-09-05 20:37:06 +0100326 mp_uint_t jl = jlen;
Damien George438c88d2014-02-22 19:25:23 +0000327 for (mpz_dig_t *jd = jdig; jl > 0; --jl, ++jd, ++id) {
Damien George9a21d2e2014-09-06 17:15:34 +0100328 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 +0000329 *id = carry & DIG_MASK;
330 carry >>= DIG_SIZE;
331 }
332
333 if (carry != 0) {
334 *id++ = carry;
335 }
336
337 ilen = id - oidig;
338 }
339
340 return ilen;
341}
342
343/* natural_div - quo * den + new_num = old_num (ie num is replaced with rem)
344 assumes den != 0
345 assumes num_dig has enough memory to be extended by 1 digit
346 assumes quo_dig has enough memory (as many digits as num)
347 assumes quo_dig is filled with zeros
348 modifies den_dig memory, but restors it to original state at end
349*/
350
Damien George40f3c022014-07-03 13:25:24 +0100351STATIC 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 +0000352 mpz_dig_t *orig_num_dig = num_dig;
353 mpz_dig_t *orig_quo_dig = quo_dig;
354 mpz_dig_t norm_shift = 0;
355 mpz_dbl_dig_t lead_den_digit;
356
357 // handle simple cases
358 {
Damien George42f3de92014-10-03 17:44:14 +0000359 int cmp = mpn_cmp(num_dig, *num_len, den_dig, den_len);
Damien George438c88d2014-02-22 19:25:23 +0000360 if (cmp == 0) {
361 *num_len = 0;
362 quo_dig[0] = 1;
363 *quo_len = 1;
364 return;
365 } else if (cmp < 0) {
366 // numerator remains the same
367 *quo_len = 0;
368 return;
369 }
370 }
371
372 // count number of leading zeros in leading digit of denominator
373 {
374 mpz_dig_t d = den_dig[den_len - 1];
Damien George9a21d2e2014-09-06 17:15:34 +0100375 while ((d & DIG_MSB) == 0) {
Damien George438c88d2014-02-22 19:25:23 +0000376 d <<= 1;
377 ++norm_shift;
378 }
379 }
380
381 // normalise denomenator (leading bit of leading digit is 1)
382 for (mpz_dig_t *den = den_dig, carry = 0; den < den_dig + den_len; ++den) {
383 mpz_dig_t d = *den;
384 *den = ((d << norm_shift) | carry) & DIG_MASK;
385 carry = d >> (DIG_SIZE - norm_shift);
386 }
387
388 // now need to shift numerator by same amount as denominator
389 // first, increase length of numerator in case we need more room to shift
390 num_dig[*num_len] = 0;
391 ++(*num_len);
392 for (mpz_dig_t *num = num_dig, carry = 0; num < num_dig + *num_len; ++num) {
393 mpz_dig_t n = *num;
394 *num = ((n << norm_shift) | carry) & DIG_MASK;
395 carry = n >> (DIG_SIZE - norm_shift);
396 }
397
398 // cache the leading digit of the denominator
399 lead_den_digit = den_dig[den_len - 1];
400
401 // point num_dig to last digit in numerator
402 num_dig += *num_len - 1;
403
404 // calculate number of digits in quotient
405 *quo_len = *num_len - den_len;
406
407 // point to last digit to store for quotient
408 quo_dig += *quo_len - 1;
409
410 // keep going while we have enough digits to divide
411 while (*num_len > den_len) {
Damien George9a21d2e2014-09-06 17:15:34 +0100412 mpz_dbl_dig_t quo = ((mpz_dbl_dig_t)*num_dig << DIG_SIZE) | num_dig[-1];
Damien George438c88d2014-02-22 19:25:23 +0000413
414 // get approximate quotient
415 quo /= lead_den_digit;
416
Damien George9a21d2e2014-09-06 17:15:34 +0100417 // Multiply quo by den and subtract from num to get remainder.
418 // We have different code here to handle different compile-time
419 // configurations of mpz:
420 //
421 // 1. DIG_SIZE is stricly less than half the number of bits
422 // available in mpz_dbl_dig_t. In this case we can use a
423 // slightly more optimal (in time and space) routine that
424 // uses the extra bits in mpz_dbl_dig_signed_t to store a
425 // sign bit.
426 //
427 // 2. DIG_SIZE is exactly half the number of bits available in
428 // mpz_dbl_dig_t. In this (common) case we need to be careful
429 // not to overflow the borrow variable. And the shifting of
430 // borrow needs some special logic (it's a shift right with
431 // round up).
432
433 if (DIG_SIZE < 8 * sizeof(mpz_dbl_dig_t) / 2) {
Damien George438c88d2014-02-22 19:25:23 +0000434 mpz_dbl_dig_signed_t borrow = 0;
435
436 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
Damien George9a21d2e2014-09-06 17:15:34 +0100437 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 +0000438 *n = borrow & DIG_MASK;
439 borrow >>= DIG_SIZE;
440 }
Damien George9a21d2e2014-09-06 17:15:34 +0100441 borrow += *num_dig; // will overflow if DIG_SIZE >= 8*sizeof(mpz_dbl_dig_t)/2
Damien George438c88d2014-02-22 19:25:23 +0000442 *num_dig = borrow & DIG_MASK;
443 borrow >>= DIG_SIZE;
444
445 // adjust quotient if it is too big
446 for (; borrow != 0; --quo) {
447 mpz_dbl_dig_t carry = 0;
448 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
Damien George9a21d2e2014-09-06 17:15:34 +0100449 carry += (mpz_dbl_dig_t)*n + (mpz_dbl_dig_t)*d;
Damien George438c88d2014-02-22 19:25:23 +0000450 *n = carry & DIG_MASK;
451 carry >>= DIG_SIZE;
452 }
453 carry += *num_dig;
454 *num_dig = carry & DIG_MASK;
455 carry >>= DIG_SIZE;
456
457 borrow += carry;
458 }
Damien George9a21d2e2014-09-06 17:15:34 +0100459 } else { // DIG_SIZE == 8 * sizeof(mpz_dbl_dig_t) / 2
460 mpz_dbl_dig_t borrow = 0;
461
462 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
463 mpz_dbl_dig_t x = (mpz_dbl_dig_t)quo * (mpz_dbl_dig_t)(*d);
464 if (x >= *n || *n - x <= borrow) {
465 borrow += (mpz_dbl_dig_t)x - (mpz_dbl_dig_t)*n;
466 *n = (-borrow) & DIG_MASK;
467 borrow = (borrow >> DIG_SIZE) + ((borrow & DIG_MASK) == 0 ? 0 : 1); // shift-right with round-up
468 } else {
469 *n = ((mpz_dbl_dig_t)*n - (mpz_dbl_dig_t)x - (mpz_dbl_dig_t)borrow) & DIG_MASK;
470 borrow = 0;
471 }
472 }
473 if (borrow >= *num_dig) {
474 borrow -= (mpz_dbl_dig_t)*num_dig;
475 *num_dig = (-borrow) & DIG_MASK;
476 borrow = (borrow >> DIG_SIZE) + ((borrow & DIG_MASK) == 0 ? 0 : 1); // shift-right with round-up
477 } else {
478 *num_dig = (*num_dig - borrow) & DIG_MASK;
479 borrow = 0;
480 }
481
482 // adjust quotient if it is too big
483 for (; borrow != 0; --quo) {
484 mpz_dbl_dig_t carry = 0;
485 for (mpz_dig_t *n = num_dig - den_len, *d = den_dig; n < num_dig; ++n, ++d) {
486 carry += (mpz_dbl_dig_t)*n + (mpz_dbl_dig_t)*d;
487 *n = carry & DIG_MASK;
488 carry >>= DIG_SIZE;
489 }
490 carry += (mpz_dbl_dig_t)*num_dig;
491 *num_dig = carry & DIG_MASK;
492 carry >>= DIG_SIZE;
493
494 //assert(borrow >= carry); // enable this to check the logic
495 borrow -= carry;
496 }
Damien George438c88d2014-02-22 19:25:23 +0000497 }
498
499 // store this digit of the quotient
500 *quo_dig = quo & DIG_MASK;
501 --quo_dig;
502
503 // move down to next digit of numerator
504 --num_dig;
505 --(*num_len);
506 }
507
508 // unnormalise denomenator
509 for (mpz_dig_t *den = den_dig + den_len - 1, carry = 0; den >= den_dig; --den) {
510 mpz_dig_t d = *den;
511 *den = ((d >> norm_shift) | carry) & DIG_MASK;
512 carry = d << (DIG_SIZE - norm_shift);
513 }
514
515 // unnormalise numerator (remainder now)
516 for (mpz_dig_t *num = orig_num_dig + *num_len - 1, carry = 0; num >= orig_num_dig; --num) {
517 mpz_dig_t n = *num;
518 *num = ((n >> norm_shift) | carry) & DIG_MASK;
519 carry = n << (DIG_SIZE - norm_shift);
520 }
521
522 // strip trailing zeros
523
524 while (*quo_len > 0 && orig_quo_dig[*quo_len - 1] == 0) {
525 --(*quo_len);
526 }
527
528 while (*num_len > 0 && orig_num_dig[*num_len - 1] == 0) {
529 --(*num_len);
530 }
531}
532
Damien George06201ff2014-03-01 19:50:50 +0000533#define MIN_ALLOC (2)
Damien George438c88d2014-02-22 19:25:23 +0000534
Damien George1c70cbf2014-08-30 00:38:16 +0100535STATIC const uint8_t log_base2_floor[] = {
Damien George438c88d2014-02-22 19:25:23 +0000536 0,
537 0, 1, 1, 2,
538 2, 2, 2, 3,
539 3, 3, 3, 3,
540 3, 3, 3, 4,
541 4, 4, 4, 4,
542 4, 4, 4, 4,
543 4, 4, 4, 4,
544 4, 4, 4, 5
545};
546
Damien George438c88d2014-02-22 19:25:23 +0000547void mpz_init_zero(mpz_t *z) {
Damien George438c88d2014-02-22 19:25:23 +0000548 z->neg = 0;
Damien George06201ff2014-03-01 19:50:50 +0000549 z->fixed_dig = 0;
550 z->alloc = 0;
Damien George438c88d2014-02-22 19:25:23 +0000551 z->len = 0;
552 z->dig = NULL;
553}
554
Damien George40f3c022014-07-03 13:25:24 +0100555void mpz_init_from_int(mpz_t *z, mp_int_t val) {
Damien George438c88d2014-02-22 19:25:23 +0000556 mpz_init_zero(z);
557 mpz_set_from_int(z, val);
558}
559
Damien Georgeafb1cf72014-09-05 20:37:06 +0100560void 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 +0000561 z->neg = 0;
562 z->fixed_dig = 1;
563 z->alloc = alloc;
564 z->len = 0;
565 z->dig = dig;
566 mpz_set_from_int(z, val);
567}
568
Damien George438c88d2014-02-22 19:25:23 +0000569void mpz_deinit(mpz_t *z) {
Damien George06201ff2014-03-01 19:50:50 +0000570 if (z != NULL && !z->fixed_dig) {
Damien George438c88d2014-02-22 19:25:23 +0000571 m_del(mpz_dig_t, z->dig, z->alloc);
572 }
573}
574
575mpz_t *mpz_zero(void) {
576 mpz_t *z = m_new_obj(mpz_t);
577 mpz_init_zero(z);
578 return z;
579}
580
Damien George40f3c022014-07-03 13:25:24 +0100581mpz_t *mpz_from_int(mp_int_t val) {
Damien George438c88d2014-02-22 19:25:23 +0000582 mpz_t *z = mpz_zero();
583 mpz_set_from_int(z, val);
584 return z;
585}
586
Damien George95307432014-09-10 22:10:33 +0100587mpz_t *mpz_from_ll(long long val, bool is_signed) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000588 mpz_t *z = mpz_zero();
Damien George95307432014-09-10 22:10:33 +0100589 mpz_set_from_ll(z, val, is_signed);
Damien Georgebb4a43f2014-03-12 15:36:06 +0000590 return z;
591}
592
David Steinberg6e0b6d02015-01-02 12:39:22 +0000593#if MICROPY_PY_BUILTINS_FLOAT
594mpz_t *mpz_from_float(mp_float_t val) {
595 mpz_t *z = mpz_zero();
596 mpz_set_from_float(z, val);
597 return z;
598}
599#endif
600
Damien Georgeafb1cf72014-09-05 20:37:06 +0100601mpz_t *mpz_from_str(const char *str, mp_uint_t len, bool neg, mp_uint_t base) {
Damien George438c88d2014-02-22 19:25:23 +0000602 mpz_t *z = mpz_zero();
603 mpz_set_from_str(z, str, len, neg, base);
604 return z;
605}
606
607void mpz_free(mpz_t *z) {
608 if (z != NULL) {
609 m_del(mpz_dig_t, z->dig, z->alloc);
610 m_del_obj(mpz_t, z);
611 }
612}
613
Damien Georgeafb1cf72014-09-05 20:37:06 +0100614STATIC void mpz_need_dig(mpz_t *z, mp_uint_t need) {
Damien George438c88d2014-02-22 19:25:23 +0000615 if (need < MIN_ALLOC) {
Damien George06201ff2014-03-01 19:50:50 +0000616 need = MIN_ALLOC;
Damien George438c88d2014-02-22 19:25:23 +0000617 }
618
Damien George06201ff2014-03-01 19:50:50 +0000619 if (z->dig == NULL || z->alloc < need) {
620 if (z->fixed_dig) {
621 // cannot reallocate fixed buffers
622 assert(0);
623 return;
624 }
625 z->dig = m_renew(mpz_dig_t, z->dig, z->alloc, need);
626 z->alloc = need;
Damien George438c88d2014-02-22 19:25:23 +0000627 }
628}
629
630mpz_t *mpz_clone(const mpz_t *src) {
631 mpz_t *z = m_new_obj(mpz_t);
Damien George438c88d2014-02-22 19:25:23 +0000632 z->neg = src->neg;
Damien George06201ff2014-03-01 19:50:50 +0000633 z->fixed_dig = 0;
634 z->alloc = src->alloc;
Damien George438c88d2014-02-22 19:25:23 +0000635 z->len = src->len;
636 if (src->dig == NULL) {
637 z->dig = NULL;
638 } else {
639 z->dig = m_new(mpz_dig_t, z->alloc);
640 memcpy(z->dig, src->dig, src->alloc * sizeof(mpz_dig_t));
641 }
642 return z;
643}
644
Damien George06201ff2014-03-01 19:50:50 +0000645/* sets dest = src
646 can have dest, src the same
647*/
Damien George438c88d2014-02-22 19:25:23 +0000648void mpz_set(mpz_t *dest, const mpz_t *src) {
649 mpz_need_dig(dest, src->len);
650 dest->neg = src->neg;
651 dest->len = src->len;
652 memcpy(dest->dig, src->dig, src->len * sizeof(mpz_dig_t));
653}
654
Damien George40f3c022014-07-03 13:25:24 +0100655void mpz_set_from_int(mpz_t *z, mp_int_t val) {
Damien George58056b02015-01-09 20:58:58 +0000656 if (val == 0) {
657 z->len = 0;
658 return;
659 }
660
Damien George06201ff2014-03-01 19:50:50 +0000661 mpz_need_dig(z, MPZ_NUM_DIG_FOR_INT);
Damien George438c88d2014-02-22 19:25:23 +0000662
Damien George40f3c022014-07-03 13:25:24 +0100663 mp_uint_t uval;
Damien George438c88d2014-02-22 19:25:23 +0000664 if (val < 0) {
665 z->neg = 1;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000666 uval = -val;
Damien George438c88d2014-02-22 19:25:23 +0000667 } else {
668 z->neg = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000669 uval = val;
Damien George438c88d2014-02-22 19:25:23 +0000670 }
671
672 z->len = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000673 while (uval > 0) {
674 z->dig[z->len++] = uval & DIG_MASK;
675 uval >>= DIG_SIZE;
676 }
677}
678
Damien George95307432014-09-10 22:10:33 +0100679void mpz_set_from_ll(mpz_t *z, long long val, bool is_signed) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000680 mpz_need_dig(z, MPZ_NUM_DIG_FOR_LL);
681
682 unsigned long long uval;
Damien George95307432014-09-10 22:10:33 +0100683 if (is_signed && val < 0) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000684 z->neg = 1;
685 uval = -val;
686 } else {
687 z->neg = 0;
688 uval = val;
689 }
690
691 z->len = 0;
692 while (uval > 0) {
693 z->dig[z->len++] = uval & DIG_MASK;
694 uval >>= DIG_SIZE;
Damien George438c88d2014-02-22 19:25:23 +0000695 }
696}
697
David Steinberg6e0b6d02015-01-02 12:39:22 +0000698#if MICROPY_PY_BUILTINS_FLOAT
699void mpz_set_from_float(mpz_t *z, mp_float_t src) {
700#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
David Steinberg6e0b6d02015-01-02 12:39:22 +0000701typedef uint64_t mp_float_int_t;
David Steinbergc585ad12015-01-13 15:19:37 +0000702#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
David Steinberg6e0b6d02015-01-02 12:39:22 +0000703typedef uint32_t mp_float_int_t;
704#endif
705 union {
706 mp_float_t f;
David Steinbergc585ad12015-01-13 15:19:37 +0000707 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 +0000708 } u = {src};
709
710 z->neg = u.p.sgn;
711 if (u.p.exp == 0) {
712 // value == 0 || value < 1
Damien George58056b02015-01-09 20:58:58 +0000713 mpz_set_from_int(z, 0);
David Steinbergc585ad12015-01-13 15:19:37 +0000714 } else if (u.p.exp == ((1 << MP_FLOAT_EXP_BITS) - 1)) {
Damien George6fd4b362015-01-02 23:04:09 +0000715 // u.p.frc == 0 indicates inf, else NaN
716 // should be handled by caller
Damien George58056b02015-01-09 20:58:58 +0000717 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000718 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000719 const int adj_exp = (int)u.p.exp - MP_FLOAT_EXP_BIAS;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000720 if (adj_exp < 0) {
721 // value < 1 , truncates to 0
Damien George58056b02015-01-09 20:58:58 +0000722 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000723 } else if (adj_exp == 0) {
724 // 1 <= value < 2 , so truncates to 1
Damien George58056b02015-01-09 20:58:58 +0000725 mpz_set_from_int(z, 1);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000726 } else {
727 // 2 <= value
728 const int dig_cnt = (adj_exp + 1 + (DIG_SIZE - 1)) / DIG_SIZE;
729 const unsigned int rem = adj_exp % DIG_SIZE;
730 int dig_ind, shft;
David Steinbergc585ad12015-01-13 15:19:37 +0000731 mp_float_int_t frc = u.p.frc | ((mp_float_int_t)1 << MP_FLOAT_FRAC_BITS);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000732
David Steinbergc585ad12015-01-13 15:19:37 +0000733 if (adj_exp < MP_FLOAT_FRAC_BITS) {
David Steinberg6e0b6d02015-01-02 12:39:22 +0000734 shft = 0;
735 dig_ind = 0;
David Steinbergc585ad12015-01-13 15:19:37 +0000736 frc >>= MP_FLOAT_FRAC_BITS - adj_exp;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000737 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000738 shft = (rem - MP_FLOAT_FRAC_BITS) % DIG_SIZE;
739 dig_ind = (adj_exp - MP_FLOAT_FRAC_BITS) / DIG_SIZE;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000740 }
741 mpz_need_dig(z, dig_cnt);
742 z->len = dig_cnt;
743 if (dig_ind != 0) {
744 memset(z->dig, 0, dig_ind * sizeof(mpz_dig_t));
745 }
746 if (shft != 0) {
747 z->dig[dig_ind++] = (frc << shft) & DIG_MASK;
748 frc >>= DIG_SIZE - shft;
749 }
David Steinberg8d427b72015-01-13 15:20:32 +0000750#if DIG_SIZE < (MP_FLOAT_FRAC_BITS + 1)
David Steinberg6e0b6d02015-01-02 12:39:22 +0000751 while (dig_ind != dig_cnt) {
752 z->dig[dig_ind++] = frc & DIG_MASK;
753 frc >>= DIG_SIZE;
754 }
David Steinberg8d427b72015-01-13 15:20:32 +0000755#else
756 if (dig_ind != dig_cnt) {
757 z->dig[dig_ind] = frc;
758 }
759#endif
David Steinberg6e0b6d02015-01-02 12:39:22 +0000760 }
761 }
762}
David Steinberg6e0b6d02015-01-02 12:39:22 +0000763#endif
764
Damien George438c88d2014-02-22 19:25:23 +0000765// returns number of bytes from str that were processed
Damien Georgeafb1cf72014-09-05 20:37:06 +0100766mp_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 +0000767 assert(base < 36);
768
769 const char *cur = str;
770 const char *top = str + len;
771
772 mpz_need_dig(z, len * 8 / DIG_SIZE + 1);
773
774 if (neg) {
775 z->neg = 1;
776 } else {
777 z->neg = 0;
778 }
779
780 z->len = 0;
781 for (; cur < top; ++cur) { // XXX UTF8 next char
Damien Georgeafb1cf72014-09-05 20:37:06 +0100782 //mp_uint_t v = char_to_numeric(cur#); // XXX UTF8 get char
783 mp_uint_t v = *cur;
Damien George438c88d2014-02-22 19:25:23 +0000784 if ('0' <= v && v <= '9') {
785 v -= '0';
786 } else if ('A' <= v && v <= 'Z') {
787 v -= 'A' - 10;
788 } else if ('a' <= v && v <= 'z') {
789 v -= 'a' - 10;
790 } else {
791 break;
792 }
793 if (v >= base) {
794 break;
795 }
796 z->len = mpn_mul_dig_add_dig(z->dig, z->len, base, v);
797 }
798
799 return cur - str;
800}
801
802bool mpz_is_zero(const mpz_t *z) {
803 return z->len == 0;
804}
805
806bool mpz_is_pos(const mpz_t *z) {
807 return z->len > 0 && z->neg == 0;
808}
809
810bool mpz_is_neg(const mpz_t *z) {
811 return z->len > 0 && z->neg != 0;
812}
813
814bool mpz_is_odd(const mpz_t *z) {
815 return z->len > 0 && (z->dig[0] & 1) != 0;
816}
817
818bool mpz_is_even(const mpz_t *z) {
819 return z->len == 0 || (z->dig[0] & 1) == 0;
820}
821
Damien George42f3de92014-10-03 17:44:14 +0000822int mpz_cmp(const mpz_t *z1, const mpz_t *z2) {
Damien Georgea9dc9b82015-01-27 17:47:38 +0000823 // to catch comparison of -0 with +0
824 if (z1->len == 0 && z2->len == 0) {
825 return 0;
826 }
Damien George42f3de92014-10-03 17:44:14 +0000827 int cmp = (int)z2->neg - (int)z1->neg;
Damien George438c88d2014-02-22 19:25:23 +0000828 if (cmp != 0) {
829 return cmp;
830 }
831 cmp = mpn_cmp(z1->dig, z1->len, z2->dig, z2->len);
832 if (z1->neg != 0) {
833 cmp = -cmp;
834 }
835 return cmp;
836}
837
Damien George06201ff2014-03-01 19:50:50 +0000838#if 0
839// obsolete
840// compares mpz with an integer that fits within DIG_SIZE bits
Damien Georgeafb1cf72014-09-05 20:37:06 +0100841mp_int_t mpz_cmp_sml_int(const mpz_t *z, mp_int_t sml_int) {
842 mp_int_t cmp;
Damien George438c88d2014-02-22 19:25:23 +0000843 if (z->neg == 0) {
844 if (sml_int < 0) return 1;
845 if (sml_int == 0) {
846 if (z->len == 0) return 0;
847 return 1;
848 }
849 if (z->len == 0) return -1;
850 assert(sml_int < (1 << DIG_SIZE));
851 if (z->len != 1) return 1;
852 cmp = z->dig[0] - sml_int;
853 } else {
854 if (sml_int > 0) return -1;
855 if (sml_int == 0) {
856 if (z->len == 0) return 0;
857 return -1;
858 }
859 if (z->len == 0) return 1;
860 assert(sml_int > -(1 << DIG_SIZE));
861 if (z->len != 1) return -1;
862 cmp = -z->dig[0] - sml_int;
863 }
864 if (cmp < 0) return -1;
865 if (cmp > 0) return 1;
866 return 0;
867}
Damien George06201ff2014-03-01 19:50:50 +0000868#endif
Damien George438c88d2014-02-22 19:25:23 +0000869
Damien George438c88d2014-02-22 19:25:23 +0000870#if 0
871these functions are unused
872
873/* returns abs(z)
874*/
875mpz_t *mpz_abs(const mpz_t *z) {
876 mpz_t *z2 = mpz_clone(z);
877 z2->neg = 0;
878 return z2;
879}
880
881/* returns -z
882*/
883mpz_t *mpz_neg(const mpz_t *z) {
884 mpz_t *z2 = mpz_clone(z);
885 z2->neg = 1 - z2->neg;
886 return z2;
887}
888
889/* returns lhs + rhs
890 can have lhs, rhs the same
891*/
892mpz_t *mpz_add(const mpz_t *lhs, const mpz_t *rhs) {
893 mpz_t *z = mpz_zero();
894 mpz_add_inpl(z, lhs, rhs);
895 return z;
896}
897
898/* returns lhs - rhs
899 can have lhs, rhs the same
900*/
901mpz_t *mpz_sub(const mpz_t *lhs, const mpz_t *rhs) {
902 mpz_t *z = mpz_zero();
903 mpz_sub_inpl(z, lhs, rhs);
904 return z;
905}
906
907/* returns lhs * rhs
908 can have lhs, rhs the same
909*/
910mpz_t *mpz_mul(const mpz_t *lhs, const mpz_t *rhs) {
911 mpz_t *z = mpz_zero();
912 mpz_mul_inpl(z, lhs, rhs);
913 return z;
914}
915
916/* returns lhs ** rhs
917 can have lhs, rhs the same
918*/
919mpz_t *mpz_pow(const mpz_t *lhs, const mpz_t *rhs) {
920 mpz_t *z = mpz_zero();
921 mpz_pow_inpl(z, lhs, rhs);
922 return z;
923}
924#endif
925
926/* computes dest = abs(z)
927 can have dest, z the same
928*/
929void mpz_abs_inpl(mpz_t *dest, const mpz_t *z) {
930 if (dest != z) {
931 mpz_set(dest, z);
932 }
933 dest->neg = 0;
934}
935
936/* computes dest = -z
937 can have dest, z the same
938*/
939void mpz_neg_inpl(mpz_t *dest, const mpz_t *z) {
940 if (dest != z) {
941 mpz_set(dest, z);
942 }
943 dest->neg = 1 - dest->neg;
944}
945
Damien George06201ff2014-03-01 19:50:50 +0000946/* computes dest = ~z (= -z - 1)
947 can have dest, z the same
948*/
949void mpz_not_inpl(mpz_t *dest, const mpz_t *z) {
950 if (dest != z) {
951 mpz_set(dest, z);
952 }
Damien Georgee0ac1942014-12-31 19:35:01 +0000953 if (dest->len == 0) {
954 mpz_need_dig(dest, 1);
955 dest->dig[0] = 1;
956 dest->len = 1;
957 dest->neg = 1;
958 } else if (dest->neg) {
Damien George06201ff2014-03-01 19:50:50 +0000959 dest->neg = 0;
960 mpz_dig_t k = 1;
961 dest->len = mpn_sub(dest->dig, dest->dig, dest->len, &k, 1);
962 } else {
Damien Georgee0ac1942014-12-31 19:35:01 +0000963 mpz_need_dig(dest, dest->len + 1);
Damien George06201ff2014-03-01 19:50:50 +0000964 mpz_dig_t k = 1;
965 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &k, 1);
966 dest->neg = 1;
967 }
968}
969
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000970/* computes dest = lhs << rhs
971 can have dest, lhs the same
972*/
Damien George40f3c022014-07-03 13:25:24 +0100973void mpz_shl_inpl(mpz_t *dest, const mpz_t *lhs, mp_int_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +0000974 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000975 mpz_set(dest, lhs);
Damien George06201ff2014-03-01 19:50:50 +0000976 } else if (rhs < 0) {
977 mpz_shr_inpl(dest, lhs, -rhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000978 } else {
Damien George06201ff2014-03-01 19:50:50 +0000979 mpz_need_dig(dest, lhs->len + (rhs + DIG_SIZE - 1) / DIG_SIZE);
980 dest->len = mpn_shl(dest->dig, lhs->dig, lhs->len, rhs);
981 dest->neg = lhs->neg;
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000982 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000983}
984
985/* computes dest = lhs >> rhs
986 can have dest, lhs the same
987*/
Damien George40f3c022014-07-03 13:25:24 +0100988void mpz_shr_inpl(mpz_t *dest, const mpz_t *lhs, mp_int_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +0000989 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000990 mpz_set(dest, lhs);
Damien George06201ff2014-03-01 19:50:50 +0000991 } else if (rhs < 0) {
992 mpz_shl_inpl(dest, lhs, -rhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000993 } else {
Damien George06201ff2014-03-01 19:50:50 +0000994 mpz_need_dig(dest, lhs->len);
995 dest->len = mpn_shr(dest->dig, lhs->dig, lhs->len, rhs);
996 dest->neg = lhs->neg;
997 if (dest->neg) {
998 // arithmetic shift right, rounding to negative infinity
Damien Georgeafb1cf72014-09-05 20:37:06 +0100999 mp_uint_t n_whole = rhs / DIG_SIZE;
1000 mp_uint_t n_part = rhs % DIG_SIZE;
Damien George06201ff2014-03-01 19:50:50 +00001001 mpz_dig_t round_up = 0;
Damien Georgeafb1cf72014-09-05 20:37:06 +01001002 for (mp_uint_t i = 0; i < lhs->len && i < n_whole; i++) {
Damien George06201ff2014-03-01 19:50:50 +00001003 if (lhs->dig[i] != 0) {
1004 round_up = 1;
1005 break;
1006 }
1007 }
1008 if (n_whole < lhs->len && (lhs->dig[n_whole] & ((1 << n_part) - 1)) != 0) {
1009 round_up = 1;
1010 }
1011 if (round_up) {
Damien Georgee0ac1942014-12-31 19:35:01 +00001012 if (dest->len == 0) {
1013 // dest == 0, so need to add 1 by hand (answer will be -1)
1014 dest->dig[0] = 1;
1015 dest->len = 1;
1016 } else {
1017 // dest > 0, so can use mpn_add to add 1
1018 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &round_up, 1);
1019 }
Damien George06201ff2014-03-01 19:50:50 +00001020 }
1021 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001022 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001023}
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001024
Damien George438c88d2014-02-22 19:25:23 +00001025/* computes dest = lhs + rhs
1026 can have dest, lhs, rhs the same
1027*/
1028void mpz_add_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1029 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1030 const mpz_t *temp = lhs;
1031 lhs = rhs;
1032 rhs = temp;
1033 }
1034
1035 if (lhs->neg == rhs->neg) {
1036 mpz_need_dig(dest, lhs->len + 1);
1037 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1038 } else {
1039 mpz_need_dig(dest, lhs->len);
1040 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1041 }
1042
1043 dest->neg = lhs->neg;
1044}
1045
1046/* computes dest = lhs - rhs
1047 can have dest, lhs, rhs the same
1048*/
1049void mpz_sub_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1050 bool neg = false;
1051
1052 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1053 const mpz_t *temp = lhs;
1054 lhs = rhs;
1055 rhs = temp;
1056 neg = true;
1057 }
1058
1059 if (lhs->neg != rhs->neg) {
1060 mpz_need_dig(dest, lhs->len + 1);
1061 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1062 } else {
1063 mpz_need_dig(dest, lhs->len);
1064 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1065 }
1066
1067 if (neg) {
1068 dest->neg = 1 - lhs->neg;
1069 } else {
1070 dest->neg = lhs->neg;
1071 }
1072}
1073
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001074/* computes dest = lhs & rhs
1075 can have dest, lhs, rhs the same
1076*/
1077void mpz_and_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001078 if (lhs->neg == rhs->neg) {
Damien Georgef55cf102014-05-29 15:01:49 +01001079 if (lhs->neg == 0) {
1080 // make sure lhs has the most digits
1081 if (lhs->len < rhs->len) {
1082 const mpz_t *temp = lhs;
1083 lhs = rhs;
1084 rhs = temp;
1085 }
1086 // do the and'ing
1087 mpz_need_dig(dest, rhs->len);
Damien Georgeff8dd3f2015-01-20 12:47:20 +00001088 dest->len = mpn_and(dest->dig, lhs->dig, rhs->dig, rhs->len);
Damien Georgef55cf102014-05-29 15:01:49 +01001089 dest->neg = 0;
1090 } else {
1091 // TODO both args are negative
1092 assert(0);
1093 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001094 } else {
Damien Georgef55cf102014-05-29 15:01:49 +01001095 // args have different sign
1096 // make sure lhs is the positive arg
1097 if (rhs->neg == 0) {
1098 const mpz_t *temp = lhs;
1099 lhs = rhs;
1100 rhs = temp;
1101 }
1102 mpz_need_dig(dest, lhs->len + 1);
1103 dest->len = mpn_and_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1104 assert(dest->len <= dest->alloc);
1105 dest->neg = 0;
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001106 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001107}
1108
1109/* computes dest = lhs | rhs
1110 can have dest, lhs, rhs the same
1111*/
1112void mpz_or_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1113 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1114 const mpz_t *temp = lhs;
1115 lhs = rhs;
1116 rhs = temp;
1117 }
1118
1119 if (lhs->neg == rhs->neg) {
1120 mpz_need_dig(dest, lhs->len);
1121 dest->len = mpn_or(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1122 } else {
1123 mpz_need_dig(dest, lhs->len);
1124 // TODO
1125 assert(0);
1126// dest->len = mpn_or_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1127 }
1128
1129 dest->neg = lhs->neg;
1130}
1131
1132/* computes dest = lhs ^ rhs
1133 can have dest, lhs, rhs the same
1134*/
1135void mpz_xor_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1136 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1137 const mpz_t *temp = lhs;
1138 lhs = rhs;
1139 rhs = temp;
1140 }
1141
1142 if (lhs->neg == rhs->neg) {
1143 mpz_need_dig(dest, lhs->len);
1144 dest->len = mpn_xor(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1145 } else {
1146 mpz_need_dig(dest, lhs->len);
1147 // TODO
1148 assert(0);
1149// dest->len = mpn_xor_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1150 }
1151
1152 dest->neg = 0;
1153}
1154
Damien George438c88d2014-02-22 19:25:23 +00001155/* computes dest = lhs * rhs
1156 can have dest, lhs, rhs the same
1157*/
1158void mpz_mul_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs)
1159{
1160 if (lhs->len == 0 || rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001161 mpz_set_from_int(dest, 0);
1162 return;
Damien George438c88d2014-02-22 19:25:23 +00001163 }
1164
1165 mpz_t *temp = NULL;
1166 if (lhs == dest) {
1167 lhs = temp = mpz_clone(lhs);
1168 if (rhs == dest) {
1169 rhs = lhs;
1170 }
1171 } else if (rhs == dest) {
1172 rhs = temp = mpz_clone(rhs);
1173 }
1174
1175 mpz_need_dig(dest, lhs->len + rhs->len); // min mem l+r-1, max mem l+r
1176 memset(dest->dig, 0, dest->alloc * sizeof(mpz_dig_t));
1177 dest->len = mpn_mul(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1178
1179 if (lhs->neg == rhs->neg) {
1180 dest->neg = 0;
1181 } else {
1182 dest->neg = 1;
1183 }
1184
1185 mpz_free(temp);
1186}
1187
1188/* computes dest = lhs ** rhs
1189 can have dest, lhs, rhs the same
1190*/
1191void mpz_pow_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1192 if (lhs->len == 0 || rhs->neg != 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001193 mpz_set_from_int(dest, 0);
1194 return;
Damien George438c88d2014-02-22 19:25:23 +00001195 }
1196
1197 if (rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001198 mpz_set_from_int(dest, 1);
1199 return;
Damien George438c88d2014-02-22 19:25:23 +00001200 }
1201
1202 mpz_t *x = mpz_clone(lhs);
1203 mpz_t *n = mpz_clone(rhs);
1204
1205 mpz_set_from_int(dest, 1);
1206
1207 while (n->len > 0) {
1208 if (mpz_is_odd(n)) {
1209 mpz_mul_inpl(dest, dest, x);
1210 }
Damien George438c88d2014-02-22 19:25:23 +00001211 n->len = mpn_shr(n->dig, n->dig, n->len, 1);
Damien George5bf565e2014-04-04 00:16:32 +01001212 if (n->len == 0) {
1213 break;
1214 }
1215 mpz_mul_inpl(x, x, x);
Damien George438c88d2014-02-22 19:25:23 +00001216 }
1217
1218 mpz_free(x);
1219 mpz_free(n);
1220}
1221
1222/* computes gcd(z1, z2)
1223 based on Knuth's modified gcd algorithm (I think?)
1224 gcd(z1, z2) >= 0
1225 gcd(0, 0) = 0
1226 gcd(z, 0) = abs(z)
1227*/
1228mpz_t *mpz_gcd(const mpz_t *z1, const mpz_t *z2) {
1229 if (z1->len == 0) {
1230 mpz_t *a = mpz_clone(z2);
1231 a->neg = 0;
1232 return a;
1233 } else if (z2->len == 0) {
1234 mpz_t *a = mpz_clone(z1);
1235 a->neg = 0;
1236 return a;
1237 }
1238
1239 mpz_t *a = mpz_clone(z1);
1240 mpz_t *b = mpz_clone(z2);
1241 mpz_t c; mpz_init_zero(&c);
1242 a->neg = 0;
1243 b->neg = 0;
1244
1245 for (;;) {
1246 if (mpz_cmp(a, b) < 0) {
1247 if (a->len == 0) {
1248 mpz_free(a);
1249 mpz_deinit(&c);
1250 return b;
1251 }
1252 mpz_t *t = a; a = b; b = t;
1253 }
1254 if (!(b->len >= 2 || (b->len == 1 && b->dig[0] > 1))) { // compute b > 0; could be mpz_cmp_small_int(b, 1) > 0
1255 break;
1256 }
1257 mpz_set(&c, b);
1258 do {
1259 mpz_add_inpl(&c, &c, &c);
1260 } while (mpz_cmp(&c, a) <= 0);
1261 c.len = mpn_shr(c.dig, c.dig, c.len, 1);
1262 mpz_sub_inpl(a, a, &c);
1263 }
1264
1265 mpz_deinit(&c);
1266
1267 if (b->len == 1 && b->dig[0] == 1) { // compute b == 1; could be mpz_cmp_small_int(b, 1) == 0
1268 mpz_free(a);
1269 return b;
1270 } else {
1271 mpz_free(b);
1272 return a;
1273 }
1274}
1275
1276/* computes lcm(z1, z2)
1277 = abs(z1) / gcd(z1, z2) * abs(z2)
1278 lcm(z1, z1) >= 0
1279 lcm(0, 0) = 0
1280 lcm(z, 0) = 0
1281*/
Damien George5d9b8162014-08-07 14:27:48 +00001282mpz_t *mpz_lcm(const mpz_t *z1, const mpz_t *z2) {
1283 if (z1->len == 0 || z2->len == 0) {
1284 return mpz_zero();
stijn01d6be42014-05-05 12:18:27 +02001285 }
Damien George438c88d2014-02-22 19:25:23 +00001286
1287 mpz_t *gcd = mpz_gcd(z1, z2);
1288 mpz_t *quo = mpz_zero();
1289 mpz_t *rem = mpz_zero();
1290 mpz_divmod_inpl(quo, rem, z1, gcd);
1291 mpz_mul_inpl(rem, quo, z2);
1292 mpz_free(gcd);
1293 mpz_free(quo);
1294 rem->neg = 0;
1295 return rem;
1296}
1297
1298/* computes new integers in quo and rem such that:
1299 quo * rhs + rem = lhs
1300 0 <= rem < rhs
1301 can have lhs, rhs the same
1302*/
1303void mpz_divmod(const mpz_t *lhs, const mpz_t *rhs, mpz_t **quo, mpz_t **rem) {
1304 *quo = mpz_zero();
1305 *rem = mpz_zero();
1306 mpz_divmod_inpl(*quo, *rem, lhs, rhs);
1307}
1308
1309/* computes new integers in quo and rem such that:
1310 quo * rhs + rem = lhs
1311 0 <= rem < rhs
1312 can have lhs, rhs the same
1313*/
1314void mpz_divmod_inpl(mpz_t *dest_quo, mpz_t *dest_rem, const mpz_t *lhs, const mpz_t *rhs) {
1315 if (rhs->len == 0) {
1316 mpz_set_from_int(dest_quo, 0);
1317 mpz_set_from_int(dest_rem, 0);
1318 return;
1319 }
1320
1321 mpz_need_dig(dest_quo, lhs->len + 1); // +1 necessary?
1322 memset(dest_quo->dig, 0, (lhs->len + 1) * sizeof(mpz_dig_t));
1323 dest_quo->len = 0;
1324 mpz_need_dig(dest_rem, lhs->len + 1); // +1 necessary?
1325 mpz_set(dest_rem, lhs);
1326 //rhs->dig[rhs->len] = 0;
1327 mpn_div(dest_rem->dig, &dest_rem->len, rhs->dig, rhs->len, dest_quo->dig, &dest_quo->len);
1328
1329 if (lhs->neg != rhs->neg) {
1330 dest_quo->neg = 1;
1331 }
1332}
1333
1334#if 0
1335these functions are unused
1336
1337/* computes floor(lhs / rhs)
1338 can have lhs, rhs the same
1339*/
1340mpz_t *mpz_div(const mpz_t *lhs, const mpz_t *rhs) {
1341 mpz_t *quo = mpz_zero();
1342 mpz_t rem; mpz_init_zero(&rem);
1343 mpz_divmod_inpl(quo, &rem, lhs, rhs);
1344 mpz_deinit(&rem);
1345 return quo;
1346}
1347
1348/* computes lhs % rhs ( >= 0)
1349 can have lhs, rhs the same
1350*/
1351mpz_t *mpz_mod(const mpz_t *lhs, const mpz_t *rhs) {
1352 mpz_t quo; mpz_init_zero(&quo);
1353 mpz_t *rem = mpz_zero();
1354 mpz_divmod_inpl(&quo, rem, lhs, rhs);
1355 mpz_deinit(&quo);
1356 return rem;
1357}
1358#endif
1359
Damien Georgeffe911d2014-07-24 14:21:37 +01001360// must return actual int value if it fits in mp_int_t
1361mp_int_t mpz_hash(const mpz_t *z) {
1362 mp_int_t val = 0;
1363 mpz_dig_t *d = z->dig + z->len;
1364
Damien George58056b02015-01-09 20:58:58 +00001365 while (d-- > z->dig) {
Damien Georgeffe911d2014-07-24 14:21:37 +01001366 val = (val << DIG_SIZE) | *d;
1367 }
1368
1369 if (z->neg != 0) {
1370 val = -val;
1371 }
1372
1373 return val;
1374}
1375
Damien George40f3c022014-07-03 13:25:24 +01001376bool mpz_as_int_checked(const mpz_t *i, mp_int_t *value) {
Damien George963a5a32015-01-16 17:47:07 +00001377 mp_uint_t val = 0;
Damien George8270e382014-04-03 11:00:54 +00001378 mpz_dig_t *d = i->dig + i->len;
1379
Damien George58056b02015-01-09 20:58:58 +00001380 while (d-- > i->dig) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001381 if (val > (~(WORD_MSBIT_HIGH) >> DIG_SIZE)) {
1382 // will overflow
Damien George8270e382014-04-03 11:00:54 +00001383 return false;
1384 }
Damien Georgec9aa58e2014-07-31 13:41:43 +00001385 val = (val << DIG_SIZE) | *d;
Damien George8270e382014-04-03 11:00:54 +00001386 }
1387
1388 if (i->neg != 0) {
1389 val = -val;
1390 }
1391
1392 *value = val;
1393 return true;
1394}
1395
Damien Georgec9aa58e2014-07-31 13:41:43 +00001396bool mpz_as_uint_checked(const mpz_t *i, mp_uint_t *value) {
1397 if (i->neg != 0) {
1398 // can't represent signed values
1399 return false;
1400 }
1401
1402 mp_uint_t val = 0;
1403 mpz_dig_t *d = i->dig + i->len;
1404
Damien George58056b02015-01-09 20:58:58 +00001405 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001406 if (val > (~(WORD_MSBIT_HIGH) >> (DIG_SIZE - 1))) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001407 // will overflow
1408 return false;
1409 }
1410 val = (val << DIG_SIZE) | *d;
1411 }
1412
1413 *value = val;
1414 return true;
1415}
1416
Damien Georgefb510b32014-06-01 13:32:54 +01001417#if MICROPY_PY_BUILTINS_FLOAT
Damien George52608102014-03-08 15:04:54 +00001418mp_float_t mpz_as_float(const mpz_t *i) {
1419 mp_float_t val = 0;
Damien George438c88d2014-02-22 19:25:23 +00001420 mpz_dig_t *d = i->dig + i->len;
1421
Damien George58056b02015-01-09 20:58:58 +00001422 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001423 val = val * DIG_BASE + *d;
Damien George438c88d2014-02-22 19:25:23 +00001424 }
1425
1426 if (i->neg != 0) {
1427 val = -val;
1428 }
1429
1430 return val;
1431}
Damien George52608102014-03-08 15:04:54 +00001432#endif
Damien George438c88d2014-02-22 19:25:23 +00001433
Damien Georgeafb1cf72014-09-05 20:37:06 +01001434mp_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 +00001435 if (base < 2 || base > 32) {
1436 return 0;
1437 }
1438
Damien Georgeafb1cf72014-09-05 20:37:06 +01001439 mp_uint_t num_digits = i->len * DIG_SIZE / log_base2_floor[base] + 1;
1440 mp_uint_t num_commas = comma ? num_digits / 3: 0;
1441 mp_uint_t prefix_len = prefix ? strlen(prefix) : 0;
Dave Hylandsc4029e52014-04-07 11:19:51 -07001442
1443 return num_digits + num_commas + prefix_len + 2; // +1 for sign, +1 for null byte
1444}
1445
Damien Georgeafb1cf72014-09-05 20:37:06 +01001446#if 0
1447this function is unused
1448char *mpz_as_str(const mpz_t *i, mp_uint_t base) {
1449 char *s = m_new(char, mpz_as_str_size(i, base, NULL, '\0'));
1450 mpz_as_str_inpl(i, base, NULL, 'a', '\0', s);
Damien George438c88d2014-02-22 19:25:23 +00001451 return s;
1452}
Damien Georgeafb1cf72014-09-05 20:37:06 +01001453#endif
Damien George438c88d2014-02-22 19:25:23 +00001454
1455// assumes enough space as calculated by mpz_as_str_size
1456// returns length of string, not including null byte
Damien Georgeafb1cf72014-09-05 20:37:06 +01001457mp_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 +00001458 if (str == NULL || base < 2 || base > 32) {
1459 str[0] = 0;
1460 return 0;
1461 }
1462
Damien Georgeafb1cf72014-09-05 20:37:06 +01001463 mp_uint_t ilen = i->len;
Damien George438c88d2014-02-22 19:25:23 +00001464
Dave Hylandsc4029e52014-04-07 11:19:51 -07001465 char *s = str;
Damien George438c88d2014-02-22 19:25:23 +00001466 if (ilen == 0) {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001467 if (prefix) {
1468 while (*prefix)
1469 *s++ = *prefix++;
1470 }
1471 *s++ = '0';
1472 *s = '\0';
1473 return s - str;
Damien George438c88d2014-02-22 19:25:23 +00001474 }
1475
Damien Georgeeec91052014-04-08 23:11:00 +01001476 // make a copy of mpz digits, so we can do the div/mod calculation
Damien George438c88d2014-02-22 19:25:23 +00001477 mpz_dig_t *dig = m_new(mpz_dig_t, ilen);
1478 memcpy(dig, i->dig, ilen * sizeof(mpz_dig_t));
1479
1480 // convert
Dave Hylandsc4029e52014-04-07 11:19:51 -07001481 char *last_comma = str;
Damien George438c88d2014-02-22 19:25:23 +00001482 bool done;
1483 do {
1484 mpz_dig_t *d = dig + ilen;
1485 mpz_dbl_dig_t a = 0;
1486
1487 // compute next remainder
1488 while (--d >= dig) {
1489 a = (a << DIG_SIZE) | *d;
1490 *d = a / base;
1491 a %= base;
1492 }
1493
1494 // convert to character
1495 a += '0';
1496 if (a > '9') {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001497 a += base_char - '9' - 1;
Damien George438c88d2014-02-22 19:25:23 +00001498 }
1499 *s++ = a;
1500
1501 // check if number is zero
1502 done = true;
1503 for (d = dig; d < dig + ilen; ++d) {
1504 if (*d != 0) {
1505 done = false;
1506 break;
1507 }
1508 }
Dave Hylandsc4029e52014-04-07 11:19:51 -07001509 if (comma && (s - last_comma) == 3) {
1510 *s++ = comma;
1511 last_comma = s;
1512 }
1513 }
1514 while (!done);
Damien George438c88d2014-02-22 19:25:23 +00001515
Damien Georgeeec91052014-04-08 23:11:00 +01001516 // free the copy of the digits array
1517 m_del(mpz_dig_t, dig, ilen);
1518
Dave Hylandsc4029e52014-04-07 11:19:51 -07001519 if (prefix) {
1520 const char *p = &prefix[strlen(prefix)];
1521 while (p > prefix) {
1522 *s++ = *--p;
1523 }
1524 }
Damien George438c88d2014-02-22 19:25:23 +00001525 if (i->neg != 0) {
1526 *s++ = '-';
1527 }
1528
1529 // reverse string
1530 for (char *u = str, *v = s - 1; u < v; ++u, --v) {
1531 char temp = *u;
1532 *u = *v;
1533 *v = temp;
1534 }
1535
Dave Hylandsc4029e52014-04-07 11:19:51 -07001536 *s = '\0'; // null termination
Damien George438c88d2014-02-22 19:25:23 +00001537
1538 return s - str;
1539}
1540
1541#endif // MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ