blob: 3bd3090411deeda467d4547c89b7f51163351cb0 [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;
Damien Georgef66ee4d2015-04-22 23:16:49 +0100100 while (jlen != 0 && idig[jlen - 1] == 0) {
Damien George06201ff2014-03-01 19:50:50 +0000101 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
Damien George848dd0e2015-03-12 15:59:40 +0000575#if 0
576these functions are unused
577
Damien George438c88d2014-02-22 19:25:23 +0000578mpz_t *mpz_zero(void) {
579 mpz_t *z = m_new_obj(mpz_t);
580 mpz_init_zero(z);
581 return z;
582}
583
Damien George40f3c022014-07-03 13:25:24 +0100584mpz_t *mpz_from_int(mp_int_t val) {
Damien George438c88d2014-02-22 19:25:23 +0000585 mpz_t *z = mpz_zero();
586 mpz_set_from_int(z, val);
587 return z;
588}
589
Damien George95307432014-09-10 22:10:33 +0100590mpz_t *mpz_from_ll(long long val, bool is_signed) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000591 mpz_t *z = mpz_zero();
Damien George95307432014-09-10 22:10:33 +0100592 mpz_set_from_ll(z, val, is_signed);
Damien Georgebb4a43f2014-03-12 15:36:06 +0000593 return z;
594}
595
David Steinberg6e0b6d02015-01-02 12:39:22 +0000596#if MICROPY_PY_BUILTINS_FLOAT
597mpz_t *mpz_from_float(mp_float_t val) {
598 mpz_t *z = mpz_zero();
599 mpz_set_from_float(z, val);
600 return z;
601}
602#endif
603
Damien Georgeafb1cf72014-09-05 20:37:06 +0100604mpz_t *mpz_from_str(const char *str, mp_uint_t len, bool neg, mp_uint_t base) {
Damien George438c88d2014-02-22 19:25:23 +0000605 mpz_t *z = mpz_zero();
606 mpz_set_from_str(z, str, len, neg, base);
607 return z;
608}
Damien George848dd0e2015-03-12 15:59:40 +0000609#endif
Damien George438c88d2014-02-22 19:25:23 +0000610
Damien George848dd0e2015-03-12 15:59:40 +0000611STATIC void mpz_free(mpz_t *z) {
Damien George438c88d2014-02-22 19:25:23 +0000612 if (z != NULL) {
613 m_del(mpz_dig_t, z->dig, z->alloc);
614 m_del_obj(mpz_t, z);
615 }
616}
617
Damien Georgeafb1cf72014-09-05 20:37:06 +0100618STATIC void mpz_need_dig(mpz_t *z, mp_uint_t need) {
Damien George438c88d2014-02-22 19:25:23 +0000619 if (need < MIN_ALLOC) {
Damien George06201ff2014-03-01 19:50:50 +0000620 need = MIN_ALLOC;
Damien George438c88d2014-02-22 19:25:23 +0000621 }
622
Damien George06201ff2014-03-01 19:50:50 +0000623 if (z->dig == NULL || z->alloc < need) {
624 if (z->fixed_dig) {
625 // cannot reallocate fixed buffers
626 assert(0);
627 return;
628 }
629 z->dig = m_renew(mpz_dig_t, z->dig, z->alloc, need);
630 z->alloc = need;
Damien George438c88d2014-02-22 19:25:23 +0000631 }
632}
633
Damien George848dd0e2015-03-12 15:59:40 +0000634STATIC mpz_t *mpz_clone(const mpz_t *src) {
Damien George438c88d2014-02-22 19:25:23 +0000635 mpz_t *z = m_new_obj(mpz_t);
Damien George438c88d2014-02-22 19:25:23 +0000636 z->neg = src->neg;
Damien George06201ff2014-03-01 19:50:50 +0000637 z->fixed_dig = 0;
638 z->alloc = src->alloc;
Damien George438c88d2014-02-22 19:25:23 +0000639 z->len = src->len;
640 if (src->dig == NULL) {
641 z->dig = NULL;
642 } else {
643 z->dig = m_new(mpz_dig_t, z->alloc);
644 memcpy(z->dig, src->dig, src->alloc * sizeof(mpz_dig_t));
645 }
646 return z;
647}
648
Damien George06201ff2014-03-01 19:50:50 +0000649/* sets dest = src
650 can have dest, src the same
651*/
Damien George438c88d2014-02-22 19:25:23 +0000652void mpz_set(mpz_t *dest, const mpz_t *src) {
653 mpz_need_dig(dest, src->len);
654 dest->neg = src->neg;
655 dest->len = src->len;
656 memcpy(dest->dig, src->dig, src->len * sizeof(mpz_dig_t));
657}
658
Damien George40f3c022014-07-03 13:25:24 +0100659void mpz_set_from_int(mpz_t *z, mp_int_t val) {
Damien George58056b02015-01-09 20:58:58 +0000660 if (val == 0) {
661 z->len = 0;
662 return;
663 }
664
Damien George06201ff2014-03-01 19:50:50 +0000665 mpz_need_dig(z, MPZ_NUM_DIG_FOR_INT);
Damien George438c88d2014-02-22 19:25:23 +0000666
Damien George40f3c022014-07-03 13:25:24 +0100667 mp_uint_t uval;
Damien George438c88d2014-02-22 19:25:23 +0000668 if (val < 0) {
669 z->neg = 1;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000670 uval = -val;
Damien George438c88d2014-02-22 19:25:23 +0000671 } else {
672 z->neg = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000673 uval = val;
Damien George438c88d2014-02-22 19:25:23 +0000674 }
675
676 z->len = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000677 while (uval > 0) {
678 z->dig[z->len++] = uval & DIG_MASK;
679 uval >>= DIG_SIZE;
680 }
681}
682
Damien George95307432014-09-10 22:10:33 +0100683void mpz_set_from_ll(mpz_t *z, long long val, bool is_signed) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000684 mpz_need_dig(z, MPZ_NUM_DIG_FOR_LL);
685
686 unsigned long long uval;
Damien George95307432014-09-10 22:10:33 +0100687 if (is_signed && val < 0) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000688 z->neg = 1;
689 uval = -val;
690 } else {
691 z->neg = 0;
692 uval = val;
693 }
694
695 z->len = 0;
696 while (uval > 0) {
697 z->dig[z->len++] = uval & DIG_MASK;
698 uval >>= DIG_SIZE;
Damien George438c88d2014-02-22 19:25:23 +0000699 }
700}
701
David Steinberg6e0b6d02015-01-02 12:39:22 +0000702#if MICROPY_PY_BUILTINS_FLOAT
703void mpz_set_from_float(mpz_t *z, mp_float_t src) {
704#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
David Steinberg6e0b6d02015-01-02 12:39:22 +0000705typedef uint64_t mp_float_int_t;
David Steinbergc585ad12015-01-13 15:19:37 +0000706#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
David Steinberg6e0b6d02015-01-02 12:39:22 +0000707typedef uint32_t mp_float_int_t;
708#endif
709 union {
710 mp_float_t f;
David Steinbergc585ad12015-01-13 15:19:37 +0000711 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 +0000712 } u = {src};
713
714 z->neg = u.p.sgn;
715 if (u.p.exp == 0) {
716 // value == 0 || value < 1
Damien George58056b02015-01-09 20:58:58 +0000717 mpz_set_from_int(z, 0);
David Steinbergc585ad12015-01-13 15:19:37 +0000718 } else if (u.p.exp == ((1 << MP_FLOAT_EXP_BITS) - 1)) {
Damien George6fd4b362015-01-02 23:04:09 +0000719 // u.p.frc == 0 indicates inf, else NaN
720 // should be handled by caller
Damien George58056b02015-01-09 20:58:58 +0000721 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000722 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000723 const int adj_exp = (int)u.p.exp - MP_FLOAT_EXP_BIAS;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000724 if (adj_exp < 0) {
725 // value < 1 , truncates to 0
Damien George58056b02015-01-09 20:58:58 +0000726 mpz_set_from_int(z, 0);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000727 } else if (adj_exp == 0) {
728 // 1 <= value < 2 , so truncates to 1
Damien George58056b02015-01-09 20:58:58 +0000729 mpz_set_from_int(z, 1);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000730 } else {
731 // 2 <= value
732 const int dig_cnt = (adj_exp + 1 + (DIG_SIZE - 1)) / DIG_SIZE;
733 const unsigned int rem = adj_exp % DIG_SIZE;
734 int dig_ind, shft;
David Steinbergc585ad12015-01-13 15:19:37 +0000735 mp_float_int_t frc = u.p.frc | ((mp_float_int_t)1 << MP_FLOAT_FRAC_BITS);
David Steinberg6e0b6d02015-01-02 12:39:22 +0000736
David Steinbergc585ad12015-01-13 15:19:37 +0000737 if (adj_exp < MP_FLOAT_FRAC_BITS) {
David Steinberg6e0b6d02015-01-02 12:39:22 +0000738 shft = 0;
739 dig_ind = 0;
David Steinbergc585ad12015-01-13 15:19:37 +0000740 frc >>= MP_FLOAT_FRAC_BITS - adj_exp;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000741 } else {
David Steinbergc585ad12015-01-13 15:19:37 +0000742 shft = (rem - MP_FLOAT_FRAC_BITS) % DIG_SIZE;
743 dig_ind = (adj_exp - MP_FLOAT_FRAC_BITS) / DIG_SIZE;
David Steinberg6e0b6d02015-01-02 12:39:22 +0000744 }
745 mpz_need_dig(z, dig_cnt);
746 z->len = dig_cnt;
747 if (dig_ind != 0) {
748 memset(z->dig, 0, dig_ind * sizeof(mpz_dig_t));
749 }
750 if (shft != 0) {
751 z->dig[dig_ind++] = (frc << shft) & DIG_MASK;
752 frc >>= DIG_SIZE - shft;
753 }
David Steinberg8d427b72015-01-13 15:20:32 +0000754#if DIG_SIZE < (MP_FLOAT_FRAC_BITS + 1)
David Steinberg6e0b6d02015-01-02 12:39:22 +0000755 while (dig_ind != dig_cnt) {
756 z->dig[dig_ind++] = frc & DIG_MASK;
757 frc >>= DIG_SIZE;
758 }
David Steinberg8d427b72015-01-13 15:20:32 +0000759#else
760 if (dig_ind != dig_cnt) {
761 z->dig[dig_ind] = frc;
762 }
763#endif
David Steinberg6e0b6d02015-01-02 12:39:22 +0000764 }
765 }
766}
David Steinberg6e0b6d02015-01-02 12:39:22 +0000767#endif
768
Damien George438c88d2014-02-22 19:25:23 +0000769// returns number of bytes from str that were processed
Damien Georgeafb1cf72014-09-05 20:37:06 +0100770mp_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 +0000771 assert(base < 36);
772
773 const char *cur = str;
774 const char *top = str + len;
775
776 mpz_need_dig(z, len * 8 / DIG_SIZE + 1);
777
778 if (neg) {
779 z->neg = 1;
780 } else {
781 z->neg = 0;
782 }
783
784 z->len = 0;
785 for (; cur < top; ++cur) { // XXX UTF8 next char
Damien Georgeafb1cf72014-09-05 20:37:06 +0100786 //mp_uint_t v = char_to_numeric(cur#); // XXX UTF8 get char
787 mp_uint_t v = *cur;
Damien George438c88d2014-02-22 19:25:23 +0000788 if ('0' <= v && v <= '9') {
789 v -= '0';
790 } else if ('A' <= v && v <= 'Z') {
791 v -= 'A' - 10;
792 } else if ('a' <= v && v <= 'z') {
793 v -= 'a' - 10;
794 } else {
795 break;
796 }
797 if (v >= base) {
798 break;
799 }
800 z->len = mpn_mul_dig_add_dig(z->dig, z->len, base, v);
801 }
802
803 return cur - str;
804}
805
806bool mpz_is_zero(const mpz_t *z) {
807 return z->len == 0;
808}
809
Damien Georgea2e38382015-03-02 12:58:06 +0000810#if 0
811these functions are unused
812
Damien George438c88d2014-02-22 19:25:23 +0000813bool mpz_is_pos(const mpz_t *z) {
814 return z->len > 0 && z->neg == 0;
815}
816
817bool mpz_is_neg(const mpz_t *z) {
818 return z->len > 0 && z->neg != 0;
819}
820
821bool mpz_is_odd(const mpz_t *z) {
822 return z->len > 0 && (z->dig[0] & 1) != 0;
823}
824
825bool mpz_is_even(const mpz_t *z) {
826 return z->len == 0 || (z->dig[0] & 1) == 0;
827}
Damien Georgea2e38382015-03-02 12:58:06 +0000828#endif
Damien George438c88d2014-02-22 19:25:23 +0000829
Damien George42f3de92014-10-03 17:44:14 +0000830int mpz_cmp(const mpz_t *z1, const mpz_t *z2) {
Damien Georgea9dc9b82015-01-27 17:47:38 +0000831 // to catch comparison of -0 with +0
832 if (z1->len == 0 && z2->len == 0) {
833 return 0;
834 }
Damien George42f3de92014-10-03 17:44:14 +0000835 int cmp = (int)z2->neg - (int)z1->neg;
Damien George438c88d2014-02-22 19:25:23 +0000836 if (cmp != 0) {
837 return cmp;
838 }
839 cmp = mpn_cmp(z1->dig, z1->len, z2->dig, z2->len);
840 if (z1->neg != 0) {
841 cmp = -cmp;
842 }
843 return cmp;
844}
845
Damien George06201ff2014-03-01 19:50:50 +0000846#if 0
847// obsolete
848// compares mpz with an integer that fits within DIG_SIZE bits
Damien Georgeafb1cf72014-09-05 20:37:06 +0100849mp_int_t mpz_cmp_sml_int(const mpz_t *z, mp_int_t sml_int) {
850 mp_int_t cmp;
Damien George438c88d2014-02-22 19:25:23 +0000851 if (z->neg == 0) {
852 if (sml_int < 0) return 1;
853 if (sml_int == 0) {
854 if (z->len == 0) return 0;
855 return 1;
856 }
857 if (z->len == 0) return -1;
858 assert(sml_int < (1 << DIG_SIZE));
859 if (z->len != 1) return 1;
860 cmp = z->dig[0] - sml_int;
861 } else {
862 if (sml_int > 0) return -1;
863 if (sml_int == 0) {
864 if (z->len == 0) return 0;
865 return -1;
866 }
867 if (z->len == 0) return 1;
868 assert(sml_int > -(1 << DIG_SIZE));
869 if (z->len != 1) return -1;
870 cmp = -z->dig[0] - sml_int;
871 }
872 if (cmp < 0) return -1;
873 if (cmp > 0) return 1;
874 return 0;
875}
Damien George06201ff2014-03-01 19:50:50 +0000876#endif
Damien George438c88d2014-02-22 19:25:23 +0000877
Damien George438c88d2014-02-22 19:25:23 +0000878#if 0
879these functions are unused
880
881/* returns abs(z)
882*/
883mpz_t *mpz_abs(const mpz_t *z) {
884 mpz_t *z2 = mpz_clone(z);
885 z2->neg = 0;
886 return z2;
887}
888
889/* returns -z
890*/
891mpz_t *mpz_neg(const mpz_t *z) {
892 mpz_t *z2 = mpz_clone(z);
893 z2->neg = 1 - z2->neg;
894 return z2;
895}
896
897/* returns lhs + rhs
898 can have lhs, rhs the same
899*/
900mpz_t *mpz_add(const mpz_t *lhs, const mpz_t *rhs) {
901 mpz_t *z = mpz_zero();
902 mpz_add_inpl(z, lhs, rhs);
903 return z;
904}
905
906/* returns lhs - rhs
907 can have lhs, rhs the same
908*/
909mpz_t *mpz_sub(const mpz_t *lhs, const mpz_t *rhs) {
910 mpz_t *z = mpz_zero();
911 mpz_sub_inpl(z, lhs, rhs);
912 return z;
913}
914
915/* returns lhs * rhs
916 can have lhs, rhs the same
917*/
918mpz_t *mpz_mul(const mpz_t *lhs, const mpz_t *rhs) {
919 mpz_t *z = mpz_zero();
920 mpz_mul_inpl(z, lhs, rhs);
921 return z;
922}
923
924/* returns lhs ** rhs
925 can have lhs, rhs the same
926*/
927mpz_t *mpz_pow(const mpz_t *lhs, const mpz_t *rhs) {
928 mpz_t *z = mpz_zero();
929 mpz_pow_inpl(z, lhs, rhs);
930 return z;
931}
Damien Georgea2e38382015-03-02 12:58:06 +0000932
933/* computes new integers in quo and rem such that:
934 quo * rhs + rem = lhs
935 0 <= rem < rhs
936 can have lhs, rhs the same
937*/
938void mpz_divmod(const mpz_t *lhs, const mpz_t *rhs, mpz_t **quo, mpz_t **rem) {
939 *quo = mpz_zero();
940 *rem = mpz_zero();
941 mpz_divmod_inpl(*quo, *rem, lhs, rhs);
942}
Damien George438c88d2014-02-22 19:25:23 +0000943#endif
944
945/* computes dest = abs(z)
946 can have dest, z the same
947*/
948void mpz_abs_inpl(mpz_t *dest, const mpz_t *z) {
949 if (dest != z) {
950 mpz_set(dest, z);
951 }
952 dest->neg = 0;
953}
954
955/* computes dest = -z
956 can have dest, z the same
957*/
958void mpz_neg_inpl(mpz_t *dest, const mpz_t *z) {
959 if (dest != z) {
960 mpz_set(dest, z);
961 }
962 dest->neg = 1 - dest->neg;
963}
964
Damien George06201ff2014-03-01 19:50:50 +0000965/* computes dest = ~z (= -z - 1)
966 can have dest, z the same
967*/
968void mpz_not_inpl(mpz_t *dest, const mpz_t *z) {
969 if (dest != z) {
970 mpz_set(dest, z);
971 }
Damien Georgee0ac1942014-12-31 19:35:01 +0000972 if (dest->len == 0) {
973 mpz_need_dig(dest, 1);
974 dest->dig[0] = 1;
975 dest->len = 1;
976 dest->neg = 1;
977 } else if (dest->neg) {
Damien George06201ff2014-03-01 19:50:50 +0000978 dest->neg = 0;
979 mpz_dig_t k = 1;
980 dest->len = mpn_sub(dest->dig, dest->dig, dest->len, &k, 1);
981 } else {
Damien Georgee0ac1942014-12-31 19:35:01 +0000982 mpz_need_dig(dest, dest->len + 1);
Damien George06201ff2014-03-01 19:50:50 +0000983 mpz_dig_t k = 1;
984 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &k, 1);
985 dest->neg = 1;
986 }
987}
988
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000989/* computes dest = lhs << rhs
990 can have dest, lhs the same
991*/
Damien George40f3c022014-07-03 13:25:24 +0100992void mpz_shl_inpl(mpz_t *dest, const mpz_t *lhs, mp_int_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +0000993 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000994 mpz_set(dest, lhs);
Damien George06201ff2014-03-01 19:50:50 +0000995 } else if (rhs < 0) {
996 mpz_shr_inpl(dest, lhs, -rhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000997 } else {
Damien George06201ff2014-03-01 19:50:50 +0000998 mpz_need_dig(dest, lhs->len + (rhs + DIG_SIZE - 1) / DIG_SIZE);
999 dest->len = mpn_shl(dest->dig, lhs->dig, lhs->len, rhs);
1000 dest->neg = lhs->neg;
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001001 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001002}
1003
1004/* computes dest = lhs >> rhs
1005 can have dest, lhs the same
1006*/
Damien George40f3c022014-07-03 13:25:24 +01001007void mpz_shr_inpl(mpz_t *dest, const mpz_t *lhs, mp_int_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +00001008 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001009 mpz_set(dest, lhs);
Damien George06201ff2014-03-01 19:50:50 +00001010 } else if (rhs < 0) {
1011 mpz_shl_inpl(dest, lhs, -rhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001012 } else {
Damien George06201ff2014-03-01 19:50:50 +00001013 mpz_need_dig(dest, lhs->len);
1014 dest->len = mpn_shr(dest->dig, lhs->dig, lhs->len, rhs);
1015 dest->neg = lhs->neg;
1016 if (dest->neg) {
1017 // arithmetic shift right, rounding to negative infinity
Damien Georgeafb1cf72014-09-05 20:37:06 +01001018 mp_uint_t n_whole = rhs / DIG_SIZE;
1019 mp_uint_t n_part = rhs % DIG_SIZE;
Damien George06201ff2014-03-01 19:50:50 +00001020 mpz_dig_t round_up = 0;
Damien Georgeafb1cf72014-09-05 20:37:06 +01001021 for (mp_uint_t i = 0; i < lhs->len && i < n_whole; i++) {
Damien George06201ff2014-03-01 19:50:50 +00001022 if (lhs->dig[i] != 0) {
1023 round_up = 1;
1024 break;
1025 }
1026 }
1027 if (n_whole < lhs->len && (lhs->dig[n_whole] & ((1 << n_part) - 1)) != 0) {
1028 round_up = 1;
1029 }
1030 if (round_up) {
Damien Georgee0ac1942014-12-31 19:35:01 +00001031 if (dest->len == 0) {
1032 // dest == 0, so need to add 1 by hand (answer will be -1)
1033 dest->dig[0] = 1;
1034 dest->len = 1;
1035 } else {
1036 // dest > 0, so can use mpn_add to add 1
1037 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &round_up, 1);
1038 }
Damien George06201ff2014-03-01 19:50:50 +00001039 }
1040 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001041 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001042}
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001043
Damien George438c88d2014-02-22 19:25:23 +00001044/* computes dest = lhs + rhs
1045 can have dest, lhs, rhs the same
1046*/
1047void mpz_add_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1048 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1049 const mpz_t *temp = lhs;
1050 lhs = rhs;
1051 rhs = temp;
1052 }
1053
1054 if (lhs->neg == rhs->neg) {
1055 mpz_need_dig(dest, lhs->len + 1);
1056 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1057 } else {
1058 mpz_need_dig(dest, lhs->len);
1059 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1060 }
1061
1062 dest->neg = lhs->neg;
1063}
1064
1065/* computes dest = lhs - rhs
1066 can have dest, lhs, rhs the same
1067*/
1068void mpz_sub_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1069 bool neg = false;
1070
1071 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1072 const mpz_t *temp = lhs;
1073 lhs = rhs;
1074 rhs = temp;
1075 neg = true;
1076 }
1077
1078 if (lhs->neg != rhs->neg) {
1079 mpz_need_dig(dest, lhs->len + 1);
1080 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1081 } else {
1082 mpz_need_dig(dest, lhs->len);
1083 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1084 }
1085
1086 if (neg) {
1087 dest->neg = 1 - lhs->neg;
1088 } else {
1089 dest->neg = lhs->neg;
1090 }
1091}
1092
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001093/* computes dest = lhs & rhs
1094 can have dest, lhs, rhs the same
1095*/
1096void mpz_and_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001097 if (lhs->neg == rhs->neg) {
Damien Georgef55cf102014-05-29 15:01:49 +01001098 if (lhs->neg == 0) {
1099 // make sure lhs has the most digits
1100 if (lhs->len < rhs->len) {
1101 const mpz_t *temp = lhs;
1102 lhs = rhs;
1103 rhs = temp;
1104 }
1105 // do the and'ing
1106 mpz_need_dig(dest, rhs->len);
Damien Georgeff8dd3f2015-01-20 12:47:20 +00001107 dest->len = mpn_and(dest->dig, lhs->dig, rhs->dig, rhs->len);
Damien Georgef55cf102014-05-29 15:01:49 +01001108 dest->neg = 0;
1109 } else {
1110 // TODO both args are negative
1111 assert(0);
1112 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001113 } else {
Damien Georgef55cf102014-05-29 15:01:49 +01001114 // args have different sign
1115 // make sure lhs is the positive arg
1116 if (rhs->neg == 0) {
1117 const mpz_t *temp = lhs;
1118 lhs = rhs;
1119 rhs = temp;
1120 }
1121 mpz_need_dig(dest, lhs->len + 1);
1122 dest->len = mpn_and_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1123 assert(dest->len <= dest->alloc);
1124 dest->neg = 0;
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001125 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001126}
1127
1128/* computes dest = lhs | rhs
1129 can have dest, lhs, rhs the same
1130*/
1131void mpz_or_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1132 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1133 const mpz_t *temp = lhs;
1134 lhs = rhs;
1135 rhs = temp;
1136 }
1137
1138 if (lhs->neg == rhs->neg) {
1139 mpz_need_dig(dest, lhs->len);
1140 dest->len = mpn_or(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1141 } else {
1142 mpz_need_dig(dest, lhs->len);
1143 // TODO
1144 assert(0);
1145// dest->len = mpn_or_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1146 }
1147
1148 dest->neg = lhs->neg;
1149}
1150
1151/* computes dest = lhs ^ rhs
1152 can have dest, lhs, rhs the same
1153*/
1154void mpz_xor_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1155 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1156 const mpz_t *temp = lhs;
1157 lhs = rhs;
1158 rhs = temp;
1159 }
1160
1161 if (lhs->neg == rhs->neg) {
1162 mpz_need_dig(dest, lhs->len);
1163 dest->len = mpn_xor(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1164 } else {
1165 mpz_need_dig(dest, lhs->len);
1166 // TODO
1167 assert(0);
1168// dest->len = mpn_xor_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1169 }
1170
1171 dest->neg = 0;
1172}
1173
Damien George438c88d2014-02-22 19:25:23 +00001174/* computes dest = lhs * rhs
1175 can have dest, lhs, rhs the same
1176*/
Damien George4dea9222015-04-09 15:29:54 +00001177void mpz_mul_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Damien George438c88d2014-02-22 19:25:23 +00001178 if (lhs->len == 0 || rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001179 mpz_set_from_int(dest, 0);
1180 return;
Damien George438c88d2014-02-22 19:25:23 +00001181 }
1182
1183 mpz_t *temp = NULL;
1184 if (lhs == dest) {
1185 lhs = temp = mpz_clone(lhs);
1186 if (rhs == dest) {
1187 rhs = lhs;
1188 }
1189 } else if (rhs == dest) {
1190 rhs = temp = mpz_clone(rhs);
1191 }
1192
1193 mpz_need_dig(dest, lhs->len + rhs->len); // min mem l+r-1, max mem l+r
1194 memset(dest->dig, 0, dest->alloc * sizeof(mpz_dig_t));
1195 dest->len = mpn_mul(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1196
1197 if (lhs->neg == rhs->neg) {
1198 dest->neg = 0;
1199 } else {
1200 dest->neg = 1;
1201 }
1202
1203 mpz_free(temp);
1204}
1205
1206/* computes dest = lhs ** rhs
1207 can have dest, lhs, rhs the same
1208*/
1209void mpz_pow_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1210 if (lhs->len == 0 || rhs->neg != 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001211 mpz_set_from_int(dest, 0);
1212 return;
Damien George438c88d2014-02-22 19:25:23 +00001213 }
1214
1215 if (rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001216 mpz_set_from_int(dest, 1);
1217 return;
Damien George438c88d2014-02-22 19:25:23 +00001218 }
1219
1220 mpz_t *x = mpz_clone(lhs);
1221 mpz_t *n = mpz_clone(rhs);
1222
1223 mpz_set_from_int(dest, 1);
1224
1225 while (n->len > 0) {
Damien Georgea2e38382015-03-02 12:58:06 +00001226 if ((n->dig[0] & 1) != 0) {
Damien George438c88d2014-02-22 19:25:23 +00001227 mpz_mul_inpl(dest, dest, x);
1228 }
Damien George438c88d2014-02-22 19:25:23 +00001229 n->len = mpn_shr(n->dig, n->dig, n->len, 1);
Damien George5bf565e2014-04-04 00:16:32 +01001230 if (n->len == 0) {
1231 break;
1232 }
1233 mpz_mul_inpl(x, x, x);
Damien George438c88d2014-02-22 19:25:23 +00001234 }
1235
1236 mpz_free(x);
1237 mpz_free(n);
1238}
1239
Damien Georgea2e38382015-03-02 12:58:06 +00001240#if 0
1241these functions are unused
1242
Damien George438c88d2014-02-22 19:25:23 +00001243/* computes gcd(z1, z2)
1244 based on Knuth's modified gcd algorithm (I think?)
1245 gcd(z1, z2) >= 0
1246 gcd(0, 0) = 0
1247 gcd(z, 0) = abs(z)
1248*/
1249mpz_t *mpz_gcd(const mpz_t *z1, const mpz_t *z2) {
1250 if (z1->len == 0) {
1251 mpz_t *a = mpz_clone(z2);
1252 a->neg = 0;
1253 return a;
1254 } else if (z2->len == 0) {
1255 mpz_t *a = mpz_clone(z1);
1256 a->neg = 0;
1257 return a;
1258 }
1259
1260 mpz_t *a = mpz_clone(z1);
1261 mpz_t *b = mpz_clone(z2);
1262 mpz_t c; mpz_init_zero(&c);
1263 a->neg = 0;
1264 b->neg = 0;
1265
1266 for (;;) {
1267 if (mpz_cmp(a, b) < 0) {
1268 if (a->len == 0) {
1269 mpz_free(a);
1270 mpz_deinit(&c);
1271 return b;
1272 }
1273 mpz_t *t = a; a = b; b = t;
1274 }
1275 if (!(b->len >= 2 || (b->len == 1 && b->dig[0] > 1))) { // compute b > 0; could be mpz_cmp_small_int(b, 1) > 0
1276 break;
1277 }
1278 mpz_set(&c, b);
1279 do {
1280 mpz_add_inpl(&c, &c, &c);
1281 } while (mpz_cmp(&c, a) <= 0);
1282 c.len = mpn_shr(c.dig, c.dig, c.len, 1);
1283 mpz_sub_inpl(a, a, &c);
1284 }
1285
1286 mpz_deinit(&c);
1287
1288 if (b->len == 1 && b->dig[0] == 1) { // compute b == 1; could be mpz_cmp_small_int(b, 1) == 0
1289 mpz_free(a);
1290 return b;
1291 } else {
1292 mpz_free(b);
1293 return a;
1294 }
1295}
1296
1297/* computes lcm(z1, z2)
1298 = abs(z1) / gcd(z1, z2) * abs(z2)
1299 lcm(z1, z1) >= 0
1300 lcm(0, 0) = 0
1301 lcm(z, 0) = 0
1302*/
Damien George5d9b8162014-08-07 14:27:48 +00001303mpz_t *mpz_lcm(const mpz_t *z1, const mpz_t *z2) {
1304 if (z1->len == 0 || z2->len == 0) {
1305 return mpz_zero();
stijn01d6be42014-05-05 12:18:27 +02001306 }
Damien George438c88d2014-02-22 19:25:23 +00001307
1308 mpz_t *gcd = mpz_gcd(z1, z2);
1309 mpz_t *quo = mpz_zero();
1310 mpz_t *rem = mpz_zero();
1311 mpz_divmod_inpl(quo, rem, z1, gcd);
1312 mpz_mul_inpl(rem, quo, z2);
1313 mpz_free(gcd);
1314 mpz_free(quo);
1315 rem->neg = 0;
1316 return rem;
1317}
Damien Georgea2e38382015-03-02 12:58:06 +00001318#endif
Damien George438c88d2014-02-22 19:25:23 +00001319
1320/* computes new integers in quo and rem such that:
1321 quo * rhs + rem = lhs
1322 0 <= rem < rhs
1323 can have lhs, rhs the same
1324*/
1325void mpz_divmod_inpl(mpz_t *dest_quo, mpz_t *dest_rem, const mpz_t *lhs, const mpz_t *rhs) {
1326 if (rhs->len == 0) {
1327 mpz_set_from_int(dest_quo, 0);
1328 mpz_set_from_int(dest_rem, 0);
1329 return;
1330 }
1331
1332 mpz_need_dig(dest_quo, lhs->len + 1); // +1 necessary?
1333 memset(dest_quo->dig, 0, (lhs->len + 1) * sizeof(mpz_dig_t));
1334 dest_quo->len = 0;
1335 mpz_need_dig(dest_rem, lhs->len + 1); // +1 necessary?
1336 mpz_set(dest_rem, lhs);
1337 //rhs->dig[rhs->len] = 0;
1338 mpn_div(dest_rem->dig, &dest_rem->len, rhs->dig, rhs->len, dest_quo->dig, &dest_quo->len);
1339
1340 if (lhs->neg != rhs->neg) {
1341 dest_quo->neg = 1;
1342 }
1343}
1344
1345#if 0
1346these functions are unused
1347
1348/* computes floor(lhs / rhs)
1349 can have lhs, rhs the same
1350*/
1351mpz_t *mpz_div(const mpz_t *lhs, const mpz_t *rhs) {
1352 mpz_t *quo = mpz_zero();
1353 mpz_t rem; mpz_init_zero(&rem);
1354 mpz_divmod_inpl(quo, &rem, lhs, rhs);
1355 mpz_deinit(&rem);
1356 return quo;
1357}
1358
1359/* computes lhs % rhs ( >= 0)
1360 can have lhs, rhs the same
1361*/
1362mpz_t *mpz_mod(const mpz_t *lhs, const mpz_t *rhs) {
1363 mpz_t quo; mpz_init_zero(&quo);
1364 mpz_t *rem = mpz_zero();
1365 mpz_divmod_inpl(&quo, rem, lhs, rhs);
1366 mpz_deinit(&quo);
1367 return rem;
1368}
1369#endif
1370
Damien Georgeffe911d2014-07-24 14:21:37 +01001371// must return actual int value if it fits in mp_int_t
1372mp_int_t mpz_hash(const mpz_t *z) {
1373 mp_int_t val = 0;
1374 mpz_dig_t *d = z->dig + z->len;
1375
Damien George58056b02015-01-09 20:58:58 +00001376 while (d-- > z->dig) {
Damien Georgeffe911d2014-07-24 14:21:37 +01001377 val = (val << DIG_SIZE) | *d;
1378 }
1379
1380 if (z->neg != 0) {
1381 val = -val;
1382 }
1383
1384 return val;
1385}
1386
Damien George40f3c022014-07-03 13:25:24 +01001387bool mpz_as_int_checked(const mpz_t *i, mp_int_t *value) {
Damien George963a5a32015-01-16 17:47:07 +00001388 mp_uint_t val = 0;
Damien George8270e382014-04-03 11:00:54 +00001389 mpz_dig_t *d = i->dig + i->len;
1390
Damien George58056b02015-01-09 20:58:58 +00001391 while (d-- > i->dig) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001392 if (val > (~(WORD_MSBIT_HIGH) >> DIG_SIZE)) {
1393 // will overflow
Damien George8270e382014-04-03 11:00:54 +00001394 return false;
1395 }
Damien Georgec9aa58e2014-07-31 13:41:43 +00001396 val = (val << DIG_SIZE) | *d;
Damien George8270e382014-04-03 11:00:54 +00001397 }
1398
1399 if (i->neg != 0) {
1400 val = -val;
1401 }
1402
1403 *value = val;
1404 return true;
1405}
1406
Damien Georgec9aa58e2014-07-31 13:41:43 +00001407bool mpz_as_uint_checked(const mpz_t *i, mp_uint_t *value) {
1408 if (i->neg != 0) {
1409 // can't represent signed values
1410 return false;
1411 }
1412
1413 mp_uint_t val = 0;
1414 mpz_dig_t *d = i->dig + i->len;
1415
Damien George58056b02015-01-09 20:58:58 +00001416 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001417 if (val > (~(WORD_MSBIT_HIGH) >> (DIG_SIZE - 1))) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001418 // will overflow
1419 return false;
1420 }
1421 val = (val << DIG_SIZE) | *d;
1422 }
1423
1424 *value = val;
1425 return true;
1426}
1427
Damien George271d18e2015-04-25 23:16:39 +01001428// writes at most len bytes to buf (so buf should be zeroed before calling)
1429void mpz_as_bytes(const mpz_t *z, bool big_endian, mp_uint_t len, byte *buf) {
1430 byte *b = buf;
1431 if (big_endian) {
1432 b += len;
1433 }
1434 mpz_dig_t *zdig = z->dig;
1435 int bits = 0;
1436 mpz_dbl_dig_t d = 0;
1437 mpz_dbl_dig_t carry = 1;
1438 for (mp_uint_t zlen = z->len; zlen > 0; --zlen) {
1439 bits += DIG_SIZE;
1440 d = (d << DIG_SIZE) | *zdig++;
1441 for (; bits >= 8; bits -= 8, d >>= 8) {
1442 mpz_dig_t val = d;
1443 if (z->neg) {
Damien George94729072015-04-25 23:51:14 +01001444 val = (~val & 0xff) + carry;
1445 carry = val >> 8;
Damien George271d18e2015-04-25 23:16:39 +01001446 }
1447 if (big_endian) {
1448 *--b = val;
1449 if (b == buf) {
1450 return;
1451 }
1452 } else {
1453 *b++ = val;
1454 if (b == buf + len) {
1455 return;
1456 }
1457 }
1458 }
1459 }
1460}
1461
Damien Georgefb510b32014-06-01 13:32:54 +01001462#if MICROPY_PY_BUILTINS_FLOAT
Damien George52608102014-03-08 15:04:54 +00001463mp_float_t mpz_as_float(const mpz_t *i) {
1464 mp_float_t val = 0;
Damien George438c88d2014-02-22 19:25:23 +00001465 mpz_dig_t *d = i->dig + i->len;
1466
Damien George58056b02015-01-09 20:58:58 +00001467 while (d-- > i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001468 val = val * DIG_BASE + *d;
Damien George438c88d2014-02-22 19:25:23 +00001469 }
1470
1471 if (i->neg != 0) {
1472 val = -val;
1473 }
1474
1475 return val;
1476}
Damien George52608102014-03-08 15:04:54 +00001477#endif
Damien George438c88d2014-02-22 19:25:23 +00001478
Damien Georgeafb1cf72014-09-05 20:37:06 +01001479mp_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 +00001480 if (base < 2 || base > 32) {
1481 return 0;
1482 }
1483
Damien Georgeafb1cf72014-09-05 20:37:06 +01001484 mp_uint_t num_digits = i->len * DIG_SIZE / log_base2_floor[base] + 1;
1485 mp_uint_t num_commas = comma ? num_digits / 3: 0;
1486 mp_uint_t prefix_len = prefix ? strlen(prefix) : 0;
Dave Hylandsc4029e52014-04-07 11:19:51 -07001487
1488 return num_digits + num_commas + prefix_len + 2; // +1 for sign, +1 for null byte
1489}
1490
Damien Georgeafb1cf72014-09-05 20:37:06 +01001491#if 0
1492this function is unused
1493char *mpz_as_str(const mpz_t *i, mp_uint_t base) {
1494 char *s = m_new(char, mpz_as_str_size(i, base, NULL, '\0'));
1495 mpz_as_str_inpl(i, base, NULL, 'a', '\0', s);
Damien George438c88d2014-02-22 19:25:23 +00001496 return s;
1497}
Damien Georgeafb1cf72014-09-05 20:37:06 +01001498#endif
Damien George438c88d2014-02-22 19:25:23 +00001499
1500// assumes enough space as calculated by mpz_as_str_size
1501// returns length of string, not including null byte
Damien Georgeafb1cf72014-09-05 20:37:06 +01001502mp_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 +00001503 if (str == NULL || base < 2 || base > 32) {
1504 str[0] = 0;
1505 return 0;
1506 }
1507
Damien Georgeafb1cf72014-09-05 20:37:06 +01001508 mp_uint_t ilen = i->len;
Damien George438c88d2014-02-22 19:25:23 +00001509
Dave Hylandsc4029e52014-04-07 11:19:51 -07001510 char *s = str;
Damien George438c88d2014-02-22 19:25:23 +00001511 if (ilen == 0) {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001512 if (prefix) {
1513 while (*prefix)
1514 *s++ = *prefix++;
1515 }
1516 *s++ = '0';
1517 *s = '\0';
1518 return s - str;
Damien George438c88d2014-02-22 19:25:23 +00001519 }
1520
Damien Georgeeec91052014-04-08 23:11:00 +01001521 // make a copy of mpz digits, so we can do the div/mod calculation
Damien George438c88d2014-02-22 19:25:23 +00001522 mpz_dig_t *dig = m_new(mpz_dig_t, ilen);
1523 memcpy(dig, i->dig, ilen * sizeof(mpz_dig_t));
1524
1525 // convert
Dave Hylandsc4029e52014-04-07 11:19:51 -07001526 char *last_comma = str;
Damien George438c88d2014-02-22 19:25:23 +00001527 bool done;
1528 do {
1529 mpz_dig_t *d = dig + ilen;
1530 mpz_dbl_dig_t a = 0;
1531
1532 // compute next remainder
1533 while (--d >= dig) {
1534 a = (a << DIG_SIZE) | *d;
1535 *d = a / base;
1536 a %= base;
1537 }
1538
1539 // convert to character
1540 a += '0';
1541 if (a > '9') {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001542 a += base_char - '9' - 1;
Damien George438c88d2014-02-22 19:25:23 +00001543 }
1544 *s++ = a;
1545
1546 // check if number is zero
1547 done = true;
1548 for (d = dig; d < dig + ilen; ++d) {
1549 if (*d != 0) {
1550 done = false;
1551 break;
1552 }
1553 }
Dave Hylandsc4029e52014-04-07 11:19:51 -07001554 if (comma && (s - last_comma) == 3) {
1555 *s++ = comma;
1556 last_comma = s;
1557 }
1558 }
1559 while (!done);
Damien George438c88d2014-02-22 19:25:23 +00001560
Damien Georgeeec91052014-04-08 23:11:00 +01001561 // free the copy of the digits array
1562 m_del(mpz_dig_t, dig, ilen);
1563
Dave Hylandsc4029e52014-04-07 11:19:51 -07001564 if (prefix) {
1565 const char *p = &prefix[strlen(prefix)];
1566 while (p > prefix) {
1567 *s++ = *--p;
1568 }
1569 }
Damien George438c88d2014-02-22 19:25:23 +00001570 if (i->neg != 0) {
1571 *s++ = '-';
1572 }
1573
1574 // reverse string
1575 for (char *u = str, *v = s - 1; u < v; ++u, --v) {
1576 char temp = *u;
1577 *u = *v;
1578 *v = temp;
1579 }
1580
Dave Hylandsc4029e52014-04-07 11:19:51 -07001581 *s = '\0'; // null termination
Damien George438c88d2014-02-22 19:25:23 +00001582
1583 return s - str;
1584}
1585
1586#endif // MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ