Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mini-gmp/mini-gmp.c
Line
Count
Source (jump to first uncovered line)
1
/* mini-gmp, a minimalistic implementation of a GNU GMP subset.
2
3
   Contributed to the GNU project by Niels Möller
4
   Additional functionalities and improvements by Marco Bodrato.
5
6
Copyright 1991-1997, 1999-2022 Free Software Foundation, Inc.
7
8
This file is part of the GNU MP Library.
9
10
The GNU MP Library is free software; you can redistribute it and/or modify
11
it under the terms of either:
12
13
  * the GNU Lesser General Public License as published by the Free
14
    Software Foundation; either version 3 of the License, or (at your
15
    option) any later version.
16
17
or
18
19
  * the GNU General Public License as published by the Free Software
20
    Foundation; either version 2 of the License, or (at your option) any
21
    later version.
22
23
or both in parallel, as here.
24
25
The GNU MP Library is distributed in the hope that it will be useful, but
26
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
27
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
28
for more details.
29
30
You should have received copies of the GNU General Public License and the
31
GNU Lesser General Public License along with the GNU MP Library.  If not,
32
see https://www.gnu.org/licenses/.  */
33
34
/* NOTE: All functions in this file which are not declared in
35
   mini-gmp.h are internal, and are not intended to be compatible
36
   with GMP or with future versions of mini-gmp. */
37
38
/* Much of the material copied from GMP files, including: gmp-impl.h,
39
   longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c,
40
   mpn/generic/lshift.c, mpn/generic/mul_1.c,
41
   mpn/generic/mul_basecase.c, mpn/generic/rshift.c,
42
   mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c,
43
   mpn/generic/submul_1.c. */
44
45
#include <assert.h>
46
#include <ctype.h>
47
#include <limits.h>
48
#include <stdio.h>
49
#include <stdlib.h>
50
#include <string.h>
51
52
#include "mini-gmp.h"
53
54
#if !defined(MINI_GMP_DONT_USE_FLOAT_H)
55
#include <float.h>
56
#endif
57
58

59
/* Macros */
60
331M
#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
61
62
148k
#define GMP_LIMB_MAX ((mp_limb_t) ~ (mp_limb_t) 0)
63
685k
#define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
64
65
89.7M
#define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
66
88.8M
#define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)
67
68
29.5M
#define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
69
142
#define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))
70
71
1.79M
#define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
72
1.56k
#define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
73
74
133
#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
75
33.0k
#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
76
77
538
#define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b)))
78
79
#if defined(DBL_MANT_DIG) && FLT_RADIX == 2
80
0
#define GMP_DBL_MANT_BITS DBL_MANT_DIG
81
#else
82
#define GMP_DBL_MANT_BITS (53)
83
#endif
84
85
/* Return non-zero if xp,xsize and yp,ysize overlap.
86
   If xp+xsize<=yp there's no overlap, or if yp+ysize<=xp there's no
87
   overlap.  If both these are false, there's an overlap. */
88
#define GMP_MPN_OVERLAP_P(xp, xsize, yp, ysize)       \
89
  ((xp) + (xsize) > (yp) && (yp) + (ysize) > (xp))
90
91
135k
#define gmp_assert_nocarry(x) do { \
92
135k
    mp_limb_t __cy = (x);    \
93
135k
    assert (__cy == 0);      \
94
135k
    (void) (__cy);       \
95
135k
  } while (0)
96
97
115k
#define gmp_clz(count, x) do {           \
98
115k
    mp_limb_t __clz_x = (x);            \
99
115k
    unsigned __clz_c = 0;           \
100
115k
    int LOCAL_SHIFT_BITS = 8;           \
101
115k
    if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS)       \
102
115k
      for (;                \
103
894k
     (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \
104
778k
     __clz_c += 8)           \
105
778k
  { __clz_x <<= LOCAL_SHIFT_BITS; }        \
106
685k
    for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++)   \
107
570k
      __clz_x <<= 1;             \
108
115k
    (count) = __clz_c;              \
109
115k
  } while (0)
110
111
100k
#define gmp_ctz(count, x) do {           \
112
100k
    mp_limb_t __ctz_x = (x);            \
113
100k
    unsigned __ctz_c = 0;           \
114
100k
    gmp_clz (__ctz_c, __ctz_x & - __ctz_x);        \
115
100k
    (count) = GMP_LIMB_BITS - 1 - __ctz_c;       \
116
100k
  } while (0)
117
118
#define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
119
12.0M
  do {                 \
120
12.0M
    mp_limb_t __x;              \
121
12.0M
    __x = (al) + (bl);              \
122
12.0M
    (sh) = (ah) + (bh) + (__x < (al));          \
123
12.0M
    (sl) = __x;               \
124
12.0M
  } while (0)
125
126
#define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
127
609k
  do {                 \
128
609k
    mp_limb_t __x;              \
129
609k
    __x = (al) - (bl);              \
130
609k
    (sh) = (ah) - (bh) - ((al) < (bl));         \
131
609k
    (sl) = __x;               \
132
609k
  } while (0)
133
134
#define gmp_umul_ppmm(w1, w0, u, v)         \
135
29.5M
  do {                 \
136
29.5M
    int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS;       \
137
29.5M
    if (sizeof(unsigned int) * CHAR_BIT >= 2 * GMP_LIMB_BITS)   \
138
29.5M
      {                 \
139
0
  unsigned int __ww = (unsigned int) (u) * (v);     \
140
0
  w0 = (mp_limb_t) __ww;            \
141
0
  w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS);     \
142
0
      }                  \
143
29.5M
    else if (GMP_ULONG_BITS >= 2 * GMP_LIMB_BITS)     \
144
29.5M
      {                 \
145
0
  unsigned long int __ww = (unsigned long int) (u) * (v);   \
146
0
  w0 = (mp_limb_t) __ww;            \
147
0
  w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS);     \
148
0
      }                  \
149
29.5M
    else {               \
150
29.5M
      mp_limb_t __x0, __x1, __x2, __x3;         \
151
29.5M
      unsigned __ul, __vl, __uh, __vh;          \
152
29.5M
      mp_limb_t __u = (u), __v = (v);         \
153
29.5M
      assert (sizeof (unsigned) * 2 >= sizeof (mp_limb_t));   \
154
29.5M
                  \
155
29.5M
      __ul = __u & GMP_LLIMB_MASK;         \
156
29.5M
      __uh = __u >> (GMP_LIMB_BITS / 2);        \
157
29.5M
      __vl = __v & GMP_LLIMB_MASK;         \
158
29.5M
      __vh = __v >> (GMP_LIMB_BITS / 2);        \
159
29.5M
                  \
160
29.5M
      __x0 = (mp_limb_t) __ul * __vl;         \
161
29.5M
      __x1 = (mp_limb_t) __ul * __vh;         \
162
29.5M
      __x2 = (mp_limb_t) __uh * __vl;         \
163
29.5M
      __x3 = (mp_limb_t) __uh * __vh;         \
164
29.5M
                  \
165
29.5M
      __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \
166
29.5M
      __x1 += __x2;   /* but this indeed can */   \
167
29.5M
      if (__x1 < __x2)   /* did we get it? */      \
168
29.5M
  __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */  \
169
29.5M
                  \
170
29.5M
      (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2));     \
171
29.5M
      (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK);  \
172
29.5M
    }                 \
173
29.5M
  } while (0)
174
175
/* If mp_limb_t is of size smaller than int, plain u*v implies
176
   automatic promotion to *signed* int, and then multiply may overflow
177
   and cause undefined behavior. Explicitly cast to unsigned int for
178
   that case. */
179
#define gmp_umullo_limb(u, v) \
180
11.7M
  ((sizeof(mp_limb_t) >= sizeof(int)) ? (u)*(v) : (unsigned int)(u) * (v))
181
182
#define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di)      \
183
11.4M
  do {                 \
184
11.4M
    mp_limb_t _qh, _ql, _r, _mask;          \
185
11.4M
    gmp_umul_ppmm (_qh, _ql, (nh), (di));        \
186
11.4M
    gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl));    \
187
11.4M
    _r = (nl) - gmp_umullo_limb (_qh, (d));        \
188
11.4M
    _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */   \
189
11.4M
    _qh += _mask;             \
190
11.4M
    _r += _mask & (d);              \
191
11.4M
    if (_r >= (d))             \
192
11.4M
      {                 \
193
25
  _r -= (d);              \
194
25
  _qh++;                \
195
25
      }                  \
196
11.4M
                  \
197
11.4M
    (r) = _r;               \
198
11.4M
    (q) = _qh;                \
199
11.4M
  } while (0)
200
201
#define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv)   \
202
303k
  do {                 \
203
303k
    mp_limb_t _q0, _t1, _t0, _mask;         \
204
303k
    gmp_umul_ppmm ((q), _q0, (n2), (dinv));        \
205
303k
    gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1));      \
206
303k
                  \
207
303k
    /* Compute the two most significant limbs of n - q'd */   \
208
303k
    (r1) = (n1) - gmp_umullo_limb ((d1), (q));       \
209
303k
    gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0));    \
210
303k
    gmp_umul_ppmm (_t1, _t0, (d0), (q));       \
211
303k
    gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0);      \
212
303k
    (q)++;                \
213
303k
                  \
214
303k
    /* Conditionally adjust q and the remainders */     \
215
303k
    _mask = - (mp_limb_t) ((r1) >= _q0);        \
216
303k
    (q) += _mask;             \
217
303k
    gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
218
303k
    if ((r1) >= (d1))             \
219
303k
      {                 \
220
1.99k
  if ((r1) > (d1) || (r0) >= (d0))       \
221
1.99k
    {               \
222
1.98k
      (q)++;              \
223
1.98k
      gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0));  \
224
1.98k
    }                \
225
1.99k
      }                  \
226
303k
  } while (0)
227
228
/* Swap macros. */
229
#define MP_LIMB_T_SWAP(x, y)            \
230
33
  do {                 \
231
33
    mp_limb_t __mp_limb_t_swap__tmp = (x);        \
232
33
    (x) = (y);                \
233
33
    (y) = __mp_limb_t_swap__tmp;          \
234
33
  } while (0)
235
#define MP_SIZE_T_SWAP(x, y)            \
236
132k
  do {                 \
237
132k
    mp_size_t __mp_size_t_swap__tmp = (x);        \
238
132k
    (x) = (y);                \
239
132k
    (y) = __mp_size_t_swap__tmp;          \
240
132k
  } while (0)
241
#define MP_BITCNT_T_SWAP(x,y)     \
242
45
  do {           \
243
45
    mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x);  \
244
45
    (x) = (y);          \
245
45
    (y) = __mp_bitcnt_t_swap__tmp;    \
246
45
  } while (0)
247
#define MP_PTR_SWAP(x, y)           \
248
69.1k
  do {                 \
249
69.1k
    mp_ptr __mp_ptr_swap__tmp = (x);          \
250
69.1k
    (x) = (y);                \
251
69.1k
    (y) = __mp_ptr_swap__tmp;           \
252
69.1k
  } while (0)
253
#define MP_SRCPTR_SWAP(x, y)            \
254
6
  do {                 \
255
6
    mp_srcptr __mp_srcptr_swap__tmp = (x);        \
256
6
    (x) = (y);                \
257
6
    (y) = __mp_srcptr_swap__tmp;          \
258
6
  } while (0)
259
260
#define MPN_PTR_SWAP(xp,xs, yp,ys)          \
261
65.0k
  do {                 \
262
65.0k
    MP_PTR_SWAP (xp, yp);            \
263
65.0k
    MP_SIZE_T_SWAP (xs, ys);            \
264
65.0k
  } while(0)
265
#define MPN_SRCPTR_SWAP(xp,xs, yp,ys)         \
266
6
  do {                 \
267
6
    MP_SRCPTR_SWAP (xp, yp);            \
268
6
    MP_SIZE_T_SWAP (xs, ys);            \
269
6
  } while(0)
270
271
#define MPZ_PTR_SWAP(x, y)            \
272
45
  do {                 \
273
45
    mpz_ptr __mpz_ptr_swap__tmp = (x);          \
274
45
    (x) = (y);                \
275
45
    (y) = __mpz_ptr_swap__tmp;            \
276
45
  } while (0)
277
#define MPZ_SRCPTR_SWAP(x, y)           \
278
2.67k
  do {                 \
279
2.67k
    mpz_srcptr __mpz_srcptr_swap__tmp = (x);        \
280
2.67k
    (x) = (y);                \
281
2.67k
    (y) = __mpz_srcptr_swap__tmp;         \
282
2.67k
  } while (0)
283
284
const int mp_bits_per_limb = GMP_LIMB_BITS;
285
286

287
/* Memory allocation and other helper functions. */
288
static void
289
gmp_die (const char *msg)
290
0
{
291
0
  fprintf (stderr, "%s\n", msg);
292
0
  abort();
293
0
}
294
295
static void *
296
gmp_default_alloc (size_t size)
297
92.0k
{
298
92.0k
  void *p;
299
300
92.0k
  assert (size > 0);
301
302
92.0k
  p = malloc (size);
303
92.0k
  if (!p)
304
0
    gmp_die("gmp_default_alloc: Virtual memory exhausted.");
305
306
92.0k
  return p;
307
92.0k
}
308
309
static void *
310
gmp_default_realloc (void *old, size_t unused_old_size, size_t new_size)
311
9.50k
{
312
9.50k
  void * p;
313
314
9.50k
  p = realloc (old, new_size);
315
316
9.50k
  if (!p)
317
0
    gmp_die("gmp_default_realloc: Virtual memory exhausted.");
318
319
9.50k
  return p;
320
9.50k
}
321
322
static void
323
gmp_default_free (void *p, size_t unused_size)
324
91.4k
{
325
91.4k
  free (p);
326
91.4k
}
327
328
static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc;
329
static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc;
330
static void (*gmp_free_func) (void *, size_t) = gmp_default_free;
331
332
void
333
mp_get_memory_functions (void *(**alloc_func) (size_t),
334
       void *(**realloc_func) (void *, size_t, size_t),
335
       void (**free_func) (void *, size_t))
336
0
{
337
0
  if (alloc_func)
338
0
    *alloc_func = gmp_allocate_func;
339
340
0
  if (realloc_func)
341
0
    *realloc_func = gmp_reallocate_func;
342
343
0
  if (free_func)
344
0
    *free_func = gmp_free_func;
345
0
}
346
347
void
348
mp_set_memory_functions (void *(*alloc_func) (size_t),
349
       void *(*realloc_func) (void *, size_t, size_t),
350
       void (*free_func) (void *, size_t))
351
0
{
352
0
  if (!alloc_func)
353
0
    alloc_func = gmp_default_alloc;
354
0
  if (!realloc_func)
355
0
    realloc_func = gmp_default_realloc;
356
0
  if (!free_func)
357
0
    free_func = gmp_default_free;
358
359
0
  gmp_allocate_func = alloc_func;
360
0
  gmp_reallocate_func = realloc_func;
361
0
  gmp_free_func = free_func;
362
0
}
363
364
92.0k
#define gmp_alloc(size) ((*gmp_allocate_func)((size)))
365
91.4k
#define gmp_free(p, size) ((*gmp_free_func) ((p), (size)))
366
9.50k
#define gmp_realloc(ptr, old_size, size) ((*gmp_reallocate_func)(ptr, old_size, size))
367
368
static mp_ptr
369
gmp_alloc_limbs (mp_size_t size)
370
86.6k
{
371
86.6k
  return (mp_ptr) gmp_alloc (size * sizeof (mp_limb_t));
372
86.6k
}
373
374
static mp_ptr
375
gmp_realloc_limbs (mp_ptr old, mp_size_t old_size, mp_size_t size)
376
8.98k
{
377
8.98k
  assert (size > 0);
378
8.98k
  return (mp_ptr) gmp_realloc (old, old_size * sizeof (mp_limb_t), size * sizeof (mp_limb_t));
379
8.98k
}
380
381
static void
382
gmp_free_limbs (mp_ptr old, mp_size_t size)
383
86.6k
{
384
86.6k
  gmp_free (old, size * sizeof (mp_limb_t));
385
86.6k
}
386
387

388
/* MPN interface */
389
390
void
391
mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
392
13.2k
{
393
13.2k
  mp_size_t i;
394
210k
  for (i = 0; i < n; i++)
395
197k
    d[i] = s[i];
396
13.2k
}
397
398
void
399
mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
400
156
{
401
1.06k
  while (--n >= 0)
402
908
    d[n] = s[n];
403
156
}
404
405
int
406
mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
407
188k
{
408
195k
  while (--n >= 0)
409
195k
    {
410
195k
      if (ap[n] != bp[n])
411
187k
  return ap[n] > bp[n] ? 1 : -1;
412
195k
    }
413
382
  return 0;
414
188k
}
415
416
static int
417
mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
418
100k
{
419
100k
  if (an != bn)
420
8.98k
    return an < bn ? -1 : 1;
421
91.2k
  else
422
91.2k
    return mpn_cmp (ap, bp, an);
423
100k
}
424
425
static mp_size_t
426
mpn_normalized_size (mp_srcptr xp, mp_size_t n)
427
132k
{
428
138k
  while (n > 0 && xp[n-1] == 0)
429
6.83k
    --n;
430
132k
  return n;
431
132k
}
432
433
int
434
mpn_zero_p(mp_srcptr rp, mp_size_t n)
435
0
{
436
0
  return mpn_normalized_size (rp, n) == 0;
437
0
}
438
439
void
440
mpn_zero (mp_ptr rp, mp_size_t n)
441
185k
{
442
185k
  while (--n >= 0)
443
38
    rp[n] = 0;
444
185k
}
445
446
mp_limb_t
447
mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
448
109k
{
449
109k
  mp_size_t i;
450
451
109k
  assert (n > 0);
452
109k
  i = 0;
453
109k
  do
454
5.09M
    {
455
5.09M
      mp_limb_t r = ap[i] + b;
456
      /* Carry out */
457
5.09M
      b = (r < b);
458
5.09M
      rp[i] = r;
459
5.09M
    }
460
5.09M
  while (++i < n);
461
462
109k
  return b;
463
109k
}
464
465
mp_limb_t
466
mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
467
377k
{
468
377k
  mp_size_t i;
469
377k
  mp_limb_t cy;
470
471
15.8M
  for (i = 0, cy = 0; i < n; i++)
472
15.4M
    {
473
15.4M
      mp_limb_t a, b, r;
474
15.4M
      a = ap[i]; b = bp[i];
475
15.4M
      r = a + cy;
476
15.4M
      cy = (r < cy);
477
15.4M
      r += b;
478
15.4M
      cy += (r < b);
479
15.4M
      rp[i] = r;
480
15.4M
    }
481
377k
  return cy;
482
377k
}
483
484
mp_limb_t
485
mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
486
377k
{
487
377k
  mp_limb_t cy;
488
489
377k
  assert (an >= bn);
490
491
377k
  cy = mpn_add_n (rp, ap, bp, bn);
492
377k
  if (an > bn)
493
4.68k
    cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy);
494
377k
  return cy;
495
377k
}
496
497
mp_limb_t
498
mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
499
9.29k
{
500
9.29k
  mp_size_t i;
501
502
9.29k
  assert (n > 0);
503
504
9.29k
  i = 0;
505
9.29k
  do
506
31.7k
    {
507
31.7k
      mp_limb_t a = ap[i];
508
      /* Carry out */
509
31.7k
      mp_limb_t cy = a < b;
510
31.7k
      rp[i] = a - b;
511
31.7k
      b = cy;
512
31.7k
    }
513
31.7k
  while (++i < n);
514
515
9.29k
  return b;
516
9.29k
}
517
518
mp_limb_t
519
mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
520
107k
{
521
107k
  mp_size_t i;
522
107k
  mp_limb_t cy;
523
524
2.63M
  for (i = 0, cy = 0; i < n; i++)
525
2.52M
    {
526
2.52M
      mp_limb_t a, b;
527
2.52M
      a = ap[i]; b = bp[i];
528
2.52M
      b += cy;
529
2.52M
      cy = (b < cy);
530
2.52M
      cy += (a < b);
531
2.52M
      rp[i] = a - b;
532
2.52M
    }
533
107k
  return cy;
534
107k
}
535
536
mp_limb_t
537
mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
538
107k
{
539
107k
  mp_limb_t cy;
540
541
107k
  assert (an >= bn);
542
543
107k
  cy = mpn_sub_n (rp, ap, bp, bn);
544
107k
  if (an > bn)
545
9.29k
    cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy);
546
107k
  return cy;
547
107k
}
548
549
mp_limb_t
550
mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
551
155k
{
552
155k
  mp_limb_t ul, cl, hpl, lpl;
553
554
155k
  assert (n >= 1);
555
556
155k
  cl = 0;
557
155k
  do
558
5.60M
    {
559
5.60M
      ul = *up++;
560
5.60M
      gmp_umul_ppmm (hpl, lpl, ul, vl);
561
562
5.60M
      lpl += cl;
563
5.60M
      cl = (lpl < cl) + hpl;
564
565
5.60M
      *rp++ = lpl;
566
5.60M
    }
567
5.60M
  while (--n != 0);
568
569
155k
  return cl;
570
155k
}
571
572
mp_limb_t
573
mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
574
243k
{
575
243k
  mp_limb_t ul, cl, hpl, lpl, rl;
576
577
243k
  assert (n >= 1);
578
579
243k
  cl = 0;
580
243k
  do
581
5.78M
    {
582
5.78M
      ul = *up++;
583
5.78M
      gmp_umul_ppmm (hpl, lpl, ul, vl);
584
585
5.78M
      lpl += cl;
586
5.78M
      cl = (lpl < cl) + hpl;
587
588
5.78M
      rl = *rp;
589
5.78M
      lpl = rl + lpl;
590
5.78M
      cl += lpl < rl;
591
5.78M
      *rp++ = lpl;
592
5.78M
    }
593
5.78M
  while (--n != 0);
594
595
243k
  return cl;
596
243k
}
597
598
mp_limb_t
599
mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
600
297k
{
601
297k
  mp_limb_t ul, cl, hpl, lpl, rl;
602
603
297k
  assert (n >= 1);
604
605
297k
  cl = 0;
606
297k
  do
607
6.17M
    {
608
6.17M
      ul = *up++;
609
6.17M
      gmp_umul_ppmm (hpl, lpl, ul, vl);
610
611
6.17M
      lpl += cl;
612
6.17M
      cl = (lpl < cl) + hpl;
613
614
6.17M
      rl = *rp;
615
6.17M
      lpl = rl - lpl;
616
6.17M
      cl += lpl > rl;
617
6.17M
      *rp++ = lpl;
618
6.17M
    }
619
6.17M
  while (--n != 0);
620
621
297k
  return cl;
622
297k
}
623
624
mp_limb_t
625
mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
626
50.4k
{
627
50.4k
  assert (un >= vn);
628
50.4k
  assert (vn >= 1);
629
50.4k
  assert (!GMP_MPN_OVERLAP_P(rp, un + vn, up, un));
630
50.4k
  assert (!GMP_MPN_OVERLAP_P(rp, un + vn, vp, vn));
631
632
  /* We first multiply by the low order limb. This result can be
633
     stored, not added, to rp. We also avoid a loop for zeroing this
634
     way. */
635
636
50.4k
  rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
637
638
  /* Now accumulate the product of up[] and the next higher limb from
639
     vp[]. */
640
641
294k
  while (--vn >= 1)
642
243k
    {
643
243k
      rp += 1, vp += 1;
644
243k
      rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
645
243k
    }
646
50.4k
  return rp[un];
647
50.4k
}
648
649
void
650
mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
651
0
{
652
0
  mpn_mul (rp, ap, n, bp, n);
653
0
}
654
655
void
656
mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n)
657
0
{
658
0
  mpn_mul (rp, ap, n, ap, n);
659
0
}
660
661
mp_limb_t
662
mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
663
438k
{
664
438k
  mp_limb_t high_limb, low_limb;
665
438k
  unsigned int tnc;
666
438k
  mp_limb_t retval;
667
668
438k
  assert (n >= 1);
669
438k
  assert (cnt >= 1);
670
438k
  assert (cnt < GMP_LIMB_BITS);
671
672
438k
  up += n;
673
438k
  rp += n;
674
675
438k
  tnc = GMP_LIMB_BITS - cnt;
676
438k
  low_limb = *--up;
677
438k
  retval = low_limb >> tnc;
678
438k
  high_limb = (low_limb << cnt);
679
680
16.2M
  while (--n != 0)
681
15.7M
    {
682
15.7M
      low_limb = *--up;
683
15.7M
      *--rp = high_limb | (low_limb >> tnc);
684
15.7M
      high_limb = (low_limb << cnt);
685
15.7M
    }
686
438k
  *--rp = high_limb;
687
688
438k
  return retval;
689
438k
}
690
691
mp_limb_t
692
mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
693
493k
{
694
493k
  mp_limb_t high_limb, low_limb;
695
493k
  unsigned int tnc;
696
493k
  mp_limb_t retval;
697
698
493k
  assert (n >= 1);
699
493k
  assert (cnt >= 1);
700
493k
  assert (cnt < GMP_LIMB_BITS);
701
702
493k
  tnc = GMP_LIMB_BITS - cnt;
703
493k
  high_limb = *up++;
704
493k
  retval = (high_limb << tnc);
705
493k
  low_limb = high_limb >> cnt;
706
707
22.7M
  while (--n != 0)
708
22.2M
    {
709
22.2M
      high_limb = *up++;
710
22.2M
      *rp++ = low_limb | (high_limb << tnc);
711
22.2M
      low_limb = high_limb >> cnt;
712
22.2M
    }
713
493k
  *rp = low_limb;
714
715
493k
  return retval;
716
493k
}
717
718
static mp_bitcnt_t
719
mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un,
720
     mp_limb_t ux)
721
92.9k
{
722
92.9k
  unsigned cnt;
723
724
92.9k
  assert (ux == 0 || ux == GMP_LIMB_MAX);
725
92.9k
  assert (0 <= i && i <= un );
726
727
93.0k
  while (limb == 0)
728
102
    {
729
102
      i++;
730
102
      if (i == un)
731
0
  return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
732
102
      limb = ux ^ up[i];
733
102
    }
734
92.9k
  gmp_ctz (cnt, limb);
735
92.9k
  return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
736
92.9k
}
737
738
mp_bitcnt_t
739
mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit)
740
92.9k
{
741
92.9k
  mp_size_t i;
742
92.9k
  i = bit / GMP_LIMB_BITS;
743
744
92.9k
  return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
745
92.9k
        i, ptr, i, 0);
746
92.9k
}
747
748
mp_bitcnt_t
749
mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
750
0
{
751
0
  mp_size_t i;
752
0
  i = bit / GMP_LIMB_BITS;
753
754
0
  return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
755
0
        i, ptr, i, GMP_LIMB_MAX);
756
0
}
757
758
void
759
mpn_com (mp_ptr rp, mp_srcptr up, mp_size_t n)
760
0
{
761
0
  while (--n >= 0)
762
0
    *rp++ = ~ *up++;
763
0
}
764
765
mp_limb_t
766
mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n)
767
0
{
768
0
  while (*up == 0)
769
0
    {
770
0
      *rp = 0;
771
0
      if (!--n)
772
0
  return 0;
773
0
      ++up; ++rp;
774
0
    }
775
0
  *rp = - *up;
776
0
  mpn_com (++rp, ++up, --n);
777
0
  return 1;
778
0
}
779
780

781
/* MPN division interface. */
782
783
/* The 3/2 inverse is defined as
784
785
     m = floor( (B^3-1) / (B u1 + u0)) - B
786
*/
787
mp_limb_t
788
mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
789
12.8k
{
790
12.8k
  mp_limb_t r, m;
791
792
12.8k
  {
793
12.8k
    mp_limb_t p, ql;
794
12.8k
    unsigned ul, uh, qh;
795
796
12.8k
    assert (sizeof (unsigned) * 2 >= sizeof (mp_limb_t));
797
    /* For notation, let b denote the half-limb base, so that B = b^2.
798
       Split u1 = b uh + ul. */
799
12.8k
    ul = u1 & GMP_LLIMB_MASK;
800
12.8k
    uh = u1 >> (GMP_LIMB_BITS / 2);
801
802
    /* Approximation of the high half of quotient. Differs from the 2/1
803
       inverse of the half limb uh, since we have already subtracted
804
       u0. */
805
12.8k
    qh = (u1 ^ GMP_LIMB_MAX) / uh;
806
807
    /* Adjust to get a half-limb 3/2 inverse, i.e., we want
808
809
       qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u
810
     = floor( (b (~u) + b-1) / u),
811
812
       and the remainder
813
814
       r = b (~u) + b-1 - qh (b uh + ul)
815
       = b (~u - qh uh) + b-1 - qh ul
816
817
       Subtraction of qh ul may underflow, which implies adjustments.
818
       But by normalization, 2 u >= B > qh ul, so we need to adjust by
819
       at most 2.
820
    */
821
822
12.8k
    r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
823
824
12.8k
    p = (mp_limb_t) qh * ul;
825
    /* Adjustment steps taken from udiv_qrnnd_c */
826
12.8k
    if (r < p)
827
3.55k
      {
828
3.55k
  qh--;
829
3.55k
  r += u1;
830
3.55k
  if (r >= u1) /* i.e. we didn't get carry when adding to r */
831
3.50k
    if (r < p)
832
51
      {
833
51
        qh--;
834
51
        r += u1;
835
51
      }
836
3.55k
      }
837
12.8k
    r -= p;
838
839
    /* Low half of the quotient is
840
841
       ql = floor ( (b r + b-1) / u1).
842
843
       This is a 3/2 division (on half-limbs), for which qh is a
844
       suitable inverse. */
845
846
12.8k
    p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
847
    /* Unlike full-limb 3/2, we can add 1 without overflow. For this to
848
       work, it is essential that ql is a full mp_limb_t. */
849
12.8k
    ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
850
851
    /* By the 3/2 trick, we don't need the high half limb. */
852
12.8k
    r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
853
854
12.8k
    if (r >= (GMP_LIMB_MAX & (p << (GMP_LIMB_BITS / 2))))
855
3.71k
      {
856
3.71k
  ql--;
857
3.71k
  r += u1;
858
3.71k
      }
859
12.8k
    m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql;
860
12.8k
    if (r >= u1)
861
13
      {
862
13
  m++;
863
13
  r -= u1;
864
13
      }
865
12.8k
  }
866
867
  /* Now m is the 2/1 inverse of u1. If u0 > 0, adjust it to become a
868
     3/2 inverse. */
869
12.8k
  if (u0 > 0)
870
10.8k
    {
871
10.8k
      mp_limb_t th, tl;
872
10.8k
      r = ~r;
873
10.8k
      r += u0;
874
10.8k
      if (r < u0)
875
9.28k
  {
876
9.28k
    m--;
877
9.28k
    if (r >= u1)
878
499
      {
879
499
        m--;
880
499
        r -= u1;
881
499
      }
882
9.28k
    r -= u1;
883
9.28k
  }
884
10.8k
      gmp_umul_ppmm (th, tl, u0, m);
885
10.8k
      r += th;
886
10.8k
      if (r < th)
887
3.79k
  {
888
3.79k
    m--;
889
3.79k
    m -= ((r > u1) | ((r == u1) & (tl > u0)));
890
3.79k
  }
891
10.8k
    }
892
893
12.8k
  return m;
894
12.8k
}
895
896
struct gmp_div_inverse
897
{
898
  /* Normalization shift count. */
899
  unsigned shift;
900
  /* Normalized divisor (d0 unused for mpn_div_qr_1) */
901
  mp_limb_t d1, d0;
902
  /* Inverse, for 2/1 or 3/2. */
903
  mp_limb_t di;
904
};
905
906
static void
907
mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
908
1.92k
{
909
1.92k
  unsigned shift;
910
911
1.92k
  assert (d > 0);
912
1.92k
  gmp_clz (shift, d);
913
1.92k
  inv->shift = shift;
914
1.92k
  inv->d1 = d << shift;
915
1.92k
  inv->di = mpn_invert_limb (inv->d1);
916
1.92k
}
917
918
static void
919
mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
920
         mp_limb_t d1, mp_limb_t d0)
921
998
{
922
998
  unsigned shift;
923
924
998
  assert (d1 > 0);
925
998
  gmp_clz (shift, d1);
926
998
  inv->shift = shift;
927
998
  if (shift > 0)
928
576
    {
929
576
      d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
930
576
      d0 <<= shift;
931
576
    }
932
998
  inv->d1 = d1;
933
998
  inv->d0 = d0;
934
998
  inv->di = mpn_invert_3by2 (d1, d0);
935
998
}
936
937
static void
938
mpn_div_qr_invert (struct gmp_div_inverse *inv,
939
       mp_srcptr dp, mp_size_t dn)
940
11.6k
{
941
11.6k
  assert (dn > 0);
942
943
11.6k
  if (dn == 1)
944
666
    mpn_div_qr_1_invert (inv, dp[0]);
945
10.9k
  else if (dn == 2)
946
998
    mpn_div_qr_2_invert (inv, dp[1], dp[0]);
947
9.96k
  else
948
9.96k
    {
949
9.96k
      unsigned shift;
950
9.96k
      mp_limb_t d1, d0;
951
952
9.96k
      d1 = dp[dn-1];
953
9.96k
      d0 = dp[dn-2];
954
9.96k
      assert (d1 > 0);
955
9.96k
      gmp_clz (shift, d1);
956
9.96k
      inv->shift = shift;
957
9.96k
      if (shift > 0)
958
9.92k
  {
959
9.92k
    d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
960
9.92k
    d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
961
9.92k
  }
962
9.96k
      inv->d1 = d1;
963
9.96k
      inv->d0 = d0;
964
9.96k
      inv->di = mpn_invert_3by2 (d1, d0);
965
9.96k
    }
966
11.6k
}
967
968
/* Not matching current public gmp interface, rather corresponding to
969
   the sbpi1_div_* functions. */
970
static mp_limb_t
971
mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
972
         const struct gmp_div_inverse *inv)
973
245k
{
974
245k
  mp_limb_t d, di;
975
245k
  mp_limb_t r;
976
245k
  mp_ptr tp = NULL;
977
245k
  mp_size_t tn = 0;
978
979
245k
  if (inv->shift > 0)
980
233k
    {
981
      /* Shift, reusing qp area if possible. In-place shift if qp == np. */
982
233k
      tp = qp;
983
233k
      if (!tp)
984
376
        {
985
376
     tn = nn;
986
376
     tp = gmp_alloc_limbs (tn);
987
376
        }
988
233k
      r = mpn_lshift (tp, np, nn, inv->shift);
989
233k
      np = tp;
990
233k
    }
991
12.1k
  else
992
12.1k
    r = 0;
993
994
245k
  d = inv->d1;
995
245k
  di = inv->di;
996
11.4M
  while (--nn >= 0)
997
11.2M
    {
998
11.2M
      mp_limb_t q;
999
1000
11.2M
      gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di);
1001
11.2M
      if (qp)
1002
11.2M
  qp[nn] = q;
1003
11.2M
    }
1004
245k
  if (tn)
1005
376
    gmp_free_limbs (tp, tn);
1006
1007
245k
  return r >> inv->shift;
1008
245k
}
1009
1010
static void
1011
mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
1012
         const struct gmp_div_inverse *inv)
1013
1.91k
{
1014
1.91k
  unsigned shift;
1015
1.91k
  mp_size_t i;
1016
1.91k
  mp_limb_t d1, d0, di, r1, r0;
1017
1018
1.91k
  assert (nn >= 2);
1019
1.91k
  shift = inv->shift;
1020
1.91k
  d1 = inv->d1;
1021
1.91k
  d0 = inv->d0;
1022
1.91k
  di = inv->di;
1023
1024
1.91k
  if (shift > 0)
1025
576
    r1 = mpn_lshift (np, np, nn, shift);
1026
1.33k
  else
1027
1.33k
    r1 = 0;
1028
1029
1.91k
  r0 = np[nn - 1];
1030
1031
1.91k
  i = nn - 2;
1032
1.91k
  do
1033
6.16k
    {
1034
6.16k
      mp_limb_t n0, q;
1035
6.16k
      n0 = np[i];
1036
6.16k
      gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di);
1037
1038
6.16k
      if (qp)
1039
511
  qp[i] = q;
1040
6.16k
    }
1041
6.16k
  while (--i >= 0);
1042
1043
1.91k
  if (shift > 0)
1044
576
    {
1045
576
      assert ((r0 & (GMP_LIMB_MAX >> (GMP_LIMB_BITS - shift))) == 0);
1046
576
      r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift));
1047
576
      r1 >>= shift;
1048
576
    }
1049
1050
1.91k
  np[1] = r1;
1051
1.91k
  np[0] = r0;
1052
1.91k
}
1053
1054
static void
1055
mpn_div_qr_pi1 (mp_ptr qp,
1056
    mp_ptr np, mp_size_t nn, mp_limb_t n1,
1057
    mp_srcptr dp, mp_size_t dn,
1058
    mp_limb_t dinv)
1059
21.9k
{
1060
21.9k
  mp_size_t i;
1061
1062
21.9k
  mp_limb_t d1, d0;
1063
21.9k
  mp_limb_t cy, cy1;
1064
21.9k
  mp_limb_t q;
1065
1066
21.9k
  assert (dn > 2);
1067
21.9k
  assert (nn >= dn);
1068
1069
21.9k
  d1 = dp[dn - 1];
1070
21.9k
  d0 = dp[dn - 2];
1071
1072
21.9k
  assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
1073
  /* Iteration variable is the index of the q limb.
1074
   *
1075
   * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
1076
   * by            <d1,          d0,        dp[dn-3],  ..., dp[0] >
1077
   */
1078
1079
21.9k
  i = nn - dn;
1080
21.9k
  do
1081
297k
    {
1082
297k
      mp_limb_t n0 = np[dn-1+i];
1083
1084
297k
      if (n1 == d1 && n0 == d0)
1085
116
  {
1086
116
    q = GMP_LIMB_MAX;
1087
116
    mpn_submul_1 (np+i, dp, dn, q);
1088
116
    n1 = np[dn-1+i];  /* update n1, last loop's value will now be invalid */
1089
116
  }
1090
297k
      else
1091
297k
  {
1092
297k
    gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv);
1093
1094
297k
    cy = mpn_submul_1 (np + i, dp, dn-2, q);
1095
1096
297k
    cy1 = n0 < cy;
1097
297k
    n0 = n0 - cy;
1098
297k
    cy = n1 < cy1;
1099
297k
    n1 = n1 - cy1;
1100
297k
    np[dn-2+i] = n0;
1101
1102
297k
    if (cy != 0)
1103
41
      {
1104
41
        n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1);
1105
41
        q--;
1106
41
      }
1107
297k
  }
1108
1109
297k
      if (qp)
1110
13.2k
  qp[i] = q;
1111
297k
    }
1112
297k
  while (--i >= 0);
1113
1114
21.9k
  np[dn - 1] = n1;
1115
21.9k
}
1116
1117
static void
1118
mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
1119
       mp_srcptr dp, mp_size_t dn,
1120
       const struct gmp_div_inverse *inv)
1121
24.6k
{
1122
24.6k
  assert (dn > 0);
1123
24.6k
  assert (nn >= dn);
1124
1125
24.6k
  if (dn == 1)
1126
783
    np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
1127
23.8k
  else if (dn == 2)
1128
1.91k
    mpn_div_qr_2_preinv (qp, np, nn, inv);
1129
21.9k
  else
1130
21.9k
    {
1131
21.9k
      mp_limb_t nh;
1132
21.9k
      unsigned shift;
1133
1134
21.9k
      assert (inv->d1 == dp[dn-1]);
1135
21.9k
      assert (inv->d0 == dp[dn-2]);
1136
21.9k
      assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);
1137
1138
21.9k
      shift = inv->shift;
1139
21.9k
      if (shift > 0)
1140
9.92k
  nh = mpn_lshift (np, np, nn, shift);
1141
12.0k
      else
1142
12.0k
  nh = 0;
1143
1144
21.9k
      mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
1145
1146
21.9k
      if (shift > 0)
1147
9.92k
  gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
1148
21.9k
    }
1149
24.6k
}
1150
1151
static void
1152
mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
1153
11.4k
{
1154
11.4k
  struct gmp_div_inverse inv;
1155
11.4k
  mp_ptr tp = NULL;
1156
1157
11.4k
  assert (dn > 0);
1158
11.4k
  assert (nn >= dn);
1159
1160
11.4k
  mpn_div_qr_invert (&inv, dp, dn);
1161
11.4k
  if (dn > 2 && inv.shift > 0)
1162
9.80k
    {
1163
9.80k
      tp = gmp_alloc_limbs (dn);
1164
9.80k
      gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
1165
9.80k
      dp = tp;
1166
9.80k
    }
1167
11.4k
  mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
1168
11.4k
  if (tp)
1169
9.80k
    gmp_free_limbs (tp, dn);
1170
11.4k
}
1171
1172

1173
/* MPN base conversion. */
1174
static unsigned
1175
mpn_base_power_of_two_p (unsigned b)
1176
5.27k
{
1177
5.27k
  switch (b)
1178
5.27k
    {
1179
0
    case 2: return 1;
1180
0
    case 4: return 2;
1181
0
    case 8: return 3;
1182
0
    case 16: return 4;
1183
0
    case 32: return 5;
1184
0
    case 64: return 6;
1185
0
    case 128: return 7;
1186
0
    case 256: return 8;
1187
5.27k
    default: return 0;
1188
5.27k
    }
1189
5.27k
}
1190
1191
struct mpn_base_info
1192
{
1193
  /* bb is the largest power of the base which fits in one limb, and
1194
     exp is the corresponding exponent. */
1195
  unsigned exp;
1196
  mp_limb_t bb;
1197
};
1198
1199
static void
1200
mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b)
1201
5.27k
{
1202
5.27k
  mp_limb_t m;
1203
5.27k
  mp_limb_t p;
1204
5.27k
  unsigned exp;
1205
1206
5.27k
  m = GMP_LIMB_MAX / b;
1207
100k
  for (exp = 1, p = b; p <= m; exp++)
1208
94.9k
    p *= b;
1209
1210
5.27k
  info->exp = exp;
1211
5.27k
  info->bb = p;
1212
5.27k
}
1213
1214
static mp_bitcnt_t
1215
mpn_limb_size_in_base_2 (mp_limb_t u)
1216
1.42k
{
1217
1.42k
  unsigned shift;
1218
1219
1.42k
  assert (u > 0);
1220
1.42k
  gmp_clz (shift, u);
1221
1.42k
  return GMP_LIMB_BITS - shift;
1222
1.42k
}
1223
1224
static size_t
1225
mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
1226
0
{
1227
0
  unsigned char mask;
1228
0
  size_t sn, j;
1229
0
  mp_size_t i;
1230
0
  unsigned shift;
1231
1232
0
  sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1233
0
  + bits - 1) / bits;
1234
1235
0
  mask = (1U << bits) - 1;
1236
1237
0
  for (i = 0, j = sn, shift = 0; j-- > 0;)
1238
0
    {
1239
0
      unsigned char digit = up[i] >> shift;
1240
1241
0
      shift += bits;
1242
1243
0
      if (shift >= GMP_LIMB_BITS && ++i < un)
1244
0
  {
1245
0
    shift -= GMP_LIMB_BITS;
1246
0
    digit |= up[i] << (bits - shift);
1247
0
  }
1248
0
      sp[j] = digit & mask;
1249
0
    }
1250
0
  return sn;
1251
0
}
1252
1253
/* We generate digits from the least significant end, and reverse at
1254
   the end. */
1255
static size_t
1256
mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
1257
      const struct gmp_div_inverse *binv)
1258
12.4k
{
1259
12.4k
  mp_size_t i;
1260
199k
  for (i = 0; w > 0; i++)
1261
186k
    {
1262
186k
      mp_limb_t h, l, r;
1263
1264
186k
      h = w >> (GMP_LIMB_BITS - binv->shift);
1265
186k
      l = w << binv->shift;
1266
1267
186k
      gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
1268
186k
      assert ((r & (GMP_LIMB_MAX >> (GMP_LIMB_BITS - binv->shift))) == 0);
1269
186k
      r >>= binv->shift;
1270
1271
186k
      sp[i] = r;
1272
186k
    }
1273
12.4k
  return i;
1274
12.4k
}
1275
1276
static size_t
1277
mpn_get_str_other (unsigned char *sp,
1278
       int base, const struct mpn_base_info *info,
1279
       mp_ptr up, mp_size_t un)
1280
445
{
1281
445
  struct gmp_div_inverse binv;
1282
445
  size_t sn;
1283
445
  size_t i;
1284
1285
445
  mpn_div_qr_1_invert (&binv, base);
1286
1287
445
  sn = 0;
1288
1289
445
  if (un > 1)
1290
365
    {
1291
365
      struct gmp_div_inverse bbinv;
1292
365
      mpn_div_qr_1_invert (&bbinv, info->bb);
1293
1294
365
      do
1295
12.0k
  {
1296
12.0k
    mp_limb_t w;
1297
12.0k
    size_t done;
1298
12.0k
    w = mpn_div_qr_1_preinv (up, up, un, &bbinv);
1299
12.0k
    un -= (up[un-1] == 0);
1300
12.0k
    done = mpn_limb_get_str (sp + sn, w, &binv);
1301
1302
57.6k
    for (sn += done; done < info->exp; done++)
1303
45.6k
      sp[sn++] = 0;
1304
12.0k
  }
1305
12.0k
      while (un > 1);
1306
365
    }
1307
445
  sn += mpn_limb_get_str (sp + sn, up[0], &binv);
1308
1309
  /* Reverse order */
1310
116k
  for (i = 0; 2*i + 1 < sn; i++)
1311
116k
    {
1312
116k
      unsigned char t = sp[i];
1313
116k
      sp[i] = sp[sn - i - 1];
1314
116k
      sp[sn - i - 1] = t;
1315
116k
    }
1316
1317
445
  return sn;
1318
445
}
1319
1320
size_t
1321
mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un)
1322
0
{
1323
0
  unsigned bits;
1324
1325
0
  assert (un > 0);
1326
0
  assert (up[un-1] > 0);
1327
1328
0
  bits = mpn_base_power_of_two_p (base);
1329
0
  if (bits)
1330
0
    return mpn_get_str_bits (sp, bits, up, un);
1331
0
  else
1332
0
    {
1333
0
      struct mpn_base_info info;
1334
1335
0
      mpn_get_base_info (&info, base);
1336
0
      return mpn_get_str_other (sp, base, &info, up, un);
1337
0
    }
1338
0
}
1339
1340
static mp_size_t
1341
mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
1342
      unsigned bits)
1343
0
{
1344
0
  mp_size_t rn;
1345
0
  mp_limb_t limb;
1346
0
  unsigned shift;
1347
1348
0
  for (limb = 0, rn = 0, shift = 0; sn-- > 0; )
1349
0
    {
1350
0
      limb |= (mp_limb_t) sp[sn] << shift;
1351
0
      shift += bits;
1352
0
      if (shift >= GMP_LIMB_BITS)
1353
0
  {
1354
0
    shift -= GMP_LIMB_BITS;
1355
0
    rp[rn++] = limb;
1356
    /* Next line is correct also if shift == 0,
1357
       bits == 8, and mp_limb_t == unsigned char. */
1358
0
    limb = (unsigned int) sp[sn] >> (bits - shift);
1359
0
  }
1360
0
    }
1361
0
  if (limb != 0)
1362
0
    rp[rn++] = limb;
1363
0
  else
1364
0
    rn = mpn_normalized_size (rp, rn);
1365
0
  return rn;
1366
0
}
1367
1368
/* Result is usually normalized, except for all-zero input, in which
1369
   case a single zero limb is written at *RP, and 1 is returned. */
1370
static mp_size_t
1371
mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1372
       mp_limb_t b, const struct mpn_base_info *info)
1373
4.83k
{
1374
4.83k
  mp_size_t rn;
1375
4.83k
  mp_limb_t w;
1376
4.83k
  unsigned k;
1377
4.83k
  size_t j;
1378
1379
4.83k
  assert (sn > 0);
1380
1381
4.83k
  k = 1 + (sn - 1) % info->exp;
1382
1383
4.83k
  j = 0;
1384
4.83k
  w = sp[j++];
1385
31.2k
  while (--k != 0)
1386
26.4k
    w = w * b + sp[j++];
1387
1388
4.83k
  rp[0] = w;
1389
1390
109k
  for (rn = 1; j < sn;)
1391
104k
    {
1392
104k
      mp_limb_t cy;
1393
1394
104k
      w = sp[j++];
1395
1.98M
      for (k = 1; k < info->exp; k++)
1396
1.88M
  w = w * b + sp[j++];
1397
1398
104k
      cy = mpn_mul_1 (rp, rp, rn, info->bb);
1399
104k
      cy += mpn_add_1 (rp, rp, rn, w);
1400
104k
      if (cy > 0)
1401
102k
  rp[rn++] = cy;
1402
104k
    }
1403
4.83k
  assert (j == sn);
1404
1405
4.83k
  return rn;
1406
4.83k
}
1407
1408
mp_size_t
1409
mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
1410
0
{
1411
0
  unsigned bits;
1412
1413
0
  if (sn == 0)
1414
0
    return 0;
1415
1416
0
  bits = mpn_base_power_of_two_p (base);
1417
0
  if (bits)
1418
0
    return mpn_set_str_bits (rp, sp, sn, bits);
1419
0
  else
1420
0
    {
1421
0
      struct mpn_base_info info;
1422
1423
0
      mpn_get_base_info (&info, base);
1424
0
      return mpn_set_str_other (rp, sp, sn, base, &info);
1425
0
    }
1426
0
}
1427
1428

1429
/* MPZ interface */
1430
void
1431
mpz_init (mpz_t r)
1432
24.8k
{
1433
24.8k
  static const mp_limb_t dummy_limb = GMP_LIMB_MAX & 0xc1a0;
1434
1435
24.8k
  r->_mp_alloc = 0;
1436
24.8k
  r->_mp_size = 0;
1437
24.8k
  r->_mp_d = (mp_ptr) &dummy_limb;
1438
24.8k
}
1439
1440
/* The utility of this function is a bit limited, since many functions
1441
   assigns the result variable using mpz_swap. */
1442
void
1443
mpz_init2 (mpz_t r, mp_bitcnt_t bits)
1444
51.3k
{
1445
51.3k
  mp_size_t rn;
1446
1447
51.3k
  bits -= (bits != 0);    /* Round down, except if 0 */
1448
51.3k
  rn = 1 + bits / GMP_LIMB_BITS;
1449
1450
51.3k
  r->_mp_alloc = rn;
1451
51.3k
  r->_mp_size = 0;
1452
51.3k
  r->_mp_d = gmp_alloc_limbs (rn);
1453
51.3k
}
1454
1455
void
1456
mpz_clear (mpz_t r)
1457
80.3k
{
1458
80.3k
  if (r->_mp_alloc)
1459
75.4k
    gmp_free_limbs (r->_mp_d, r->_mp_alloc);
1460
80.3k
}
1461
1462
static mp_ptr
1463
mpz_realloc (mpz_t r, mp_size_t size)
1464
33.0k
{
1465
33.0k
  size = GMP_MAX (size, 1);
1466
1467
33.0k
  if (r->_mp_alloc)
1468
8.98k
    r->_mp_d = gmp_realloc_limbs (r->_mp_d, r->_mp_alloc, size);
1469
24.0k
  else
1470
24.0k
    r->_mp_d = gmp_alloc_limbs (size);
1471
33.0k
  r->_mp_alloc = size;
1472
1473
33.0k
  if (GMP_ABS (r->_mp_size) > size)
1474
0
    r->_mp_size = 0;
1475
1476
33.0k
  return r->_mp_d;
1477
33.0k
}
1478
1479
/* Realloc for an mpz_t WHAT if it has less than NEEDED limbs.  */
1480
1.16M
#define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc      \
1481
1.16M
        ? mpz_realloc(z,n)      \
1482
1.16M
        : (z)->_mp_d)
1483

1484
/* MPZ assignment and basic conversions. */
1485
void
1486
mpz_set_si (mpz_t r, signed long int x)
1487
21
{
1488
21
  if (x >= 0)
1489
18
    mpz_set_ui (r, x);
1490
3
  else /* (x < 0) */
1491
3
    if (GMP_LIMB_BITS < GMP_ULONG_BITS)
1492
0
      {
1493
0
  mpz_set_ui (r, GMP_NEG_CAST (unsigned long int, x));
1494
0
  mpz_neg (r, r);
1495
0
      }
1496
3
  else
1497
3
    {
1498
3
      r->_mp_size = -1;
1499
3
      MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x);
1500
3
    }
1501
21
}
1502
1503
void
1504
mpz_set_ui (mpz_t r, unsigned long int x)
1505
6.90k
{
1506
6.90k
  if (x > 0)
1507
6.89k
    {
1508
6.89k
      r->_mp_size = 1;
1509
6.89k
      MPZ_REALLOC (r, 1)[0] = x;
1510
6.89k
      if (GMP_LIMB_BITS < GMP_ULONG_BITS)
1511
0
  {
1512
0
    int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS;
1513
0
    while (x >>= LOCAL_GMP_LIMB_BITS)
1514
0
      {
1515
0
        ++ r->_mp_size;
1516
0
        MPZ_REALLOC (r, r->_mp_size)[r->_mp_size - 1] = x;
1517
0
      }
1518
0
  }
1519
6.89k
    }
1520
8
  else
1521
8
    r->_mp_size = 0;
1522
6.90k
}
1523
1524
void
1525
mpz_set (mpz_t r, const mpz_t x)
1526
206k
{
1527
  /* Allow the NOP r == x */
1528
206k
  if (r != x)
1529
12.1k
    {
1530
12.1k
      mp_size_t n;
1531
12.1k
      mp_ptr rp;
1532
1533
12.1k
      n = GMP_ABS (x->_mp_size);
1534
12.1k
      rp = MPZ_REALLOC (r, n);
1535
1536
12.1k
      mpn_copyi (rp, x->_mp_d, n);
1537
12.1k
      r->_mp_size = x->_mp_size;
1538
12.1k
    }
1539
206k
}
1540
1541
void
1542
mpz_init_set_si (mpz_t r, signed long int x)
1543
2
{
1544
2
  mpz_init (r);
1545
2
  mpz_set_si (r, x);
1546
2
}
1547
1548
void
1549
mpz_init_set_ui (mpz_t r, unsigned long int x)
1550
6.86k
{
1551
6.86k
  mpz_init (r);
1552
6.86k
  mpz_set_ui (r, x);
1553
6.86k
}
1554
1555
void
1556
mpz_init_set (mpz_t r, const mpz_t x)
1557
11.4k
{
1558
11.4k
  mpz_init (r);
1559
11.4k
  mpz_set (r, x);
1560
11.4k
}
1561
1562
int
1563
mpz_fits_slong_p (const mpz_t u)
1564
21
{
1565
21
  return mpz_cmp_si (u, LONG_MAX) <= 0 && mpz_cmp_si (u, LONG_MIN) >= 0;
1566
21
}
1567
1568
static int
1569
mpn_absfits_ulong_p (mp_srcptr up, mp_size_t un)
1570
885
{
1571
885
  int ulongsize = GMP_ULONG_BITS / GMP_LIMB_BITS;
1572
885
  mp_limb_t ulongrem = 0;
1573
1574
885
  if (GMP_ULONG_BITS % GMP_LIMB_BITS != 0)
1575
0
    ulongrem = (mp_limb_t) (ULONG_MAX >> GMP_LIMB_BITS * ulongsize) + 1;
1576
1577
885
  return un <= ulongsize || (up[ulongsize] < ulongrem && un == ulongsize + 1);
1578
885
}
1579
1580
int
1581
mpz_fits_ulong_p (const mpz_t u)
1582
89
{
1583
89
  mp_size_t us = u->_mp_size;
1584
1585
89
  return us >= 0 && mpn_absfits_ulong_p (u->_mp_d, us);
1586
89
}
1587
1588
int
1589
mpz_fits_sint_p (const mpz_t u)
1590
0
{
1591
0
  return mpz_cmp_si (u, INT_MAX) <= 0 && mpz_cmp_si (u, INT_MIN) >= 0;
1592
0
}
1593
1594
int
1595
mpz_fits_uint_p (const mpz_t u)
1596
0
{
1597
0
  return u->_mp_size >= 0 && mpz_cmpabs_ui (u, UINT_MAX) <= 0;
1598
0
}
1599
1600
int
1601
mpz_fits_sshort_p (const mpz_t u)
1602
0
{
1603
0
  return mpz_cmp_si (u, SHRT_MAX) <= 0 && mpz_cmp_si (u, SHRT_MIN) >= 0;
1604
0
}
1605
1606
int
1607
mpz_fits_ushort_p (const mpz_t u)
1608
0
{
1609
0
  return u->_mp_size >= 0 && mpz_cmpabs_ui (u, USHRT_MAX) <= 0;
1610
0
}
1611
1612
long int
1613
mpz_get_si (const mpz_t u)
1614
11
{
1615
11
  unsigned long r = mpz_get_ui (u);
1616
11
  unsigned long c = -LONG_MAX - LONG_MIN;
1617
1618
11
  if (u->_mp_size < 0)
1619
    /* This expression is necessary to properly handle -LONG_MIN */
1620
0
    return -(long) c - (long) ((r - c) & LONG_MAX);
1621
11
  else
1622
11
    return (long) (r & LONG_MAX);
1623
11
}
1624
1625
unsigned long int
1626
mpz_get_ui (const mpz_t u)
1627
315
{
1628
315
  if (GMP_LIMB_BITS < GMP_ULONG_BITS)
1629
0
    {
1630
0
      int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS;
1631
0
      unsigned long r = 0;
1632
0
      mp_size_t n = GMP_ABS (u->_mp_size);
1633
0
      n = GMP_MIN (n, 1 + (mp_size_t) (GMP_ULONG_BITS - 1) / GMP_LIMB_BITS);
1634
0
      while (--n >= 0)
1635
0
  r = (r << LOCAL_GMP_LIMB_BITS) + u->_mp_d[n];
1636
0
      return r;
1637
0
    }
1638
1639
315
  return u->_mp_size == 0 ? 0 : u->_mp_d[0];
1640
315
}
1641
1642
size_t
1643
mpz_size (const mpz_t u)
1644
0
{
1645
0
  return GMP_ABS (u->_mp_size);
1646
0
}
1647
1648
mp_limb_t
1649
mpz_getlimbn (const mpz_t u, mp_size_t n)
1650
0
{
1651
0
  if (n >= 0 && n < GMP_ABS (u->_mp_size))
1652
0
    return u->_mp_d[n];
1653
0
  else
1654
0
    return 0;
1655
0
}
1656
1657
void
1658
mpz_realloc2 (mpz_t x, mp_bitcnt_t n)
1659
0
{
1660
0
  mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS);
1661
0
}
1662
1663
mp_srcptr
1664
mpz_limbs_read (mpz_srcptr x)
1665
0
{
1666
0
  return x->_mp_d;
1667
0
}
1668
1669
mp_ptr
1670
mpz_limbs_modify (mpz_t x, mp_size_t n)
1671
0
{
1672
0
  assert (n > 0);
1673
0
  return MPZ_REALLOC (x, n);
1674
0
}
1675
1676
mp_ptr
1677
mpz_limbs_write (mpz_t x, mp_size_t n)
1678
0
{
1679
0
  return mpz_limbs_modify (x, n);
1680
0
}
1681
1682
void
1683
mpz_limbs_finish (mpz_t x, mp_size_t xs)
1684
0
{
1685
0
  mp_size_t xn;
1686
0
  xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs));
1687
0
  x->_mp_size = xs < 0 ? -xn : xn;
1688
0
}
1689
1690
static mpz_srcptr
1691
mpz_roinit_normal_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
1692
17
{
1693
17
  x->_mp_alloc = 0;
1694
17
  x->_mp_d = (mp_ptr) xp;
1695
17
  x->_mp_size = xs;
1696
17
  return x;
1697
17
}
1698
1699
mpz_srcptr
1700
mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
1701
0
{
1702
0
  mpz_roinit_normal_n (x, xp, xs);
1703
0
  mpz_limbs_finish (x, xs);
1704
0
  return x;
1705
0
}
1706
1707

1708
/* Conversions and comparison to double. */
1709
void
1710
mpz_set_d (mpz_t r, double x)
1711
0
{
1712
0
  int sign;
1713
0
  mp_ptr rp;
1714
0
  mp_size_t rn, i;
1715
0
  double B;
1716
0
  double Bi;
1717
0
  mp_limb_t f;
1718
1719
  /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
1720
     zero or infinity. */
1721
0
  if (x != x || x == x * 0.5)
1722
0
    {
1723
0
      r->_mp_size = 0;
1724
0
      return;
1725
0
    }
1726
1727
0
  sign = x < 0.0 ;
1728
0
  if (sign)
1729
0
    x = - x;
1730
1731
0
  if (x < 1.0)
1732
0
    {
1733
0
      r->_mp_size = 0;
1734
0
      return;
1735
0
    }
1736
0
  B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
1737
0
  Bi = 1.0 / B;
1738
0
  for (rn = 1; x >= B; rn++)
1739
0
    x *= Bi;
1740
1741
0
  rp = MPZ_REALLOC (r, rn);
1742
1743
0
  f = (mp_limb_t) x;
1744
0
  x -= f;
1745
0
  assert (x < 1.0);
1746
0
  i = rn-1;
1747
0
  rp[i] = f;
1748
0
  while (--i >= 0)
1749
0
    {
1750
0
      x = B * x;
1751
0
      f = (mp_limb_t) x;
1752
0
      x -= f;
1753
0
      assert (x < 1.0);
1754
0
      rp[i] = f;
1755
0
    }
1756
1757
0
  r->_mp_size = sign ? - rn : rn;
1758
0
}
1759
1760
void
1761
mpz_init_set_d (mpz_t r, double x)
1762
0
{
1763
0
  mpz_init (r);
1764
0
  mpz_set_d (r, x);
1765
0
}
1766
1767
double
1768
mpz_get_d (const mpz_t u)
1769
0
{
1770
0
  int m;
1771
0
  mp_limb_t l;
1772
0
  mp_size_t un;
1773
0
  double x;
1774
0
  double B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
1775
1776
0
  un = GMP_ABS (u->_mp_size);
1777
1778
0
  if (un == 0)
1779
0
    return 0.0;
1780
1781
0
  l = u->_mp_d[--un];
1782
0
  gmp_clz (m, l);
1783
0
  m = m + GMP_DBL_MANT_BITS - GMP_LIMB_BITS;
1784
0
  if (m < 0)
1785
0
    l &= GMP_LIMB_MAX << -m;
1786
1787
0
  for (x = l; --un >= 0;)
1788
0
    {
1789
0
      x = B*x;
1790
0
      if (m > 0) {
1791
0
  l = u->_mp_d[un];
1792
0
  m -= GMP_LIMB_BITS;
1793
0
  if (m < 0)
1794
0
    l &= GMP_LIMB_MAX << -m;
1795
0
  x += l;
1796
0
      }
1797
0
    }
1798
1799
0
  if (u->_mp_size < 0)
1800
0
    x = -x;
1801
1802
0
  return x;
1803
0
}
1804
1805
int
1806
mpz_cmpabs_d (const mpz_t x, double d)
1807
0
{
1808
0
  mp_size_t xn;
1809
0
  double B, Bi;
1810
0
  mp_size_t i;
1811
1812
0
  xn = x->_mp_size;
1813
0
  d = GMP_ABS (d);
1814
1815
0
  if (xn != 0)
1816
0
    {
1817
0
      xn = GMP_ABS (xn);
1818
1819
0
      B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1);
1820
0
      Bi = 1.0 / B;
1821
1822
      /* Scale d so it can be compared with the top limb. */
1823
0
      for (i = 1; i < xn; i++)
1824
0
  d *= Bi;
1825
1826
0
      if (d >= B)
1827
0
  return -1;
1828
1829
      /* Compare floor(d) to top limb, subtract and cancel when equal. */
1830
0
      for (i = xn; i-- > 0;)
1831
0
  {
1832
0
    mp_limb_t f, xl;
1833
1834
0
    f = (mp_limb_t) d;
1835
0
    xl = x->_mp_d[i];
1836
0
    if (xl > f)
1837
0
      return 1;
1838
0
    else if (xl < f)
1839
0
      return -1;
1840
0
    d = B * (d - f);
1841
0
  }
1842
0
    }
1843
0
  return - (d > 0.0);
1844
0
}
1845
1846
int
1847
mpz_cmp_d (const mpz_t x, double d)
1848
0
{
1849
0
  if (x->_mp_size < 0)
1850
0
    {
1851
0
      if (d >= 0.0)
1852
0
  return -1;
1853
0
      else
1854
0
  return -mpz_cmpabs_d (x, d);
1855
0
    }
1856
0
  else
1857
0
    {
1858
0
      if (d < 0.0)
1859
0
  return 1;
1860
0
      else
1861
0
  return mpz_cmpabs_d (x, d);
1862
0
    }
1863
0
}
1864
1865

1866
/* MPZ comparisons and the like. */
1867
int
1868
mpz_sgn (const mpz_t u)
1869
400
{
1870
400
  return GMP_CMP (u->_mp_size, 0);
1871
400
}
1872
1873
int
1874
mpz_cmp_si (const mpz_t u, long v)
1875
36
{
1876
36
  mp_size_t usize = u->_mp_size;
1877
1878
36
  if (v >= 0)
1879
25
    return mpz_cmp_ui (u, v);
1880
11
  else if (usize >= 0)
1881
11
    return 1;
1882
0
  else
1883
0
    return - mpz_cmpabs_ui (u, GMP_NEG_CAST (unsigned long int, v));
1884
36
}
1885
1886
int
1887
mpz_cmp_ui (const mpz_t u, unsigned long v)
1888
570
{
1889
570
  mp_size_t usize = u->_mp_size;
1890
1891
570
  if (usize < 0)
1892
0
    return -1;
1893
570
  else
1894
570
    return mpz_cmpabs_ui (u, v);
1895
570
}
1896
1897
int
1898
mpz_cmp (const mpz_t a, const mpz_t b)
1899
93.3k
{
1900
93.3k
  mp_size_t asize = a->_mp_size;
1901
93.3k
  mp_size_t bsize = b->_mp_size;
1902
1903
93.3k
  if (asize != bsize)
1904
3.98k
    return (asize < bsize) ? -1 : 1;
1905
89.3k
  else if (asize >= 0)
1906
89.3k
    return mpn_cmp (a->_mp_d, b->_mp_d, asize);
1907
0
  else
1908
0
    return mpn_cmp (b->_mp_d, a->_mp_d, -asize);
1909
93.3k
}
1910
1911
int
1912
mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1913
796
{
1914
796
  mp_size_t un = GMP_ABS (u->_mp_size);
1915
1916
796
  if (! mpn_absfits_ulong_p (u->_mp_d, un))
1917
658
    return 1;
1918
138
  else
1919
138
    {
1920
138
      unsigned long uu = mpz_get_ui (u);
1921
138
      return GMP_CMP(uu, v);
1922
138
    }
1923
796
}
1924
1925
int
1926
mpz_cmpabs (const mpz_t u, const mpz_t v)
1927
707
{
1928
707
  return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size),
1929
707
       v->_mp_d, GMP_ABS (v->_mp_size));
1930
707
}
1931
1932
void
1933
mpz_abs (mpz_t r, const mpz_t u)
1934
606
{
1935
606
  mpz_set (r, u);
1936
606
  r->_mp_size = GMP_ABS (r->_mp_size);
1937
606
}
1938
1939
void
1940
mpz_neg (mpz_t r, const mpz_t u)
1941
193k
{
1942
193k
  mpz_set (r, u);
1943
193k
  r->_mp_size = -r->_mp_size;
1944
193k
}
1945
1946
void
1947
mpz_swap (mpz_t u, mpz_t v)
1948
64.9k
{
1949
64.9k
  MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc);
1950
64.9k
  MPN_PTR_SWAP (u->_mp_d, u->_mp_size, v->_mp_d, v->_mp_size);
1951
64.9k
}
1952
1953

1954
/* MPZ addition and subtraction */
1955
1956
1957
void
1958
mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1959
23
{
1960
23
  mpz_t bb;
1961
23
  mpz_init_set_ui (bb, b);
1962
23
  mpz_add (r, a, bb);
1963
23
  mpz_clear (bb);
1964
23
}
1965
1966
void
1967
mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1968
2
{
1969
2
  mpz_ui_sub (r, b, a);
1970
2
  mpz_neg (r, r);
1971
2
}
1972
1973
void
1974
mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b)
1975
5
{
1976
5
  mpz_neg (r, b);
1977
5
  mpz_add_ui (r, r, a);
1978
5
}
1979
1980
static mp_size_t
1981
mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
1982
377k
{
1983
377k
  mp_size_t an = GMP_ABS (a->_mp_size);
1984
377k
  mp_size_t bn = GMP_ABS (b->_mp_size);
1985
377k
  mp_ptr rp;
1986
377k
  mp_limb_t cy;
1987
1988
377k
  if (an < bn)
1989
2.61k
    {
1990
2.61k
      MPZ_SRCPTR_SWAP (a, b);
1991
2.61k
      MP_SIZE_T_SWAP (an, bn);
1992
2.61k
    }
1993
1994
377k
  rp = MPZ_REALLOC (r, an + 1);
1995
377k
  cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
1996
1997
377k
  rp[an] = cy;
1998
1999
377k
  return an + cy;
2000
377k
}
2001
2002
static mp_size_t
2003
mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
2004
99.5k
{
2005
99.5k
  mp_size_t an = GMP_ABS (a->_mp_size);
2006
99.5k
  mp_size_t bn = GMP_ABS (b->_mp_size);
2007
99.5k
  int cmp;
2008
99.5k
  mp_ptr rp;
2009
2010
99.5k
  cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
2011
99.5k
  if (cmp > 0)
2012
96.3k
    {
2013
96.3k
      rp = MPZ_REALLOC (r, an);
2014
96.3k
      gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn));
2015
96.3k
      return mpn_normalized_size (rp, an);
2016
96.3k
    }
2017
3.15k
  else if (cmp < 0)
2018
3.15k
    {
2019
3.15k
      rp = MPZ_REALLOC (r, bn);
2020
3.15k
      gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an));
2021
3.15k
      return -mpn_normalized_size (rp, bn);
2022
3.15k
    }
2023
1
  else
2024
1
    return 0;
2025
99.5k
}
2026
2027
void
2028
mpz_add (mpz_t r, const mpz_t a, const mpz_t b)
2029
285k
{
2030
285k
  mp_size_t rn;
2031
2032
285k
  if ( (a->_mp_size ^ b->_mp_size) >= 0)
2033
281k
    rn = mpz_abs_add (r, a, b);
2034
3.67k
  else
2035
3.67k
    rn = mpz_abs_sub (r, a, b);
2036
2037
285k
  r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
2038
285k
}
2039
2040
void
2041
mpz_sub (mpz_t r, const mpz_t a, const mpz_t b)
2042
191k
{
2043
191k
  mp_size_t rn;
2044
2045
191k
  if ( (a->_mp_size ^ b->_mp_size) >= 0)
2046
95.8k
    rn = mpz_abs_sub (r, a, b);
2047
95.5k
  else
2048
95.5k
    rn = mpz_abs_add (r, a, b);
2049
2050
191k
  r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
2051
191k
}
2052
2053

2054
/* MPZ multiplication */
2055
void
2056
mpz_mul_si (mpz_t r, const mpz_t u, long int v)
2057
3.12k
{
2058
3.12k
  if (v < 0)
2059
1.56k
    {
2060
1.56k
      mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v));
2061
1.56k
      mpz_neg (r, r);
2062
1.56k
    }
2063
1.56k
  else
2064
1.56k
    mpz_mul_ui (r, u, v);
2065
3.12k
}
2066
2067
void
2068
mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2069
3.19k
{
2070
3.19k
  mpz_t vv;
2071
3.19k
  mpz_init_set_ui (vv, v);
2072
3.19k
  mpz_mul (r, u, vv);
2073
3.19k
  mpz_clear (vv);
2074
3.19k
  return;
2075
3.19k
}
2076
2077
void
2078
mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
2079
50.4k
{
2080
50.4k
  int sign;
2081
50.4k
  mp_size_t un, vn, rn;
2082
50.4k
  mpz_t t;
2083
50.4k
  mp_ptr tp;
2084
2085
50.4k
  un = u->_mp_size;
2086
50.4k
  vn = v->_mp_size;
2087
2088
50.4k
  if (un == 0 || vn == 0)
2089
3
    {
2090
3
      r->_mp_size = 0;
2091
3
      return;
2092
3
    }
2093
2094
50.4k
  sign = (un ^ vn) < 0;
2095
2096
50.4k
  un = GMP_ABS (un);
2097
50.4k
  vn = GMP_ABS (vn);
2098
2099
50.4k
  mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
2100
2101
50.4k
  tp = t->_mp_d;
2102
50.4k
  if (un >= vn)
2103
50.0k
    mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn);
2104
401
  else
2105
401
    mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un);
2106
2107
50.4k
  rn = un + vn;
2108
50.4k
  rn -= tp[rn-1] == 0;
2109
2110
50.4k
  t->_mp_size = sign ? - rn : rn;
2111
50.4k
  mpz_swap (r, t);
2112
50.4k
  mpz_clear (t);
2113
50.4k
}
2114
2115
void
2116
mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
2117
185k
{
2118
185k
  mp_size_t un, rn;
2119
185k
  mp_size_t limbs;
2120
185k
  unsigned shift;
2121
185k
  mp_ptr rp;
2122
2123
185k
  un = GMP_ABS (u->_mp_size);
2124
185k
  if (un == 0)
2125
2
    {
2126
2
      r->_mp_size = 0;
2127
2
      return;
2128
2
    }
2129
2130
185k
  limbs = bits / GMP_LIMB_BITS;
2131
185k
  shift = bits % GMP_LIMB_BITS;
2132
2133
185k
  rn = un + limbs + (shift > 0);
2134
185k
  rp = MPZ_REALLOC (r, rn);
2135
185k
  if (shift > 0)
2136
185k
    {
2137
185k
      mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
2138
185k
      rp[rn-1] = cy;
2139
185k
      rn -= (cy == 0);
2140
185k
    }
2141
156
  else
2142
156
    mpn_copyd (rp + limbs, u->_mp_d, un);
2143
2144
185k
  mpn_zero (rp, limbs);
2145
2146
185k
  r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
2147
185k
}
2148
2149
void
2150
mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2151
0
{
2152
0
  mpz_t t;
2153
0
  mpz_init_set_ui (t, v);
2154
0
  mpz_mul (t, u, t);
2155
0
  mpz_add (r, r, t);
2156
0
  mpz_clear (t);
2157
0
}
2158
2159
void
2160
mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2161
3.11k
{
2162
3.11k
  mpz_t t;
2163
3.11k
  mpz_init_set_ui (t, v);
2164
3.11k
  mpz_mul (t, u, t);
2165
3.11k
  mpz_sub (r, r, t);
2166
3.11k
  mpz_clear (t);
2167
3.11k
}
2168
2169
void
2170
mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v)
2171
2
{
2172
2
  mpz_t t;
2173
2
  mpz_init (t);
2174
2
  mpz_mul (t, u, v);
2175
2
  mpz_add (r, r, t);
2176
2
  mpz_clear (t);
2177
2
}
2178
2179
void
2180
mpz_submul (mpz_t r, const mpz_t u, const mpz_t v)
2181
2
{
2182
2
  mpz_t t;
2183
2
  mpz_init (t);
2184
2
  mpz_mul (t, u, v);
2185
2
  mpz_sub (r, r, t);
2186
2
  mpz_clear (t);
2187
2
}
2188
2189

2190
/* MPZ division */
2191
enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
2192
2193
/* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
2194
static int
2195
mpz_div_qr (mpz_t q, mpz_t r,
2196
      const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode)
2197
11.9k
{
2198
11.9k
  mp_size_t ns, ds, nn, dn, qs;
2199
11.9k
  ns = n->_mp_size;
2200
11.9k
  ds = d->_mp_size;
2201
2202
11.9k
  if (ds == 0)
2203
0
    gmp_die("mpz_div_qr: Divide by zero.");
2204
2205
11.9k
  if (ns == 0)
2206
9
    {
2207
9
      if (q)
2208
0
  q->_mp_size = 0;
2209
9
      if (r)
2210
9
  r->_mp_size = 0;
2211
9
      return 0;
2212
9
    }
2213
2214
11.9k
  nn = GMP_ABS (ns);
2215
11.9k
  dn = GMP_ABS (ds);
2216
2217
11.9k
  qs = ds ^ ns;
2218
2219
11.9k
  if (nn < dn)
2220
551
    {
2221
551
      if (mode == GMP_DIV_CEIL && qs >= 0)
2222
0
  {
2223
    /* q = 1, r = n - d */
2224
0
    if (r)
2225
0
      mpz_sub (r, n, d);
2226
0
    if (q)
2227
0
      mpz_set_ui (q, 1);
2228
0
  }
2229
551
      else if (mode == GMP_DIV_FLOOR && qs < 0)
2230
91
  {
2231
    /* q = -1, r = n + d */
2232
91
    if (r)
2233
91
      mpz_add (r, n, d);
2234
91
    if (q)
2235
0
      mpz_set_si (q, -1);
2236
91
  }
2237
460
      else
2238
460
  {
2239
    /* q = 0, r = d */
2240
460
    if (r)
2241
457
      mpz_set (r, n);
2242
460
    if (q)
2243
0
      q->_mp_size = 0;
2244
460
  }
2245
551
      return 1;
2246
551
    }
2247
11.4k
  else
2248
11.4k
    {
2249
11.4k
      mp_ptr np, qp;
2250
11.4k
      mp_size_t qn, rn;
2251
11.4k
      mpz_t tq, tr;
2252
2253
11.4k
      mpz_init_set (tr, n);
2254
11.4k
      np = tr->_mp_d;
2255
2256
11.4k
      qn = nn - dn + 1;
2257
2258
11.4k
      if (q)
2259
926
  {
2260
926
    mpz_init2 (tq, qn * GMP_LIMB_BITS);
2261
926
    qp = tq->_mp_d;
2262
926
  }
2263
10.4k
      else
2264
10.4k
  qp = NULL;
2265
2266
11.4k
      mpn_div_qr (qp, np, nn, d->_mp_d, dn);
2267
2268
11.4k
      if (qp)
2269
926
  {
2270
926
    qn -= (qp[qn-1] == 0);
2271
2272
926
    tq->_mp_size = qs < 0 ? -qn : qn;
2273
926
  }
2274
11.4k
      rn = mpn_normalized_size (np, dn);
2275
11.4k
      tr->_mp_size = ns < 0 ? - rn : rn;
2276
2277
11.4k
      if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0)
2278
125
  {
2279
125
    if (q)
2280
0
      mpz_sub_ui (tq, tq, 1);
2281
125
    if (r)
2282
125
      mpz_add (tr, tr, d);
2283
125
  }
2284
11.2k
      else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0)
2285
0
  {
2286
0
    if (q)
2287
0
      mpz_add_ui (tq, tq, 1);
2288
0
    if (r)
2289
0
      mpz_sub (tr, tr, d);
2290
0
  }
2291
2292
11.4k
      if (q)
2293
926
  {
2294
926
    mpz_swap (tq, q);
2295
926
    mpz_clear (tq);
2296
926
  }
2297
11.4k
      if (r)
2298
10.6k
  mpz_swap (tr, r);
2299
2300
11.4k
      mpz_clear (tr);
2301
2302
11.4k
      return rn != 0;
2303
11.4k
    }
2304
11.9k
}
2305
2306
void
2307
mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2308
0
{
2309
0
  mpz_div_qr (q, r, n, d, GMP_DIV_CEIL);
2310
0
}
2311
2312
void
2313
mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2314
0
{
2315
0
  mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR);
2316
0
}
2317
2318
void
2319
mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2320
75
{
2321
75
  mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC);
2322
75
}
2323
2324
void
2325
mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2326
0
{
2327
0
  mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL);
2328
0
}
2329
2330
void
2331
mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2332
0
{
2333
0
  mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR);
2334
0
}
2335
2336
void
2337
mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2338
630
{
2339
630
  mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC);
2340
630
}
2341
2342
void
2343
mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2344
0
{
2345
0
  mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2346
0
}
2347
2348
void
2349
mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2350
0
{
2351
0
  mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2352
0
}
2353
2354
void
2355
mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2356
7.79k
{
2357
7.79k
  mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC);
2358
7.79k
}
2359
2360
void
2361
mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
2362
3.22k
{
2363
3.22k
  mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL);
2364
3.22k
}
2365
2366
static void
2367
mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
2368
    enum mpz_div_round_mode mode)
2369
475k
{
2370
475k
  mp_size_t un, qn;
2371
475k
  mp_size_t limb_cnt;
2372
475k
  mp_ptr qp;
2373
475k
  int adjust;
2374
2375
475k
  un = u->_mp_size;
2376
475k
  if (un == 0)
2377
112
    {
2378
112
      q->_mp_size = 0;
2379
112
      return;
2380
112
    }
2381
475k
  limb_cnt = bit_index / GMP_LIMB_BITS;
2382
475k
  qn = GMP_ABS (un) - limb_cnt;
2383
475k
  bit_index %= GMP_LIMB_BITS;
2384
2385
475k
  if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */
2386
    /* Note: Below, the final indexing at limb_cnt is valid because at
2387
       that point we have qn > 0. */
2388
0
    adjust = (qn <= 0
2389
0
        || !mpn_zero_p (u->_mp_d, limb_cnt)
2390
0
        || (u->_mp_d[limb_cnt]
2391
0
      & (((mp_limb_t) 1 << bit_index) - 1)));
2392
475k
  else
2393
475k
    adjust = 0;
2394
2395
475k
  if (qn <= 0)
2396
0
    qn = 0;
2397
475k
  else
2398
475k
    {
2399
475k
      qp = MPZ_REALLOC (q, qn);
2400
2401
475k
      if (bit_index != 0)
2402
475k
  {
2403
475k
    mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
2404
475k
    qn -= qp[qn - 1] == 0;
2405
475k
  }
2406
216
      else
2407
216
  {
2408
216
    mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
2409
216
  }
2410
475k
    }
2411
2412
475k
  q->_mp_size = qn;
2413
2414
475k
  if (adjust)
2415
0
    mpz_add_ui (q, q, 1);
2416
475k
  if (un < 0)
2417
191k
    mpz_neg (q, q);
2418
475k
}
2419
2420
static void
2421
mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
2422
    enum mpz_div_round_mode mode)
2423
0
{
2424
0
  mp_size_t us, un, rn;
2425
0
  mp_ptr rp;
2426
0
  mp_limb_t mask;
2427
2428
0
  us = u->_mp_size;
2429
0
  if (us == 0 || bit_index == 0)
2430
0
    {
2431
0
      r->_mp_size = 0;
2432
0
      return;
2433
0
    }
2434
0
  rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
2435
0
  assert (rn > 0);
2436
2437
0
  rp = MPZ_REALLOC (r, rn);
2438
0
  un = GMP_ABS (us);
2439
2440
0
  mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);
2441
2442
0
  if (rn > un)
2443
0
    {
2444
      /* Quotient (with truncation) is zero, and remainder is
2445
   non-zero */
2446
0
      if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2447
0
  {
2448
    /* Have to negate and sign extend. */
2449
0
    mp_size_t i;
2450
2451
0
    gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un));
2452
0
    for (i = un; i < rn - 1; i++)
2453
0
      rp[i] = GMP_LIMB_MAX;
2454
2455
0
    rp[rn-1] = mask;
2456
0
    us = -us;
2457
0
  }
2458
0
      else
2459
0
  {
2460
    /* Just copy */
2461
0
    if (r != u)
2462
0
      mpn_copyi (rp, u->_mp_d, un);
2463
2464
0
    rn = un;
2465
0
  }
2466
0
    }
2467
0
  else
2468
0
    {
2469
0
      if (r != u)
2470
0
  mpn_copyi (rp, u->_mp_d, rn - 1);
2471
2472
0
      rp[rn-1] = u->_mp_d[rn-1] & mask;
2473
2474
0
      if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2475
0
  {
2476
    /* If r != 0, compute 2^{bit_count} - r. */
2477
0
    mpn_neg (rp, rp, rn);
2478
2479
0
    rp[rn-1] &= mask;
2480
2481
    /* us is not used for anything else, so we can modify it
2482
       here to indicate flipped sign. */
2483
0
    us = -us;
2484
0
  }
2485
0
    }
2486
0
  rn = mpn_normalized_size (rp, rn);
2487
0
  r->_mp_size = us < 0 ? -rn : rn;
2488
0
}
2489
2490
void
2491
mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2492
0
{
2493
0
  mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL);
2494
0
}
2495
2496
void
2497
mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2498
0
{
2499
0
  mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR);
2500
0
}
2501
2502
void
2503
mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2504
475k
{
2505
475k
  mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC);
2506
475k
}
2507
2508
void
2509
mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2510
0
{
2511
0
  mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL);
2512
0
}
2513
2514
void
2515
mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2516
0
{
2517
0
  mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR);
2518
0
}
2519
2520
void
2521
mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2522
0
{
2523
0
  mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC);
2524
0
}
2525
2526
void
2527
mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d)
2528
156
{
2529
156
  gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC));
2530
156
}
2531
2532
int
2533
mpz_divisible_p (const mpz_t n, const mpz_t d)
2534
8
{
2535
8
  return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2536
8
}
2537
2538
int
2539
mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m)
2540
0
{
2541
0
  mpz_t t;
2542
0
  int res;
2543
2544
  /* a == b (mod 0) iff a == b */
2545
0
  if (mpz_sgn (m) == 0)
2546
0
    return (mpz_cmp (a, b) == 0);
2547
2548
0
  mpz_init (t);
2549
0
  mpz_sub (t, a, b);
2550
0
  res = mpz_divisible_p (t, m);
2551
0
  mpz_clear (t);
2552
2553
0
  return res;
2554
0
}
2555
2556
static unsigned long
2557
mpz_div_qr_ui (mpz_t q, mpz_t r,
2558
         const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
2559
100
{
2560
100
  unsigned long ret;
2561
100
  mpz_t rr, dd;
2562
2563
100
  mpz_init (rr);
2564
100
  mpz_init_set_ui (dd, d);
2565
100
  mpz_div_qr (q, rr, n, dd, mode);
2566
100
  mpz_clear (dd);
2567
100
  ret = mpz_get_ui (rr);
2568
2569
100
  if (r)
2570
0
    mpz_swap (r, rr);
2571
100
  mpz_clear (rr);
2572
2573
100
  return ret;
2574
100
}
2575
2576
unsigned long
2577
mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2578
0
{
2579
0
  return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL);
2580
0
}
2581
2582
unsigned long
2583
mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2584
0
{
2585
0
  return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR);
2586
0
}
2587
2588
unsigned long
2589
mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2590
0
{
2591
0
  return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC);
2592
0
}
2593
2594
unsigned long
2595
mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2596
0
{
2597
0
  return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL);
2598
0
}
2599
2600
unsigned long
2601
mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2602
0
{
2603
0
  return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR);
2604
0
}
2605
2606
unsigned long
2607
mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2608
65
{
2609
65
  return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC);
2610
65
}
2611
2612
unsigned long
2613
mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2614
0
{
2615
0
  return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL);
2616
0
}
2617
unsigned long
2618
mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2619
0
{
2620
0
  return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2621
0
}
2622
unsigned long
2623
mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2624
0
{
2625
0
  return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC);
2626
0
}
2627
2628
unsigned long
2629
mpz_cdiv_ui (const mpz_t n, unsigned long d)
2630
0
{
2631
0
  return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL);
2632
0
}
2633
2634
unsigned long
2635
mpz_fdiv_ui (const mpz_t n, unsigned long d)
2636
0
{
2637
0
  return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR);
2638
0
}
2639
2640
unsigned long
2641
mpz_tdiv_ui (const mpz_t n, unsigned long d)
2642
35
{
2643
35
  return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC);
2644
35
}
2645
2646
unsigned long
2647
mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d)
2648
0
{
2649
0
  return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2650
0
}
2651
2652
void
2653
mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d)
2654
0
{
2655
0
  gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC));
2656
0
}
2657
2658
int
2659
mpz_divisible_ui_p (const mpz_t n, unsigned long d)
2660
0
{
2661
0
  return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2662
0
}
2663
2664

2665
/* GCD */
2666
static mp_limb_t
2667
mpn_gcd_11 (mp_limb_t u, mp_limb_t v)
2668
58
{
2669
58
  unsigned shift;
2670
2671
58
  assert ( (u | v) > 0);
2672
2673
58
  if (u == 0)
2674
0
    return v;
2675
58
  else if (v == 0)
2676
0
    return u;
2677
2678
58
  gmp_ctz (shift, u | v);
2679
2680
58
  u >>= shift;
2681
58
  v >>= shift;
2682
2683
58
  if ( (u & 1) == 0)
2684
33
    MP_LIMB_T_SWAP (u, v);
2685
2686
110
  while ( (v & 1) == 0)
2687
52
    v >>= 1;
2688
2689
1.50k
  while (u != v)
2690
1.44k
    {
2691
1.44k
      if (u > v)
2692
754
  {
2693
754
    u -= v;
2694
754
    do
2695
1.50k
      u >>= 1;
2696
1.50k
    while ( (u & 1) == 0);
2697
754
  }
2698
690
      else
2699
690
  {
2700
690
    v -= u;
2701
690
    do
2702
1.39k
      v >>= 1;
2703
1.39k
    while ( (v & 1) == 0);
2704
690
  }
2705
1.44k
    }
2706
58
  return u << shift;
2707
58
}
2708
2709
mp_size_t
2710
mpn_gcd (mp_ptr rp, mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t vn)
2711
58
{
2712
58
  assert (un >= vn);
2713
58
  assert (vn > 0);
2714
58
  assert (!GMP_MPN_OVERLAP_P (up, un, vp, vn));
2715
58
  assert (vp[vn-1] > 0);
2716
58
  assert ((up[0] | vp[0]) & 1);
2717
2718
58
  if (un > vn)
2719
53
    mpn_div_qr (NULL, up, un, vp, vn);
2720
2721
58
  un = mpn_normalized_size (up, vn);
2722
58
  if (un == 0)
2723
0
    {
2724
0
      mpn_copyi (rp, vp, vn);
2725
0
      return vn;
2726
0
    }
2727
2728
58
  if (!(vp[0] & 1))
2729
0
    MPN_PTR_SWAP (up, un, vp, vn);
2730
2731
7.83k
  while (un > 1 || vn > 1)
2732
7.77k
    {
2733
7.77k
      int shift;
2734
7.77k
      assert (vp[0] & 1);
2735
2736
7.77k
      while (up[0] == 0)
2737
0
  {
2738
0
    up++;
2739
0
    un--;
2740
0
  }
2741
7.77k
      gmp_ctz (shift, up[0]);
2742
7.77k
      if (shift > 0)
2743
7.77k
  {
2744
7.77k
    gmp_assert_nocarry (mpn_rshift(up, up, un, shift));
2745
7.77k
    un -= (up[un-1] == 0);
2746
7.77k
  }
2747
2748
7.77k
      if (un < vn)
2749
175
  MPN_PTR_SWAP (up, un, vp, vn);
2750
7.60k
      else if (un == vn)
2751
7.46k
  {
2752
7.46k
    int c = mpn_cmp (up, vp, un);
2753
7.46k
    if (c == 0)
2754
0
      {
2755
0
        mpn_copyi (rp, up, un);
2756
0
        return un;
2757
0
      }
2758
7.46k
    else if (c < 0)
2759
4.10k
      MP_PTR_SWAP (up, vp);
2760
7.46k
  }
2761
2762
7.77k
      gmp_assert_nocarry (mpn_sub (up, up, un, vp, vn));
2763
7.77k
      un = mpn_normalized_size (up, un);
2764
7.77k
    }
2765
58
  rp[0] = mpn_gcd_11 (up[0], vp[0]);
2766
58
  return 1;
2767
58
}
2768
2769
unsigned long
2770
mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
2771
48
{
2772
48
  mpz_t t;
2773
48
  mpz_init_set_ui(t, v);
2774
48
  mpz_gcd (t, u, t);
2775
48
  if (v > 0)
2776
48
    v = mpz_get_ui (t);
2777
2778
48
  if (g)
2779
1
    mpz_swap (t, g);
2780
2781
48
  mpz_clear (t);
2782
2783
48
  return v;
2784
48
}
2785
2786
static mp_bitcnt_t
2787
mpz_make_odd (mpz_t r)
2788
92.8k
{
2789
92.8k
  mp_bitcnt_t shift;
2790
2791
92.8k
  assert (r->_mp_size > 0);
2792
  /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
2793
92.8k
  shift = mpn_scan1 (r->_mp_d, 0);
2794
92.8k
  mpz_tdiv_q_2exp (r, r, shift);
2795
2796
92.8k
  return shift;
2797
92.8k
}
2798
2799
void
2800
mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
2801
62
{
2802
62
  mpz_t tu, tv;
2803
62
  mp_bitcnt_t uz, vz, gz;
2804
2805
62
  if (u->_mp_size == 0)
2806
0
    {
2807
0
      mpz_abs (g, v);
2808
0
      return;
2809
0
    }
2810
62
  if (v->_mp_size == 0)
2811
4
    {
2812
4
      mpz_abs (g, u);
2813
4
      return;
2814
4
    }
2815
2816
58
  mpz_init (tu);
2817
58
  mpz_init (tv);
2818
2819
58
  mpz_abs (tu, u);
2820
58
  uz = mpz_make_odd (tu);
2821
58
  mpz_abs (tv, v);
2822
58
  vz = mpz_make_odd (tv);
2823
58
  gz = GMP_MIN (uz, vz);
2824
2825
58
  if (tu->_mp_size < tv->_mp_size)
2826
6
    mpz_swap (tu, tv);
2827
2828
58
  tu->_mp_size = mpn_gcd (tu->_mp_d, tu->_mp_d, tu->_mp_size, tv->_mp_d, tv->_mp_size);
2829
58
  mpz_mul_2exp (g, tu, gz);
2830
2831
58
  mpz_clear (tu);
2832
58
  mpz_clear (tv);
2833
58
}
2834
2835
void
2836
mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
2837
77
{
2838
77
  mpz_t tu, tv, s0, s1, t0, t1;
2839
77
  mp_bitcnt_t uz, vz, gz;
2840
77
  mp_bitcnt_t power;
2841
77
  int cmp;
2842
2843
77
  if (u->_mp_size == 0)
2844
0
    {
2845
      /* g = 0 u + sgn(v) v */
2846
0
      signed long sign = mpz_sgn (v);
2847
0
      mpz_abs (g, v);
2848
0
      if (s)
2849
0
  s->_mp_size = 0;
2850
0
      if (t)
2851
0
  mpz_set_si (t, sign);
2852
0
      return;
2853
0
    }
2854
2855
77
  if (v->_mp_size == 0)
2856
2
    {
2857
      /* g = sgn(u) u + 0 v */
2858
2
      signed long sign = mpz_sgn (u);
2859
2
      mpz_abs (g, u);
2860
2
      if (s)
2861
2
  mpz_set_si (s, sign);
2862
2
      if (t)
2863
2
  t->_mp_size = 0;
2864
2
      return;
2865
2
    }
2866
2867
75
  mpz_init (tu);
2868
75
  mpz_init (tv);
2869
75
  mpz_init (s0);
2870
75
  mpz_init (s1);
2871
75
  mpz_init (t0);
2872
75
  mpz_init (t1);
2873
2874
75
  mpz_abs (tu, u);
2875
75
  uz = mpz_make_odd (tu);
2876
75
  mpz_abs (tv, v);
2877
75
  vz = mpz_make_odd (tv);
2878
75
  gz = GMP_MIN (uz, vz);
2879
2880
75
  uz -= gz;
2881
75
  vz -= gz;
2882
2883
  /* Cofactors corresponding to odd gcd. gz handled later. */
2884
75
  if (tu->_mp_size < tv->_mp_size)
2885
45
    {
2886
45
      mpz_swap (tu, tv);
2887
45
      MPZ_SRCPTR_SWAP (u, v);
2888
45
      MPZ_PTR_SWAP (s, t);
2889
45
      MP_BITCNT_T_SWAP (uz, vz);
2890
45
    }
2891
2892
  /* Maintain
2893
   *
2894
   * u = t0 tu + t1 tv
2895
   * v = s0 tu + s1 tv
2896
   *
2897
   * where u and v denote the inputs with common factors of two
2898
   * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
2899
   *
2900
   * 2^p tu =  s1 u - t1 v
2901
   * 2^p tv = -s0 u + t0 v
2902
   */
2903
2904
  /* After initial division, tu = q tv + tu', we have
2905
   *
2906
   * u = 2^uz (tu' + q tv)
2907
   * v = 2^vz tv
2908
   *
2909
   * or
2910
   *
2911
   * t0 = 2^uz, t1 = 2^uz q
2912
   * s0 = 0,    s1 = 2^vz
2913
   */
2914
2915
75
  mpz_tdiv_qr (t1, tu, tu, tv);
2916
75
  mpz_mul_2exp (t1, t1, uz);
2917
2918
75
  mpz_setbit (s1, vz);
2919
75
  power = uz + vz;
2920
2921
75
  if (tu->_mp_size > 0)
2922
71
    {
2923
71
      mp_bitcnt_t shift;
2924
71
      shift = mpz_make_odd (tu);
2925
71
      mpz_setbit (t0, uz + shift);
2926
71
      power += shift;
2927
2928
71
      for (;;)
2929
92.6k
  {
2930
92.6k
    int c;
2931
92.6k
    c = mpz_cmp (tu, tv);
2932
92.6k
    if (c == 0)
2933
71
      break;
2934
2935
92.5k
    if (c < 0)
2936
46.3k
      {
2937
        /* tv = tv' + tu
2938
         *
2939
         * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
2940
         * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */
2941
2942
46.3k
        mpz_sub (tv, tv, tu);
2943
46.3k
        mpz_add (t0, t0, t1);
2944
46.3k
        mpz_add (s0, s0, s1);
2945
2946
46.3k
        shift = mpz_make_odd (tv);
2947
46.3k
        mpz_mul_2exp (t1, t1, shift);
2948
46.3k
        mpz_mul_2exp (s1, s1, shift);
2949
46.3k
      }
2950
46.1k
    else
2951
46.1k
      {
2952
46.1k
        mpz_sub (tu, tu, tv);
2953
46.1k
        mpz_add (t1, t0, t1);
2954
46.1k
        mpz_add (s1, s0, s1);
2955
2956
46.1k
        shift = mpz_make_odd (tu);
2957
46.1k
        mpz_mul_2exp (t0, t0, shift);
2958
46.1k
        mpz_mul_2exp (s0, s0, shift);
2959
46.1k
      }
2960
92.5k
    power += shift;
2961
92.5k
  }
2962
71
    }
2963
4
  else
2964
4
    mpz_setbit (t0, uz);
2965
2966
  /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
2967
     cofactors. */
2968
2969
75
  mpz_mul_2exp (tv, tv, gz);
2970
75
  mpz_neg (s0, s0);
2971
2972
  /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
2973
     adjust cofactors, we need u / g and v / g */
2974
2975
75
  mpz_divexact (s1, v, tv);
2976
75
  mpz_abs (s1, s1);
2977
75
  mpz_divexact (t1, u, tv);
2978
75
  mpz_abs (t1, t1);
2979
2980
190k
  while (power-- > 0)
2981
190k
    {
2982
      /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
2983
190k
      if (mpz_odd_p (s0) || mpz_odd_p (t0))
2984
95.5k
  {
2985
95.5k
    mpz_sub (s0, s0, s1);
2986
95.5k
    mpz_add (t0, t0, t1);
2987
95.5k
  }
2988
190k
      assert (mpz_even_p (t0) && mpz_even_p (s0));
2989
190k
      mpz_tdiv_q_2exp (s0, s0, 1);
2990
190k
      mpz_tdiv_q_2exp (t0, t0, 1);
2991
190k
    }
2992
2993
  /* Choose small cofactors (they should generally satify
2994
2995
       |s| < |u| / 2g and |t| < |v| / 2g,
2996
2997
     with some documented exceptions). Always choose the smallest s,
2998
     if there are two choices for s with same absolute value, choose
2999
     the one with smallest corresponding t (this asymmetric condition
3000
     is needed to prefer s = 0, |t| = 1 when g = |a| = |b|). */
3001
75
  mpz_add (s1, s0, s1);
3002
75
  mpz_sub (t1, t0, t1);
3003
75
  cmp = mpz_cmpabs (s0, s1);
3004
75
  if (cmp > 0 || (cmp == 0 && mpz_cmpabs (t0, t1) > 0))
3005
38
    {
3006
38
      mpz_swap (s0, s1);
3007
38
      mpz_swap (t0, t1);
3008
38
    }
3009
75
  if (u->_mp_size < 0)
3010
0
    mpz_neg (s0, s0);
3011
75
  if (v->_mp_size < 0)
3012
0
    mpz_neg (t0, t0);
3013
3014
75
  mpz_swap (g, tv);
3015
75
  if (s)
3016
34
    mpz_swap (s, s0);
3017
75
  if (t)
3018
55
    mpz_swap (t, t0);
3019
3020
75
  mpz_clear (tu);
3021
75
  mpz_clear (tv);
3022
75
  mpz_clear (s0);
3023
75
  mpz_clear (s1);
3024
75
  mpz_clear (t0);
3025
75
  mpz_clear (t1);
3026
75
}
3027
3028
void
3029
mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v)
3030
6
{
3031
6
  mpz_t g;
3032
3033
6
  if (u->_mp_size == 0 || v->_mp_size == 0)
3034
0
    {
3035
0
      r->_mp_size = 0;
3036
0
      return;
3037
0
    }
3038
3039
6
  mpz_init (g);
3040
3041
6
  mpz_gcd (g, u, v);
3042
6
  mpz_divexact (g, u, g);
3043
6
  mpz_mul (r, g, v);
3044
3045
6
  mpz_clear (g);
3046
6
  mpz_abs (r, r);
3047
6
}
3048
3049
void
3050
mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v)
3051
1
{
3052
1
  if (v == 0 || u->_mp_size == 0)
3053
0
    {
3054
0
      r->_mp_size = 0;
3055
0
      return;
3056
0
    }
3057
3058
1
  v /= mpz_gcd_ui (NULL, u, v);
3059
1
  mpz_mul_ui (r, u, v);
3060
3061
1
  mpz_abs (r, r);
3062
1
}
3063
3064
int
3065
mpz_invert (mpz_t r, const mpz_t u, const mpz_t m)
3066
63
{
3067
63
  mpz_t g, tr;
3068
63
  int invertible;
3069
3070
63
  if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0)
3071
2
    return 0;
3072
3073
61
  mpz_init (g);
3074
61
  mpz_init (tr);
3075
3076
61
  mpz_gcdext (g, tr, NULL, u, m);
3077
61
  invertible = (mpz_cmp_ui (g, 1) == 0);
3078
3079
61
  if (invertible)
3080
45
    {
3081
45
      if (tr->_mp_size < 0)
3082
23
  {
3083
23
    if (m->_mp_size >= 0)
3084
23
      mpz_add (tr, tr, m);
3085
0
    else
3086
0
      mpz_sub (tr, tr, m);
3087
23
  }
3088
45
      mpz_swap (r, tr);
3089
45
    }
3090
3091
61
  mpz_clear (g);
3092
61
  mpz_clear (tr);
3093
61
  return invertible;
3094
63
}
3095
3096

3097
/* Higher level operations (sqrt, pow and root) */
3098
3099
void
3100
mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
3101
142
{
3102
142
  unsigned long bit;
3103
142
  mpz_t tr;
3104
142
  mpz_init_set_ui (tr, 1);
3105
3106
142
  bit = GMP_ULONG_HIGHBIT;
3107
142
  do
3108
9.08k
    {
3109
9.08k
      mpz_mul (tr, tr, tr);
3110
9.08k
      if (e & bit)
3111
150
  mpz_mul (tr, tr, b);
3112
9.08k
      bit >>= 1;
3113
9.08k
    }
3114
9.08k
  while (bit > 0);
3115
3116
142
  mpz_swap (r, tr);
3117
142
  mpz_clear (tr);
3118
142
}
3119
3120
void
3121
mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
3122
0
{
3123
0
  mpz_t b;
3124
3125
0
  mpz_init_set_ui (b, blimb);
3126
0
  mpz_pow_ui (r, b, e);
3127
0
  mpz_clear (b);
3128
0
}
3129
3130
void
3131
mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
3132
148
{
3133
148
  mpz_t tr;
3134
148
  mpz_t base;
3135
148
  mp_size_t en, mn;
3136
148
  mp_srcptr mp;
3137
148
  struct gmp_div_inverse minv;
3138
148
  unsigned shift;
3139
148
  mp_ptr tp = NULL;
3140
3141
148
  en = GMP_ABS (e->_mp_size);
3142
148
  mn = GMP_ABS (m->_mp_size);
3143
148
  if (mn == 0)
3144
0
    gmp_die ("mpz_powm: Zero modulo.");
3145
3146
148
  if (en == 0)
3147
0
    {
3148
0
      mpz_set_ui (r, mpz_cmpabs_ui (m, 1));
3149
0
      return;
3150
0
    }
3151
3152
148
  mp = m->_mp_d;
3153
148
  mpn_div_qr_invert (&minv, mp, mn);
3154
148
  shift = minv.shift;
3155
3156
148
  if (shift > 0)
3157
142
    {
3158
      /* To avoid shifts, we do all our reductions, except the final
3159
   one, using a *normalized* m. */
3160
142
      minv.shift = 0;
3161
3162
142
      tp = gmp_alloc_limbs (mn);
3163
142
      gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift));
3164
142
      mp = tp;
3165
142
    }
3166
3167
148
  mpz_init (base);
3168
3169
148
  if (e->_mp_size < 0)
3170
0
    {
3171
0
      if (!mpz_invert (base, b, m))
3172
0
  gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
3173
0
    }
3174
148
  else
3175
148
    {
3176
148
      mp_size_t bn;
3177
148
      mpz_abs (base, b);
3178
3179
148
      bn = base->_mp_size;
3180
148
      if (bn >= mn)
3181
25
  {
3182
25
    mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv);
3183
25
    bn = mn;
3184
25
  }
3185
3186
      /* We have reduced the absolute value. Now take care of the
3187
   sign. Note that we get zero represented non-canonically as
3188
   m. */
3189
148
      if (b->_mp_size < 0)
3190
0
  {
3191
0
    mp_ptr bp = MPZ_REALLOC (base, mn);
3192
0
    gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
3193
0
    bn = mn;
3194
0
  }
3195
148
      base->_mp_size = mpn_normalized_size (base->_mp_d, bn);
3196
148
    }
3197
148
  mpz_init_set_ui (tr, 1);
3198
3199
442
  while (--en >= 0)
3200
294
    {
3201
294
      mp_limb_t w = e->_mp_d[en];
3202
294
      mp_limb_t bit;
3203
3204
294
      bit = GMP_LIMB_HIGHBIT;
3205
294
      do
3206
18.8k
  {
3207
18.8k
    mpz_mul (tr, tr, tr);
3208
18.8k
    if (w & bit)
3209
6.70k
      mpz_mul (tr, tr, base);
3210
18.8k
    if (tr->_mp_size > mn)
3211
13.0k
      {
3212
13.0k
        mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3213
13.0k
        tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3214
13.0k
      }
3215
18.8k
    bit >>= 1;
3216
18.8k
  }
3217
18.8k
      while (bit > 0);
3218
294
    }
3219
3220
  /* Final reduction */
3221
148
  if (tr->_mp_size >= mn)
3222
147
    {
3223
147
      minv.shift = shift;
3224
147
      mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3225
147
      tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3226
147
    }
3227
148
  if (tp)
3228
142
    gmp_free_limbs (tp, mn);
3229
3230
148
  mpz_swap (r, tr);
3231
148
  mpz_clear (tr);
3232
148
  mpz_clear (base);
3233
148
}
3234
3235
void
3236
mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3237
51
{
3238
51
  mpz_t e;
3239
3240
51
  mpz_init_set_ui (e, elimb);
3241
51
  mpz_powm (r, b, e, m);
3242
51
  mpz_clear (e);
3243
51
}
3244
3245
/* x=trunc(y^(1/z)), r=y-x^z */
3246
void
3247
mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
3248
85
{
3249
85
  int sgn;
3250
85
  mp_bitcnt_t bc;
3251
85
  mpz_t t, u;
3252
3253
85
  sgn = y->_mp_size < 0;
3254
85
  if ((~z & sgn) != 0)
3255
0
    gmp_die ("mpz_rootrem: Negative argument, with even root.");
3256
85
  if (z == 0)
3257
0
    gmp_die ("mpz_rootrem: Zeroth root.");
3258
3259
85
  if (mpz_cmpabs_ui (y, 1) <= 0) {
3260
4
    if (x)
3261
4
      mpz_set (x, y);
3262
4
    if (r)
3263
1
      r->_mp_size = 0;
3264
4
    return;
3265
4
  }
3266
3267
81
  mpz_init (u);
3268
81
  mpz_init (t);
3269
81
  bc = (mpz_sizeinbase (y, 2) - 1) / z + 1;
3270
81
  mpz_setbit (t, bc);
3271
3272
81
  if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3273
565
    do {
3274
565
      mpz_swap (u, t);      /* u = x */
3275
565
      mpz_tdiv_q (t, y, u);   /* t = y/x */
3276
565
      mpz_add (t, t, u);    /* t = y/x + x */
3277
565
      mpz_tdiv_q_2exp (t, t, 1);  /* x'= (y/x + x)/2 */
3278
565
    } while (mpz_cmpabs (t, u) < 0);  /* |x'| < |x| */
3279
8
  else /* z != 2 */ {
3280
8
    mpz_t v;
3281
3282
8
    mpz_init (v);
3283
8
    if (sgn)
3284
0
      mpz_neg (t, t);
3285
3286
65
    do {
3287
65
      mpz_swap (u, t);      /* u = x */
3288
65
      mpz_pow_ui (t, u, z - 1);   /* t = x^(z-1) */
3289
65
      mpz_tdiv_q (t, y, t);   /* t = y/x^(z-1) */
3290
65
      mpz_mul_ui (v, u, z - 1);   /* v = x*(z-1) */
3291
65
      mpz_add (t, t, v);    /* t = y/x^(z-1) + x*(z-1) */
3292
65
      mpz_tdiv_q_ui (t, t, z);    /* x'=(y/x^(z-1) + x*(z-1))/z */
3293
65
    } while (mpz_cmpabs (t, u) < 0);  /* |x'| < |x| */
3294
3295
8
    mpz_clear (v);
3296
8
  }
3297
3298
81
  if (r) {
3299
45
    mpz_pow_ui (t, u, z);
3300
45
    mpz_sub (r, y, t);
3301
45
  }
3302
81
  if (x)
3303
81
    mpz_swap (x, u);
3304
81
  mpz_clear (u);
3305
81
  mpz_clear (t);
3306
81
}
3307
3308
int
3309
mpz_root (mpz_t x, const mpz_t y, unsigned long z)
3310
22
{
3311
22
  int res;
3312
22
  mpz_t r;
3313
3314
22
  mpz_init (r);
3315
22
  mpz_rootrem (x, r, y, z);
3316
22
  res = r->_mp_size == 0;
3317
22
  mpz_clear (r);
3318
3319
22
  return res;
3320
22
}
3321
3322
/* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
3323
void
3324
mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
3325
20
{
3326
20
  mpz_rootrem (s, r, u, 2);
3327
20
}
3328
3329
void
3330
mpz_sqrt (mpz_t s, const mpz_t u)
3331
39
{
3332
39
  mpz_rootrem (s, NULL, u, 2);
3333
39
}
3334
3335
int
3336
mpz_perfect_square_p (const mpz_t u)
3337
0
{
3338
0
  if (u->_mp_size <= 0)
3339
0
    return (u->_mp_size == 0);
3340
0
  else
3341
0
    return mpz_root (NULL, u, 2);
3342
0
}
3343
3344
int
3345
mpn_perfect_square_p (mp_srcptr p, mp_size_t n)
3346
0
{
3347
0
  mpz_t t;
3348
3349
0
  assert (n > 0);
3350
0
  assert (p [n-1] != 0);
3351
0
  return mpz_root (NULL, mpz_roinit_normal_n (t, p, n), 2);
3352
0
}
3353
3354
mp_size_t
3355
mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n)
3356
0
{
3357
0
  mpz_t s, r, u;
3358
0
  mp_size_t res;
3359
3360
0
  assert (n > 0);
3361
0
  assert (p [n-1] != 0);
3362
3363
0
  mpz_init (r);
3364
0
  mpz_init (s);
3365
0
  mpz_rootrem (s, r, mpz_roinit_normal_n (u, p, n), 2);
3366
3367
0
  assert (s->_mp_size == (n+1)/2);
3368
0
  mpn_copyd (sp, s->_mp_d, s->_mp_size);
3369
0
  mpz_clear (s);
3370
0
  res = r->_mp_size;
3371
0
  if (rp)
3372
0
    mpn_copyd (rp, r->_mp_d, res);
3373
0
  mpz_clear (r);
3374
0
  return res;
3375
0
}
3376

3377
/* Combinatorics */
3378
3379
void
3380
mpz_mfac_uiui (mpz_t x, unsigned long n, unsigned long m)
3381
2
{
3382
2
  mpz_set_ui (x, n + (n == 0));
3383
2
  if (m + 1 < 2) return;
3384
2
  while (n > m + 1)
3385
0
    mpz_mul_ui (x, x, n -= m);
3386
2
}
3387
3388
void
3389
mpz_2fac_ui (mpz_t x, unsigned long n)
3390
0
{
3391
0
  mpz_mfac_uiui (x, n, 2);
3392
0
}
3393
3394
void
3395
mpz_fac_ui (mpz_t x, unsigned long n)
3396
2
{
3397
2
  mpz_mfac_uiui (x, n, 1);
3398
2
}
3399
3400
void
3401
mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
3402
0
{
3403
0
  mpz_t t;
3404
3405
0
  mpz_set_ui (r, k <= n);
3406
3407
0
  if (k > (n >> 1))
3408
0
    k = (k <= n) ? n - k : 0;
3409
3410
0
  mpz_init (t);
3411
0
  mpz_fac_ui (t, k);
3412
3413
0
  for (; k > 0; --k)
3414
0
    mpz_mul_ui (r, r, n--);
3415
3416
0
  mpz_divexact (r, r, t);
3417
0
  mpz_clear (t);
3418
0
}
3419
3420

3421
/* Primality testing */
3422
3423
/* Computes Kronecker (a/b) with odd b, a!=0 and GCD(a,b) = 1 */
3424
/* Adapted from JACOBI_BASE_METHOD==4 in mpn/generic/jacbase.c */
3425
static int
3426
gmp_jacobi_coprime (mp_limb_t a, mp_limb_t b)
3427
35
{
3428
35
  int c, bit = 0;
3429
3430
35
  assert (b & 1);
3431
35
  assert (a != 0);
3432
  /* assert (mpn_gcd_11 (a, b) == 1); */
3433
3434
  /* Below, we represent a and b shifted right so that the least
3435
     significant one bit is implicit. */
3436
35
  b >>= 1;
3437
3438
35
  gmp_ctz(c, a);
3439
35
  a >>= 1;
3440
3441
35
  for (;;)
3442
52
    {
3443
52
      a >>= c;
3444
      /* (2/b) = -1 if b = 3 or 5 mod 8 */
3445
52
      bit ^= c & (b ^ (b >> 1));
3446
52
      if (a < b)
3447
52
  {
3448
52
    if (a == 0)
3449
35
      return bit & 1 ? -1 : 1;
3450
17
    bit ^= a & b;
3451
17
    a = b - a;
3452
17
    b -= a;
3453
17
  }
3454
0
      else
3455
0
  {
3456
0
    a -= b;
3457
0
    assert (a != 0);
3458
0
  }
3459
3460
17
      gmp_ctz(c, a);
3461
17
      ++c;
3462
17
    }
3463
35
}
3464
3465
static void
3466
gmp_lucas_step_k_2k (mpz_t V, mpz_t Qk, const mpz_t n)
3467
3.11k
{
3468
3.11k
  mpz_mod (Qk, Qk, n);
3469
  /* V_{2k} <- V_k ^ 2 - 2Q^k */
3470
3.11k
  mpz_mul (V, V, V);
3471
3.11k
  mpz_submul_ui (V, Qk, 2);
3472
3.11k
  mpz_tdiv_r (V, V, n);
3473
  /* Q^{2k} = (Q^k)^2 */
3474
3.11k
  mpz_mul (Qk, Qk, Qk);
3475
3.11k
}
3476
3477
/* Computes V_k, Q^k (mod n) for the Lucas' sequence */
3478
/* with P=1, Q=Q; k = (n>>b0)|1. */
3479
/* Requires an odd n > 4; b0 > 0; -2*Q must not overflow a long */
3480
/* Returns (U_k == 0) and sets V=V_k and Qk=Q^k. */
3481
static int
3482
gmp_lucas_mod (mpz_t V, mpz_t Qk, long Q,
3483
         mp_bitcnt_t b0, const mpz_t n)
3484
17
{
3485
17
  mp_bitcnt_t bs;
3486
17
  mpz_t U;
3487
17
  int res;
3488
3489
17
  assert (b0 > 0);
3490
17
  assert (Q <= - (LONG_MIN / 2));
3491
17
  assert (Q >= - (LONG_MAX / 2));
3492
17
  assert (mpz_cmp_ui (n, 4) > 0);
3493
17
  assert (mpz_odd_p (n));
3494
3495
17
  mpz_init_set_ui (U, 1); /* U1 = 1 */
3496
17
  mpz_set_ui (V, 1); /* V1 = 1 */
3497
17
  mpz_set_si (Qk, Q);
3498
3499
3.12k
  for (bs = mpz_sizeinbase (n, 2) - 1; --bs >= b0;)
3500
3.11k
    {
3501
      /* U_{2k} <- U_k * V_k */
3502
3.11k
      mpz_mul (U, U, V);
3503
      /* V_{2k} <- V_k ^ 2 - 2Q^k */
3504
      /* Q^{2k} = (Q^k)^2 */
3505
3.11k
      gmp_lucas_step_k_2k (V, Qk, n);
3506
3507
      /* A step k->k+1 is performed if the bit in $n$ is 1  */
3508
      /* mpz_tstbit(n,bs) or the bit is 0 in $n$ but  */
3509
      /* should be 1 in $n+1$ (bs == b0)      */
3510
3.11k
      if (b0 == bs || mpz_tstbit (n, bs))
3511
1.56k
  {
3512
    /* Q^{k+1} <- Q^k * Q */
3513
1.56k
    mpz_mul_si (Qk, Qk, Q);
3514
    /* U_{k+1} <- (U_k + V_k) / 2 */
3515
1.56k
    mpz_swap (U, V); /* Keep in V the old value of U_k */
3516
1.56k
    mpz_add (U, U, V);
3517
    /* We have to compute U/2, so we need an even value, */
3518
    /* equivalent (mod n) */
3519
1.56k
    if (mpz_odd_p (U))
3520
749
      mpz_add (U, U, n);
3521
1.56k
    mpz_tdiv_q_2exp (U, U, 1);
3522
    /* V_{k+1} <-(D*U_k + V_k) / 2 =
3523
      U_{k+1} + (D-1)/2*U_k = U_{k+1} - 2Q*U_k */
3524
1.56k
    mpz_mul_si (V, V, -2*Q);
3525
1.56k
    mpz_add (V, U, V);
3526
1.56k
    mpz_tdiv_r (V, V, n);
3527
1.56k
  }
3528
3.11k
      mpz_tdiv_r (U, U, n);
3529
3.11k
    }
3530
3531
17
  res = U->_mp_size == 0;
3532
17
  mpz_clear (U);
3533
17
  return res;
3534
17
}
3535
3536
/* Performs strong Lucas' test on x, with parameters suggested */
3537
/* for the BPSW test. Qk is only passed to recycle a variable. */
3538
/* Requires GCD (x,6) = 1.*/
3539
static int
3540
gmp_stronglucas (const mpz_t x, mpz_t Qk)
3541
17
{
3542
17
  mp_bitcnt_t b0;
3543
17
  mpz_t V, n;
3544
17
  mp_limb_t maxD, D; /* The absolute value is stored. */
3545
17
  long Q;
3546
17
  mp_limb_t tl;
3547
3548
  /* Test on the absolute value. */
3549
17
  mpz_roinit_normal_n (n, x->_mp_d, GMP_ABS (x->_mp_size));
3550
3551
17
  assert (mpz_odd_p (n));
3552
  /* assert (mpz_gcd_ui (NULL, n, 6) == 1); */
3553
17
  if (mpz_root (Qk, n, 2))
3554
0
    return 0; /* A square is composite. */
3555
3556
  /* Check Ds up to square root (in case, n is prime)
3557
     or avoid overflows */
3558
17
  maxD = (Qk->_mp_size == 1) ? Qk->_mp_d [0] - 1 : GMP_LIMB_MAX;
3559
3560
17
  D = 3;
3561
  /* Search a D such that (D/n) = -1 in the sequence 5,-7,9,-11,.. */
3562
  /* For those Ds we have (D/n) = (n/|D|) */
3563
17
  do
3564
35
    {
3565
35
      if (D >= maxD)
3566
0
  return 1 + (D != GMP_LIMB_MAX); /* (1 + ! ~ D) */
3567
35
      D += 2;
3568
35
      tl = mpz_tdiv_ui (n, D);
3569
35
      if (tl == 0)
3570
0
  return 0;
3571
35
    }
3572
35
  while (gmp_jacobi_coprime (tl, D) == 1);
3573
3574
17
  mpz_init (V);
3575
3576
  /* n-(D/n) = n+1 = d*2^{b0}, with d = (n>>b0) | 1 */
3577
17
  b0 = mpn_common_scan (~ n->_mp_d[0], 0, n->_mp_d, n->_mp_size, GMP_LIMB_MAX);
3578
  /* b0 = mpz_scan0 (n, 0); */
3579
3580
  /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */
3581
17
  Q = (D & 2) ? (long) (D >> 2) + 1 : -(long) (D >> 2);
3582
3583
17
  if (! gmp_lucas_mod (V, Qk, Q, b0, n))  /* If Ud != 0 */
3584
13
    while (V->_mp_size != 0 && --b0 != 0)  /* while Vk != 0 */
3585
      /* V <- V ^ 2 - 2Q^k */
3586
      /* Q^{2k} = (Q^k)^2 */
3587
7
      gmp_lucas_step_k_2k (V, Qk, n);
3588
3589
17
  mpz_clear (V);
3590
17
  return (b0 != 0);
3591
17
}
3592
3593
static int
3594
gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
3595
     const mpz_t q, mp_bitcnt_t k)
3596
24
{
3597
24
  assert (k > 0);
3598
3599
  /* Caller must initialize y to the base. */
3600
24
  mpz_powm (y, y, q, n);
3601
3602
24
  if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
3603
10
    return 1;
3604
3605
58
  while (--k > 0)
3606
51
    {
3607
51
      mpz_powm_ui (y, y, 2, n);
3608
51
      if (mpz_cmp (y, nm1) == 0)
3609
7
  return 1;
3610
51
    }
3611
7
  return 0;
3612
14
}
3613
3614
/* This product is 0xc0cfd797, and fits in 32 bits. */
3615
#define GMP_PRIME_PRODUCT \
3616
46
  (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)
3617
3618
/* Bit (p+1)/2 is set, for each odd prime <= 61 */
3619
0
#define GMP_PRIME_MASK 0xc96996dcUL
3620
3621
int
3622
mpz_probab_prime_p (const mpz_t n, int reps)
3623
55
{
3624
55
  mpz_t nm1;
3625
55
  mpz_t q;
3626
55
  mpz_t y;
3627
55
  mp_bitcnt_t k;
3628
55
  int is_prime;
3629
55
  int j;
3630
3631
  /* Note that we use the absolute value of n only, for compatibility
3632
     with the real GMP. */
3633
55
  if (mpz_even_p (n))
3634
9
    return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0;
3635
3636
  /* Above test excludes n == 0 */
3637
46
  assert (n->_mp_size != 0);
3638
3639
46
  if (mpz_cmpabs_ui (n, 64) < 0)
3640
0
    return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2;
3641
3642
46
  if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
3643
22
    return 0;
3644
3645
  /* All prime factors are >= 31. */
3646
24
  if (mpz_cmpabs_ui (n, 31*31) < 0)
3647
0
    return 2;
3648
3649
24
  mpz_init (nm1);
3650
24
  mpz_init (q);
3651
3652
  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
3653
24
  mpz_abs (nm1, n);
3654
24
  nm1->_mp_d[0] -= 1;
3655
  /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
3656
24
  k = mpn_scan1 (nm1->_mp_d, 0);
3657
24
  mpz_tdiv_q_2exp (q, nm1, k);
3658
3659
  /* BPSW test */
3660
24
  mpz_init_set_ui (y, 2);
3661
24
  is_prime = gmp_millerrabin (n, nm1, y, q, k) && gmp_stronglucas (n, y);
3662
24
  reps -= 24; /* skip the first 24 repetitions */
3663
3664
  /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
3665
     j^2 + j + 41 using Euler's polynomial. We potentially stop early,
3666
     if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
3667
     30 (a[30] == 971 > 31*31 == 961). */
3668
3669
24
  for (j = 0; is_prime & (j < reps); j++)
3670
0
    {
3671
0
      mpz_set_ui (y, (unsigned long) j*j+j+41);
3672
0
      if (mpz_cmp (y, nm1) >= 0)
3673
0
  {
3674
    /* Don't try any further bases. This "early" break does not affect
3675
       the result for any reasonable reps value (<=5000 was tested) */
3676
0
    assert (j >= 30);
3677
0
    break;
3678
0
  }
3679
0
      is_prime = gmp_millerrabin (n, nm1, y, q, k);
3680
0
    }
3681
24
  mpz_clear (nm1);
3682
24
  mpz_clear (q);
3683
24
  mpz_clear (y);
3684
3685
24
  return is_prime;
3686
24
}
3687
3688

3689
/* Logical operations and bit manipulation. */
3690
3691
/* Numbers are treated as if represented in two's complement (and
3692
   infinitely sign extended). For a negative values we get the two's
3693
   complement from -x = ~x + 1, where ~ is bitwise complement.
3694
   Negation transforms
3695
3696
     xxxx10...0
3697
3698
   into
3699
3700
     yyyy10...0
3701
3702
   where yyyy is the bitwise complement of xxxx. So least significant
3703
   bits, up to and including the first one bit, are unchanged, and
3704
   the more significant bits are all complemented.
3705
3706
   To change a bit from zero to one in a negative number, subtract the
3707
   corresponding power of two from the absolute value. This can never
3708
   underflow. To change a bit from one to zero, add the corresponding
3709
   power of two, and this might overflow. E.g., if x = -001111, the
3710
   two's complement is 110001. Clearing the least significant bit, we
3711
   get two's complement 110000, and -010000. */
3712
3713
int
3714
mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
3715
3.32k
{
3716
3.32k
  mp_size_t limb_index;
3717
3.32k
  unsigned shift;
3718
3.32k
  mp_size_t ds;
3719
3.32k
  mp_size_t dn;
3720
3.32k
  mp_limb_t w;
3721
3.32k
  int bit;
3722
3723
3.32k
  ds = d->_mp_size;
3724
3.32k
  dn = GMP_ABS (ds);
3725
3.32k
  limb_index = bit_index / GMP_LIMB_BITS;
3726
3.32k
  if (limb_index >= dn)
3727
232
    return ds < 0;
3728
3729
3.09k
  shift = bit_index % GMP_LIMB_BITS;
3730
3.09k
  w = d->_mp_d[limb_index];
3731
3.09k
  bit = (w >> shift) & 1;
3732
3733
3.09k
  if (ds < 0)
3734
0
    {
3735
      /* d < 0. Check if any of the bits below is set: If so, our bit
3736
   must be complemented. */
3737
0
      if (shift > 0 && (mp_limb_t) (w << (GMP_LIMB_BITS - shift)) > 0)
3738
0
  return bit ^ 1;
3739
0
      while (--limb_index >= 0)
3740
0
  if (d->_mp_d[limb_index] > 0)
3741
0
    return bit ^ 1;
3742
0
    }
3743
3.09k
  return bit;
3744
3.09k
}
3745
3746
static void
3747
mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
3748
231
{
3749
231
  mp_size_t dn, limb_index;
3750
231
  mp_limb_t bit;
3751
231
  mp_ptr dp;
3752
3753
231
  dn = GMP_ABS (d->_mp_size);
3754
3755
231
  limb_index = bit_index / GMP_LIMB_BITS;
3756
231
  bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3757
3758
231
  if (limb_index >= dn)
3759
231
    {
3760
231
      mp_size_t i;
3761
      /* The bit should be set outside of the end of the number.
3762
   We have to increase the size of the number. */
3763
231
      dp = MPZ_REALLOC (d, limb_index + 1);
3764
3765
231
      dp[limb_index] = bit;
3766
1.41k
      for (i = dn; i < limb_index; i++)
3767
1.18k
  dp[i] = 0;
3768
231
      dn = limb_index + 1;
3769
231
    }
3770
0
  else
3771
0
    {
3772
0
      mp_limb_t cy;
3773
3774
0
      dp = d->_mp_d;
3775
3776
0
      cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit);
3777
0
      if (cy > 0)
3778
0
  {
3779
0
    dp = MPZ_REALLOC (d, dn + 1);
3780
0
    dp[dn++] = cy;
3781
0
  }
3782
0
    }
3783
3784
231
  d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3785
231
}
3786
3787
static void
3788
mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
3789
1
{
3790
1
  mp_size_t dn, limb_index;
3791
1
  mp_ptr dp;
3792
1
  mp_limb_t bit;
3793
3794
1
  dn = GMP_ABS (d->_mp_size);
3795
1
  dp = d->_mp_d;
3796
3797
1
  limb_index = bit_index / GMP_LIMB_BITS;
3798
1
  bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3799
3800
1
  assert (limb_index < dn);
3801
3802
1
  gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3803
1
         dn - limb_index, bit));
3804
1
  dn = mpn_normalized_size (dp, dn);
3805
1
  d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3806
1
}
3807
3808
void
3809
mpz_setbit (mpz_t d, mp_bitcnt_t bit_index)
3810
232
{
3811
232
  if (!mpz_tstbit (d, bit_index))
3812
231
    {
3813
231
      if (d->_mp_size >= 0)
3814
231
  mpz_abs_add_bit (d, bit_index);
3815
0
      else
3816
0
  mpz_abs_sub_bit (d, bit_index);
3817
231
    }
3818
232
}
3819
3820
void
3821
mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index)
3822
2
{
3823
2
  if (mpz_tstbit (d, bit_index))
3824
1
    {
3825
1
      if (d->_mp_size >= 0)
3826
1
  mpz_abs_sub_bit (d, bit_index);
3827
0
      else
3828
0
  mpz_abs_add_bit (d, bit_index);
3829
1
    }
3830
2
}
3831
3832
void
3833
mpz_combit (mpz_t d, mp_bitcnt_t bit_index)
3834
0
{
3835
0
  if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0))
3836
0
    mpz_abs_sub_bit (d, bit_index);
3837
0
  else
3838
0
    mpz_abs_add_bit (d, bit_index);
3839
0
}
3840
3841
void
3842
mpz_com (mpz_t r, const mpz_t u)
3843
0
{
3844
0
  mpz_add_ui (r, u, 1);
3845
0
  mpz_neg (r, r);
3846
0
}
3847
3848
void
3849
mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
3850
9
{
3851
9
  mp_size_t un, vn, rn, i;
3852
9
  mp_ptr up, vp, rp;
3853
3854
9
  mp_limb_t ux, vx, rx;
3855
9
  mp_limb_t uc, vc, rc;
3856
9
  mp_limb_t ul, vl, rl;
3857
3858
9
  un = GMP_ABS (u->_mp_size);
3859
9
  vn = GMP_ABS (v->_mp_size);
3860
9
  if (un < vn)
3861
2
    {
3862
2
      MPZ_SRCPTR_SWAP (u, v);
3863
2
      MP_SIZE_T_SWAP (un, vn);
3864
2
    }
3865
9
  if (vn == 0)
3866
2
    {
3867
2
      r->_mp_size = 0;
3868
2
      return;
3869
2
    }
3870
3871
7
  uc = u->_mp_size < 0;
3872
7
  vc = v->_mp_size < 0;
3873
7
  rc = uc & vc;
3874
3875
7
  ux = -uc;
3876
7
  vx = -vc;
3877
7
  rx = -rc;
3878
3879
  /* If the smaller input is positive, higher limbs don't matter. */
3880
7
  rn = vx ? un : vn;
3881
3882
7
  rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3883
3884
7
  up = u->_mp_d;
3885
7
  vp = v->_mp_d;
3886
3887
7
  i = 0;
3888
7
  do
3889
97
    {
3890
97
      ul = (up[i] ^ ux) + uc;
3891
97
      uc = ul < uc;
3892
3893
97
      vl = (vp[i] ^ vx) + vc;
3894
97
      vc = vl < vc;
3895
3896
97
      rl = ( (ul & vl) ^ rx) + rc;
3897
97
      rc = rl < rc;
3898
97
      rp[i] = rl;
3899
97
    }
3900
97
  while (++i < vn);
3901
7
  assert (vc == 0);
3902
3903
7
  for (; i < rn; i++)
3904
0
    {
3905
0
      ul = (up[i] ^ ux) + uc;
3906
0
      uc = ul < uc;
3907
3908
0
      rl = ( (ul & vx) ^ rx) + rc;
3909
0
      rc = rl < rc;
3910
0
      rp[i] = rl;
3911
0
    }
3912
7
  if (rc)
3913
0
    rp[rn++] = rc;
3914
7
  else
3915
7
    rn = mpn_normalized_size (rp, rn);
3916
3917
7
  r->_mp_size = rx ? -rn : rn;
3918
7
}
3919
3920
void
3921
mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
3922
13
{
3923
13
  mp_size_t un, vn, rn, i;
3924
13
  mp_ptr up, vp, rp;
3925
3926
13
  mp_limb_t ux, vx, rx;
3927
13
  mp_limb_t uc, vc, rc;
3928
13
  mp_limb_t ul, vl, rl;
3929
3930
13
  un = GMP_ABS (u->_mp_size);
3931
13
  vn = GMP_ABS (v->_mp_size);
3932
13
  if (un < vn)
3933
6
    {
3934
6
      MPZ_SRCPTR_SWAP (u, v);
3935
6
      MP_SIZE_T_SWAP (un, vn);
3936
6
    }
3937
13
  if (vn == 0)
3938
5
    {
3939
5
      mpz_set (r, u);
3940
5
      return;
3941
5
    }
3942
3943
8
  uc = u->_mp_size < 0;
3944
8
  vc = v->_mp_size < 0;
3945
8
  rc = uc | vc;
3946
3947
8
  ux = -uc;
3948
8
  vx = -vc;
3949
8
  rx = -rc;
3950
3951
  /* If the smaller input is negative, by sign extension higher limbs
3952
     don't matter. */
3953
8
  rn = vx ? vn : un;
3954
3955
8
  rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3956
3957
8
  up = u->_mp_d;
3958
8
  vp = v->_mp_d;
3959
3960
8
  i = 0;
3961
8
  do
3962
148
    {
3963
148
      ul = (up[i] ^ ux) + uc;
3964
148
      uc = ul < uc;
3965
3966
148
      vl = (vp[i] ^ vx) + vc;
3967
148
      vc = vl < vc;
3968
3969
148
      rl = ( (ul | vl) ^ rx) + rc;
3970
148
      rc = rl < rc;
3971
148
      rp[i] = rl;
3972
148
    }
3973
148
  while (++i < vn);
3974
8
  assert (vc == 0);
3975
3976
264
  for (; i < rn; i++)
3977
256
    {
3978
256
      ul = (up[i] ^ ux) + uc;
3979
256
      uc = ul < uc;
3980
3981
256
      rl = ( (ul | vx) ^ rx) + rc;
3982
256
      rc = rl < rc;
3983
256
      rp[i] = rl;
3984
256
    }
3985
8
  if (rc)
3986
0
    rp[rn++] = rc;
3987
8
  else
3988
8
    rn = mpn_normalized_size (rp, rn);
3989
3990
8
  r->_mp_size = rx ? -rn : rn;
3991
8
}
3992
3993
void
3994
mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
3995
15
{
3996
15
  mp_size_t un, vn, i;
3997
15
  mp_ptr up, vp, rp;
3998
3999
15
  mp_limb_t ux, vx, rx;
4000
15
  mp_limb_t uc, vc, rc;
4001
15
  mp_limb_t ul, vl, rl;
4002
4003
15
  un = GMP_ABS (u->_mp_size);
4004
15
  vn = GMP_ABS (v->_mp_size);
4005
15
  if (un < vn)
4006
9
    {
4007
9
      MPZ_SRCPTR_SWAP (u, v);
4008
9
      MP_SIZE_T_SWAP (un, vn);
4009
9
    }
4010
15
  if (vn == 0)
4011
1
    {
4012
1
      mpz_set (r, u);
4013
1
      return;
4014
1
    }
4015
4016
14
  uc = u->_mp_size < 0;
4017
14
  vc = v->_mp_size < 0;
4018
14
  rc = uc ^ vc;
4019
4020
14
  ux = -uc;
4021
14
  vx = -vc;
4022
14
  rx = -rc;
4023
4024
14
  rp = MPZ_REALLOC (r, un + (mp_size_t) rc);
4025
4026
14
  up = u->_mp_d;
4027
14
  vp = v->_mp_d;
4028
4029
14
  i = 0;
4030
14
  do
4031
244
    {
4032
244
      ul = (up[i] ^ ux) + uc;
4033
244
      uc = ul < uc;
4034
4035
244
      vl = (vp[i] ^ vx) + vc;
4036
244
      vc = vl < vc;
4037
4038
244
      rl = (ul ^ vl ^ rx) + rc;
4039
244
      rc = rl < rc;
4040
244
      rp[i] = rl;
4041
244
    }
4042
244
  while (++i < vn);
4043
14
  assert (vc == 0);
4044
4045
685
  for (; i < un; i++)
4046
671
    {
4047
671
      ul = (up[i] ^ ux) + uc;
4048
671
      uc = ul < uc;
4049
4050
671
      rl = (ul ^ ux) + rc;
4051
671
      rc = rl < rc;
4052
671
      rp[i] = rl;
4053
671
    }
4054
14
  if (rc)
4055
0
    rp[un++] = rc;
4056
14
  else
4057
14
    un = mpn_normalized_size (rp, un);
4058
4059
14
  r->_mp_size = rx ? -un : un;
4060
14
}
4061
4062
static unsigned
4063
gmp_popcount_limb (mp_limb_t x)
4064
1.25k
{
4065
1.25k
  unsigned c;
4066
4067
  /* Do 16 bits at a time, to avoid limb-sized constants. */
4068
1.25k
  int LOCAL_SHIFT_BITS = 16;
4069
6.10k
  for (c = 0; x > 0;)
4070
4.84k
    {
4071
4.84k
      unsigned w = x - ((x >> 1) & 0x5555);
4072
4.84k
      w = ((w >> 2) & 0x3333) + (w & 0x3333);
4073
4.84k
      w =  (w >> 4) + w;
4074
4.84k
      w = ((w >> 8) & 0x000f) + (w & 0x000f);
4075
4.84k
      c += w;
4076
4.84k
      if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS)
4077
4.84k
  x >>= LOCAL_SHIFT_BITS;
4078
0
      else
4079
0
  x = 0;
4080
4.84k
    }
4081
1.25k
  return c;
4082
1.25k
}
4083
4084
mp_bitcnt_t
4085
mpn_popcount (mp_srcptr p, mp_size_t n)
4086
10
{
4087
10
  mp_size_t i;
4088
10
  mp_bitcnt_t c;
4089
4090
421
  for (c = 0, i = 0; i < n; i++)
4091
411
    c += gmp_popcount_limb (p[i]);
4092
4093
10
  return c;
4094
10
}
4095
4096
mp_bitcnt_t
4097
mpz_popcount (const mpz_t u)
4098
10
{
4099
10
  mp_size_t un;
4100
4101
10
  un = u->_mp_size;
4102
4103
10
  if (un < 0)
4104
0
    return ~(mp_bitcnt_t) 0;
4105
4106
10
  return mpn_popcount (u->_mp_d, un);
4107
10
}
4108
4109
mp_bitcnt_t
4110
mpz_hamdist (const mpz_t u, const mpz_t v)
4111
16
{
4112
16
  mp_size_t un, vn, i;
4113
16
  mp_limb_t uc, vc, ul, vl, comp;
4114
16
  mp_srcptr up, vp;
4115
16
  mp_bitcnt_t c;
4116
4117
16
  un = u->_mp_size;
4118
16
  vn = v->_mp_size;
4119
4120
16
  if ( (un ^ vn) < 0)
4121
0
    return ~(mp_bitcnt_t) 0;
4122
4123
16
  comp = - (uc = vc = (un < 0));
4124
16
  if (uc)
4125
0
    {
4126
0
      assert (vn < 0);
4127
0
      un = -un;
4128
0
      vn = -vn;
4129
0
    }
4130
4131
16
  up = u->_mp_d;
4132
16
  vp = v->_mp_d;
4133
4134
16
  if (un < vn)
4135
6
    MPN_SRCPTR_SWAP (up, un, vp, vn);
4136
4137
290
  for (i = 0, c = 0; i < vn; i++)
4138
274
    {
4139
274
      ul = (up[i] ^ comp) + uc;
4140
274
      uc = ul < uc;
4141
4142
274
      vl = (vp[i] ^ comp) + vc;
4143
274
      vc = vl < vc;
4144
4145
274
      c += gmp_popcount_limb (ul ^ vl);
4146
274
    }
4147
16
  assert (vc == 0);
4148
4149
586
  for (; i < un; i++)
4150
570
    {
4151
570
      ul = (up[i] ^ comp) + uc;
4152
570
      uc = ul < uc;
4153
4154
570
      c += gmp_popcount_limb (ul ^ comp);
4155
570
    }
4156
4157
16
  return c;
4158
16
}
4159
4160
mp_bitcnt_t
4161
mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
4162
3
{
4163
3
  mp_ptr up;
4164
3
  mp_size_t us, un, i;
4165
3
  mp_limb_t limb, ux;
4166
4167
3
  us = u->_mp_size;
4168
3
  un = GMP_ABS (us);
4169
3
  i = starting_bit / GMP_LIMB_BITS;
4170
4171
  /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit
4172
     for u<0. Notice this test picks up any u==0 too. */
4173
3
  if (i >= un)
4174
1
    return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
4175
4176
2
  up = u->_mp_d;
4177
2
  ux = 0;
4178
2
  limb = up[i];
4179
4180
2
  if (starting_bit != 0)
4181
0
    {
4182
0
      if (us < 0)
4183
0
  {
4184
0
    ux = mpn_zero_p (up, i);
4185
0
    limb = ~ limb + ux;
4186
0
    ux = - (mp_limb_t) (limb >= ux);
4187
0
  }
4188
4189
      /* Mask to 0 all bits before starting_bit, thus ignoring them. */
4190
0
      limb &= GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS);
4191
0
    }
4192
4193
2
  return mpn_common_scan (limb, i, up, un, ux);
4194
3
}
4195
4196
mp_bitcnt_t
4197
mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
4198
0
{
4199
0
  mp_ptr up;
4200
0
  mp_size_t us, un, i;
4201
0
  mp_limb_t limb, ux;
4202
4203
0
  us = u->_mp_size;
4204
0
  ux = - (mp_limb_t) (us >= 0);
4205
0
  un = GMP_ABS (us);
4206
0
  i = starting_bit / GMP_LIMB_BITS;
4207
4208
  /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
4209
     u<0.  Notice this test picks up all cases of u==0 too. */
4210
0
  if (i >= un)
4211
0
    return (ux ? starting_bit : ~(mp_bitcnt_t) 0);
4212
4213
0
  up = u->_mp_d;
4214
0
  limb = up[i] ^ ux;
4215
4216
0
  if (ux == 0)
4217
0
    limb -= mpn_zero_p (up, i); /* limb = ~(~limb + zero_p) */
4218
4219
  /* Mask all bits before starting_bit, thus ignoring them. */
4220
0
  limb &= GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS);
4221
4222
0
  return mpn_common_scan (limb, i, up, un, ux);
4223
0
}
4224
4225

4226
/* MPZ base conversion. */
4227
4228
size_t
4229
mpz_sizeinbase (const mpz_t u, int base)
4230
1.53k
{
4231
1.53k
  mp_size_t un, tn;
4232
1.53k
  mp_srcptr up;
4233
1.53k
  mp_ptr tp;
4234
1.53k
  mp_bitcnt_t bits;
4235
1.53k
  struct gmp_div_inverse bi;
4236
1.53k
  size_t ndigits;
4237
4238
1.53k
  assert (base >= 2);
4239
1.53k
  assert (base <= 62);
4240
4241
1.53k
  un = GMP_ABS (u->_mp_size);
4242
1.53k
  if (un == 0)
4243
114
    return 1;
4244
4245
1.42k
  up = u->_mp_d;
4246
4247
1.42k
  bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]);
4248
1.42k
  switch (base)
4249
1.42k
    {
4250
977
    case 2:
4251
977
      return bits;
4252
0
    case 4:
4253
0
      return (bits + 1) / 2;
4254
0
    case 8:
4255
0
      return (bits + 2) / 3;
4256
0
    case 16:
4257
0
      return (bits + 3) / 4;
4258
0
    case 32:
4259
0
      return (bits + 4) / 5;
4260
      /* FIXME: Do something more clever for the common case of base
4261
   10. */
4262
1.42k
    }
4263
4264
445
  tp = gmp_alloc_limbs (un);
4265
445
  mpn_copyi (tp, up, un);
4266
445
  mpn_div_qr_1_invert (&bi, base);
4267
4268
445
  tn = un;
4269
445
  ndigits = 0;
4270
445
  do
4271
232k
    {
4272
232k
      ndigits++;
4273
232k
      mpn_div_qr_1_preinv (tp, tp, tn, &bi);
4274
232k
      tn -= (tp[tn-1] == 0);
4275
232k
    }
4276
232k
  while (tn > 0);
4277
4278
445
  gmp_free_limbs (tp, un);
4279
445
  return ndigits;
4280
1.42k
}
4281
4282
char *
4283
mpz_get_str (char *sp, int base, const mpz_t u)
4284
549
{
4285
549
  unsigned bits;
4286
549
  const char *digits;
4287
549
  mp_size_t un;
4288
549
  size_t i, sn, osn;
4289
4290
549
  digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
4291
549
  if (base > 1)
4292
549
    {
4293
549
      if (base <= 36)
4294
549
  digits = "0123456789abcdefghijklmnopqrstuvwxyz";
4295
0
      else if (base > 62)
4296
0
  return NULL;
4297
549
    }
4298
0
  else if (base >= -1)
4299
0
    base = 10;
4300
0
  else
4301
0
    {
4302
0
      base = -base;
4303
0
      if (base > 36)
4304
0
  return NULL;
4305
0
    }
4306
4307
549
  sn = 1 + mpz_sizeinbase (u, base);
4308
549
  if (!sp)
4309
549
    {
4310
549
      osn = 1 + sn;
4311
549
      sp = (char *) gmp_alloc (osn);
4312
549
    }
4313
0
  else
4314
0
    osn = 0;
4315
549
  un = GMP_ABS (u->_mp_size);
4316
4317
549
  if (un == 0)
4318
104
    {
4319
104
      sp[0] = '0';
4320
104
      sn = 1;
4321
104
      goto ret;
4322
104
    }
4323
4324
445
  i = 0;
4325
4326
445
  if (u->_mp_size < 0)
4327
27
    sp[i++] = '-';
4328
4329
445
  bits = mpn_base_power_of_two_p (base);
4330
4331
445
  if (bits)
4332
    /* Not modified in this case. */
4333
0
    sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un);
4334
445
  else
4335
445
    {
4336
445
      struct mpn_base_info info;
4337
445
      mp_ptr tp;
4338
4339
445
      mpn_get_base_info (&info, base);
4340
445
      tp = gmp_alloc_limbs (un);
4341
445
      mpn_copyi (tp, u->_mp_d, un);
4342
4343
445
      sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un);
4344
445
      gmp_free_limbs (tp, un);
4345
445
    }
4346
4347
232k
  for (; i < sn; i++)
4348
232k
    sp[i] = digits[(unsigned char) sp[i]];
4349
4350
549
ret:
4351
549
  sp[sn] = '\0';
4352
549
  if (osn && osn != sn + 1)
4353
522
    sp = (char*) gmp_realloc (sp, osn, sn + 1);
4354
549
  return sp;
4355
445
}
4356
4357
int
4358
mpz_set_str (mpz_t r, const char *sp, int base)
4359
4.83k
{
4360
4.83k
  unsigned bits, value_of_a;
4361
4.83k
  mp_size_t rn, alloc;
4362
4.83k
  mp_ptr rp;
4363
4.83k
  size_t dn, sn;
4364
4.83k
  int sign;
4365
4.83k
  unsigned char *dp;
4366
4367
4.83k
  assert (base == 0 || (base >= 2 && base <= 62));
4368
4369
4.83k
  while (isspace( (unsigned char) *sp))
4370
0
    sp++;
4371
4372
4.83k
  sign = (*sp == '-');
4373
4.83k
  sp += sign;
4374
4375
4.83k
  if (base == 0)
4376
0
    {
4377
0
      if (sp[0] == '0')
4378
0
  {
4379
0
    if (sp[1] == 'x' || sp[1] == 'X')
4380
0
      {
4381
0
        base = 16;
4382
0
        sp += 2;
4383
0
      }
4384
0
    else if (sp[1] == 'b' || sp[1] == 'B')
4385
0
      {
4386
0
        base = 2;
4387
0
        sp += 2;
4388
0
      }
4389
0
    else
4390
0
      base = 8;
4391
0
  }
4392
0
      else
4393
0
  base = 10;
4394
0
    }
4395
4396
4.83k
  if (!*sp)
4397
0
    {
4398
0
      r->_mp_size = 0;
4399
0
      return -1;
4400
0
    }
4401
4.83k
  sn = strlen(sp);
4402
4.83k
  dp = (unsigned char *) gmp_alloc (sn);
4403
4404
4.83k
  value_of_a = (base > 36) ? 36 : 10;
4405
2.02M
  for (dn = 0; *sp; sp++)
4406
2.01M
    {
4407
2.01M
      unsigned digit;
4408
4409
2.01M
      if (isspace ((unsigned char) *sp))
4410
0
  continue;
4411
2.01M
      else if (*sp >= '0' && *sp <= '9')
4412
2.01M
  digit = *sp - '0';
4413
0
      else if (*sp >= 'a' && *sp <= 'z')
4414
0
  digit = *sp - 'a' + value_of_a;
4415
0
      else if (*sp >= 'A' && *sp <= 'Z')
4416
0
  digit = *sp - 'A' + 10;
4417
0
      else
4418
0
  digit = base; /* fail */
4419
4420
2.01M
      if (digit >= (unsigned) base)
4421
0
  {
4422
0
    gmp_free (dp, sn);
4423
0
    r->_mp_size = 0;
4424
0
    return -1;
4425
0
  }
4426
4427
2.01M
      dp[dn++] = digit;
4428
2.01M
    }
4429
4430
4.83k
  if (!dn)
4431
0
    {
4432
0
      gmp_free (dp, sn);
4433
0
      r->_mp_size = 0;
4434
0
      return -1;
4435
0
    }
4436
4.83k
  bits = mpn_base_power_of_two_p (base);
4437
4438
4.83k
  if (bits > 0)
4439
0
    {
4440
0
      alloc = (dn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
4441
0
      rp = MPZ_REALLOC (r, alloc);
4442
0
      rn = mpn_set_str_bits (rp, dp, dn, bits);
4443
0
    }
4444
4.83k
  else
4445
4.83k
    {
4446
4.83k
      struct mpn_base_info info;
4447
4.83k
      mpn_get_base_info (&info, base);
4448
4.83k
      alloc = (dn + info.exp - 1) / info.exp;
4449
4.83k
      rp = MPZ_REALLOC (r, alloc);
4450
4.83k
      rn = mpn_set_str_other (rp, dp, dn, base, &info);
4451
      /* Normalization, needed for all-zero input. */
4452
4.83k
      assert (rn > 0);
4453
4.83k
      rn -= rp[rn-1] == 0;
4454
4.83k
    }
4455
4.83k
  assert (rn <= alloc);
4456
4.83k
  gmp_free (dp, sn);
4457
4458
4.83k
  r->_mp_size = sign ? - rn : rn;
4459
4460
4.83k
  return 0;
4461
4.83k
}
4462
4463
int
4464
mpz_init_set_str (mpz_t r, const char *sp, int base)
4465
0
{
4466
0
  mpz_init (r);
4467
0
  return mpz_set_str (r, sp, base);
4468
0
}
4469
4470
size_t
4471
mpz_out_str (FILE *stream, int base, const mpz_t x)
4472
0
{
4473
0
  char *str;
4474
0
  size_t len, n;
4475
4476
0
  str = mpz_get_str (NULL, base, x);
4477
0
  if (!str)
4478
0
    return 0;
4479
0
  len = strlen (str);
4480
0
  n = fwrite (str, 1, len, stream);
4481
0
  gmp_free (str, len + 1);
4482
0
  return n;
4483
0
}
4484
4485

4486
static int
4487
gmp_detect_endian (void)
4488
0
{
4489
0
  static const int i = 2;
4490
0
  const unsigned char *p = (const unsigned char *) &i;
4491
0
  return 1 - *p;
4492
0
}
4493
4494
/* Import and export. Does not support nails. */
4495
void
4496
mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
4497
      size_t nails, const void *src)
4498
0
{
4499
0
  const unsigned char *p;
4500
0
  ptrdiff_t word_step;
4501
0
  mp_ptr rp;
4502
0
  mp_size_t rn;
4503
4504
  /* The current (partial) limb. */
4505
0
  mp_limb_t limb;
4506
  /* The number of bytes already copied to this limb (starting from
4507
     the low end). */
4508
0
  size_t bytes;
4509
  /* The index where the limb should be stored, when completed. */
4510
0
  mp_size_t i;
4511
4512
0
  if (nails != 0)
4513
0
    gmp_die ("mpz_import: Nails not supported.");
4514
4515
0
  assert (order == 1 || order == -1);
4516
0
  assert (endian >= -1 && endian <= 1);
4517
4518
0
  if (endian == 0)
4519
0
    endian = gmp_detect_endian ();
4520
4521
0
  p = (unsigned char *) src;
4522
4523
0
  word_step = (order != endian) ? 2 * size : 0;
4524
4525
  /* Process bytes from the least significant end, so point p at the
4526
     least significant word. */
4527
0
  if (order == 1)
4528
0
    {
4529
0
      p += size * (count - 1);
4530
0
      word_step = - word_step;
4531
0
    }
4532
4533
  /* And at least significant byte of that word. */
4534
0
  if (endian == 1)
4535
0
    p += (size - 1);
4536
4537
0
  rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
4538
0
  rp = MPZ_REALLOC (r, rn);
4539
4540
0
  for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step)
4541
0
    {
4542
0
      size_t j;
4543
0
      for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4544
0
  {
4545
0
    limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT);
4546
0
    if (bytes == sizeof(mp_limb_t))
4547
0
      {
4548
0
        rp[i++] = limb;
4549
0
        bytes = 0;
4550
0
        limb = 0;
4551
0
      }
4552
0
  }
4553
0
    }
4554
0
  assert (i + (bytes > 0) == rn);
4555
0
  if (limb != 0)
4556
0
    rp[i++] = limb;
4557
0
  else
4558
0
    i = mpn_normalized_size (rp, i);
4559
4560
0
  r->_mp_size = i;
4561
0
}
4562
4563
void *
4564
mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
4565
      size_t nails, const mpz_t u)
4566
0
{
4567
0
  size_t count;
4568
0
  mp_size_t un;
4569
4570
0
  if (nails != 0)
4571
0
    gmp_die ("mpz_export: Nails not supported.");
4572
4573
0
  assert (order == 1 || order == -1);
4574
0
  assert (endian >= -1 && endian <= 1);
4575
0
  assert (size > 0 || u->_mp_size == 0);
4576
4577
0
  un = u->_mp_size;
4578
0
  count = 0;
4579
0
  if (un != 0)
4580
0
    {
4581
0
      size_t k;
4582
0
      unsigned char *p;
4583
0
      ptrdiff_t word_step;
4584
      /* The current (partial) limb. */
4585
0
      mp_limb_t limb;
4586
      /* The number of bytes left to do in this limb. */
4587
0
      size_t bytes;
4588
      /* The index where the limb was read. */
4589
0
      mp_size_t i;
4590
4591
0
      un = GMP_ABS (un);
4592
4593
      /* Count bytes in top limb. */
4594
0
      limb = u->_mp_d[un-1];
4595
0
      assert (limb != 0);
4596
4597
0
      k = (GMP_LIMB_BITS <= CHAR_BIT);
4598
0
      if (!k)
4599
0
  {
4600
0
    do {
4601
0
      int LOCAL_CHAR_BIT = CHAR_BIT;
4602
0
      k++; limb >>= LOCAL_CHAR_BIT;
4603
0
    } while (limb != 0);
4604
0
  }
4605
      /* else limb = 0; */
4606
4607
0
      count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
4608
4609
0
      if (!r)
4610
0
  r = gmp_alloc (count * size);
4611
4612
0
      if (endian == 0)
4613
0
  endian = gmp_detect_endian ();
4614
4615
0
      p = (unsigned char *) r;
4616
4617
0
      word_step = (order != endian) ? 2 * size : 0;
4618
4619
      /* Process bytes from the least significant end, so point p at the
4620
   least significant word. */
4621
0
      if (order == 1)
4622
0
  {
4623
0
    p += size * (count - 1);
4624
0
    word_step = - word_step;
4625
0
  }
4626
4627
      /* And at least significant byte of that word. */
4628
0
      if (endian == 1)
4629
0
  p += (size - 1);
4630
4631
0
      for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
4632
0
  {
4633
0
    size_t j;
4634
0
    for (j = 0; j < size; ++j, p -= (ptrdiff_t) endian)
4635
0
      {
4636
0
        if (sizeof (mp_limb_t) == 1)
4637
0
    {
4638
0
      if (i < un)
4639
0
        *p = u->_mp_d[i++];
4640
0
      else
4641
0
        *p = 0;
4642
0
    }
4643
0
        else
4644
0
    {
4645
0
      int LOCAL_CHAR_BIT = CHAR_BIT;
4646
0
      if (bytes == 0)
4647
0
        {
4648
0
          if (i < un)
4649
0
      limb = u->_mp_d[i++];
4650
0
          bytes = sizeof (mp_limb_t);
4651
0
        }
4652
0
      *p = limb;
4653
0
      limb >>= LOCAL_CHAR_BIT;
4654
0
      bytes--;
4655
0
    }
4656
0
      }
4657
0
  }
4658
0
      assert (i == un);
4659
0
      assert (k == count);
4660
0
    }
4661
4662
0
  if (countp)
4663
0
    *countp = count;
4664
4665
0
  return r;
4666
0
}