blob: b1d6e2b3226e519552682658b00b6a79823d7241 [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
201 assumes enough memory in i; assumes normalised j, k; assumes jlen >= klen
202 can have i, j, k pointing to same memory
203*/
Damien Georgeafb1cf72014-09-05 20:37:06 +0100204STATIC mp_uint_t mpn_and(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 +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 George06201ff2014-03-01 19:50:50 +0000656 mpz_need_dig(z, MPZ_NUM_DIG_FOR_INT);
Damien George438c88d2014-02-22 19:25:23 +0000657
Damien George40f3c022014-07-03 13:25:24 +0100658 mp_uint_t uval;
Damien George438c88d2014-02-22 19:25:23 +0000659 if (val < 0) {
660 z->neg = 1;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000661 uval = -val;
Damien George438c88d2014-02-22 19:25:23 +0000662 } else {
663 z->neg = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000664 uval = val;
Damien George438c88d2014-02-22 19:25:23 +0000665 }
666
667 z->len = 0;
Damien Georgebb4a43f2014-03-12 15:36:06 +0000668 while (uval > 0) {
669 z->dig[z->len++] = uval & DIG_MASK;
670 uval >>= DIG_SIZE;
671 }
672}
673
Damien George95307432014-09-10 22:10:33 +0100674void mpz_set_from_ll(mpz_t *z, long long val, bool is_signed) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000675 mpz_need_dig(z, MPZ_NUM_DIG_FOR_LL);
676
677 unsigned long long uval;
Damien George95307432014-09-10 22:10:33 +0100678 if (is_signed && val < 0) {
Damien Georgebb4a43f2014-03-12 15:36:06 +0000679 z->neg = 1;
680 uval = -val;
681 } else {
682 z->neg = 0;
683 uval = val;
684 }
685
686 z->len = 0;
687 while (uval > 0) {
688 z->dig[z->len++] = uval & DIG_MASK;
689 uval >>= DIG_SIZE;
Damien George438c88d2014-02-22 19:25:23 +0000690 }
691}
692
David Steinberg6e0b6d02015-01-02 12:39:22 +0000693#if MICROPY_PY_BUILTINS_FLOAT
694void mpz_set_from_float(mpz_t *z, mp_float_t src) {
695#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
696#define EXP_SZ 11
697#define FRC_SZ 52
698typedef uint64_t mp_float_int_t;
699#else
700#define EXP_SZ 8
701#define FRC_SZ 23
702typedef uint32_t mp_float_int_t;
703#endif
704 union {
705 mp_float_t f;
706 struct { mp_float_int_t frc:FRC_SZ, exp:EXP_SZ, sgn:1; } p;
707 } u = {src};
708
709 z->neg = u.p.sgn;
710 if (u.p.exp == 0) {
711 // value == 0 || value < 1
712 mpz_init_zero(z);
713 } else if (u.p.exp == ((1 << EXP_SZ) - 1)) {
714 // inf or NaN
715#if 0
716 // TODO: this probably isn't the right place to throw an exception
717 if(u.p.frc == 0)
718 nlr_raise(mp_obj_new_exception_msg_varg(&mp_type_OverflowError, "cannot convert float infinity to integer"));
719 else
720 nlr_raise(mp_obj_new_exception_msg_varg(&mp_type_ValueError, "cannot convert float NaN to integer"));
721#else
722 mpz_init_zero(z);
723#endif
724 } else {
725 const int adj_exp = (int)u.p.exp - ((1 << (EXP_SZ - 1)) - 1);
726 if (adj_exp < 0) {
727 // value < 1 , truncates to 0
728 mpz_init_zero(z);
729 } else if (adj_exp == 0) {
730 // 1 <= value < 2 , so truncates to 1
731 mpz_init_from_int(z, 1);
732 } else {
733 // 2 <= value
734 const int dig_cnt = (adj_exp + 1 + (DIG_SIZE - 1)) / DIG_SIZE;
735 const unsigned int rem = adj_exp % DIG_SIZE;
736 int dig_ind, shft;
737 mp_float_int_t frc = u.p.frc | ((mp_float_int_t)1 << FRC_SZ);
738
739 if (adj_exp < FRC_SZ) {
740 shft = 0;
741 dig_ind = 0;
742 frc >>= FRC_SZ - adj_exp;
743 } else {
744 shft = (rem - FRC_SZ) % DIG_SIZE;
745 dig_ind = (adj_exp - FRC_SZ) / DIG_SIZE;
746 }
747 mpz_need_dig(z, dig_cnt);
748 z->len = dig_cnt;
749 if (dig_ind != 0) {
750 memset(z->dig, 0, dig_ind * sizeof(mpz_dig_t));
751 }
752 if (shft != 0) {
753 z->dig[dig_ind++] = (frc << shft) & DIG_MASK;
754 frc >>= DIG_SIZE - shft;
755 }
756 while (dig_ind != dig_cnt) {
757 z->dig[dig_ind++] = frc & DIG_MASK;
758 frc >>= DIG_SIZE;
759 }
760 }
761 }
762}
763#undef EXP_SZ
764#undef FRC_SZ
765#endif
766
Damien George438c88d2014-02-22 19:25:23 +0000767// returns number of bytes from str that were processed
Damien Georgeafb1cf72014-09-05 20:37:06 +0100768mp_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 +0000769 assert(base < 36);
770
771 const char *cur = str;
772 const char *top = str + len;
773
774 mpz_need_dig(z, len * 8 / DIG_SIZE + 1);
775
776 if (neg) {
777 z->neg = 1;
778 } else {
779 z->neg = 0;
780 }
781
782 z->len = 0;
783 for (; cur < top; ++cur) { // XXX UTF8 next char
Damien Georgeafb1cf72014-09-05 20:37:06 +0100784 //mp_uint_t v = char_to_numeric(cur#); // XXX UTF8 get char
785 mp_uint_t v = *cur;
Damien George438c88d2014-02-22 19:25:23 +0000786 if ('0' <= v && v <= '9') {
787 v -= '0';
788 } else if ('A' <= v && v <= 'Z') {
789 v -= 'A' - 10;
790 } else if ('a' <= v && v <= 'z') {
791 v -= 'a' - 10;
792 } else {
793 break;
794 }
795 if (v >= base) {
796 break;
797 }
798 z->len = mpn_mul_dig_add_dig(z->dig, z->len, base, v);
799 }
800
801 return cur - str;
802}
803
804bool mpz_is_zero(const mpz_t *z) {
805 return z->len == 0;
806}
807
808bool mpz_is_pos(const mpz_t *z) {
809 return z->len > 0 && z->neg == 0;
810}
811
812bool mpz_is_neg(const mpz_t *z) {
813 return z->len > 0 && z->neg != 0;
814}
815
816bool mpz_is_odd(const mpz_t *z) {
817 return z->len > 0 && (z->dig[0] & 1) != 0;
818}
819
820bool mpz_is_even(const mpz_t *z) {
821 return z->len == 0 || (z->dig[0] & 1) == 0;
822}
823
Damien George42f3de92014-10-03 17:44:14 +0000824int mpz_cmp(const mpz_t *z1, const mpz_t *z2) {
825 int cmp = (int)z2->neg - (int)z1->neg;
Damien George438c88d2014-02-22 19:25:23 +0000826 if (cmp != 0) {
827 return cmp;
828 }
829 cmp = mpn_cmp(z1->dig, z1->len, z2->dig, z2->len);
830 if (z1->neg != 0) {
831 cmp = -cmp;
832 }
833 return cmp;
834}
835
Damien George06201ff2014-03-01 19:50:50 +0000836#if 0
837// obsolete
838// compares mpz with an integer that fits within DIG_SIZE bits
Damien Georgeafb1cf72014-09-05 20:37:06 +0100839mp_int_t mpz_cmp_sml_int(const mpz_t *z, mp_int_t sml_int) {
840 mp_int_t cmp;
Damien George438c88d2014-02-22 19:25:23 +0000841 if (z->neg == 0) {
842 if (sml_int < 0) return 1;
843 if (sml_int == 0) {
844 if (z->len == 0) return 0;
845 return 1;
846 }
847 if (z->len == 0) return -1;
848 assert(sml_int < (1 << DIG_SIZE));
849 if (z->len != 1) return 1;
850 cmp = z->dig[0] - sml_int;
851 } else {
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 }
862 if (cmp < 0) return -1;
863 if (cmp > 0) return 1;
864 return 0;
865}
Damien George06201ff2014-03-01 19:50:50 +0000866#endif
Damien George438c88d2014-02-22 19:25:23 +0000867
Damien George438c88d2014-02-22 19:25:23 +0000868#if 0
869these functions are unused
870
871/* returns abs(z)
872*/
873mpz_t *mpz_abs(const mpz_t *z) {
874 mpz_t *z2 = mpz_clone(z);
875 z2->neg = 0;
876 return z2;
877}
878
879/* returns -z
880*/
881mpz_t *mpz_neg(const mpz_t *z) {
882 mpz_t *z2 = mpz_clone(z);
883 z2->neg = 1 - z2->neg;
884 return z2;
885}
886
887/* returns lhs + rhs
888 can have lhs, rhs the same
889*/
890mpz_t *mpz_add(const mpz_t *lhs, const mpz_t *rhs) {
891 mpz_t *z = mpz_zero();
892 mpz_add_inpl(z, lhs, rhs);
893 return z;
894}
895
896/* returns lhs - rhs
897 can have lhs, rhs the same
898*/
899mpz_t *mpz_sub(const mpz_t *lhs, const mpz_t *rhs) {
900 mpz_t *z = mpz_zero();
901 mpz_sub_inpl(z, lhs, rhs);
902 return z;
903}
904
905/* returns lhs * rhs
906 can have lhs, rhs the same
907*/
908mpz_t *mpz_mul(const mpz_t *lhs, const mpz_t *rhs) {
909 mpz_t *z = mpz_zero();
910 mpz_mul_inpl(z, lhs, rhs);
911 return z;
912}
913
914/* returns lhs ** rhs
915 can have lhs, rhs the same
916*/
917mpz_t *mpz_pow(const mpz_t *lhs, const mpz_t *rhs) {
918 mpz_t *z = mpz_zero();
919 mpz_pow_inpl(z, lhs, rhs);
920 return z;
921}
922#endif
923
924/* computes dest = abs(z)
925 can have dest, z the same
926*/
927void mpz_abs_inpl(mpz_t *dest, const mpz_t *z) {
928 if (dest != z) {
929 mpz_set(dest, z);
930 }
931 dest->neg = 0;
932}
933
934/* computes dest = -z
935 can have dest, z the same
936*/
937void mpz_neg_inpl(mpz_t *dest, const mpz_t *z) {
938 if (dest != z) {
939 mpz_set(dest, z);
940 }
941 dest->neg = 1 - dest->neg;
942}
943
Damien George06201ff2014-03-01 19:50:50 +0000944/* computes dest = ~z (= -z - 1)
945 can have dest, z the same
946*/
947void mpz_not_inpl(mpz_t *dest, const mpz_t *z) {
948 if (dest != z) {
949 mpz_set(dest, z);
950 }
Damien Georgee0ac1942014-12-31 19:35:01 +0000951 if (dest->len == 0) {
952 mpz_need_dig(dest, 1);
953 dest->dig[0] = 1;
954 dest->len = 1;
955 dest->neg = 1;
956 } else if (dest->neg) {
Damien George06201ff2014-03-01 19:50:50 +0000957 dest->neg = 0;
958 mpz_dig_t k = 1;
959 dest->len = mpn_sub(dest->dig, dest->dig, dest->len, &k, 1);
960 } else {
Damien Georgee0ac1942014-12-31 19:35:01 +0000961 mpz_need_dig(dest, dest->len + 1);
Damien George06201ff2014-03-01 19:50:50 +0000962 mpz_dig_t k = 1;
963 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &k, 1);
964 dest->neg = 1;
965 }
966}
967
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000968/* computes dest = lhs << rhs
969 can have dest, lhs the same
970*/
Damien George40f3c022014-07-03 13:25:24 +0100971void mpz_shl_inpl(mpz_t *dest, const mpz_t *lhs, mp_int_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +0000972 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000973 mpz_set(dest, lhs);
Damien George06201ff2014-03-01 19:50:50 +0000974 } else if (rhs < 0) {
975 mpz_shr_inpl(dest, lhs, -rhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000976 } else {
Damien George06201ff2014-03-01 19:50:50 +0000977 mpz_need_dig(dest, lhs->len + (rhs + DIG_SIZE - 1) / DIG_SIZE);
978 dest->len = mpn_shl(dest->dig, lhs->dig, lhs->len, rhs);
979 dest->neg = lhs->neg;
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000980 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000981}
982
983/* computes dest = lhs >> rhs
984 can have dest, lhs the same
985*/
Damien George40f3c022014-07-03 13:25:24 +0100986void mpz_shr_inpl(mpz_t *dest, const mpz_t *lhs, mp_int_t rhs) {
Damien George06201ff2014-03-01 19:50:50 +0000987 if (lhs->len == 0 || rhs == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000988 mpz_set(dest, lhs);
Damien George06201ff2014-03-01 19:50:50 +0000989 } else if (rhs < 0) {
990 mpz_shl_inpl(dest, lhs, -rhs);
Damien Georgec5ac2ac2014-02-26 16:56:30 +0000991 } else {
Damien George06201ff2014-03-01 19:50:50 +0000992 mpz_need_dig(dest, lhs->len);
993 dest->len = mpn_shr(dest->dig, lhs->dig, lhs->len, rhs);
994 dest->neg = lhs->neg;
995 if (dest->neg) {
996 // arithmetic shift right, rounding to negative infinity
Damien Georgeafb1cf72014-09-05 20:37:06 +0100997 mp_uint_t n_whole = rhs / DIG_SIZE;
998 mp_uint_t n_part = rhs % DIG_SIZE;
Damien George06201ff2014-03-01 19:50:50 +0000999 mpz_dig_t round_up = 0;
Damien Georgeafb1cf72014-09-05 20:37:06 +01001000 for (mp_uint_t i = 0; i < lhs->len && i < n_whole; i++) {
Damien George06201ff2014-03-01 19:50:50 +00001001 if (lhs->dig[i] != 0) {
1002 round_up = 1;
1003 break;
1004 }
1005 }
1006 if (n_whole < lhs->len && (lhs->dig[n_whole] & ((1 << n_part) - 1)) != 0) {
1007 round_up = 1;
1008 }
1009 if (round_up) {
Damien Georgee0ac1942014-12-31 19:35:01 +00001010 if (dest->len == 0) {
1011 // dest == 0, so need to add 1 by hand (answer will be -1)
1012 dest->dig[0] = 1;
1013 dest->len = 1;
1014 } else {
1015 // dest > 0, so can use mpn_add to add 1
1016 dest->len = mpn_add(dest->dig, dest->dig, dest->len, &round_up, 1);
1017 }
Damien George06201ff2014-03-01 19:50:50 +00001018 }
1019 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001020 }
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001021}
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001022
Damien George438c88d2014-02-22 19:25:23 +00001023/* computes dest = lhs + rhs
1024 can have dest, lhs, rhs the same
1025*/
1026void mpz_add_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1027 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1028 const mpz_t *temp = lhs;
1029 lhs = rhs;
1030 rhs = temp;
1031 }
1032
1033 if (lhs->neg == rhs->neg) {
1034 mpz_need_dig(dest, lhs->len + 1);
1035 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1036 } else {
1037 mpz_need_dig(dest, lhs->len);
1038 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1039 }
1040
1041 dest->neg = lhs->neg;
1042}
1043
1044/* computes dest = lhs - rhs
1045 can have dest, lhs, rhs the same
1046*/
1047void mpz_sub_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1048 bool neg = false;
1049
1050 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1051 const mpz_t *temp = lhs;
1052 lhs = rhs;
1053 rhs = temp;
1054 neg = true;
1055 }
1056
1057 if (lhs->neg != rhs->neg) {
1058 mpz_need_dig(dest, lhs->len + 1);
1059 dest->len = mpn_add(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1060 } else {
1061 mpz_need_dig(dest, lhs->len);
1062 dest->len = mpn_sub(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1063 }
1064
1065 if (neg) {
1066 dest->neg = 1 - lhs->neg;
1067 } else {
1068 dest->neg = lhs->neg;
1069 }
1070}
1071
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001072/* computes dest = lhs & rhs
1073 can have dest, lhs, rhs the same
1074*/
1075void mpz_and_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001076 if (lhs->neg == rhs->neg) {
Damien Georgef55cf102014-05-29 15:01:49 +01001077 if (lhs->neg == 0) {
1078 // make sure lhs has the most digits
1079 if (lhs->len < rhs->len) {
1080 const mpz_t *temp = lhs;
1081 lhs = rhs;
1082 rhs = temp;
1083 }
1084 // do the and'ing
1085 mpz_need_dig(dest, rhs->len);
1086 dest->len = mpn_and(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1087 dest->neg = 0;
1088 } else {
1089 // TODO both args are negative
1090 assert(0);
1091 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001092 } else {
Damien Georgef55cf102014-05-29 15:01:49 +01001093 // args have different sign
1094 // make sure lhs is the positive arg
1095 if (rhs->neg == 0) {
1096 const mpz_t *temp = lhs;
1097 lhs = rhs;
1098 rhs = temp;
1099 }
1100 mpz_need_dig(dest, lhs->len + 1);
1101 dest->len = mpn_and_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1102 assert(dest->len <= dest->alloc);
1103 dest->neg = 0;
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001104 }
Paul Sokolovsky57207b82014-03-23 01:52:36 +02001105}
1106
1107/* computes dest = lhs | rhs
1108 can have dest, lhs, rhs the same
1109*/
1110void mpz_or_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1111 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1112 const mpz_t *temp = lhs;
1113 lhs = rhs;
1114 rhs = temp;
1115 }
1116
1117 if (lhs->neg == rhs->neg) {
1118 mpz_need_dig(dest, lhs->len);
1119 dest->len = mpn_or(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1120 } else {
1121 mpz_need_dig(dest, lhs->len);
1122 // TODO
1123 assert(0);
1124// dest->len = mpn_or_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1125 }
1126
1127 dest->neg = lhs->neg;
1128}
1129
1130/* computes dest = lhs ^ rhs
1131 can have dest, lhs, rhs the same
1132*/
1133void mpz_xor_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1134 if (mpn_cmp(lhs->dig, lhs->len, rhs->dig, rhs->len) < 0) {
1135 const mpz_t *temp = lhs;
1136 lhs = rhs;
1137 rhs = temp;
1138 }
1139
1140 if (lhs->neg == rhs->neg) {
1141 mpz_need_dig(dest, lhs->len);
1142 dest->len = mpn_xor(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1143 } else {
1144 mpz_need_dig(dest, lhs->len);
1145 // TODO
1146 assert(0);
1147// dest->len = mpn_xor_neg(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1148 }
1149
1150 dest->neg = 0;
1151}
1152
Damien George438c88d2014-02-22 19:25:23 +00001153/* computes dest = lhs * rhs
1154 can have dest, lhs, rhs the same
1155*/
1156void mpz_mul_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs)
1157{
1158 if (lhs->len == 0 || rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001159 mpz_set_from_int(dest, 0);
1160 return;
Damien George438c88d2014-02-22 19:25:23 +00001161 }
1162
1163 mpz_t *temp = NULL;
1164 if (lhs == dest) {
1165 lhs = temp = mpz_clone(lhs);
1166 if (rhs == dest) {
1167 rhs = lhs;
1168 }
1169 } else if (rhs == dest) {
1170 rhs = temp = mpz_clone(rhs);
1171 }
1172
1173 mpz_need_dig(dest, lhs->len + rhs->len); // min mem l+r-1, max mem l+r
1174 memset(dest->dig, 0, dest->alloc * sizeof(mpz_dig_t));
1175 dest->len = mpn_mul(dest->dig, lhs->dig, lhs->len, rhs->dig, rhs->len);
1176
1177 if (lhs->neg == rhs->neg) {
1178 dest->neg = 0;
1179 } else {
1180 dest->neg = 1;
1181 }
1182
1183 mpz_free(temp);
1184}
1185
1186/* computes dest = lhs ** rhs
1187 can have dest, lhs, rhs the same
1188*/
1189void mpz_pow_inpl(mpz_t *dest, const mpz_t *lhs, const mpz_t *rhs) {
1190 if (lhs->len == 0 || rhs->neg != 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001191 mpz_set_from_int(dest, 0);
1192 return;
Damien George438c88d2014-02-22 19:25:23 +00001193 }
1194
1195 if (rhs->len == 0) {
Damien Georgec5ac2ac2014-02-26 16:56:30 +00001196 mpz_set_from_int(dest, 1);
1197 return;
Damien George438c88d2014-02-22 19:25:23 +00001198 }
1199
1200 mpz_t *x = mpz_clone(lhs);
1201 mpz_t *n = mpz_clone(rhs);
1202
1203 mpz_set_from_int(dest, 1);
1204
1205 while (n->len > 0) {
1206 if (mpz_is_odd(n)) {
1207 mpz_mul_inpl(dest, dest, x);
1208 }
Damien George438c88d2014-02-22 19:25:23 +00001209 n->len = mpn_shr(n->dig, n->dig, n->len, 1);
Damien George5bf565e2014-04-04 00:16:32 +01001210 if (n->len == 0) {
1211 break;
1212 }
1213 mpz_mul_inpl(x, x, x);
Damien George438c88d2014-02-22 19:25:23 +00001214 }
1215
1216 mpz_free(x);
1217 mpz_free(n);
1218}
1219
1220/* computes gcd(z1, z2)
1221 based on Knuth's modified gcd algorithm (I think?)
1222 gcd(z1, z2) >= 0
1223 gcd(0, 0) = 0
1224 gcd(z, 0) = abs(z)
1225*/
1226mpz_t *mpz_gcd(const mpz_t *z1, const mpz_t *z2) {
1227 if (z1->len == 0) {
1228 mpz_t *a = mpz_clone(z2);
1229 a->neg = 0;
1230 return a;
1231 } else if (z2->len == 0) {
1232 mpz_t *a = mpz_clone(z1);
1233 a->neg = 0;
1234 return a;
1235 }
1236
1237 mpz_t *a = mpz_clone(z1);
1238 mpz_t *b = mpz_clone(z2);
1239 mpz_t c; mpz_init_zero(&c);
1240 a->neg = 0;
1241 b->neg = 0;
1242
1243 for (;;) {
1244 if (mpz_cmp(a, b) < 0) {
1245 if (a->len == 0) {
1246 mpz_free(a);
1247 mpz_deinit(&c);
1248 return b;
1249 }
1250 mpz_t *t = a; a = b; b = t;
1251 }
1252 if (!(b->len >= 2 || (b->len == 1 && b->dig[0] > 1))) { // compute b > 0; could be mpz_cmp_small_int(b, 1) > 0
1253 break;
1254 }
1255 mpz_set(&c, b);
1256 do {
1257 mpz_add_inpl(&c, &c, &c);
1258 } while (mpz_cmp(&c, a) <= 0);
1259 c.len = mpn_shr(c.dig, c.dig, c.len, 1);
1260 mpz_sub_inpl(a, a, &c);
1261 }
1262
1263 mpz_deinit(&c);
1264
1265 if (b->len == 1 && b->dig[0] == 1) { // compute b == 1; could be mpz_cmp_small_int(b, 1) == 0
1266 mpz_free(a);
1267 return b;
1268 } else {
1269 mpz_free(b);
1270 return a;
1271 }
1272}
1273
1274/* computes lcm(z1, z2)
1275 = abs(z1) / gcd(z1, z2) * abs(z2)
1276 lcm(z1, z1) >= 0
1277 lcm(0, 0) = 0
1278 lcm(z, 0) = 0
1279*/
Damien George5d9b8162014-08-07 14:27:48 +00001280mpz_t *mpz_lcm(const mpz_t *z1, const mpz_t *z2) {
1281 if (z1->len == 0 || z2->len == 0) {
1282 return mpz_zero();
stijn01d6be42014-05-05 12:18:27 +02001283 }
Damien George438c88d2014-02-22 19:25:23 +00001284
1285 mpz_t *gcd = mpz_gcd(z1, z2);
1286 mpz_t *quo = mpz_zero();
1287 mpz_t *rem = mpz_zero();
1288 mpz_divmod_inpl(quo, rem, z1, gcd);
1289 mpz_mul_inpl(rem, quo, z2);
1290 mpz_free(gcd);
1291 mpz_free(quo);
1292 rem->neg = 0;
1293 return rem;
1294}
1295
1296/* computes new integers in quo and rem such that:
1297 quo * rhs + rem = lhs
1298 0 <= rem < rhs
1299 can have lhs, rhs the same
1300*/
1301void mpz_divmod(const mpz_t *lhs, const mpz_t *rhs, mpz_t **quo, mpz_t **rem) {
1302 *quo = mpz_zero();
1303 *rem = mpz_zero();
1304 mpz_divmod_inpl(*quo, *rem, lhs, rhs);
1305}
1306
1307/* computes new integers in quo and rem such that:
1308 quo * rhs + rem = lhs
1309 0 <= rem < rhs
1310 can have lhs, rhs the same
1311*/
1312void mpz_divmod_inpl(mpz_t *dest_quo, mpz_t *dest_rem, const mpz_t *lhs, const mpz_t *rhs) {
1313 if (rhs->len == 0) {
1314 mpz_set_from_int(dest_quo, 0);
1315 mpz_set_from_int(dest_rem, 0);
1316 return;
1317 }
1318
1319 mpz_need_dig(dest_quo, lhs->len + 1); // +1 necessary?
1320 memset(dest_quo->dig, 0, (lhs->len + 1) * sizeof(mpz_dig_t));
1321 dest_quo->len = 0;
1322 mpz_need_dig(dest_rem, lhs->len + 1); // +1 necessary?
1323 mpz_set(dest_rem, lhs);
1324 //rhs->dig[rhs->len] = 0;
1325 mpn_div(dest_rem->dig, &dest_rem->len, rhs->dig, rhs->len, dest_quo->dig, &dest_quo->len);
1326
1327 if (lhs->neg != rhs->neg) {
1328 dest_quo->neg = 1;
1329 }
1330}
1331
1332#if 0
1333these functions are unused
1334
1335/* computes floor(lhs / rhs)
1336 can have lhs, rhs the same
1337*/
1338mpz_t *mpz_div(const mpz_t *lhs, const mpz_t *rhs) {
1339 mpz_t *quo = mpz_zero();
1340 mpz_t rem; mpz_init_zero(&rem);
1341 mpz_divmod_inpl(quo, &rem, lhs, rhs);
1342 mpz_deinit(&rem);
1343 return quo;
1344}
1345
1346/* computes lhs % rhs ( >= 0)
1347 can have lhs, rhs the same
1348*/
1349mpz_t *mpz_mod(const mpz_t *lhs, const mpz_t *rhs) {
1350 mpz_t quo; mpz_init_zero(&quo);
1351 mpz_t *rem = mpz_zero();
1352 mpz_divmod_inpl(&quo, rem, lhs, rhs);
1353 mpz_deinit(&quo);
1354 return rem;
1355}
1356#endif
1357
Damien Georgeffe911d2014-07-24 14:21:37 +01001358// must return actual int value if it fits in mp_int_t
1359mp_int_t mpz_hash(const mpz_t *z) {
1360 mp_int_t val = 0;
1361 mpz_dig_t *d = z->dig + z->len;
1362
1363 while (--d >= z->dig) {
1364 val = (val << DIG_SIZE) | *d;
1365 }
1366
1367 if (z->neg != 0) {
1368 val = -val;
1369 }
1370
1371 return val;
1372}
1373
Damien George40f3c022014-07-03 13:25:24 +01001374bool mpz_as_int_checked(const mpz_t *i, mp_int_t *value) {
1375 mp_int_t val = 0;
Damien George8270e382014-04-03 11:00:54 +00001376 mpz_dig_t *d = i->dig + i->len;
1377
1378 while (--d >= i->dig) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001379 if (val > (~(WORD_MSBIT_HIGH) >> DIG_SIZE)) {
1380 // will overflow
Damien George8270e382014-04-03 11:00:54 +00001381 return false;
1382 }
Damien Georgec9aa58e2014-07-31 13:41:43 +00001383 val = (val << DIG_SIZE) | *d;
Damien George8270e382014-04-03 11:00:54 +00001384 }
1385
1386 if (i->neg != 0) {
1387 val = -val;
1388 }
1389
1390 *value = val;
1391 return true;
1392}
1393
Damien Georgec9aa58e2014-07-31 13:41:43 +00001394bool mpz_as_uint_checked(const mpz_t *i, mp_uint_t *value) {
1395 if (i->neg != 0) {
1396 // can't represent signed values
1397 return false;
1398 }
1399
1400 mp_uint_t val = 0;
1401 mpz_dig_t *d = i->dig + i->len;
1402
1403 while (--d >= i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001404 if (val > (~(WORD_MSBIT_HIGH) >> (DIG_SIZE - 1))) {
Damien Georgec9aa58e2014-07-31 13:41:43 +00001405 // will overflow
1406 return false;
1407 }
1408 val = (val << DIG_SIZE) | *d;
1409 }
1410
1411 *value = val;
1412 return true;
1413}
1414
Damien Georgefb510b32014-06-01 13:32:54 +01001415#if MICROPY_PY_BUILTINS_FLOAT
Damien George52608102014-03-08 15:04:54 +00001416mp_float_t mpz_as_float(const mpz_t *i) {
1417 mp_float_t val = 0;
Damien George438c88d2014-02-22 19:25:23 +00001418 mpz_dig_t *d = i->dig + i->len;
1419
1420 while (--d >= i->dig) {
Damien George9a21d2e2014-09-06 17:15:34 +01001421 val = val * DIG_BASE + *d;
Damien George438c88d2014-02-22 19:25:23 +00001422 }
1423
1424 if (i->neg != 0) {
1425 val = -val;
1426 }
1427
1428 return val;
1429}
Damien George52608102014-03-08 15:04:54 +00001430#endif
Damien George438c88d2014-02-22 19:25:23 +00001431
Damien Georgeafb1cf72014-09-05 20:37:06 +01001432mp_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 +00001433 if (base < 2 || base > 32) {
1434 return 0;
1435 }
1436
Damien Georgeafb1cf72014-09-05 20:37:06 +01001437 mp_uint_t num_digits = i->len * DIG_SIZE / log_base2_floor[base] + 1;
1438 mp_uint_t num_commas = comma ? num_digits / 3: 0;
1439 mp_uint_t prefix_len = prefix ? strlen(prefix) : 0;
Dave Hylandsc4029e52014-04-07 11:19:51 -07001440
1441 return num_digits + num_commas + prefix_len + 2; // +1 for sign, +1 for null byte
1442}
1443
Damien Georgeafb1cf72014-09-05 20:37:06 +01001444#if 0
1445this function is unused
1446char *mpz_as_str(const mpz_t *i, mp_uint_t base) {
1447 char *s = m_new(char, mpz_as_str_size(i, base, NULL, '\0'));
1448 mpz_as_str_inpl(i, base, NULL, 'a', '\0', s);
Damien George438c88d2014-02-22 19:25:23 +00001449 return s;
1450}
Damien Georgeafb1cf72014-09-05 20:37:06 +01001451#endif
Damien George438c88d2014-02-22 19:25:23 +00001452
1453// assumes enough space as calculated by mpz_as_str_size
1454// returns length of string, not including null byte
Damien Georgeafb1cf72014-09-05 20:37:06 +01001455mp_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 +00001456 if (str == NULL || base < 2 || base > 32) {
1457 str[0] = 0;
1458 return 0;
1459 }
1460
Damien Georgeafb1cf72014-09-05 20:37:06 +01001461 mp_uint_t ilen = i->len;
Damien George438c88d2014-02-22 19:25:23 +00001462
Dave Hylandsc4029e52014-04-07 11:19:51 -07001463 char *s = str;
Damien George438c88d2014-02-22 19:25:23 +00001464 if (ilen == 0) {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001465 if (prefix) {
1466 while (*prefix)
1467 *s++ = *prefix++;
1468 }
1469 *s++ = '0';
1470 *s = '\0';
1471 return s - str;
Damien George438c88d2014-02-22 19:25:23 +00001472 }
1473
Damien Georgeeec91052014-04-08 23:11:00 +01001474 // make a copy of mpz digits, so we can do the div/mod calculation
Damien George438c88d2014-02-22 19:25:23 +00001475 mpz_dig_t *dig = m_new(mpz_dig_t, ilen);
1476 memcpy(dig, i->dig, ilen * sizeof(mpz_dig_t));
1477
1478 // convert
Dave Hylandsc4029e52014-04-07 11:19:51 -07001479 char *last_comma = str;
Damien George438c88d2014-02-22 19:25:23 +00001480 bool done;
1481 do {
1482 mpz_dig_t *d = dig + ilen;
1483 mpz_dbl_dig_t a = 0;
1484
1485 // compute next remainder
1486 while (--d >= dig) {
1487 a = (a << DIG_SIZE) | *d;
1488 *d = a / base;
1489 a %= base;
1490 }
1491
1492 // convert to character
1493 a += '0';
1494 if (a > '9') {
Dave Hylandsc4029e52014-04-07 11:19:51 -07001495 a += base_char - '9' - 1;
Damien George438c88d2014-02-22 19:25:23 +00001496 }
1497 *s++ = a;
1498
1499 // check if number is zero
1500 done = true;
1501 for (d = dig; d < dig + ilen; ++d) {
1502 if (*d != 0) {
1503 done = false;
1504 break;
1505 }
1506 }
Dave Hylandsc4029e52014-04-07 11:19:51 -07001507 if (comma && (s - last_comma) == 3) {
1508 *s++ = comma;
1509 last_comma = s;
1510 }
1511 }
1512 while (!done);
Damien George438c88d2014-02-22 19:25:23 +00001513
Damien Georgeeec91052014-04-08 23:11:00 +01001514 // free the copy of the digits array
1515 m_del(mpz_dig_t, dig, ilen);
1516
Dave Hylandsc4029e52014-04-07 11:19:51 -07001517 if (prefix) {
1518 const char *p = &prefix[strlen(prefix)];
1519 while (p > prefix) {
1520 *s++ = *--p;
1521 }
1522 }
Damien George438c88d2014-02-22 19:25:23 +00001523 if (i->neg != 0) {
1524 *s++ = '-';
1525 }
1526
1527 // reverse string
1528 for (char *u = str, *v = s - 1; u < v; ++u, --v) {
1529 char temp = *u;
1530 *u = *v;
1531 *v = temp;
1532 }
1533
Dave Hylandsc4029e52014-04-07 11:19:51 -07001534 *s = '\0'; // null termination
Damien George438c88d2014-02-22 19:25:23 +00001535
1536 return s - str;
1537}
1538
1539#endif // MICROPY_LONGINT_IMPL == MICROPY_LONGINT_IMPL_MPZ