commit d5eab1233f12a705a3bd3aabffd562764b8eef7b parent 816523ac5ec43d452ef3e19b1d29fbdc9ff3ff5a Author: Cedric <cedric.zwahlen@students.bfh.ch> Date: Wed, 22 Nov 2023 00:30:54 +0100 Add montgomery multiplication Not stresstested, and does not yet work with even moduli Diffstat:
14 files changed, 5558 insertions(+), 40 deletions(-)
diff --git a/.DS_Store b/.DS_Store Binary files differ. diff --git a/source/.DS_Store b/source/.DS_Store Binary files differ. diff --git a/source/gmp.c b/source/gmp.c @@ -0,0 +1,4639 @@ +// +// gmp.c +// lib-gpu-verify +// +// Created by Cedric Zwahlen on 12.11.2023. +// + +#include "gmp.h" + + + +/* mini-gmp, a minimalistic implementation of a GNU GMP subset. + + Contributed to the GNU project by Niels Möller + Additional functionalities and improvements by Marco Bodrato. + +Copyright 1991-1997, 1999-2022 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +/* NOTE: All functions in this file which are not declared in + mini-gmp.h are internal, and are not intended to be compatible + with GMP or with future versions of mini-gmp. */ + +/* Much of the material copied from GMP files, including: gmp-impl.h, + longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c, + mpn/generic/lshift.c, mpn/generic/mul_1.c, + mpn/generic/mul_basecase.c, mpn/generic/rshift.c, + mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c, + mpn/generic/submul_1.c. */ + +#include <assert.h> +#include <ctype.h> +#include <limits.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +//#include "mini-gmp.h" + +#if !defined(MINI_GMP_DONT_USE_FLOAT_H) +#include <float.h> +#endif + + +/* Macros */ +#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT) + +#define GMP_LIMB_MAX ((mp_limb_t) ~ (mp_limb_t) 0) +#define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1)) + +#define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2)) +#define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1) + +#define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT) +#define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1)) + +#define GMP_ABS(x) ((x) >= 0 ? (x) : -(x)) +#define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1)) + +#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b)) +#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b)) + +#define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b))) + +#if defined(DBL_MANT_DIG) && FLT_RADIX == 2 +#define GMP_DBL_MANT_BITS DBL_MANT_DIG +#else +#define GMP_DBL_MANT_BITS (53) +#endif + +/* Return non-zero if xp,xsize and yp,ysize overlap. + If xp+xsize<=yp there's no overlap, or if yp+ysize<=xp there's no + overlap. If both these are false, there's an overlap. */ +#define GMP_MPN_OVERLAP_P(xp, xsize, yp, ysize) \ + ((xp) + (xsize) > (yp) && (yp) + (ysize) > (xp)) + +#define gmp_assert_nocarry(x) do { \ + mp_limb_t __cy = (x); \ + assert (__cy == 0); \ + (void) (__cy); \ + } while (0) + +#define gmp_clz(count, x) do { \ + mp_limb_t __clz_x = (x); \ + unsigned __clz_c = 0; \ + int LOCAL_SHIFT_BITS = 8; \ + if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS) \ + for (; \ + (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \ + __clz_c += 8) \ + { __clz_x <<= LOCAL_SHIFT_BITS; } \ + for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \ + __clz_x <<= 1; \ + (count) = __clz_c; \ + } while (0) + +#define gmp_ctz(count, x) do { \ + mp_limb_t __ctz_x = (x); \ + unsigned __ctz_c = 0; \ + gmp_clz (__ctz_c, __ctz_x & - __ctz_x); \ + (count) = GMP_LIMB_BITS - 1 - __ctz_c; \ + } while (0) + +#define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \ + do { \ + mp_limb_t __x; \ + __x = (al) + (bl); \ + (sh) = (ah) + (bh) + (__x < (al)); \ + (sl) = __x; \ + } while (0) + +#define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \ + do { \ + mp_limb_t __x; \ + __x = (al) - (bl); \ + (sh) = (ah) - (bh) - ((al) < (bl)); \ + (sl) = __x; \ + } while (0) + +#define gmp_umul_ppmm(w1, w0, u, v) \ + do { \ + int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS; \ + if (sizeof(unsigned int) * CHAR_BIT >= 2 * GMP_LIMB_BITS) \ + { \ + unsigned int __ww = (unsigned int) (u) * (v); \ + w0 = (mp_limb_t) __ww; \ + w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS); \ + } \ + else if (GMP_ULONG_BITS >= 2 * GMP_LIMB_BITS) \ + { \ + unsigned long int __ww = (unsigned long int) (u) * (v); \ + w0 = (mp_limb_t) __ww; \ + w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS); \ + } \ + else { \ + mp_limb_t __x0, __x1, __x2, __x3; \ + unsigned __ul, __vl, __uh, __vh; \ + mp_limb_t __u = (u), __v = (v); \ + assert (sizeof (unsigned) * 2 >= sizeof (mp_limb_t)); \ + \ + __ul = __u & GMP_LLIMB_MASK; \ + __uh = __u >> (GMP_LIMB_BITS / 2); \ + __vl = __v & GMP_LLIMB_MASK; \ + __vh = __v >> (GMP_LIMB_BITS / 2); \ + \ + __x0 = (mp_limb_t) __ul * __vl; \ + __x1 = (mp_limb_t) __ul * __vh; \ + __x2 = (mp_limb_t) __uh * __vl; \ + __x3 = (mp_limb_t) __uh * __vh; \ + \ + __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \ + __x1 += __x2; /* but this indeed can */ \ + if (__x1 < __x2) /* did we get it? */ \ + __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \ + \ + (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \ + (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \ + } \ + } while (0) + +/* If mp_limb_t is of size smaller than int, plain u*v implies + automatic promotion to *signed* int, and then multiply may overflow + and cause undefined behavior. Explicitly cast to unsigned int for + that case. */ +#define gmp_umullo_limb(u, v) \ + ((sizeof(mp_limb_t) >= sizeof(int)) ? (u)*(v) : (unsigned int)(u) * (v)) + +#define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \ + do { \ + mp_limb_t _qh, _ql, _r, _mask; \ + gmp_umul_ppmm (_qh, _ql, (nh), (di)); \ + gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \ + _r = (nl) - gmp_umullo_limb (_qh, (d)); \ + _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \ + _qh += _mask; \ + _r += _mask & (d); \ + if (_r >= (d)) \ + { \ + _r -= (d); \ + _qh++; \ + } \ + \ + (r) = _r; \ + (q) = _qh; \ + } while (0) + +#define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \ + do { \ + mp_limb_t _q0, _t1, _t0, _mask; \ + gmp_umul_ppmm ((q), _q0, (n2), (dinv)); \ + gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1)); \ + \ + /* Compute the two most significant limbs of n - q'd */ \ + (r1) = (n1) - gmp_umullo_limb ((d1), (q)); \ + gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \ + gmp_umul_ppmm (_t1, _t0, (d0), (q)); \ + gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0); \ + (q)++; \ + \ + /* Conditionally adjust q and the remainders */ \ + _mask = - (mp_limb_t) ((r1) >= _q0); \ + (q) += _mask; \ + gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \ + if ((r1) >= (d1)) \ + { \ + if ((r1) > (d1) || (r0) >= (d0)) \ + { \ + (q)++; \ + gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \ + } \ + } \ + } while (0) + +/* Swap macros. */ +#define MP_LIMB_T_SWAP(x, y) \ + do { \ + mp_limb_t __mp_limb_t_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_limb_t_swap__tmp; \ + } while (0) +#define MP_SIZE_T_SWAP(x, y) \ + do { \ + mp_size_t __mp_size_t_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_size_t_swap__tmp; \ + } while (0) +#define MP_BITCNT_T_SWAP(x,y) \ + do { \ + mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_bitcnt_t_swap__tmp; \ + } while (0) +#define MP_PTR_SWAP(x, y) \ + do { \ + mp_ptr __mp_ptr_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_ptr_swap__tmp; \ + } while (0) +#define MP_SRCPTR_SWAP(x, y) \ + do { \ + mp_srcptr __mp_srcptr_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mp_srcptr_swap__tmp; \ + } while (0) + +#define MPN_PTR_SWAP(xp,xs, yp,ys) \ + do { \ + MP_PTR_SWAP (xp, yp); \ + MP_SIZE_T_SWAP (xs, ys); \ + } while(0) +#define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \ + do { \ + MP_SRCPTR_SWAP (xp, yp); \ + MP_SIZE_T_SWAP (xs, ys); \ + } while(0) + +#define MPZ_PTR_SWAP(x, y) \ + do { \ + mpz_ptr __mpz_ptr_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mpz_ptr_swap__tmp; \ + } while (0) +#define MPZ_SRCPTR_SWAP(x, y) \ + do { \ + mpz_srcptr __mpz_srcptr_swap__tmp = (x); \ + (x) = (y); \ + (y) = __mpz_srcptr_swap__tmp; \ + } while (0) + +const int mp_bits_per_limb = GMP_LIMB_BITS; + + +/* Memory allocation and other helper functions. */ +static void +gmp_die (const char *msg) +{ + fprintf (stderr, "%s\n", msg); + abort(); +} + +static void * +gmp_default_alloc (size_t size) +{ + void *p; + + assert (size > 0); + + p = malloc (size); + if (!p) + gmp_die("gmp_default_alloc: Virtual memory exhausted."); + + return p; +} + +static void * +gmp_default_realloc (void *old, size_t unused_old_size, size_t new_size) +{ + void * p; + + p = realloc (old, new_size); + + if (!p) + gmp_die("gmp_default_realloc: Virtual memory exhausted."); + + return p; +} + +static void +gmp_default_free (void *p, size_t unused_size) +{ + free (p); +} + +static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc; +static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc; +static void (*gmp_free_func) (void *, size_t) = gmp_default_free; + +void +mp_get_memory_functions (void *(**alloc_func) (size_t), + void *(**realloc_func) (void *, size_t, size_t), + void (**free_func) (void *, size_t)) +{ + if (alloc_func) + *alloc_func = gmp_allocate_func; + + if (realloc_func) + *realloc_func = gmp_reallocate_func; + + if (free_func) + *free_func = gmp_free_func; +} + +void +mp_set_memory_functions (void *(*alloc_func) (size_t), + void *(*realloc_func) (void *, size_t, size_t), + void (*free_func) (void *, size_t)) +{ + if (!alloc_func) + alloc_func = gmp_default_alloc; + if (!realloc_func) + realloc_func = gmp_default_realloc; + if (!free_func) + free_func = gmp_default_free; + + gmp_allocate_func = alloc_func; + gmp_reallocate_func = realloc_func; + gmp_free_func = free_func; +} + +#define gmp_alloc(size) ((*gmp_allocate_func)((size))) +#define gmp_free(p, size) ((*gmp_free_func) ((p), (size))) +#define gmp_realloc(ptr, old_size, size) ((*gmp_reallocate_func)(ptr, old_size, size)) + +static mp_ptr +gmp_alloc_limbs (mp_size_t size) +{ + return (mp_ptr) gmp_alloc (size * sizeof (mp_limb_t)); +} + +static mp_ptr +gmp_realloc_limbs (mp_ptr old, mp_size_t old_size, mp_size_t size) +{ + assert (size > 0); + return (mp_ptr) gmp_realloc (old, old_size * sizeof (mp_limb_t), size * sizeof (mp_limb_t)); +} + +static void +gmp_free_limbs (mp_ptr old, mp_size_t size) +{ + gmp_free (old, size * sizeof (mp_limb_t)); +} + + +/* MPN interface */ + +void +mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n) +{ + mp_size_t i; + for (i = 0; i < n; i++) + d[i] = s[i]; +} + +void +mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n) +{ + while (--n >= 0) + d[n] = s[n]; +} + +int +mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + while (--n >= 0) + { + if (ap[n] != bp[n]) + return ap[n] > bp[n] ? 1 : -1; + } + return 0; +} + +static int +mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) +{ + if (an != bn) + return an < bn ? -1 : 1; + else + return mpn_cmp (ap, bp, an); +} + +static mp_size_t +mpn_normalized_size (mp_srcptr xp, mp_size_t n) +{ + while (n > 0 && xp[n-1] == 0) + --n; + return n; +} + +int +mpn_zero_p(mp_srcptr rp, mp_size_t n) +{ + return mpn_normalized_size (rp, n) == 0; +} + +void +mpn_zero (mp_ptr rp, mp_size_t n) +{ + while (--n >= 0) + rp[n] = 0; +} + +mp_limb_t +mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b) +{ + mp_size_t i; + + assert (n > 0); + i = 0; + do + { + mp_limb_t r = ap[i] + b; + /* Carry out */ + b = (r < b); + rp[i] = r; + } + while (++i < n); + + return b; +} + +mp_limb_t +mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mp_size_t i; + mp_limb_t cy; + + for (i = 0, cy = 0; i < n; i++) + { + mp_limb_t a, b, r; + a = ap[i]; b = bp[i]; + r = a + cy; + cy = (r < cy); + r += b; + cy += (r < b); + rp[i] = r; + } + return cy; +} + +mp_limb_t +mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) +{ + mp_limb_t cy; + + assert (an >= bn); + + cy = mpn_add_n (rp, ap, bp, bn); + if (an > bn) + cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy); + return cy; +} + +mp_limb_t +mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b) +{ + mp_size_t i; + + assert (n > 0); + + i = 0; + do + { + mp_limb_t a = ap[i]; + /* Carry out */ + mp_limb_t cy = a < b; + rp[i] = a - b; + b = cy; + } + while (++i < n); + + return b; +} + +mp_limb_t +mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mp_size_t i; + mp_limb_t cy; + + for (i = 0, cy = 0; i < n; i++) + { + mp_limb_t a, b; + a = ap[i]; b = bp[i]; + b += cy; + cy = (b < cy); + cy += (a < b); + rp[i] = a - b; + } + return cy; +} + +mp_limb_t +mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) +{ + mp_limb_t cy; + + assert (an >= bn); + + cy = mpn_sub_n (rp, ap, bp, bn); + if (an > bn) + cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy); + return cy; +} + +mp_limb_t +mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) +{ + mp_limb_t ul, cl, hpl, lpl; + + assert (n >= 1); + + cl = 0; + do + { + ul = *up++; + gmp_umul_ppmm (hpl, lpl, ul, vl); + + lpl += cl; + cl = (lpl < cl) + hpl; + + *rp++ = lpl; + } + while (--n != 0); + + return cl; +} + +mp_limb_t +mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) +{ + mp_limb_t ul, cl, hpl, lpl, rl; + + assert (n >= 1); + + cl = 0; + do + { + ul = *up++; + gmp_umul_ppmm (hpl, lpl, ul, vl); + + lpl += cl; + cl = (lpl < cl) + hpl; + + rl = *rp; + lpl = rl + lpl; + cl += lpl < rl; + *rp++ = lpl; + } + while (--n != 0); + + return cl; +} + +mp_limb_t +mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) +{ + mp_limb_t ul, cl, hpl, lpl, rl; + + assert (n >= 1); + + cl = 0; + do + { + ul = *up++; + gmp_umul_ppmm (hpl, lpl, ul, vl); + + lpl += cl; + cl = (lpl < cl) + hpl; + + rl = *rp; + lpl = rl - lpl; + cl += lpl > rl; + *rp++ = lpl; + } + while (--n != 0); + + return cl; +} + +mp_limb_t +mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn) +{ + assert (un >= vn); + assert (vn >= 1); + assert (!GMP_MPN_OVERLAP_P(rp, un + vn, up, un)); + assert (!GMP_MPN_OVERLAP_P(rp, un + vn, vp, vn)); + + /* We first multiply by the low order limb. This result can be + stored, not added, to rp. We also avoid a loop for zeroing this + way. */ + + rp[un] = mpn_mul_1 (rp, up, un, vp[0]); + + /* Now accumulate the product of up[] and the next higher limb from + vp[]. */ + + while (--vn >= 1) + { + rp += 1, vp += 1; + rp[un] = mpn_addmul_1 (rp, up, un, vp[0]); + } + return rp[un]; +} + +void +mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mpn_mul (rp, ap, n, bp, n); +} + +void +mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n) +{ + mpn_mul (rp, ap, n, ap, n); +} + +mp_limb_t +mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) +{ + mp_limb_t high_limb, low_limb; + unsigned int tnc; + mp_limb_t retval; + + assert (n >= 1); + assert (cnt >= 1); + assert (cnt < GMP_LIMB_BITS); + + up += n; + rp += n; + + tnc = GMP_LIMB_BITS - cnt; + low_limb = *--up; + retval = low_limb >> tnc; + high_limb = (low_limb << cnt); + + while (--n != 0) + { + low_limb = *--up; + *--rp = high_limb | (low_limb >> tnc); + high_limb = (low_limb << cnt); + } + *--rp = high_limb; + + return retval; +} + +mp_limb_t +mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) +{ + mp_limb_t high_limb, low_limb; + unsigned int tnc; + mp_limb_t retval; + + assert (n >= 1); + assert (cnt >= 1); + assert (cnt < GMP_LIMB_BITS); + + tnc = GMP_LIMB_BITS - cnt; + high_limb = *up++; + retval = (high_limb << tnc); + low_limb = high_limb >> cnt; + + while (--n != 0) + { + high_limb = *up++; + *rp++ = low_limb | (high_limb << tnc); + low_limb = high_limb >> cnt; + } + *rp = low_limb; + + return retval; +} + +static mp_bitcnt_t +mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un, + mp_limb_t ux) +{ + unsigned cnt; + + assert (ux == 0 || ux == GMP_LIMB_MAX); + assert (0 <= i && i <= un ); + + while (limb == 0) + { + i++; + if (i == un) + return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS); + limb = ux ^ up[i]; + } + gmp_ctz (cnt, limb); + return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt; +} + +mp_bitcnt_t +mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit) +{ + mp_size_t i; + i = bit / GMP_LIMB_BITS; + + return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)), + i, ptr, i, 0); +} + +mp_bitcnt_t +mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit) +{ + mp_size_t i; + i = bit / GMP_LIMB_BITS; + + return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)), + i, ptr, i, GMP_LIMB_MAX); +} + +void +mpn_com (mp_ptr rp, mp_srcptr up, mp_size_t n) +{ + while (--n >= 0) + *rp++ = ~ *up++; +} + +mp_limb_t +mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n) +{ + while (*up == 0) + { + *rp = 0; + if (!--n) + return 0; + ++up; ++rp; + } + *rp = - *up; + mpn_com (++rp, ++up, --n); + return 1; +} + + +/* MPN division interface. */ + +/* The 3/2 inverse is defined as + + m = floor( (B^3-1) / (B u1 + u0)) - B +*/ +mp_limb_t +mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) +{ + mp_limb_t r, m; + + { + mp_limb_t p, ql; + unsigned ul, uh, qh; + + assert (sizeof (unsigned) * 2 >= sizeof (mp_limb_t)); + /* For notation, let b denote the half-limb base, so that B = b^2. + Split u1 = b uh + ul. */ + ul = u1 & GMP_LLIMB_MASK; + uh = u1 >> (GMP_LIMB_BITS / 2); + + /* Approximation of the high half of quotient. Differs from the 2/1 + inverse of the half limb uh, since we have already subtracted + u0. */ + qh = (u1 ^ GMP_LIMB_MAX) / uh; + + /* Adjust to get a half-limb 3/2 inverse, i.e., we want + + qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u + = floor( (b (~u) + b-1) / u), + + and the remainder + + r = b (~u) + b-1 - qh (b uh + ul) + = b (~u - qh uh) + b-1 - qh ul + + Subtraction of qh ul may underflow, which implies adjustments. + But by normalization, 2 u >= B > qh ul, so we need to adjust by + at most 2. + */ + + r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK; + + p = (mp_limb_t) qh * ul; + /* Adjustment steps taken from udiv_qrnnd_c */ + if (r < p) + { + qh--; + r += u1; + if (r >= u1) /* i.e. we didn't get carry when adding to r */ + if (r < p) + { + qh--; + r += u1; + } + } + r -= p; + + /* Low half of the quotient is + + ql = floor ( (b r + b-1) / u1). + + This is a 3/2 division (on half-limbs), for which qh is a + suitable inverse. */ + + p = (r >> (GMP_LIMB_BITS / 2)) * qh + r; + /* Unlike full-limb 3/2, we can add 1 without overflow. For this to + work, it is essential that ql is a full mp_limb_t. */ + ql = (p >> (GMP_LIMB_BITS / 2)) + 1; + + /* By the 3/2 trick, we don't need the high half limb. */ + r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1; + + if (r >= (GMP_LIMB_MAX & (p << (GMP_LIMB_BITS / 2)))) + { + ql--; + r += u1; + } + m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql; + if (r >= u1) + { + m++; + r -= u1; + } + } + + /* Now m is the 2/1 inverse of u1. If u0 > 0, adjust it to become a + 3/2 inverse. */ + if (u0 > 0) + { + mp_limb_t th, tl; + r = ~r; + r += u0; + if (r < u0) + { + m--; + if (r >= u1) + { + m--; + r -= u1; + } + r -= u1; + } + gmp_umul_ppmm (th, tl, u0, m); + r += th; + if (r < th) + { + m--; + m -= ((r > u1) | ((r == u1) & (tl > u0))); + } + } + + return m; +} + +struct gmp_div_inverse +{ + /* Normalization shift count. */ + unsigned shift; + /* Normalized divisor (d0 unused for mpn_div_qr_1) */ + mp_limb_t d1, d0; + /* Inverse, for 2/1 or 3/2. */ + mp_limb_t di; +}; + +static void +mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d) +{ + unsigned shift; + + assert (d > 0); + gmp_clz (shift, d); + inv->shift = shift; + inv->d1 = d << shift; + inv->di = mpn_invert_limb (inv->d1); +} + +static void +mpn_div_qr_2_invert (struct gmp_div_inverse *inv, + mp_limb_t d1, mp_limb_t d0) +{ + unsigned shift; + + assert (d1 > 0); + gmp_clz (shift, d1); + inv->shift = shift; + if (shift > 0) + { + d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift)); + d0 <<= shift; + } + inv->d1 = d1; + inv->d0 = d0; + inv->di = mpn_invert_3by2 (d1, d0); +} + +static void +mpn_div_qr_invert (struct gmp_div_inverse *inv, + mp_srcptr dp, mp_size_t dn) +{ + assert (dn > 0); + + if (dn == 1) + mpn_div_qr_1_invert (inv, dp[0]); + else if (dn == 2) + mpn_div_qr_2_invert (inv, dp[1], dp[0]); + else + { + unsigned shift; + mp_limb_t d1, d0; + + d1 = dp[dn-1]; + d0 = dp[dn-2]; + assert (d1 > 0); + gmp_clz (shift, d1); + inv->shift = shift; + if (shift > 0) + { + d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift)); + d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift)); + } + inv->d1 = d1; + inv->d0 = d0; + inv->di = mpn_invert_3by2 (d1, d0); + } +} + +/* Not matching current public gmp interface, rather corresponding to + the sbpi1_div_* functions. */ +static mp_limb_t +mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn, + const struct gmp_div_inverse *inv) +{ + mp_limb_t d, di; + mp_limb_t r; + mp_ptr tp = NULL; + mp_size_t tn = 0; + + if (inv->shift > 0) + { + /* Shift, reusing qp area if possible. In-place shift if qp == np. */ + tp = qp; + if (!tp) + { + tn = nn; + tp = gmp_alloc_limbs (tn); + } + r = mpn_lshift (tp, np, nn, inv->shift); + np = tp; + } + else + r = 0; + + d = inv->d1; + di = inv->di; + while (--nn >= 0) + { + mp_limb_t q; + + gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di); + if (qp) + qp[nn] = q; + } + if (tn) + gmp_free_limbs (tp, tn); + + return r >> inv->shift; +} + +static void +mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn, + const struct gmp_div_inverse *inv) +{ + unsigned shift; + mp_size_t i; + mp_limb_t d1, d0, di, r1, r0; + + assert (nn >= 2); + shift = inv->shift; + d1 = inv->d1; + d0 = inv->d0; + di = inv->di; + + if (shift > 0) + r1 = mpn_lshift (np, np, nn, shift); + else + r1 = 0; + + r0 = np[nn - 1]; + + i = nn - 2; + do + { + mp_limb_t n0, q; + n0 = np[i]; + gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di); + + if (qp) + qp[i] = q; + } + while (--i >= 0); + + if (shift > 0) + { + assert ((r0 & (GMP_LIMB_MAX >> (GMP_LIMB_BITS - shift))) == 0); + r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift)); + r1 >>= shift; + } + + np[1] = r1; + np[0] = r0; +} + +static void +mpn_div_qr_pi1 (mp_ptr qp, + mp_ptr np, mp_size_t nn, mp_limb_t n1, + mp_srcptr dp, mp_size_t dn, + mp_limb_t dinv) +{ + mp_size_t i; + + mp_limb_t d1, d0; + mp_limb_t cy, cy1; + mp_limb_t q; + + assert (dn > 2); + assert (nn >= dn); + + d1 = dp[dn - 1]; + d0 = dp[dn - 2]; + + assert ((d1 & GMP_LIMB_HIGHBIT) != 0); + /* Iteration variable is the index of the q limb. + * + * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]> + * by <d1, d0, dp[dn-3], ..., dp[0] > + */ + + i = nn - dn; + do + { + mp_limb_t n0 = np[dn-1+i]; + + if (n1 == d1 && n0 == d0) + { + q = GMP_LIMB_MAX; + mpn_submul_1 (np+i, dp, dn, q); + n1 = np[dn-1+i]; /* update n1, last loop's value will now be invalid */ + } + else + { + gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv); + + cy = mpn_submul_1 (np + i, dp, dn-2, q); + + cy1 = n0 < cy; + n0 = n0 - cy; + cy = n1 < cy1; + n1 = n1 - cy1; + np[dn-2+i] = n0; + + if (cy != 0) + { + n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1); + q--; + } + } + + if (qp) + qp[i] = q; + } + while (--i >= 0); + + np[dn - 1] = n1; +} + +static void +mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn, + mp_srcptr dp, mp_size_t dn, + const struct gmp_div_inverse *inv) +{ + assert (dn > 0); + assert (nn >= dn); + + if (dn == 1) + np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv); + else if (dn == 2) + mpn_div_qr_2_preinv (qp, np, nn, inv); + else + { + mp_limb_t nh; + unsigned shift; + + assert (inv->d1 == dp[dn-1]); + assert (inv->d0 == dp[dn-2]); + assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0); + + shift = inv->shift; + if (shift > 0) + nh = mpn_lshift (np, np, nn, shift); + else + nh = 0; + + mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di); + + if (shift > 0) + gmp_assert_nocarry (mpn_rshift (np, np, dn, shift)); + } +} + +static void +mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn) +{ + struct gmp_div_inverse inv; + mp_ptr tp = NULL; + + assert (dn > 0); + assert (nn >= dn); + + mpn_div_qr_invert (&inv, dp, dn); + if (dn > 2 && inv.shift > 0) + { + tp = gmp_alloc_limbs (dn); + gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift)); + dp = tp; + } + mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv); + if (tp) + gmp_free_limbs (tp, dn); +} + + +/* MPN base conversion. */ +static unsigned +mpn_base_power_of_two_p (unsigned b) +{ + switch (b) + { + case 2: return 1; + case 4: return 2; + case 8: return 3; + case 16: return 4; + case 32: return 5; + case 64: return 6; + case 128: return 7; + case 256: return 8; + default: return 0; + } +} + +struct mpn_base_info +{ + /* bb is the largest power of the base which fits in one limb, and + exp is the corresponding exponent. */ + unsigned exp; + mp_limb_t bb; +}; + +static void +mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b) +{ + mp_limb_t m; + mp_limb_t p; + unsigned exp; + + m = GMP_LIMB_MAX / b; + for (exp = 1, p = b; p <= m; exp++) + p *= b; + + info->exp = exp; + info->bb = p; +} + +static mp_bitcnt_t +mpn_limb_size_in_base_2 (mp_limb_t u) +{ + unsigned shift; + + assert (u > 0); + gmp_clz (shift, u); + return GMP_LIMB_BITS - shift; +} + +static size_t +mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un) +{ + unsigned char mask; + size_t sn, j; + mp_size_t i; + unsigned shift; + + sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]) + + bits - 1) / bits; + + mask = (1U << bits) - 1; + + for (i = 0, j = sn, shift = 0; j-- > 0;) + { + unsigned char digit = up[i] >> shift; + + shift += bits; + + if (shift >= GMP_LIMB_BITS && ++i < un) + { + shift -= GMP_LIMB_BITS; + digit |= up[i] << (bits - shift); + } + sp[j] = digit & mask; + } + return sn; +} + +/* We generate digits from the least significant end, and reverse at + the end. */ +static size_t +mpn_limb_get_str (unsigned char *sp, mp_limb_t w, + const struct gmp_div_inverse *binv) +{ + mp_size_t i; + for (i = 0; w > 0; i++) + { + mp_limb_t h, l, r; + + h = w >> (GMP_LIMB_BITS - binv->shift); + l = w << binv->shift; + + gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di); + assert ((r & (GMP_LIMB_MAX >> (GMP_LIMB_BITS - binv->shift))) == 0); + r >>= binv->shift; + + sp[i] = r; + } + return i; +} + +static size_t +mpn_get_str_other (unsigned char *sp, + int base, const struct mpn_base_info *info, + mp_ptr up, mp_size_t un) +{ + struct gmp_div_inverse binv; + size_t sn; + size_t i; + + mpn_div_qr_1_invert (&binv, base); + + sn = 0; + + if (un > 1) + { + struct gmp_div_inverse bbinv; + mpn_div_qr_1_invert (&bbinv, info->bb); + + do + { + mp_limb_t w; + size_t done; + w = mpn_div_qr_1_preinv (up, up, un, &bbinv); + un -= (up[un-1] == 0); + done = mpn_limb_get_str (sp + sn, w, &binv); + + for (sn += done; done < info->exp; done++) + sp[sn++] = 0; + } + while (un > 1); + } + sn += mpn_limb_get_str (sp + sn, up[0], &binv); + + /* Reverse order */ + for (i = 0; 2*i + 1 < sn; i++) + { + unsigned char t = sp[i]; + sp[i] = sp[sn - i - 1]; + sp[sn - i - 1] = t; + } + + return sn; +} + +size_t +mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un) +{ + unsigned bits; + + assert (un > 0); + assert (up[un-1] > 0); + + bits = mpn_base_power_of_two_p (base); + if (bits) + return mpn_get_str_bits (sp, bits, up, un); + else + { + struct mpn_base_info info; + + mpn_get_base_info (&info, base); + return mpn_get_str_other (sp, base, &info, up, un); + } +} + +static mp_size_t +mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn, + unsigned bits) +{ + mp_size_t rn; + mp_limb_t limb; + unsigned shift; + + for (limb = 0, rn = 0, shift = 0; sn-- > 0; ) + { + limb |= (mp_limb_t) sp[sn] << shift; + shift += bits; + if (shift >= GMP_LIMB_BITS) + { + shift -= GMP_LIMB_BITS; + rp[rn++] = limb; + /* Next line is correct also if shift == 0, + bits == 8, and mp_limb_t == unsigned char. */ + limb = (unsigned int) sp[sn] >> (bits - shift); + } + } + if (limb != 0) + rp[rn++] = limb; + else + rn = mpn_normalized_size (rp, rn); + return rn; +} + +/* Result is usually normalized, except for all-zero input, in which + case a single zero limb is written at *RP, and 1 is returned. */ +static mp_size_t +mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn, + mp_limb_t b, const struct mpn_base_info *info) +{ + mp_size_t rn; + mp_limb_t w; + unsigned k; + size_t j; + + assert (sn > 0); + + k = 1 + (sn - 1) % info->exp; + + j = 0; + w = sp[j++]; + while (--k != 0) + w = w * b + sp[j++]; + + rp[0] = w; + + for (rn = 1; j < sn;) + { + mp_limb_t cy; + + w = sp[j++]; + for (k = 1; k < info->exp; k++) + w = w * b + sp[j++]; + + cy = mpn_mul_1 (rp, rp, rn, info->bb); + cy += mpn_add_1 (rp, rp, rn, w); + if (cy > 0) + rp[rn++] = cy; + } + assert (j == sn); + + return rn; +} + +mp_size_t +mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base) +{ + unsigned bits; + + if (sn == 0) + return 0; + + bits = mpn_base_power_of_two_p (base); + if (bits) + return mpn_set_str_bits (rp, sp, sn, bits); + else + { + struct mpn_base_info info; + + mpn_get_base_info (&info, base); + return mpn_set_str_other (rp, sp, sn, base, &info); + } +} + + +/* MPZ interface */ +void +mpz_init (mpz_t r) +{ + static const mp_limb_t dummy_limb = GMP_LIMB_MAX & 0xc1a0; + + r->_mp_alloc = 0; + r->_mp_size = 0; + r->_mp_d = (mp_ptr) &dummy_limb; +} + +/* The utility of this function is a bit limited, since many functions + assigns the result variable using mpz_swap. */ +void +mpz_init2 (mpz_t r, mp_bitcnt_t bits) +{ + mp_size_t rn; + + bits -= (bits != 0); /* Round down, except if 0 */ + rn = 1 + bits / GMP_LIMB_BITS; + + r->_mp_alloc = rn; + r->_mp_size = 0; + r->_mp_d = gmp_alloc_limbs (rn); +} + +void +mpz_clear (mpz_t r) +{ + if (r->_mp_alloc) + gmp_free_limbs (r->_mp_d, r->_mp_alloc); +} + +static mp_ptr +mpz_realloc (mpz_t r, mp_size_t size) +{ + size = GMP_MAX (size, 1); + + if (r->_mp_alloc) + r->_mp_d = gmp_realloc_limbs (r->_mp_d, r->_mp_alloc, size); + else + r->_mp_d = gmp_alloc_limbs (size); + r->_mp_alloc = size; + + if (GMP_ABS (r->_mp_size) > size) + r->_mp_size = 0; + + return r->_mp_d; +} + +/* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */ +#define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \ + ? mpz_realloc(z,n) \ + : (z)->_mp_d) + +/* MPZ assignment and basic conversions. */ +void +mpz_set_si (mpz_t r, signed long int x) +{ + if (x >= 0) + mpz_set_ui (r, x); + else /* (x < 0) */ + if (GMP_LIMB_BITS < GMP_ULONG_BITS) + { + mpz_set_ui (r, GMP_NEG_CAST (unsigned long int, x)); + mpz_neg (r, r); + } + else + { + r->_mp_size = -1; + MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x); + } +} + +void +mpz_set_ui (mpz_t r, unsigned long int x) +{ + if (x > 0) + { + r->_mp_size = 1; + MPZ_REALLOC (r, 1)[0] = x; + if (GMP_LIMB_BITS < GMP_ULONG_BITS) + { + int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS; + while (x >>= LOCAL_GMP_LIMB_BITS) + { + ++ r->_mp_size; + MPZ_REALLOC (r, r->_mp_size)[r->_mp_size - 1] = x; + } + } + } + else + r->_mp_size = 0; +} + +void +mpz_set (mpz_t r, const mpz_t x) +{ + /* Allow the NOP r == x */ + if (r != x) + { + mp_size_t n; + mp_ptr rp; + + n = GMP_ABS (x->_mp_size); + rp = MPZ_REALLOC (r, n); + + mpn_copyi (rp, x->_mp_d, n); + r->_mp_size = x->_mp_size; + } +} + +void +mpz_init_set_si (mpz_t r, signed long int x) +{ + mpz_init (r); + mpz_set_si (r, x); +} + +void +mpz_init_set_ui (mpz_t r, unsigned long int x) +{ + mpz_init (r); + mpz_set_ui (r, x); +} + +void +mpz_init_set (mpz_t r, const mpz_t x) +{ + mpz_init (r); + mpz_set (r, x); +} + +int +mpz_fits_slong_p (const mpz_t u) +{ + return mpz_cmp_si (u, LONG_MAX) <= 0 && mpz_cmp_si (u, LONG_MIN) >= 0; +} + +static int +mpn_absfits_ulong_p (mp_srcptr up, mp_size_t un) +{ + int ulongsize = GMP_ULONG_BITS / GMP_LIMB_BITS; + mp_limb_t ulongrem = 0; + + if (GMP_ULONG_BITS % GMP_LIMB_BITS != 0) + ulongrem = (mp_limb_t) (ULONG_MAX >> GMP_LIMB_BITS * ulongsize) + 1; + + return un <= ulongsize || (up[ulongsize] < ulongrem && un == ulongsize + 1); +} + +int +mpz_fits_ulong_p (const mpz_t u) +{ + mp_size_t us = u->_mp_size; + + return us >= 0 && mpn_absfits_ulong_p (u->_mp_d, us); +} + +int +mpz_fits_sint_p (const mpz_t u) +{ + return mpz_cmp_si (u, INT_MAX) <= 0 && mpz_cmp_si (u, INT_MIN) >= 0; +} + +int +mpz_fits_uint_p (const mpz_t u) +{ + return u->_mp_size >= 0 && mpz_cmpabs_ui (u, UINT_MAX) <= 0; +} + +int +mpz_fits_sshort_p (const mpz_t u) +{ + return mpz_cmp_si (u, SHRT_MAX) <= 0 && mpz_cmp_si (u, SHRT_MIN) >= 0; +} + +int +mpz_fits_ushort_p (const mpz_t u) +{ + return u->_mp_size >= 0 && mpz_cmpabs_ui (u, USHRT_MAX) <= 0; +} + +long int +mpz_get_si (const mpz_t u) +{ + unsigned long r = mpz_get_ui (u); + unsigned long c = -LONG_MAX - LONG_MIN; + + if (u->_mp_size < 0) + /* This expression is necessary to properly handle -LONG_MIN */ + return -(long) c - (long) ((r - c) & LONG_MAX); + else + return (long) (r & LONG_MAX); +} + +unsigned long int +mpz_get_ui (const mpz_t u) +{ + if (GMP_LIMB_BITS < GMP_ULONG_BITS) + { + int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS; + unsigned long r = 0; + mp_size_t n = GMP_ABS (u->_mp_size); + n = GMP_MIN (n, 1 + (mp_size_t) (GMP_ULONG_BITS - 1) / GMP_LIMB_BITS); + while (--n >= 0) + r = (r << LOCAL_GMP_LIMB_BITS) + u->_mp_d[n]; + return r; + } + + return u->_mp_size == 0 ? 0 : u->_mp_d[0]; +} + +size_t +mpz_size (const mpz_t u) +{ + return GMP_ABS (u->_mp_size); +} + +mp_limb_t +mpz_getlimbn (const mpz_t u, mp_size_t n) +{ + if (n >= 0 && n < GMP_ABS (u->_mp_size)) + return u->_mp_d[n]; + else + return 0; +} + +void +mpz_realloc2 (mpz_t x, mp_bitcnt_t n) +{ + mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS); +} + +mp_srcptr +mpz_limbs_read (mpz_srcptr x) +{ + return x->_mp_d; +} + +mp_ptr +mpz_limbs_modify (mpz_t x, mp_size_t n) +{ + assert (n > 0); + return MPZ_REALLOC (x, n); +} + +mp_ptr +mpz_limbs_write (mpz_t x, mp_size_t n) +{ + return mpz_limbs_modify (x, n); +} + +void +mpz_limbs_finish (mpz_t x, mp_size_t xs) +{ + mp_size_t xn; + xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs)); + x->_mp_size = xs < 0 ? -xn : xn; +} + +static mpz_srcptr +mpz_roinit_normal_n (mpz_t x, mp_srcptr xp, mp_size_t xs) +{ + x->_mp_alloc = 0; + x->_mp_d = (mp_ptr) xp; + x->_mp_size = xs; + return x; +} + +mpz_srcptr +mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs) +{ + mpz_roinit_normal_n (x, xp, xs); + mpz_limbs_finish (x, xs); + return x; +} + + +/* Conversions and comparison to double. */ +void +mpz_set_d (mpz_t r, double x) +{ + int sign; + mp_ptr rp; + mp_size_t rn, i; + double B; + double Bi; + mp_limb_t f; + + /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is + zero or infinity. */ + if (x != x || x == x * 0.5) + { + r->_mp_size = 0; + return; + } + + sign = x < 0.0 ; + if (sign) + x = - x; + + if (x < 1.0) + { + r->_mp_size = 0; + return; + } + B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1); + Bi = 1.0 / B; + for (rn = 1; x >= B; rn++) + x *= Bi; + + rp = MPZ_REALLOC (r, rn); + + f = (mp_limb_t) x; + x -= f; + assert (x < 1.0); + i = rn-1; + rp[i] = f; + while (--i >= 0) + { + x = B * x; + f = (mp_limb_t) x; + x -= f; + assert (x < 1.0); + rp[i] = f; + } + + r->_mp_size = sign ? - rn : rn; +} + +void +mpz_init_set_d (mpz_t r, double x) +{ + mpz_init (r); + mpz_set_d (r, x); +} + +double +mpz_get_d (const mpz_t u) +{ + int m; + mp_limb_t l; + mp_size_t un; + double x; + double B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1); + + un = GMP_ABS (u->_mp_size); + + if (un == 0) + return 0.0; + + l = u->_mp_d[--un]; + gmp_clz (m, l); + m = m + GMP_DBL_MANT_BITS - GMP_LIMB_BITS; + if (m < 0) + l &= GMP_LIMB_MAX << -m; + + for (x = l; --un >= 0;) + { + x = B*x; + if (m > 0) { + l = u->_mp_d[un]; + m -= GMP_LIMB_BITS; + if (m < 0) + l &= GMP_LIMB_MAX << -m; + x += l; + } + } + + if (u->_mp_size < 0) + x = -x; + + return x; +} + +int +mpz_cmpabs_d (const mpz_t x, double d) +{ + mp_size_t xn; + double B, Bi; + mp_size_t i; + + xn = x->_mp_size; + d = GMP_ABS (d); + + if (xn != 0) + { + xn = GMP_ABS (xn); + + B = 4.0 * (double) (GMP_LIMB_HIGHBIT >> 1); + Bi = 1.0 / B; + + /* Scale d so it can be compared with the top limb. */ + for (i = 1; i < xn; i++) + d *= Bi; + + if (d >= B) + return -1; + + /* Compare floor(d) to top limb, subtract and cancel when equal. */ + for (i = xn; i-- > 0;) + { + mp_limb_t f, xl; + + f = (mp_limb_t) d; + xl = x->_mp_d[i]; + if (xl > f) + return 1; + else if (xl < f) + return -1; + d = B * (d - f); + } + } + return - (d > 0.0); +} + +int +mpz_cmp_d (const mpz_t x, double d) +{ + if (x->_mp_size < 0) + { + if (d >= 0.0) + return -1; + else + return -mpz_cmpabs_d (x, d); + } + else + { + if (d < 0.0) + return 1; + else + return mpz_cmpabs_d (x, d); + } +} + + +/* MPZ comparisons and the like. */ +int +mpz_sgn (const mpz_t u) +{ + return GMP_CMP (u->_mp_size, 0); +} + +int +mpz_cmp_si (const mpz_t u, long v) +{ + mp_size_t usize = u->_mp_size; + + if (v >= 0) + return mpz_cmp_ui (u, v); + else if (usize >= 0) + return 1; + else + return - mpz_cmpabs_ui (u, GMP_NEG_CAST (unsigned long int, v)); +} + +int +mpz_cmp_ui (const mpz_t u, unsigned long v) +{ + mp_size_t usize = u->_mp_size; + + if (usize < 0) + return -1; + else + return mpz_cmpabs_ui (u, v); +} + +int +mpz_cmp (const mpz_t a, const mpz_t b) +{ + mp_size_t asize = a->_mp_size; + mp_size_t bsize = b->_mp_size; + + if (asize != bsize) + return (asize < bsize) ? -1 : 1; + else if (asize >= 0) + return mpn_cmp (a->_mp_d, b->_mp_d, asize); + else + return mpn_cmp (b->_mp_d, a->_mp_d, -asize); +} + +int +mpz_cmpabs_ui (const mpz_t u, unsigned long v) +{ + mp_size_t un = GMP_ABS (u->_mp_size); + + if (! mpn_absfits_ulong_p (u->_mp_d, un)) + return 1; + else + { + unsigned long uu = mpz_get_ui (u); + return GMP_CMP(uu, v); + } +} + +int +mpz_cmpabs (const mpz_t u, const mpz_t v) +{ + return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size), + v->_mp_d, GMP_ABS (v->_mp_size)); +} + +void +mpz_abs (mpz_t r, const mpz_t u) +{ + mpz_set (r, u); + r->_mp_size = GMP_ABS (r->_mp_size); +} + +void +mpz_neg (mpz_t r, const mpz_t u) +{ + mpz_set (r, u); + r->_mp_size = -r->_mp_size; +} + +void +mpz_swap (mpz_t u, mpz_t v) +{ + MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc); + MPN_PTR_SWAP (u->_mp_d, u->_mp_size, v->_mp_d, v->_mp_size); +} + + +/* MPZ addition and subtraction */ + + +void +mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + mpz_t bb; + mpz_init_set_ui (bb, b); + mpz_add (r, a, bb); + mpz_clear (bb); +} + +void +mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + mpz_ui_sub (r, b, a); + mpz_neg (r, r); +} + +void +mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b) +{ + mpz_neg (r, b); + mpz_add_ui (r, r, a); +} + +static mp_size_t +mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t an = GMP_ABS (a->_mp_size); + mp_size_t bn = GMP_ABS (b->_mp_size); + mp_ptr rp; + mp_limb_t cy; + + if (an < bn) + { + MPZ_SRCPTR_SWAP (a, b); + MP_SIZE_T_SWAP (an, bn); + } + + rp = MPZ_REALLOC (r, an + 1); + cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn); + + rp[an] = cy; + + return an + cy; +} + +static mp_size_t +mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t an = GMP_ABS (a->_mp_size); + mp_size_t bn = GMP_ABS (b->_mp_size); + int cmp; + mp_ptr rp; + + cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn); + if (cmp > 0) + { + rp = MPZ_REALLOC (r, an); + gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn)); + return mpn_normalized_size (rp, an); + } + else if (cmp < 0) + { + rp = MPZ_REALLOC (r, bn); + gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an)); + return -mpn_normalized_size (rp, bn); + } + else + return 0; +} + +void +mpz_add (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t rn; + + if ( (a->_mp_size ^ b->_mp_size) >= 0) + rn = mpz_abs_add (r, a, b); + else + rn = mpz_abs_sub (r, a, b); + + r->_mp_size = a->_mp_size >= 0 ? rn : - rn; +} + +void +mpz_sub (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t rn; + + if ( (a->_mp_size ^ b->_mp_size) >= 0) + rn = mpz_abs_sub (r, a, b); + else + rn = mpz_abs_add (r, a, b); + + r->_mp_size = a->_mp_size >= 0 ? rn : - rn; +} + + +/* MPZ multiplication */ +void +mpz_mul_si (mpz_t r, const mpz_t u, long int v) +{ + if (v < 0) + { + mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v)); + mpz_neg (r, r); + } + else + mpz_mul_ui (r, u, v); +} + +void +mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v) +{ + mpz_t vv; + mpz_init_set_ui (vv, v); + mpz_mul (r, u, vv); + mpz_clear (vv); + return; +} + +void +mpz_mul (mpz_t r, const mpz_t u, const mpz_t v) +{ + int sign; + mp_size_t un, vn, rn; + mpz_t t; + mp_ptr tp; + + un = u->_mp_size; + vn = v->_mp_size; + + if (un == 0 || vn == 0) + { + r->_mp_size = 0; + return; + } + + sign = (un ^ vn) < 0; + + un = GMP_ABS (un); + vn = GMP_ABS (vn); + + mpz_init2 (t, (un + vn) * GMP_LIMB_BITS); + + tp = t->_mp_d; + if (un >= vn) + mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn); + else + mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un); + + rn = un + vn; + rn -= tp[rn-1] == 0; + + t->_mp_size = sign ? - rn : rn; + mpz_swap (r, t); + mpz_clear (t); +} + +void +mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits) +{ + mp_size_t un, rn; + mp_size_t limbs; + unsigned shift; + mp_ptr rp; + + un = GMP_ABS (u->_mp_size); + if (un == 0) + { + r->_mp_size = 0; + return; + } + + limbs = bits / GMP_LIMB_BITS; + shift = bits % GMP_LIMB_BITS; + + rn = un + limbs + (shift > 0); + rp = MPZ_REALLOC (r, rn); + if (shift > 0) + { + mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift); + rp[rn-1] = cy; + rn -= (cy == 0); + } + else + mpn_copyd (rp + limbs, u->_mp_d, un); + + mpn_zero (rp, limbs); + + r->_mp_size = (u->_mp_size < 0) ? - rn : rn; +} + +void +mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v) +{ + mpz_t t; + mpz_init_set_ui (t, v); + mpz_mul (t, u, t); + mpz_add (r, r, t); + mpz_clear (t); +} + +void +mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v) +{ + mpz_t t; + mpz_init_set_ui (t, v); + mpz_mul (t, u, t); + mpz_sub (r, r, t); + mpz_clear (t); +} + +void +mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v) +{ + mpz_t t; + mpz_init (t); + mpz_mul (t, u, v); + mpz_add (r, r, t); + mpz_clear (t); +} + +void +mpz_submul (mpz_t r, const mpz_t u, const mpz_t v) +{ + mpz_t t; + mpz_init (t); + mpz_mul (t, u, v); + mpz_sub (r, r, t); + mpz_clear (t); +} + + +/* MPZ division */ +enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC }; + +/* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */ +static int +mpz_div_qr (mpz_t q, mpz_t r, + const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode) +{ + mp_size_t ns, ds, nn, dn, qs; + ns = n->_mp_size; + ds = d->_mp_size; + + if (ds == 0) + gmp_die("mpz_div_qr: Divide by zero."); + + if (ns == 0) + { + if (q) + q->_mp_size = 0; + if (r) + r->_mp_size = 0; + return 0; + } + + nn = GMP_ABS (ns); + dn = GMP_ABS (ds); + + qs = ds ^ ns; + + if (nn < dn) + { + if (mode == GMP_DIV_CEIL && qs >= 0) + { + /* q = 1, r = n - d */ + if (r) + mpz_sub (r, n, d); + if (q) + mpz_set_ui (q, 1); + } + else if (mode == GMP_DIV_FLOOR && qs < 0) + { + /* q = -1, r = n + d */ + if (r) + mpz_add (r, n, d); + if (q) + mpz_set_si (q, -1); + } + else + { + /* q = 0, r = d */ + if (r) + mpz_set (r, n); + if (q) + q->_mp_size = 0; + } + return 1; + } + else + { + mp_ptr np, qp; + mp_size_t qn, rn; + mpz_t tq, tr; + + mpz_init_set (tr, n); + np = tr->_mp_d; + + qn = nn - dn + 1; + + if (q) + { + mpz_init2 (tq, qn * GMP_LIMB_BITS); + qp = tq->_mp_d; + } + else + qp = NULL; + + mpn_div_qr (qp, np, nn, d->_mp_d, dn); + + if (qp) + { + qn -= (qp[qn-1] == 0); + + tq->_mp_size = qs < 0 ? -qn : qn; + } + rn = mpn_normalized_size (np, dn); + tr->_mp_size = ns < 0 ? - rn : rn; + + if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0) + { + if (q) + mpz_sub_ui (tq, tq, 1); + if (r) + mpz_add (tr, tr, d); + } + else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0) + { + if (q) + mpz_add_ui (tq, tq, 1); + if (r) + mpz_sub (tr, tr, d); + } + + if (q) + { + mpz_swap (tq, q); + mpz_clear (tq); + } + if (r) + mpz_swap (tr, r); + + mpz_clear (tr); + + return rn != 0; + } +} + +void +mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, r, n, d, GMP_DIV_CEIL); +} + +void +mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC); +} + +void +mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL); +} + +void +mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC); +} + +void +mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL); +} + +void +mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC); +} + +void +mpz_mod (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL); +} + +static void +mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index, + enum mpz_div_round_mode mode) +{ + mp_size_t un, qn; + mp_size_t limb_cnt; + mp_ptr qp; + int adjust; + + un = u->_mp_size; + if (un == 0) + { + q->_mp_size = 0; + return; + } + limb_cnt = bit_index / GMP_LIMB_BITS; + qn = GMP_ABS (un) - limb_cnt; + bit_index %= GMP_LIMB_BITS; + + if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */ + /* Note: Below, the final indexing at limb_cnt is valid because at + that point we have qn > 0. */ + adjust = (qn <= 0 + || !mpn_zero_p (u->_mp_d, limb_cnt) + || (u->_mp_d[limb_cnt] + & (((mp_limb_t) 1 << bit_index) - 1))); + else + adjust = 0; + + if (qn <= 0) + qn = 0; + else + { + qp = MPZ_REALLOC (q, qn); + + if (bit_index != 0) + { + mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index); + qn -= qp[qn - 1] == 0; + } + else + { + mpn_copyi (qp, u->_mp_d + limb_cnt, qn); + } + } + + q->_mp_size = qn; + + if (adjust) + mpz_add_ui (q, q, 1); + if (un < 0) + mpz_neg (q, q); +} + +static void +mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index, + enum mpz_div_round_mode mode) +{ + mp_size_t us, un, rn; + mp_ptr rp; + mp_limb_t mask; + + us = u->_mp_size; + if (us == 0 || bit_index == 0) + { + r->_mp_size = 0; + return; + } + rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; + assert (rn > 0); + + rp = MPZ_REALLOC (r, rn); + un = GMP_ABS (us); + + mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index); + + if (rn > un) + { + /* Quotient (with truncation) is zero, and remainder is + non-zero */ + if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */ + { + /* Have to negate and sign extend. */ + mp_size_t i; + + gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un)); + for (i = un; i < rn - 1; i++) + rp[i] = GMP_LIMB_MAX; + + rp[rn-1] = mask; + us = -us; + } + else + { + /* Just copy */ + if (r != u) + mpn_copyi (rp, u->_mp_d, un); + + rn = un; + } + } + else + { + if (r != u) + mpn_copyi (rp, u->_mp_d, rn - 1); + + rp[rn-1] = u->_mp_d[rn-1] & mask; + + if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */ + { + /* If r != 0, compute 2^{bit_count} - r. */ + mpn_neg (rp, rp, rn); + + rp[rn-1] &= mask; + + /* us is not used for anything else, so we can modify it + here to indicate flipped sign. */ + us = -us; + } + } + rn = mpn_normalized_size (rp, rn); + r->_mp_size = us < 0 ? -rn : rn; +} + +void +mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL); +} + +void +mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC); +} + +void +mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL); +} + +void +mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC); +} + +void +mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d) +{ + gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC)); +} + +int +mpz_divisible_p (const mpz_t n, const mpz_t d) +{ + return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0; +} + +int +mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m) +{ + mpz_t t; + int res; + + /* a == b (mod 0) iff a == b */ + if (mpz_sgn (m) == 0) + return (mpz_cmp (a, b) == 0); + + mpz_init (t); + mpz_sub (t, a, b); + res = mpz_divisible_p (t, m); + mpz_clear (t); + + return res; +} + +static unsigned long +mpz_div_qr_ui (mpz_t q, mpz_t r, + const mpz_t n, unsigned long d, enum mpz_div_round_mode mode) +{ + unsigned long ret; + mpz_t rr, dd; + + mpz_init (rr); + mpz_init_set_ui (dd, d); + mpz_div_qr (q, rr, n, dd, mode); + mpz_clear (dd); + ret = mpz_get_ui (rr); + + if (r) + mpz_swap (r, rr); + mpz_clear (rr); + + return ret; +} + +unsigned long +mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL); +} + +unsigned long +mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR); +} + +unsigned long +mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC); +} + +unsigned long +mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL); +} + +unsigned long +mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR); +} + +unsigned long +mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC); +} + +unsigned long +mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL); +} +unsigned long +mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR); +} +unsigned long +mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC); +} + +unsigned long +mpz_cdiv_ui (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL); +} + +unsigned long +mpz_fdiv_ui (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR); +} + +unsigned long +mpz_tdiv_ui (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC); +} + +unsigned long +mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR); +} + +void +mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC)); +} + +int +mpz_divisible_ui_p (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0; +} + + +/* GCD */ +static mp_limb_t +mpn_gcd_11 (mp_limb_t u, mp_limb_t v) +{ + unsigned shift; + + assert ( (u | v) > 0); + + if (u == 0) + return v; + else if (v == 0) + return u; + + gmp_ctz (shift, u | v); + + u >>= shift; + v >>= shift; + + if ( (u & 1) == 0) + MP_LIMB_T_SWAP (u, v); + + while ( (v & 1) == 0) + v >>= 1; + + while (u != v) + { + if (u > v) + { + u -= v; + do + u >>= 1; + while ( (u & 1) == 0); + } + else + { + v -= u; + do + v >>= 1; + while ( (v & 1) == 0); + } + } + return u << shift; +} + +unsigned long +mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v) +{ + mpz_t t; + mpz_init_set_ui(t, v); + mpz_gcd (t, u, t); + if (v > 0) + v = mpz_get_ui (t); + + if (g) + mpz_swap (t, g); + + mpz_clear (t); + + return v; +} + +static mp_bitcnt_t +mpz_make_odd (mpz_t r) +{ + mp_bitcnt_t shift; + + assert (r->_mp_size > 0); + /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */ + shift = mpn_scan1 (r->_mp_d, 0); + mpz_tdiv_q_2exp (r, r, shift); + + return shift; +} + +void +mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v) +{ + mpz_t tu, tv; + mp_bitcnt_t uz, vz, gz; + + if (u->_mp_size == 0) + { + mpz_abs (g, v); + return; + } + if (v->_mp_size == 0) + { + mpz_abs (g, u); + return; + } + + mpz_init (tu); + mpz_init (tv); + + mpz_abs (tu, u); + uz = mpz_make_odd (tu); + mpz_abs (tv, v); + vz = mpz_make_odd (tv); + gz = GMP_MIN (uz, vz); + + if (tu->_mp_size < tv->_mp_size) + mpz_swap (tu, tv); + + mpz_tdiv_r (tu, tu, tv); + if (tu->_mp_size == 0) + { + mpz_swap (g, tv); + } + else + for (;;) + { + int c; + + mpz_make_odd (tu); + c = mpz_cmp (tu, tv); + if (c == 0) + { + mpz_swap (g, tu); + break; + } + if (c < 0) + mpz_swap (tu, tv); + + if (tv->_mp_size == 1) + { + mp_limb_t *gp; + + mpz_tdiv_r (tu, tu, tv); + gp = MPZ_REALLOC (g, 1); /* gp = mpz_limbs_modify (g, 1); */ + *gp = mpn_gcd_11 (tu->_mp_d[0], tv->_mp_d[0]); + + g->_mp_size = *gp != 0; /* mpz_limbs_finish (g, 1); */ + break; + } + mpz_sub (tu, tu, tv); + } + mpz_clear (tu); + mpz_clear (tv); + mpz_mul_2exp (g, g, gz); +} + +void +mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v) +{ + mpz_t tu, tv, s0, s1, t0, t1; + mp_bitcnt_t uz, vz, gz; + mp_bitcnt_t power; + + if (u->_mp_size == 0) + { + /* g = 0 u + sgn(v) v */ + signed long sign = mpz_sgn (v); + mpz_abs (g, v); + if (s) + s->_mp_size = 0; + if (t) + mpz_set_si (t, sign); + return; + } + + if (v->_mp_size == 0) + { + /* g = sgn(u) u + 0 v */ + signed long sign = mpz_sgn (u); + mpz_abs (g, u); + if (s) + mpz_set_si (s, sign); + if (t) + t->_mp_size = 0; + return; + } + + mpz_init (tu); + mpz_init (tv); + mpz_init (s0); + mpz_init (s1); + mpz_init (t0); + mpz_init (t1); + + mpz_abs (tu, u); + uz = mpz_make_odd (tu); + mpz_abs (tv, v); + vz = mpz_make_odd (tv); + gz = GMP_MIN (uz, vz); + + uz -= gz; + vz -= gz; + + /* Cofactors corresponding to odd gcd. gz handled later. */ + if (tu->_mp_size < tv->_mp_size) + { + mpz_swap (tu, tv); + MPZ_SRCPTR_SWAP (u, v); + MPZ_PTR_SWAP (s, t); + MP_BITCNT_T_SWAP (uz, vz); + } + + /* Maintain + * + * u = t0 tu + t1 tv + * v = s0 tu + s1 tv + * + * where u and v denote the inputs with common factors of two + * eliminated, and det (s0, t0; s1, t1) = 2^p. Then + * + * 2^p tu = s1 u - t1 v + * 2^p tv = -s0 u + t0 v + */ + + /* After initial division, tu = q tv + tu', we have + * + * u = 2^uz (tu' + q tv) + * v = 2^vz tv + * + * or + * + * t0 = 2^uz, t1 = 2^uz q + * s0 = 0, s1 = 2^vz + */ + + mpz_tdiv_qr (t1, tu, tu, tv); + mpz_mul_2exp (t1, t1, uz); + + mpz_setbit (s1, vz); + power = uz + vz; + + if (tu->_mp_size > 0) + { + mp_bitcnt_t shift; + shift = mpz_make_odd (tu); + mpz_setbit (t0, uz + shift); + power += shift; + + for (;;) + { + int c; + c = mpz_cmp (tu, tv); + if (c == 0) + break; + + if (c < 0) + { + /* tv = tv' + tu + * + * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv' + * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */ + + mpz_sub (tv, tv, tu); + mpz_add (t0, t0, t1); + mpz_add (s0, s0, s1); + + shift = mpz_make_odd (tv); + mpz_mul_2exp (t1, t1, shift); + mpz_mul_2exp (s1, s1, shift); + } + else + { + mpz_sub (tu, tu, tv); + mpz_add (t1, t0, t1); + mpz_add (s1, s0, s1); + + shift = mpz_make_odd (tu); + mpz_mul_2exp (t0, t0, shift); + mpz_mul_2exp (s0, s0, shift); + } + power += shift; + } + } + else + mpz_setbit (t0, uz); + + /* Now tv = odd part of gcd, and -s0 and t0 are corresponding + cofactors. */ + + mpz_mul_2exp (tv, tv, gz); + mpz_neg (s0, s0); + + /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To + adjust cofactors, we need u / g and v / g */ + + mpz_divexact (s1, v, tv); + mpz_abs (s1, s1); + mpz_divexact (t1, u, tv); + mpz_abs (t1, t1); + + while (power-- > 0) + { + /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */ + if (mpz_odd_p (s0) || mpz_odd_p (t0)) + { + mpz_sub (s0, s0, s1); + mpz_add (t0, t0, t1); + } + assert (mpz_even_p (t0) && mpz_even_p (s0)); + mpz_tdiv_q_2exp (s0, s0, 1); + mpz_tdiv_q_2exp (t0, t0, 1); + } + + /* Arrange so that |s| < |u| / 2g */ + mpz_add (s1, s0, s1); + if (mpz_cmpabs (s0, s1) > 0) + { + mpz_swap (s0, s1); + mpz_sub (t0, t0, t1); + } + if (u->_mp_size < 0) + mpz_neg (s0, s0); + if (v->_mp_size < 0) + mpz_neg (t0, t0); + + mpz_swap (g, tv); + if (s) + mpz_swap (s, s0); + if (t) + mpz_swap (t, t0); + + mpz_clear (tu); + mpz_clear (tv); + mpz_clear (s0); + mpz_clear (s1); + mpz_clear (t0); + mpz_clear (t1); +} + +void +mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v) +{ + mpz_t g; + + if (u->_mp_size == 0 || v->_mp_size == 0) + { + r->_mp_size = 0; + return; + } + + mpz_init (g); + + mpz_gcd (g, u, v); + mpz_divexact (g, u, g); + mpz_mul (r, g, v); + + mpz_clear (g); + mpz_abs (r, r); +} + +void +mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v) +{ + if (v == 0 || u->_mp_size == 0) + { + r->_mp_size = 0; + return; + } + + v /= mpz_gcd_ui (NULL, u, v); + mpz_mul_ui (r, u, v); + + mpz_abs (r, r); +} + +int +mpz_invert (mpz_t r, const mpz_t u, const mpz_t m) +{ + mpz_t g, tr; + int invertible; + + if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0) + return 0; + + mpz_init (g); + mpz_init (tr); + + mpz_gcdext (g, tr, NULL, u, m); + invertible = (mpz_cmp_ui (g, 1) == 0); + + if (invertible) + { + if (tr->_mp_size < 0) + { + if (m->_mp_size >= 0) + mpz_add (tr, tr, m); + else + mpz_sub (tr, tr, m); + } + mpz_swap (r, tr); + } + + mpz_clear (g); + mpz_clear (tr); + return invertible; +} + + +/* Higher level operations (sqrt, pow and root) */ + +void +mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e) +{ + unsigned long bit; + mpz_t tr; + mpz_init_set_ui (tr, 1); + + bit = GMP_ULONG_HIGHBIT; + do + { + mpz_mul (tr, tr, tr); + if (e & bit) + mpz_mul (tr, tr, b); + bit >>= 1; + } + while (bit > 0); + + mpz_swap (r, tr); + mpz_clear (tr); +} + +void +mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e) +{ + mpz_t b; + + mpz_init_set_ui (b, blimb); + mpz_pow_ui (r, b, e); + mpz_clear (b); +} + +void +mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m) +{ + mpz_t tr; + mpz_t base; + mp_size_t en, mn; + mp_srcptr mp; + struct gmp_div_inverse minv; + unsigned shift; + mp_ptr tp = NULL; + + en = GMP_ABS (e->_mp_size); + mn = GMP_ABS (m->_mp_size); + if (mn == 0) + gmp_die ("mpz_powm: Zero modulo."); + + if (en == 0) + { + mpz_set_ui (r, mpz_cmpabs_ui (m, 1)); + return; + } + + mp = m->_mp_d; + mpn_div_qr_invert (&minv, mp, mn); + shift = minv.shift; + + if (shift > 0) + { + /* To avoid shifts, we do all our reductions, except the final + one, using a *normalized* m. */ + minv.shift = 0; + + tp = gmp_alloc_limbs (mn); + gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift)); + mp = tp; + } + + mpz_init (base); + + if (e->_mp_size < 0) + { + if (!mpz_invert (base, b, m)) + gmp_die ("mpz_powm: Negative exponent and non-invertible base."); + } + else + { + mp_size_t bn; + mpz_abs (base, b); + + bn = base->_mp_size; + if (bn >= mn) + { + mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv); + bn = mn; + } + + /* We have reduced the absolute value. Now take care of the + sign. Note that we get zero represented non-canonically as + m. */ + if (b->_mp_size < 0) + { + mp_ptr bp = MPZ_REALLOC (base, mn); + gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn)); + bn = mn; + } + base->_mp_size = mpn_normalized_size (base->_mp_d, bn); + } + mpz_init_set_ui (tr, 1); + + while (--en >= 0) + { + mp_limb_t w = e->_mp_d[en]; + mp_limb_t bit; + + bit = GMP_LIMB_HIGHBIT; + do + { + mpz_mul (tr, tr, tr); + if (w & bit) + mpz_mul (tr, tr, base); + if (tr->_mp_size > mn) + { + mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv); + tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn); + } + bit >>= 1; + } + while (bit > 0); + } + + /* Final reduction */ + if (tr->_mp_size >= mn) + { + minv.shift = shift; + mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv); + tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn); + } + if (tp) + gmp_free_limbs (tp, mn); + + mpz_swap (r, tr); + mpz_clear (tr); + mpz_clear (base); +} + +void +mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m) +{ + mpz_t e; + + mpz_init_set_ui (e, elimb); + mpz_powm (r, b, e, m); + mpz_clear (e); +} + +/* x=trunc(y^(1/z)), r=y-x^z */ +void +mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z) +{ + int sgn; + mp_bitcnt_t bc; + mpz_t t, u; + + sgn = y->_mp_size < 0; + if ((~z & sgn) != 0) + gmp_die ("mpz_rootrem: Negative argument, with even root."); + if (z == 0) + gmp_die ("mpz_rootrem: Zeroth root."); + + if (mpz_cmpabs_ui (y, 1) <= 0) { + if (x) + mpz_set (x, y); + if (r) + r->_mp_size = 0; + return; + } + + mpz_init (u); + mpz_init (t); + bc = (mpz_sizeinbase (y, 2) - 1) / z + 1; + mpz_setbit (t, bc); + + if (z == 2) /* simplify sqrt loop: z-1 == 1 */ + do { + mpz_swap (u, t); /* u = x */ + mpz_tdiv_q (t, y, u); /* t = y/x */ + mpz_add (t, t, u); /* t = y/x + x */ + mpz_tdiv_q_2exp (t, t, 1); /* x'= (y/x + x)/2 */ + } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */ + else /* z != 2 */ { + mpz_t v; + + mpz_init (v); + if (sgn) + mpz_neg (t, t); + + do { + mpz_swap (u, t); /* u = x */ + mpz_pow_ui (t, u, z - 1); /* t = x^(z-1) */ + mpz_tdiv_q (t, y, t); /* t = y/x^(z-1) */ + mpz_mul_ui (v, u, z - 1); /* v = x*(z-1) */ + mpz_add (t, t, v); /* t = y/x^(z-1) + x*(z-1) */ + mpz_tdiv_q_ui (t, t, z); /* x'=(y/x^(z-1) + x*(z-1))/z */ + } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */ + + mpz_clear (v); + } + + if (r) { + mpz_pow_ui (t, u, z); + mpz_sub (r, y, t); + } + if (x) + mpz_swap (x, u); + mpz_clear (u); + mpz_clear (t); +} + +int +mpz_root (mpz_t x, const mpz_t y, unsigned long z) +{ + int res; + mpz_t r; + + mpz_init (r); + mpz_rootrem (x, r, y, z); + res = r->_mp_size == 0; + mpz_clear (r); + + return res; +} + +/* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */ +void +mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u) +{ + mpz_rootrem (s, r, u, 2); +} + +void +mpz_sqrt (mpz_t s, const mpz_t u) +{ + mpz_rootrem (s, NULL, u, 2); +} + +int +mpz_perfect_square_p (const mpz_t u) +{ + if (u->_mp_size <= 0) + return (u->_mp_size == 0); + else + return mpz_root (NULL, u, 2); +} + +int +mpn_perfect_square_p (mp_srcptr p, mp_size_t n) +{ + mpz_t t; + + assert (n > 0); + assert (p [n-1] != 0); + return mpz_root (NULL, mpz_roinit_normal_n (t, p, n), 2); +} + +mp_size_t +mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n) +{ + mpz_t s, r, u; + mp_size_t res; + + assert (n > 0); + assert (p [n-1] != 0); + + mpz_init (r); + mpz_init (s); + mpz_rootrem (s, r, mpz_roinit_normal_n (u, p, n), 2); + + assert (s->_mp_size == (n+1)/2); + mpn_copyd (sp, s->_mp_d, s->_mp_size); + mpz_clear (s); + res = r->_mp_size; + if (rp) + mpn_copyd (rp, r->_mp_d, res); + mpz_clear (r); + return res; +} + +/* Combinatorics */ + +void +mpz_mfac_uiui (mpz_t x, unsigned long n, unsigned long m) +{ + mpz_set_ui (x, n + (n == 0)); + if (m + 1 < 2) return; + while (n > m + 1) + mpz_mul_ui (x, x, n -= m); +} + +void +mpz_2fac_ui (mpz_t x, unsigned long n) +{ + mpz_mfac_uiui (x, n, 2); +} + +void +mpz_fac_ui (mpz_t x, unsigned long n) +{ + mpz_mfac_uiui (x, n, 1); +} + +void +mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k) +{ + mpz_t t; + + mpz_set_ui (r, k <= n); + + if (k > (n >> 1)) + k = (k <= n) ? n - k : 0; + + mpz_init (t); + mpz_fac_ui (t, k); + + for (; k > 0; --k) + mpz_mul_ui (r, r, n--); + + mpz_divexact (r, r, t); + mpz_clear (t); +} + + +/* Primality testing */ + +/* Computes Kronecker (a/b) with odd b, a!=0 and GCD(a,b) = 1 */ +/* Adapted from JACOBI_BASE_METHOD==4 in mpn/generic/jacbase.c */ +static int +gmp_jacobi_coprime (mp_limb_t a, mp_limb_t b) +{ + int c, bit = 0; + + assert (b & 1); + assert (a != 0); + /* assert (mpn_gcd_11 (a, b) == 1); */ + + /* Below, we represent a and b shifted right so that the least + significant one bit is implicit. */ + b >>= 1; + + gmp_ctz(c, a); + a >>= 1; + + for (;;) + { + a >>= c; + /* (2/b) = -1 if b = 3 or 5 mod 8 */ + bit ^= c & (b ^ (b >> 1)); + if (a < b) + { + if (a == 0) + return bit & 1 ? -1 : 1; + bit ^= a & b; + a = b - a; + b -= a; + } + else + { + a -= b; + assert (a != 0); + } + + gmp_ctz(c, a); + ++c; + } +} + +static void +gmp_lucas_step_k_2k (mpz_t V, mpz_t Qk, const mpz_t n) +{ + mpz_mod (Qk, Qk, n); + /* V_{2k} <- V_k ^ 2 - 2Q^k */ + mpz_mul (V, V, V); + mpz_submul_ui (V, Qk, 2); + mpz_tdiv_r (V, V, n); + /* Q^{2k} = (Q^k)^2 */ + mpz_mul (Qk, Qk, Qk); +} + +/* Computes V_k, Q^k (mod n) for the Lucas' sequence */ +/* with P=1, Q=Q; k = (n>>b0)|1. */ +/* Requires an odd n > 4; b0 > 0; -2*Q must not overflow a long */ +/* Returns (U_k == 0) and sets V=V_k and Qk=Q^k. */ +static int +gmp_lucas_mod (mpz_t V, mpz_t Qk, long Q, + mp_bitcnt_t b0, const mpz_t n) +{ + mp_bitcnt_t bs; + mpz_t U; + int res; + + assert (b0 > 0); + assert (Q <= - (LONG_MIN / 2)); + assert (Q >= - (LONG_MAX / 2)); + assert (mpz_cmp_ui (n, 4) > 0); + assert (mpz_odd_p (n)); + + mpz_init_set_ui (U, 1); /* U1 = 1 */ + mpz_set_ui (V, 1); /* V1 = 1 */ + mpz_set_si (Qk, Q); + + for (bs = mpz_sizeinbase (n, 2) - 1; --bs >= b0;) + { + /* U_{2k} <- U_k * V_k */ + mpz_mul (U, U, V); + /* V_{2k} <- V_k ^ 2 - 2Q^k */ + /* Q^{2k} = (Q^k)^2 */ + gmp_lucas_step_k_2k (V, Qk, n); + + /* A step k->k+1 is performed if the bit in $n$ is 1 */ + /* mpz_tstbit(n,bs) or the bit is 0 in $n$ but */ + /* should be 1 in $n+1$ (bs == b0) */ + if (b0 == bs || mpz_tstbit (n, bs)) + { + /* Q^{k+1} <- Q^k * Q */ + mpz_mul_si (Qk, Qk, Q); + /* U_{k+1} <- (U_k + V_k) / 2 */ + mpz_swap (U, V); /* Keep in V the old value of U_k */ + mpz_add (U, U, V); + /* We have to compute U/2, so we need an even value, */ + /* equivalent (mod n) */ + if (mpz_odd_p (U)) + mpz_add (U, U, n); + mpz_tdiv_q_2exp (U, U, 1); + /* V_{k+1} <-(D*U_k + V_k) / 2 = + U_{k+1} + (D-1)/2*U_k = U_{k+1} - 2Q*U_k */ + mpz_mul_si (V, V, -2*Q); + mpz_add (V, U, V); + mpz_tdiv_r (V, V, n); + } + mpz_tdiv_r (U, U, n); + } + + res = U->_mp_size == 0; + mpz_clear (U); + return res; +} + +/* Performs strong Lucas' test on x, with parameters suggested */ +/* for the BPSW test. Qk is only passed to recycle a variable. */ +/* Requires GCD (x,6) = 1.*/ +static int +gmp_stronglucas (const mpz_t x, mpz_t Qk) +{ + mp_bitcnt_t b0; + mpz_t V, n; + mp_limb_t maxD, D; /* The absolute value is stored. */ + long Q; + mp_limb_t tl; + + /* Test on the absolute value. */ + mpz_roinit_normal_n (n, x->_mp_d, GMP_ABS (x->_mp_size)); + + assert (mpz_odd_p (n)); + /* assert (mpz_gcd_ui (NULL, n, 6) == 1); */ + if (mpz_root (Qk, n, 2)) + return 0; /* A square is composite. */ + + /* Check Ds up to square root (in case, n is prime) + or avoid overflows */ + maxD = (Qk->_mp_size == 1) ? Qk->_mp_d [0] - 1 : GMP_LIMB_MAX; + + D = 3; + /* Search a D such that (D/n) = -1 in the sequence 5,-7,9,-11,.. */ + /* For those Ds we have (D/n) = (n/|D|) */ + do + { + if (D >= maxD) + return 1 + (D != GMP_LIMB_MAX); /* (1 + ! ~ D) */ + D += 2; + tl = mpz_tdiv_ui (n, D); + if (tl == 0) + return 0; + } + while (gmp_jacobi_coprime (tl, D) == 1); + + mpz_init (V); + + /* n-(D/n) = n+1 = d*2^{b0}, with d = (n>>b0) | 1 */ + b0 = mpn_common_scan (~ n->_mp_d[0], 0, n->_mp_d, n->_mp_size, GMP_LIMB_MAX); + /* b0 = mpz_scan0 (n, 0); */ + + /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */ + Q = (D & 2) ? (long) (D >> 2) + 1 : -(long) (D >> 2); + + if (! gmp_lucas_mod (V, Qk, Q, b0, n)) /* If Ud != 0 */ + while (V->_mp_size != 0 && --b0 != 0) /* while Vk != 0 */ + /* V <- V ^ 2 - 2Q^k */ + /* Q^{2k} = (Q^k)^2 */ + gmp_lucas_step_k_2k (V, Qk, n); + + mpz_clear (V); + return (b0 != 0); +} + +static int +gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y, + const mpz_t q, mp_bitcnt_t k) +{ + assert (k > 0); + + /* Caller must initialize y to the base. */ + mpz_powm (y, y, q, n); + + if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0) + return 1; + + while (--k > 0) + { + mpz_powm_ui (y, y, 2, n); + if (mpz_cmp (y, nm1) == 0) + return 1; + } + return 0; +} + +/* This product is 0xc0cfd797, and fits in 32 bits. */ +#define GMP_PRIME_PRODUCT \ + (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL) + +/* Bit (p+1)/2 is set, for each odd prime <= 61 */ +#define GMP_PRIME_MASK 0xc96996dcUL + +int +mpz_probab_prime_p (const mpz_t n, int reps) +{ + mpz_t nm1; + mpz_t q; + mpz_t y; + mp_bitcnt_t k; + int is_prime; + int j; + + /* Note that we use the absolute value of n only, for compatibility + with the real GMP. */ + if (mpz_even_p (n)) + return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0; + + /* Above test excludes n == 0 */ + assert (n->_mp_size != 0); + + if (mpz_cmpabs_ui (n, 64) < 0) + return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2; + + if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1) + return 0; + + /* All prime factors are >= 31. */ + if (mpz_cmpabs_ui (n, 31*31) < 0) + return 2; + + mpz_init (nm1); + mpz_init (q); + + /* Find q and k, where q is odd and n = 1 + 2**k * q. */ + mpz_abs (nm1, n); + nm1->_mp_d[0] -= 1; + /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */ + k = mpn_scan1 (nm1->_mp_d, 0); + mpz_tdiv_q_2exp (q, nm1, k); + + /* BPSW test */ + mpz_init_set_ui (y, 2); + is_prime = gmp_millerrabin (n, nm1, y, q, k) && gmp_stronglucas (n, y); + reps -= 24; /* skip the first 24 repetitions */ + + /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] = + j^2 + j + 41 using Euler's polynomial. We potentially stop early, + if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps > + 30 (a[30] == 971 > 31*31 == 961). */ + + for (j = 0; is_prime & (j < reps); j++) + { + mpz_set_ui (y, (unsigned long) j*j+j+41); + if (mpz_cmp (y, nm1) >= 0) + { + /* Don't try any further bases. This "early" break does not affect + the result for any reasonable reps value (<=5000 was tested) */ + assert (j >= 30); + break; + } + is_prime = gmp_millerrabin (n, nm1, y, q, k); + } + mpz_clear (nm1); + mpz_clear (q); + mpz_clear (y); + + return is_prime; +} + + +/* Logical operations and bit manipulation. */ + +/* Numbers are treated as if represented in two's complement (and + infinitely sign extended). For a negative values we get the two's + complement from -x = ~x + 1, where ~ is bitwise complement. + Negation transforms + + xxxx10...0 + + into + + yyyy10...0 + + where yyyy is the bitwise complement of xxxx. So least significant + bits, up to and including the first one bit, are unchanged, and + the more significant bits are all complemented. + + To change a bit from zero to one in a negative number, subtract the + corresponding power of two from the absolute value. This can never + underflow. To change a bit from one to zero, add the corresponding + power of two, and this might overflow. E.g., if x = -001111, the + two's complement is 110001. Clearing the least significant bit, we + get two's complement 110000, and -010000. */ + +int +mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index) +{ + mp_size_t limb_index; + unsigned shift; + mp_size_t ds; + mp_size_t dn; + mp_limb_t w; + int bit; + + ds = d->_mp_size; + dn = GMP_ABS (ds); + limb_index = bit_index / GMP_LIMB_BITS; + if (limb_index >= dn) + return ds < 0; + + shift = bit_index % GMP_LIMB_BITS; + w = d->_mp_d[limb_index]; + bit = (w >> shift) & 1; + + if (ds < 0) + { + /* d < 0. Check if any of the bits below is set: If so, our bit + must be complemented. */ + if (shift > 0 && (mp_limb_t) (w << (GMP_LIMB_BITS - shift)) > 0) + return bit ^ 1; + while (--limb_index >= 0) + if (d->_mp_d[limb_index] > 0) + return bit ^ 1; + } + return bit; +} + +static void +mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index) +{ + mp_size_t dn, limb_index; + mp_limb_t bit; + mp_ptr dp; + + dn = GMP_ABS (d->_mp_size); + + limb_index = bit_index / GMP_LIMB_BITS; + bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS); + + if (limb_index >= dn) + { + mp_size_t i; + /* The bit should be set outside of the end of the number. + We have to increase the size of the number. */ + dp = MPZ_REALLOC (d, limb_index + 1); + + dp[limb_index] = bit; + for (i = dn; i < limb_index; i++) + dp[i] = 0; + dn = limb_index + 1; + } + else + { + mp_limb_t cy; + + dp = d->_mp_d; + + cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit); + if (cy > 0) + { + dp = MPZ_REALLOC (d, dn + 1); + dp[dn++] = cy; + } + } + + d->_mp_size = (d->_mp_size < 0) ? - dn : dn; +} + +static void +mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index) +{ + mp_size_t dn, limb_index; + mp_ptr dp; + mp_limb_t bit; + + dn = GMP_ABS (d->_mp_size); + dp = d->_mp_d; + + limb_index = bit_index / GMP_LIMB_BITS; + bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS); + + assert (limb_index < dn); + + gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index, + dn - limb_index, bit)); + dn = mpn_normalized_size (dp, dn); + d->_mp_size = (d->_mp_size < 0) ? - dn : dn; +} + +void +mpz_setbit (mpz_t d, mp_bitcnt_t bit_index) +{ + if (!mpz_tstbit (d, bit_index)) + { + if (d->_mp_size >= 0) + mpz_abs_add_bit (d, bit_index); + else + mpz_abs_sub_bit (d, bit_index); + } +} + +void +mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index) +{ + if (mpz_tstbit (d, bit_index)) + { + if (d->_mp_size >= 0) + mpz_abs_sub_bit (d, bit_index); + else + mpz_abs_add_bit (d, bit_index); + } +} + +void +mpz_combit (mpz_t d, mp_bitcnt_t bit_index) +{ + if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0)) + mpz_abs_sub_bit (d, bit_index); + else + mpz_abs_add_bit (d, bit_index); +} + +void +mpz_com (mpz_t r, const mpz_t u) +{ + mpz_add_ui (r, u, 1); + mpz_neg (r, r); +} + +void +mpz_and (mpz_t r, const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, rn, i; + mp_ptr up, vp, rp; + + mp_limb_t ux, vx, rx; + mp_limb_t uc, vc, rc; + mp_limb_t ul, vl, rl; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + if (un < vn) + { + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (un, vn); + } + if (vn == 0) + { + r->_mp_size = 0; + return; + } + + uc = u->_mp_size < 0; + vc = v->_mp_size < 0; + rc = uc & vc; + + ux = -uc; + vx = -vc; + rx = -rc; + + /* If the smaller input is positive, higher limbs don't matter. */ + rn = vx ? un : vn; + + rp = MPZ_REALLOC (r, rn + (mp_size_t) rc); + + up = u->_mp_d; + vp = v->_mp_d; + + i = 0; + do + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + vl = (vp[i] ^ vx) + vc; + vc = vl < vc; + + rl = ( (ul & vl) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + while (++i < vn); + assert (vc == 0); + + for (; i < rn; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + rl = ( (ul & vx) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + if (rc) + rp[rn++] = rc; + else + rn = mpn_normalized_size (rp, rn); + + r->_mp_size = rx ? -rn : rn; +} + +void +mpz_ior (mpz_t r, const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, rn, i; + mp_ptr up, vp, rp; + + mp_limb_t ux, vx, rx; + mp_limb_t uc, vc, rc; + mp_limb_t ul, vl, rl; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + if (un < vn) + { + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (un, vn); + } + if (vn == 0) + { + mpz_set (r, u); + return; + } + + uc = u->_mp_size < 0; + vc = v->_mp_size < 0; + rc = uc | vc; + + ux = -uc; + vx = -vc; + rx = -rc; + + /* If the smaller input is negative, by sign extension higher limbs + don't matter. */ + rn = vx ? vn : un; + + rp = MPZ_REALLOC (r, rn + (mp_size_t) rc); + + up = u->_mp_d; + vp = v->_mp_d; + + i = 0; + do + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + vl = (vp[i] ^ vx) + vc; + vc = vl < vc; + + rl = ( (ul | vl) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + while (++i < vn); + assert (vc == 0); + + for (; i < rn; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + rl = ( (ul | vx) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + if (rc) + rp[rn++] = rc; + else + rn = mpn_normalized_size (rp, rn); + + r->_mp_size = rx ? -rn : rn; +} + +void +mpz_xor (mpz_t r, const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, i; + mp_ptr up, vp, rp; + + mp_limb_t ux, vx, rx; + mp_limb_t uc, vc, rc; + mp_limb_t ul, vl, rl; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + if (un < vn) + { + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (un, vn); + } + if (vn == 0) + { + mpz_set (r, u); + return; + } + + uc = u->_mp_size < 0; + vc = v->_mp_size < 0; + rc = uc ^ vc; + + ux = -uc; + vx = -vc; + rx = -rc; + + rp = MPZ_REALLOC (r, un + (mp_size_t) rc); + + up = u->_mp_d; + vp = v->_mp_d; + + i = 0; + do + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + vl = (vp[i] ^ vx) + vc; + vc = vl < vc; + + rl = (ul ^ vl ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + while (++i < vn); + assert (vc == 0); + + for (; i < un; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + rl = (ul ^ ux) + rc; + rc = rl < rc; + rp[i] = rl; + } + if (rc) + rp[un++] = rc; + else + un = mpn_normalized_size (rp, un); + + r->_mp_size = rx ? -un : un; +} + +static unsigned +gmp_popcount_limb (mp_limb_t x) +{ + unsigned c; + + /* Do 16 bits at a time, to avoid limb-sized constants. */ + int LOCAL_SHIFT_BITS = 16; + for (c = 0; x > 0;) + { + unsigned w = x - ((x >> 1) & 0x5555); + w = ((w >> 2) & 0x3333) + (w & 0x3333); + w = (w >> 4) + w; + w = ((w >> 8) & 0x000f) + (w & 0x000f); + c += w; + if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS) + x >>= LOCAL_SHIFT_BITS; + else + x = 0; + } + return c; +} + +mp_bitcnt_t +mpn_popcount (mp_srcptr p, mp_size_t n) +{ + mp_size_t i; + mp_bitcnt_t c; + + for (c = 0, i = 0; i < n; i++) + c += gmp_popcount_limb (p[i]); + + return c; +} + +mp_bitcnt_t +mpz_popcount (const mpz_t u) +{ + mp_size_t un; + + un = u->_mp_size; + + if (un < 0) + return ~(mp_bitcnt_t) 0; + + return mpn_popcount (u->_mp_d, un); +} + +mp_bitcnt_t +mpz_hamdist (const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, i; + mp_limb_t uc, vc, ul, vl, comp; + mp_srcptr up, vp; + mp_bitcnt_t c; + + un = u->_mp_size; + vn = v->_mp_size; + + if ( (un ^ vn) < 0) + return ~(mp_bitcnt_t) 0; + + comp = - (uc = vc = (un < 0)); + if (uc) + { + assert (vn < 0); + un = -un; + vn = -vn; + } + + up = u->_mp_d; + vp = v->_mp_d; + + if (un < vn) + MPN_SRCPTR_SWAP (up, un, vp, vn); + + for (i = 0, c = 0; i < vn; i++) + { + ul = (up[i] ^ comp) + uc; + uc = ul < uc; + + vl = (vp[i] ^ comp) + vc; + vc = vl < vc; + + c += gmp_popcount_limb (ul ^ vl); + } + assert (vc == 0); + + for (; i < un; i++) + { + ul = (up[i] ^ comp) + uc; + uc = ul < uc; + + c += gmp_popcount_limb (ul ^ comp); + } + + return c; +} + +mp_bitcnt_t +mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit) +{ + mp_ptr up; + mp_size_t us, un, i; + mp_limb_t limb, ux; + + us = u->_mp_size; + un = GMP_ABS (us); + i = starting_bit / GMP_LIMB_BITS; + + /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit + for u<0. Notice this test picks up any u==0 too. */ + if (i >= un) + return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit); + + up = u->_mp_d; + ux = 0; + limb = up[i]; + + if (starting_bit != 0) + { + if (us < 0) + { + ux = mpn_zero_p (up, i); + limb = ~ limb + ux; + ux = - (mp_limb_t) (limb >= ux); + } + + /* Mask to 0 all bits before starting_bit, thus ignoring them. */ + limb &= GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS); + } + + return mpn_common_scan (limb, i, up, un, ux); +} + +mp_bitcnt_t +mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit) +{ + mp_ptr up; + mp_size_t us, un, i; + mp_limb_t limb, ux; + + us = u->_mp_size; + ux = - (mp_limb_t) (us >= 0); + un = GMP_ABS (us); + i = starting_bit / GMP_LIMB_BITS; + + /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for + u<0. Notice this test picks up all cases of u==0 too. */ + if (i >= un) + return (ux ? starting_bit : ~(mp_bitcnt_t) 0); + + up = u->_mp_d; + limb = up[i] ^ ux; + + if (ux == 0) + limb -= mpn_zero_p (up, i); /* limb = ~(~limb + zero_p) */ + + /* Mask all bits before starting_bit, thus ignoring them. */ + limb &= GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS); + + return mpn_common_scan (limb, i, up, un, ux); +} + + +/* MPZ base conversion. */ + +size_t +mpz_sizeinbase (const mpz_t u, int base) +{ + mp_size_t un, tn; + mp_srcptr up; + mp_ptr tp; + mp_bitcnt_t bits; + struct gmp_div_inverse bi; + size_t ndigits; + + assert (base >= 2); + assert (base <= 62); + + un = GMP_ABS (u->_mp_size); + if (un == 0) + return 1; + + up = u->_mp_d; + + bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]); + switch (base) + { + case 2: + return bits; + case 4: + return (bits + 1) / 2; + case 8: + return (bits + 2) / 3; + case 16: + return (bits + 3) / 4; + case 32: + return (bits + 4) / 5; + /* FIXME: Do something more clever for the common case of base + 10. */ + } + + tp = gmp_alloc_limbs (un); + mpn_copyi (tp, up, un); + mpn_div_qr_1_invert (&bi, base); + + tn = un; + ndigits = 0; + do + { + ndigits++; + mpn_div_qr_1_preinv (tp, tp, tn, &bi); + tn -= (tp[tn-1] == 0); + } + while (tn > 0); + + gmp_free_limbs (tp, un); + return ndigits; +} + +char * +mpz_get_str (char *sp, int base, const mpz_t u) +{ + unsigned bits; + const char *digits; + mp_size_t un; + size_t i, sn, osn; + + digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"; + if (base > 1) + { + if (base <= 36) + digits = "0123456789abcdefghijklmnopqrstuvwxyz"; + else if (base > 62) + return NULL; + } + else if (base >= -1) + base = 10; + else + { + base = -base; + if (base > 36) + return NULL; + } + + sn = 1 + mpz_sizeinbase (u, base); + if (!sp) + { + osn = 1 + sn; + sp = (char *) gmp_alloc (osn); + } + else + osn = 0; + un = GMP_ABS (u->_mp_size); + + if (un == 0) + { + sp[0] = '0'; + sn = 1; + goto ret; + } + + i = 0; + + if (u->_mp_size < 0) + sp[i++] = '-'; + + bits = mpn_base_power_of_two_p (base); + + if (bits) + /* Not modified in this case. */ + sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un); + else + { + struct mpn_base_info info; + mp_ptr tp; + + mpn_get_base_info (&info, base); + tp = gmp_alloc_limbs (un); + mpn_copyi (tp, u->_mp_d, un); + + sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un); + gmp_free_limbs (tp, un); + } + + for (; i < sn; i++) + sp[i] = digits[(unsigned char) sp[i]]; + +ret: + sp[sn] = '\0'; + if (osn && osn != sn + 1) + sp = (char*) gmp_realloc (sp, osn, sn + 1); + return sp; +} + +int +mpz_set_str (mpz_t r, const char *sp, int base) +{ + unsigned bits, value_of_a; + mp_size_t rn, alloc; + mp_ptr rp; + size_t dn, sn; + int sign; + unsigned char *dp; + + assert (base == 0 || (base >= 2 && base <= 62)); + + while (isspace( (unsigned char) *sp)) + sp++; + + sign = (*sp == '-'); + sp += sign; + + if (base == 0) + { + if (sp[0] == '0') + { + if (sp[1] == 'x' || sp[1] == 'X') + { + base = 16; + sp += 2; + } + else if (sp[1] == 'b' || sp[1] == 'B') + { + base = 2; + sp += 2; + } + else + base = 8; + } + else + base = 10; + } + + if (!*sp) + { + r->_mp_size = 0; + return -1; + } + sn = strlen(sp); + dp = (unsigned char *) gmp_alloc (sn); + + value_of_a = (base > 36) ? 36 : 10; + for (dn = 0; *sp; sp++) + { + unsigned digit; + + if (isspace ((unsigned char) *sp)) + continue; + else if (*sp >= '0' && *sp <= '9') + digit = *sp - '0'; + else if (*sp >= 'a' && *sp <= 'z') + digit = *sp - 'a' + value_of_a; + else if (*sp >= 'A' && *sp <= 'Z') + digit = *sp - 'A' + 10; + else + digit = base; /* fail */ + + if (digit >= (unsigned) base) + { + gmp_free (dp, sn); + r->_mp_size = 0; + return -1; + } + + dp[dn++] = digit; + } + + if (!dn) + { + gmp_free (dp, sn); + r->_mp_size = 0; + return -1; + } + bits = mpn_base_power_of_two_p (base); + + if (bits > 0) + { + alloc = (dn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; + rp = MPZ_REALLOC (r, alloc); + rn = mpn_set_str_bits (rp, dp, dn, bits); + } + else + { + struct mpn_base_info info; + mpn_get_base_info (&info, base); + alloc = (dn + info.exp - 1) / info.exp; + rp = MPZ_REALLOC (r, alloc); + rn = mpn_set_str_other (rp, dp, dn, base, &info); + /* Normalization, needed for all-zero input. */ + assert (rn > 0); + rn -= rp[rn-1] == 0; + } + assert (rn <= alloc); + gmp_free (dp, sn); + + r->_mp_size = sign ? - rn : rn; + + return 0; +} + +int +mpz_init_set_str (mpz_t r, const char *sp, int base) +{ + mpz_init (r); + return mpz_set_str (r, sp, base); +} + +size_t +mpz_out_str (FILE *stream, int base, const mpz_t x) +{ + char *str; + size_t len, n; + + str = mpz_get_str (NULL, base, x); + if (!str) + return 0; + len = strlen (str); + n = fwrite (str, 1, len, stream); + gmp_free (str, len + 1); + return n; +} + + +static int +gmp_detect_endian (void) +{ + static const int i = 2; + const unsigned char *p = (const unsigned char *) &i; + return 1 - *p; +} + +/* Import and export. Does not support nails. */ +void +mpz_import (mpz_t r, size_t count, int order, size_t size, int endian, + size_t nails, const void *src) +{ + const unsigned char *p; + ptrdiff_t word_step; + mp_ptr rp; + mp_size_t rn; + + /* The current (partial) limb. */ + mp_limb_t limb; + /* The number of bytes already copied to this limb (starting from + the low end). */ + size_t bytes; + /* The index where the limb should be stored, when completed. */ + mp_size_t i; + + if (nails != 0) + gmp_die ("mpz_import: Nails not supported."); + + assert (order == 1 || order == -1); + assert (endian >= -1 && endian <= 1); + + if (endian == 0) + endian = gmp_detect_endian (); + + p = (unsigned char *) src; + + word_step = (order != endian) ? 2 * size : 0; + + /* Process bytes from the least significant end, so point p at the + least significant word. */ + if (order == 1) + { + p += size * (count - 1); + word_step = - word_step; + } + + /* And at least significant byte of that word. */ + if (endian == 1) + p += (size - 1); + + rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t); + rp = MPZ_REALLOC (r, rn); + + for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step) + { + size_t j; + for (j = 0; j < size; j++, p -= (ptrdiff_t) endian) + { + limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT); + if (bytes == sizeof(mp_limb_t)) + { + rp[i++] = limb; + bytes = 0; + limb = 0; + } + } + } + assert (i + (bytes > 0) == rn); + if (limb != 0) + rp[i++] = limb; + else + i = mpn_normalized_size (rp, i); + + r->_mp_size = i; +} + +void * +mpz_export (void *r, size_t *countp, int order, size_t size, int endian, + size_t nails, const mpz_t u) +{ + size_t count; + mp_size_t un; + + if (nails != 0) + gmp_die ("mpz_export: Nails not supported."); + + assert (order == 1 || order == -1); + assert (endian >= -1 && endian <= 1); + assert (size > 0 || u->_mp_size == 0); + + un = u->_mp_size; + count = 0; + if (un != 0) + { + size_t k; + unsigned char *p; + ptrdiff_t word_step; + /* The current (partial) limb. */ + mp_limb_t limb; + /* The number of bytes left to do in this limb. */ + size_t bytes; + /* The index where the limb was read. */ + mp_size_t i; + + un = GMP_ABS (un); + + /* Count bytes in top limb. */ + limb = u->_mp_d[un-1]; + assert (limb != 0); + + k = (GMP_LIMB_BITS <= CHAR_BIT); + if (!k) + { + do { + int LOCAL_CHAR_BIT = CHAR_BIT; + k++; limb >>= LOCAL_CHAR_BIT; + } while (limb != 0); + } + /* else limb = 0; */ + + count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size; + + if (!r) + r = gmp_alloc (count * size); + + if (endian == 0) + endian = gmp_detect_endian (); + + p = (unsigned char *) r; + + word_step = (order != endian) ? 2 * size : 0; + + /* Process bytes from the least significant end, so point p at the + least significant word. */ + if (order == 1) + { + p += size * (count - 1); + word_step = - word_step; + } + + /* And at least significant byte of that word. */ + if (endian == 1) + p += (size - 1); + + for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step) + { + size_t j; + for (j = 0; j < size; ++j, p -= (ptrdiff_t) endian) + { + if (sizeof (mp_limb_t) == 1) + { + if (i < un) + *p = u->_mp_d[i++]; + else + *p = 0; + } + else + { + int LOCAL_CHAR_BIT = CHAR_BIT; + if (bytes == 0) + { + if (i < un) + limb = u->_mp_d[i++]; + bytes = sizeof (mp_limb_t); + } + *p = limb; + limb >>= LOCAL_CHAR_BIT; + bytes--; + } + } + } + assert (i == un); + assert (k == count); + } + + if (countp) + *countp = count; + + return r; +} + diff --git a/source/gmp.h b/source/gmp.h @@ -0,0 +1,316 @@ +// +// gmp.h +// lib-gpu-verify +// +// Created by Cedric Zwahlen on 12.11.2023. +// + +#ifndef gmp_h +#define gmp_h + +#include <stdio.h> + + +/* mini-gmp, a minimalistic implementation of a GNU GMP subset. + +Copyright 2011-2015, 2017, 2019-2021 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +/* About mini-gmp: This is a minimal implementation of a subset of the + GMP interface. It is intended for inclusion into applications which + have modest bignums needs, as a fallback when the real GMP library + is not installed. + + This file defines the public interface. */ + + +/* For size_t */ +#include <stddef.h> + + +void mp_set_memory_functions (void *(*) (size_t), + void *(*) (void *, size_t, size_t), + void (*) (void *, size_t)); + +void mp_get_memory_functions (void *(**) (size_t), + void *(**) (void *, size_t, size_t), + void (**) (void *, size_t)); + +#ifndef MINI_GMP_LIMB_TYPE +#define MINI_GMP_LIMB_TYPE long +#endif + +typedef unsigned MINI_GMP_LIMB_TYPE mp_limb_t; +typedef long mp_size_t; +typedef unsigned long mp_bitcnt_t; + +typedef mp_limb_t *mp_ptr; +typedef const mp_limb_t *mp_srcptr; + +typedef struct +{ + int _mp_alloc; /* Number of *limbs* allocated and pointed + to by the _mp_d field. */ + int _mp_size; /* abs(_mp_size) is the number of limbs the + last field points to. If _mp_size is + negative this is a negative number. */ + mp_limb_t *_mp_d; /* Pointer to the limbs. */ +} __mpz_struct; + +typedef __mpz_struct mpz_t[1]; + +typedef __mpz_struct *mpz_ptr; +typedef const __mpz_struct *mpz_srcptr; + +extern const int mp_bits_per_limb; + +void mpn_copyi (mp_ptr, mp_srcptr, mp_size_t); +void mpn_copyd (mp_ptr, mp_srcptr, mp_size_t); +void mpn_zero (mp_ptr, mp_size_t); + +int mpn_cmp (mp_srcptr, mp_srcptr, mp_size_t); +int mpn_zero_p (mp_srcptr, mp_size_t); + +mp_limb_t mpn_add_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_add_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); +mp_limb_t mpn_add (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t); + +mp_limb_t mpn_sub_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_sub_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); +mp_limb_t mpn_sub (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t); + +mp_limb_t mpn_mul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_addmul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_submul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); + +mp_limb_t mpn_mul (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t); +void mpn_mul_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); +void mpn_sqr (mp_ptr, mp_srcptr, mp_size_t); +int mpn_perfect_square_p (mp_srcptr, mp_size_t); +mp_size_t mpn_sqrtrem (mp_ptr, mp_ptr, mp_srcptr, mp_size_t); + +mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int); +mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int); + +mp_bitcnt_t mpn_scan0 (mp_srcptr, mp_bitcnt_t); +mp_bitcnt_t mpn_scan1 (mp_srcptr, mp_bitcnt_t); + +void mpn_com (mp_ptr, mp_srcptr, mp_size_t); +mp_limb_t mpn_neg (mp_ptr, mp_srcptr, mp_size_t); + +mp_bitcnt_t mpn_popcount (mp_srcptr, mp_size_t); + +mp_limb_t mpn_invert_3by2 (mp_limb_t, mp_limb_t); +#define mpn_invert_limb(x) mpn_invert_3by2 ((x), 0) + +size_t mpn_get_str (unsigned char *, int, mp_ptr, mp_size_t); +mp_size_t mpn_set_str (mp_ptr, const unsigned char *, size_t, int); + +void mpz_init (mpz_t); +void mpz_init2 (mpz_t, mp_bitcnt_t); +void mpz_clear (mpz_t); + +#define mpz_odd_p(z) (((z)->_mp_size != 0) & (int) (z)->_mp_d[0]) +#define mpz_even_p(z) (! mpz_odd_p (z)) + +int mpz_sgn (const mpz_t); +int mpz_cmp_si (const mpz_t, long); +int mpz_cmp_ui (const mpz_t, unsigned long); +int mpz_cmp (const mpz_t, const mpz_t); +int mpz_cmpabs_ui (const mpz_t, unsigned long); +int mpz_cmpabs (const mpz_t, const mpz_t); +int mpz_cmp_d (const mpz_t, double); +int mpz_cmpabs_d (const mpz_t, double); + +void mpz_abs (mpz_t, const mpz_t); +void mpz_neg (mpz_t, const mpz_t); +void mpz_swap (mpz_t, mpz_t); + +void mpz_add_ui (mpz_t, const mpz_t, unsigned long); +void mpz_add (mpz_t, const mpz_t, const mpz_t); +void mpz_sub_ui (mpz_t, const mpz_t, unsigned long); +void mpz_ui_sub (mpz_t, unsigned long, const mpz_t); +void mpz_sub (mpz_t, const mpz_t, const mpz_t); + +void mpz_mul_si (mpz_t, const mpz_t, long int); +void mpz_mul_ui (mpz_t, const mpz_t, unsigned long int); +void mpz_mul (mpz_t, const mpz_t, const mpz_t); +void mpz_mul_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_addmul_ui (mpz_t, const mpz_t, unsigned long int); +void mpz_addmul (mpz_t, const mpz_t, const mpz_t); +void mpz_submul_ui (mpz_t, const mpz_t, unsigned long int); +void mpz_submul (mpz_t, const mpz_t, const mpz_t); + +void mpz_cdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_fdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_tdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_cdiv_q (mpz_t, const mpz_t, const mpz_t); +void mpz_fdiv_q (mpz_t, const mpz_t, const mpz_t); +void mpz_tdiv_q (mpz_t, const mpz_t, const mpz_t); +void mpz_cdiv_r (mpz_t, const mpz_t, const mpz_t); +void mpz_fdiv_r (mpz_t, const mpz_t, const mpz_t); +void mpz_tdiv_r (mpz_t, const mpz_t, const mpz_t); + +void mpz_cdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_fdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_tdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_cdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_fdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_tdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t); + +void mpz_mod (mpz_t, const mpz_t, const mpz_t); + +void mpz_divexact (mpz_t, const mpz_t, const mpz_t); + +int mpz_divisible_p (const mpz_t, const mpz_t); +int mpz_congruent_p (const mpz_t, const mpz_t, const mpz_t); + +unsigned long mpz_cdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long); +unsigned long mpz_fdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long); +unsigned long mpz_tdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long); +unsigned long mpz_cdiv_q_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_fdiv_q_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_tdiv_q_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_cdiv_r_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_fdiv_r_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_tdiv_r_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_cdiv_ui (const mpz_t, unsigned long); +unsigned long mpz_fdiv_ui (const mpz_t, unsigned long); +unsigned long mpz_tdiv_ui (const mpz_t, unsigned long); + +unsigned long mpz_mod_ui (mpz_t, const mpz_t, unsigned long); + +void mpz_divexact_ui (mpz_t, const mpz_t, unsigned long); + +int mpz_divisible_ui_p (const mpz_t, unsigned long); + +unsigned long mpz_gcd_ui (mpz_t, const mpz_t, unsigned long); +void mpz_gcd (mpz_t, const mpz_t, const mpz_t); +void mpz_gcdext (mpz_t, mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_lcm_ui (mpz_t, const mpz_t, unsigned long); +void mpz_lcm (mpz_t, const mpz_t, const mpz_t); +int mpz_invert (mpz_t, const mpz_t, const mpz_t); + +void mpz_sqrtrem (mpz_t, mpz_t, const mpz_t); +void mpz_sqrt (mpz_t, const mpz_t); +int mpz_perfect_square_p (const mpz_t); + +void mpz_pow_ui (mpz_t, const mpz_t, unsigned long); +void mpz_ui_pow_ui (mpz_t, unsigned long, unsigned long); +void mpz_powm (mpz_t, const mpz_t, const mpz_t, const mpz_t); +void mpz_powm_ui (mpz_t, const mpz_t, unsigned long, const mpz_t); + +void mpz_rootrem (mpz_t, mpz_t, const mpz_t, unsigned long); +int mpz_root (mpz_t, const mpz_t, unsigned long); + +void mpz_fac_ui (mpz_t, unsigned long); +void mpz_2fac_ui (mpz_t, unsigned long); +void mpz_mfac_uiui (mpz_t, unsigned long, unsigned long); +void mpz_bin_uiui (mpz_t, unsigned long, unsigned long); + +int mpz_probab_prime_p (const mpz_t, int); + +int mpz_tstbit (const mpz_t, mp_bitcnt_t); +void mpz_setbit (mpz_t, mp_bitcnt_t); +void mpz_clrbit (mpz_t, mp_bitcnt_t); +void mpz_combit (mpz_t, mp_bitcnt_t); + +void mpz_com (mpz_t, const mpz_t); +void mpz_and (mpz_t, const mpz_t, const mpz_t); +void mpz_ior (mpz_t, const mpz_t, const mpz_t); +void mpz_xor (mpz_t, const mpz_t, const mpz_t); + +mp_bitcnt_t mpz_popcount (const mpz_t); +mp_bitcnt_t mpz_hamdist (const mpz_t, const mpz_t); +mp_bitcnt_t mpz_scan0 (const mpz_t, mp_bitcnt_t); +mp_bitcnt_t mpz_scan1 (const mpz_t, mp_bitcnt_t); + +int mpz_fits_slong_p (const mpz_t); +int mpz_fits_ulong_p (const mpz_t); +int mpz_fits_sint_p (const mpz_t); +int mpz_fits_uint_p (const mpz_t); +int mpz_fits_sshort_p (const mpz_t); +int mpz_fits_ushort_p (const mpz_t); +long int mpz_get_si (const mpz_t); +unsigned long int mpz_get_ui (const mpz_t); +double mpz_get_d (const mpz_t); +size_t mpz_size (const mpz_t); +mp_limb_t mpz_getlimbn (const mpz_t, mp_size_t); + +void mpz_realloc2 (mpz_t, mp_bitcnt_t); +mp_srcptr mpz_limbs_read (mpz_srcptr); +mp_ptr mpz_limbs_modify (mpz_t, mp_size_t); +mp_ptr mpz_limbs_write (mpz_t, mp_size_t); +void mpz_limbs_finish (mpz_t, mp_size_t); +mpz_srcptr mpz_roinit_n (mpz_t, mp_srcptr, mp_size_t); + +#define MPZ_ROINIT_N(xp, xs) {{0, (xs),(xp) }} + +void mpz_set_si (mpz_t, signed long int); +void mpz_set_ui (mpz_t, unsigned long int); +void mpz_set (mpz_t, const mpz_t); +void mpz_set_d (mpz_t, double); + +void mpz_init_set_si (mpz_t, signed long int); +void mpz_init_set_ui (mpz_t, unsigned long int); +void mpz_init_set (mpz_t, const mpz_t); +void mpz_init_set_d (mpz_t, double); + +size_t mpz_sizeinbase (const mpz_t, int); +char *mpz_get_str (char *, int, const mpz_t); +int mpz_set_str (mpz_t, const char *, int); +int mpz_init_set_str (mpz_t, const char *, int); + +/* This long list taken from gmp.h. */ +/* For reference, "defined(EOF)" cannot be used here. In g++ 2.95.4, + <iostream> defines EOF but not FILE. */ +#if defined (FILE) \ + || defined (H_STDIO) \ + || defined (_H_STDIO) /* AIX */ \ + || defined (_STDIO_H) /* glibc, Sun, SCO */ \ + || defined (_STDIO_H_) /* BSD, OSF */ \ + || defined (__STDIO_H) /* Borland */ \ + || defined (__STDIO_H__) /* IRIX */ \ + || defined (_STDIO_INCLUDED) /* HPUX */ \ + || defined (__dj_include_stdio_h_) /* DJGPP */ \ + || defined (_FILE_DEFINED) /* Microsoft */ \ + || defined (__STDIO__) /* Apple MPW MrC */ \ + || defined (_MSL_STDIO_H) /* Metrowerks */ \ + || defined (_STDIO_H_INCLUDED) /* QNX4 */ \ + || defined (_ISO_STDIO_ISO_H) /* Sun C++ */ \ + || defined (__STDIO_LOADED) /* VMS */ \ + || defined (_STDIO) /* HPE NonStop */ \ + || defined (__DEFINED_FILE) /* musl */ +size_t mpz_out_str (FILE *, int, const mpz_t); +#endif + +void mpz_import (mpz_t, size_t, int, size_t, int, size_t, const void *); +void *mpz_export (void *, size_t *, int, size_t, int, size_t, const mpz_t); + + +#endif /* gmp_h */ diff --git a/source/lib-gpu-verify.c b/source/lib-gpu-verify.c @@ -5,15 +5,16 @@ #include "rsa-test.h" #include "opencl-test.h" +#include "montgomery.h" int main(int argc, char** argv) { - //mont_prepare("07", "0A", "0D"); + mont_prepare("07", "0A", "0D"); //opencl_tests(); - rsa_tests(); + // rsa_tests(); // montgomery_test(); diff --git a/source/montgomery.c b/source/montgomery.c @@ -0,0 +1,184 @@ +// +// montgomery.c +// lib-gpu-verify +// +// Created by Cedric Zwahlen on 14.11.2023. +// + +#include "montgomery.h" + +// assume 64 bit, keep in mind that on the GPU, we can't dynamically allocate + +void mont_product(mpz_t ret, + const mpz_t a, const mpz_t b, + const mpz_t r, const mpz_t r_1, + const mpz_t n, const mpz_t ni + ); + +void mont_modexp(mpz_t ret, + mpz_t a, mpz_t e, + const mpz_t M, + const mpz_t n, const mpz_t ni, + const mpz_t r, const mpz_t r_1 + ); + +// CPU +void mont_prepare(char *base, char *exponent, char *modulus) { + + mpz_t b,e,m; + + mpz_init_set_str(b,base,16); // M + mpz_init_set_str(e,exponent,16); + mpz_init_set_str(m,modulus,16); // n + + mpz_t r, r_1, ni, M, x; + mpz_init(r); + mpz_init(r_1); + mpz_init(ni); + mpz_init(M); + mpz_init(x); + + + // r and n (modulus) must be relatively prime (this is a given if n (modulus) is odd) + + + // calculate r, which must be larger than the modulo and also a power of 2 + size_t len = mpz_sizeinbase(m,2); + mpz_init_set_si(r, 1 << len); + + mpz_t one, o, oo; // some helper variables + mpz_init(one); // must equal 1 + mpz_init_set_si(o,1); + mpz_init_set_si(oo,0); + + mpz_t temp_r_1, temp_ni; // more helper variables + mpz_init(temp_r_1); + mpz_init(temp_ni); + + mpz_gcdext(one, temp_r_1, temp_ni, r, m); // set r_1 and ni + + mpz_abs(temp_r_1, temp_r_1); + mpz_abs(temp_ni, temp_ni); + + mpz_sub(ni, r, temp_ni); + mpz_sub(r_1, m, temp_r_1); + + // test if one == 1 + + if (mpz_cmp(o, one)) + assert(1); + + mpz_mul(one, r, r_1); + mpz_mul(o,ni,m); + + mpz_sub(oo, one, o); // oo must be one + + mpz_init_set_si(one, 1); + + if (mpz_cmp(oo, one)) + assert(1); + + mpz_mul(M, b, r); + mpz_mod(M, M, m); // set M + + mpz_mod(x, r, m); // set x + + mpz_t res; + + mpz_init(res); + + //mont_product(res, x, x, r, r_1, m, ni); + + mont_modexp(res, x, e, M, m, ni, r, r_1); + + + +} + +// maybe GPU? +// MARK: n MUST be an odd number +void mont_modexp(mpz_t ret, + mpz_t a, mpz_t e, + const mpz_t M, + const mpz_t n, const mpz_t ni, + const mpz_t r, const mpz_t r_1 + ) { + + mpz_t aa,xx; + + mpz_init_set(aa, M); + mpz_init_set(xx, a); + + int k = (int)mpz_sizeinbase(e,2); + + for (int i = k - 1; i >= 0; i--) { + printf("%d",i); + + mont_product(xx, xx, xx, r, r_1, n, ni); + + if (mpz_tstbit(e, i)) + mont_product(xx, aa, xx, r, r_1, n, ni); + + + } + + // this transforms out of the montgomery form, so maybe move that to finish + + mpz_t x,one; + + mpz_init(x); + mpz_init_set_ui(one, 1); + + mont_product(x, xx, one, r, r_1, n, ni); + + mpz_set(ret, x); + +} + +void mont_finish(void) { + + + +} + + +// GPU +void mont_product(mpz_t ret, + const mpz_t a, const mpz_t b, + const mpz_t r, const mpz_t r_1, + const mpz_t n, const mpz_t ni + ) { + + mpz_t t,m,u; + + mpz_init(t); + mpz_init(m); + mpz_init(u); + + // MARK: this is the part we can accelerate... + // because r is a power of 2 we can take shortcuts + mpz_mul(t, a, b); + mpz_mod(t, t, r); + + // ditto + mpz_mul(m, t, ni); + mpz_mod(m, m, r); + + mpz_t ab,mn; + + mpz_init(ab); + mpz_init(mn); + + mpz_mul(ab, a, b); + mpz_mul(mn, m, n); + + mpz_add(ab, ab, mn); + + // this as well + mpz_divexact(u, ab, r); + + if (mpz_cmp(u, n) >= 0) + mpz_sub(u, u, n); + + mpz_set(ret, u); +} diff --git a/source/montgomery.h b/source/montgomery.h @@ -0,0 +1,19 @@ +// +// montgomery.h +// lib-gpu-verify +// +// Created by Cedric Zwahlen on 14.11.2023. +// + +#ifndef montgomery_h +#define montgomery_h + +#include <stdio.h> + +#include "gmp.h" + +#include <assert.h> + +void mont_prepare(char *base, char *exponent, char *modulus); + +#endif /* montgomery_h */ diff --git a/xcode/.DS_Store b/xcode/.DS_Store Binary files differ. diff --git a/xcode/lib-gpu-generate/msgsig.txt b/xcode/lib-gpu-generate/msgsig.txt @@ -1,32 +1,64 @@ -0090D8325E43E77A1B -009A2E7A3ABAC99E3531FE4DD97E8823EA04FAE925BE05F8DE5D87836FEFFD7D28FECFA41C573BBAE99FA26333E66B7EC7B13EA893395B831221EEE1BAD36F6D75F403FC5E9F64A440A39D34E756D33302B9BA915FAA3A6D819A40B0AEF9A5AEF34C024C45B68139BCB2D4970D538E2F1416482B406487CD49A15CE7F73C7F3E49457D0F27FCF44F896C340AE6FB18E378C54B0E99C1BB0D5B86F88C9F5739CEC50B76474EE4E8D803D1BEC51BB4CDD00ADA20DC267312BB305EF8DB7D488D32D08940853FA3CEA4C5005195465FFE411F3C08D13EE3B84B0E3F436A2162EB5E373B47CE0C4C10CA9EA19003412420E8668659D3012BAE5F373720ADCC11EA1497 -0D939385724C504B -78B5660DEE2FC8104BFD5E52EDD024EFF5C0692455BEC165BA758238BC9FA05E597AE246366825A736673AB3F0922BF8E040343512E548BEA1F4AA54F493FDC95198CC082271D986F1EE86F7B8FE38C6BE91D4E772F583374A4AB6CB60A518D0EC08A6A0C7BC43CF908D9E97D4D1DD578E369505574D9BCAEAB64030F7F52984EEBEF3F71B5D51CFA10989CB1C1E2250B6E9554A5E5F402FB845C2271662FED518FC8356B618B3BC299B0614E806DAE5BB84958186E9A1CF94E434D6F9135088037BA0FD1F81576B00F4F26648E15AA113E8A500437DEDFFE7D68719546A12A4BE94072D746CFEC4D320099992FB831B5E57699E71CD439BB06F7FA5263154D4 -329BB45CE74502D2 -1E6B8D60EEBE802A3603E0EF51622A78E7A0D0CBFAF932B953A0EC8CB0ECE830495D8A6618777B5439AA558FA2B9AB1F38638197591CADC929B892F64E129F5D584F37515BAFB0728D57B77F0136AE07F332964E6F53DA576859BCA785349291A62EDEB708FBA386E302B07C06A8DE9D83C717B14FEDA82CDB749B2E9E622EDBA67501BDF1B05D511D5EE4B95B9B764F584032ABA5685ED3680F43095436951596D4190FA93A345FCBE29599BD073990E412902B98D36122CF70BABF0D96FE35C818AD175E667CA6D0C194ADF7BCE6ACD6B948C10D215E3C9E901511D5318997D8E592B2627DCCB7D897CA54F42A1418EDE4A2C510BE4B83246A5BC02E414A5F -4690F9B10DC84BA0 -3833F1D92E17CEE8C71B29CDCD9611120CC20DAC979536075756F4B454D36E931A74DAF91B9F09DCEEB8DF1704F78122AF1BFA4E512DE3FF6D2245CD332A7713FFEB68A4150775CC42229A0EB0E24DAFB141D9698B397C0A7879C75975D2743A0EDE9F26CBBED21FFA121B179460DB81CA73EC85FA3BFAB5C0C5E058A2F7A1D202C242C6CAB84ACFEB6FC8C20F3B65BA309741127FA3D03C285FEEC3A20854F8B70A2FDEE253A5FDB2EFF61082CA0F630312B7959AC4E38AFA0BA82847450A0F7595521FAF14EA163C6BB526851922FC9959111111ECEA1C7EC13F86AB49D3521A061D69AEB40E3E5C6ED9572DFAE4E9A08135D23B4293AE11D940116A47FFC8 -17507073E32A3562 -121B51C2C8FDAA14BB2157C62F57B137B3777C7034F58B180437CD481F6183FF8BC76B631C0495C26240BBE532B6BC6A68C263B5593317AF35E45A03DAFE8CB154DD14EED8DCB6AB558BA001E5B4751383BB123578FD6926449AF09C9D69C376EFFEA35AB16D866D0A39DE16DD400859AC642BB5BADEFB9895FFA778DA0885C15FC0F1B3D6747ECE8EAE00F395CD42F96C4A7C0E757C523F3E0122FE2D1B1991B8894E16125662C51C9A75F64D88658BA6B52F25795753897555BFDAC950995E26218C884474E1AE74C9B269CA8C5BFC5E4B3FF23041626948125316AB6F0BBBD14C26D939F2395342FE1D4A3AF4857D39B9267C94E5C50F5424E17D882EB5C3 -00ACFCF8A7DC12E7FD -009149648B0B020681C7BB7966D9A5D95F9B96EDF74A9748ACA85D4D3463CAFD0B956EA115B0188C8AFCBF3A777041026C3AAF8312C65EAB1563C31F777D4A370C1E3A3690B94DD1DF44DF64E0D4957C8ADF5C8708DEE96DF601476DFE751ED209452AF612FE065749CCA16A9876B96B596CF24826C619706AE7FC31C06E7B726A3C2CA9F6CD278C68209EE9BFB40A6599DC2137099A9BBDF8E434617F1A8048FE82CC6990710B7E347C3FDEC7766134F5F2B4D1B1BF8CEFBFE8A902624A50EF90FC4915D8F443BD6094DBD367F81C3BB6043949406B783A4D23FD2732568F50490F4878AE325C60891BE2F18190B37E62BEA9C88001A3A1E66D133FE6D9E59EA4 -0098C4901E4E475AB9 -494D83F5AB6DDE76517E2A610F7F7412FB0FC44B2474F3BB404102813CDF7675CB16F4BD193B3B99EB2CAD063D381A7556C389A9147A4338185DB0A25F80B8F4C8A2247F1224AEB08CFDC3FDF47A32702FAD243719466DB3B290B4250CD98C7F855337BA8570A095556BBAF0AA9AF04F7B70ADE40A8C785C7243288A584154601417E6262276A291FAF715AD89ED01A9E102EC81C77A5770B14D731D32F796CF41ABB9D138AE44B09289AEF28B126B4D61B4BE58DAFF139859183856B6CD7741B31364CC62906023DAAD6AEE39CC310FE3CD7769E8835D015B82C9D8DB5B9294F52C88DFE0371BA3EF1A5C75574678AA66BD55A4EAB3867872751833785853C4 -00DD4ADC302D41DE2A -7671B379F7E8216D159598B28F78055D05A1313E42446F7B818E0C0E98B9B5163127AEDF2A35EEFD27FEAADA7A5816AD1468C113EB69E0F8E2C2690999C6709174A70F56D8E10307F5DBCD1A6745B057C29B9F5502FE1E5C8E0CE82EEAF9FD53A144450213BD94AEC27139D4D0D8057EA2266B8CC472E8FBE5818660905D8259609EDE3E0709243867177894C5B99208C1C9640DD77771668631446373F0F56F47C6994B27E9E9AD77EBE14FD328656F1103D915FD5A8C35CC955EE2AEE8DBB8D7C2A93440FC02AF8F908D9803EC0FF0B59A2A37E264B5C439AE34E0BE26691CD19F8E508B345843440332441327D32B9A8F984CAC48B83139144ED183C8D537 -11BFF9FB7C1681B0 -1D65893019C63343A22038380E1553698AEE82E895545AA24B8E6EB3F6FBD8ED2FC940C3EA33FCA6F990F84C14A0145D3DDD42E9E30FAC93D3F711D32AD6DFFE3DAB594E3666E3891DB9007569E998210CC7B2363338D8233004AB8D69EADC5DC68B60D398F803F6BBEC24F1C6458479DEF087E6D0CE5C659E2F94DBAD85DC6895D7ACF595B3D3F90B1186672DAB0CFAEFAFDDFA39EC649DA18830CD2EF2FC49A3B5667E462B03599D3DAA7B06BA663985C2BBEC53E9B71E546EAE199D806CA41E939B1F3C149799DB2131B5B88884F86A760429485A40C9F6C29A1ED5CB28BFA8C48E66E69A10925FEFACEF3F5A3B244DB9E24CD1D88DAA20E76A3E60DD00FA -00915C5BCB718411C6 -0087B8189FB165EE2E9C1CA6ED4A7C0D98F33A4CB6724BE52E60C53345C546B22205F462FB803DEAA817146E3060935F015509A3EEA30491A92CE6B970980A11E07E5F7ADE2FFE89F39E239F7EDABE36D6DD3D52ADB1D7595962E074B0B08E8E097D34E862D4BEFA090EF94B6ED9FB1A7815399CB2DE82EF9B18B14D07DB4AE96B3ACF29A98B477B29E267F8249F998AD07C1AA59FC67423B82F497DA2486BA02223228BDF663F9274951F7974C3AA8A3285EE5EF2D07D9A34E716A62CF111C683E8C889BDCBC5BCE685AFFEE161E36B974EA7F2D6A24B30077665962D11F23B6313EE53352384559094E60FE7D12F083E84A8C1400FD08F03DC290E095B3160BB -592A6A94F317309B -7A7BCE2335533AF60929795A3DE8A56BE1CD0615CA11005625B216014908130C9F8130202975BD5930B5CC1346ADF21522579060BDAE4FE6765FA278B469C0365A0C514772D4E2BAB4A8BF92716840FEEEB8DBDC00A90D10DE8105D1E6A5FBEEF84706295B02ED0D96D35FD80EDF6F0AE35AC58D988335D1AB229E93E502B4A9DE723661D9FFA923E48282F2B74E6115F02E6B5174B7AB79FA4B3532CBC27CC3677833ADD638E32B7CF2C6EBE104A174BA94E508E9741285DB9BE9EA70076838C3AFDCD71C81F171718055D51A714A10FBF0AC16B8A7F2D0114A9FC1DF65CA873D472779057EEFE92C192AF2728A7004910C756C9BD93028C20BB179DDF215F1 -0E9D0414A05989D9 -5A2BDC6DD3013DCE4E15FE50424CFC275BE3518F56A592D4290FDAA3058252B51010CDCF6B1455F45561CC51B916A569BBCC7EF84B973B6C1C4713E3AF0B2EC18A3A6EE6067C16131B73DD279B452F0FACB256261D9C1464E20A372048F51E881F1F1B9B828C31D9B1EBE05EE26E59DF68CF8149326A0D6F70B8A3A7BA6D993440DA3A6C74C745B59F4DB6B4542C8E51113817D7A92EE8F1E7CFF67607F97F6B9DA318247EF9EFC484BE7B8EE58B3460889BC490B0153306F69C6809792656E500311953DF90EA44DAC1E1C4E8189E83308D9443B8EE89009B0BB214B4F80E9817F217CF7AA071B875CC784B88866E040FD04A450BCF5BB0125DF185D579AAD5 -0096AA789DD034DB92 -387AA57B2027DB04CB8C7A3714C55228F209026EC3442AD6FBFBFDC953C6E1F3782108E1AB6FB92200992FFF13F5B41B5DCDE57B8E0C4D1168C65369BE641231244690CB5ED6B121CBEF3D08D7C9376938211A67B886AC35D35DF45F4300B3B090BEC933005DC7DD7B63FC2558B135033DB44DB7191791EE811F04F4C8830956DF10F9C23BB8E867FC09D62953D9804ACD95746ED37166897AF0B616B4E5F68F8A2C8D3721FD1EE2087D6DA2C294D7F76D9AE275300385C92D03B58E67DDBAFDDCC4982A2E3EFA2491BFBF110CBCB849234EBB3B8BB51BCCD2E9A9F8285BECDBB43B23FE58EA6BFDF95E4BB570A4FEDB184FC70F78A74DFE257E857C9FCE2772 -0098FF0D98581E7B2F -26F5D1CA1A9106EB39ACD225B74DCB51FCBD200237B8082D3ACBC15C5BA83021806941DCB4BD97BA58D1956046B52EFD7D8023A187F30CC46F0AF56ABDA373CC55DE92C70896A4B7EEFA97A4FC4FE56088BA4E73450616A2F234BDB3039FB5D7C0475EC6558ED9732D2CEA9C7DA5664E2218DE210CFB403AB4A74459ED2ADA3273362632BCCD5D6D3A36837510F28CD0B14C396CF3FDD1ADB0FCF6071B5C71052F19594FBDE997BA166E5C445FE2285EBEE4773B623FC300EE0C370D60E595F638813DF105481BA2009D386E59502016768916A9EB2982E50639F636D0DE7FDB5EF6613711796F62ACE42D12F0900C445280B55A81A0B61B572EA34D1C5CAFA0 -3012F2FB2F58E6EC -71BF267C57B5847244A217B02CDF8B5595574EECE917C44F0A9B6723A2945AD138644A25CFD89F9700F337AA42C56AF10CF968976AECC5D5FAD438F0B86417208511D43AB73BC45C8622AF5409B49184346A667E9303FB872F6D20C2C2A72BAD110CB3422E33782278D452A48B7985FC67CDEF2EEFC59BC25C3F4614E66E8CC82207574BFB431B7BCD17CC9D75D47DC86DA6775DAAD045B9DC6C5B80003C23FB766FF14393EC38A4E0751E0369C5B3D416E5648FCD715C2E1850974A218267410B7454AAD27E9B9090361EA3DC679A2CFA93EAC15DA4B218C6AE04D79EA8BDFD9EAFC9B08BC706F7F9F0C19E2FAE7A771F02195FECB61A678259081666CF066F -3C2BDE5A818B5439 -08191FAB1992DB18592519FD2F36BB86034953AA143816C93699A34BE405A7B0183691205636E69DFA4A872E1000878226F4E1A0CF030023E09BEF303240174EDF7CCAABE052CA0EF397CF0AABCBFD2574C93C39B21659EA0CE4BA9AC6BAFDEF0B7C3DB4AB5927F9F0C2CE55699B4E636AB484CC68BA5F4892012CC36992D82CAC33135265133F532C7330871A32DEA89E3363D5E1E65CE047CF296A5CA840C93B86074111E1CE1E273FF23C5F5372AAF75C1D448193E14DB99BB8BBFB08FED72DF19AB78B1CD030C3679E581F46DCBBF4C193BE00911C974BE03A51F2705E59E4484C7BBA118988FC5478E75FF2DBB247D49F4F9A706D5F0FD73C0AB68C03D8 +00D2165BC2B4BFDD64 +5BD6158BDE0AC0655B6FCEA57994011D18B6B3C9E5FF75C45FC1E5EC2C1F26D6AB8547A17C0BC15D40F4346CFE74CF4EB417E6850D45C3B49E9389DAF400BC5E5B3F5D8E1E45A23DD042A87E82703209F9EA9808A002FEC00C96A5F0D9B7673B4B0A224438D81C0A9CEAD0DD22802B409230072768E73688D63EAB1C9BC242FAEDCFE0C8478B38254BAAC07AD6F82A27A0C3893FCB604BB57158F9125027AECC91D55B364B5C2BB9FE07FB6AB69F5A65112A2B7D5A805CA9B2C1CB75D315DE345BA68100DD5E46FA3BA54B614C298E60EBAF95CEC738DA2513736ECE051D153CECAC29F4A432A5FEB287E2A1B8C4640C58FF9E9E7DB6889E4865D1F1C8CF4E47 +00D138A9266AD7D5BB +00B8B50049C1A29E9ADCC23209AEFE7FFB1458B554B98FD244F4A658E5D6C6F61E42C2393D0C59B5EB4AB1989FD2C93430602B61008BD0C4059131E35DC5F28099BE8CD55B1970627D8AA3070B066DDBD4E2B20584A35BF2F54E2F2C532429BF2D709CFF82DDAE94D3EE34DE4E6D348473C513C93515B90ABE0C6442C2749553A51B991AC5A70DDF21383DC2120B9990BD279CCF0EF0DAE177CA8E6F8479418337F05E7B1C9617092BF143F1A16414B1F268F3C531651E2AF6497A1D7F219BFB56AD875C77D8F6F31E27A1029F51BF819A8F4E5D22334AFB67873CC7459E0F2284490F675D9479A0C9C50E78765EA088977FB10871EC39B38CDA8B73A8A7C4B893 +00B894254930C49655 +00956E3E7B09F7FECEF26CA44FFD69F19DC8DB6C3A29A707C2CDAD56994A58D6ACB8B275678D0D8670D3C716AC5C98398C8067943C7292F787F5451E8202F4C8BAEFA6CA787BC79B73A99CC4C85743EC7320E17195D560A380356A9D32AA81EF276A9DE8B9F6728647851AAD0090A458FB928BCE86884BD7CC7AC3CF226CE546E596135A948B820E1865D6A3395DF2BD5EB26FE5259B2B950CC61F887C0D5A81F77549D8F792D32552870358EC5B2B45552C35829D732CC1A08898FD2FFDFF5EBFE0BEE7D5702FCA240B377BFE7D2821E123F2A146725D01A5CF0A6C89FB7E73CA6F3B8640C44B0FA1A51B429BB3D4668495F20A25FB4185831C3B479C5041713C +33D86BA3A9A40EF8 +00B91210473412D5D2BAED71B445E2FF12BB6477119409AE029E307038E781B942D3717DBF580C10A0BA1C6F1115E699E573C643E3FD30A3B08D02E22A20C34872E841FB242455A92178686D11D33C63BB9F6A276B006BF494852406E275A2965D00714E87E6782DB11E9EC6D5108A21678B8DBE5DAF871792E13EF0F24490864004BC76B978E1BFF178BCBFAFA6C4496DEFDA235922A689C83D84BA4C04031264C4213935099A374EE94AE649037000856989F48C8E7E70D553806CB4919E7042ED1C8A9045080A64341F2DFAE15DD5EE443AF7405D8BFF70A8D36CCD838CA5206FF6032709D990D9C151AC05D8A2C57E63934A31957D2ADA44D52BEE2D1830D9 +0098407339FD84FD9E +0087F0BD6130721104A5F15CCABED06660A657DF8EA7563DF7034E95DF721EC845D3B6B9F9BADC4DA99F0FF29C7FFED88895C1028EF5E8C456699656B50B50D7B37617E282E7D6A0C41BDC0F33D64981D8DA84AD43D2E338D21010494E9A622F611F295816EA99B0D5270A4A7A0894311C9C2F2D23E6FE85F96648641294B3A28C1291D4257E6CDAC2E4F5F37CAF7243E02CDD02A4E2B364E4BD70FDB3B102312F6C369E5E9032327E49681F871061975427FC749E73CA133112553B7B903653310439186CA8E45555806A71EC9E0761B8ABC91D8A11ECE3D97EA98E38CCF16A955E2BCF0C9F395F4F670FE3FA39F4BFE44CD9C71435E360B74485A1461C9B561C +6B5C23584BB6FB94 +30E14CCAF4E499F60E1254740D02BA84349FF9C936EA12E95964CEE90079070B966E92DE6098EE4B0075C99109EDAFA5983351FCC4F1E78897E2EDCD4BDD98D0E5FE0BB07F1B69844B0372D344878252D3EA2967F593E4F2CD71AB0F92A656DE57BE8DCD54795AD826E9357A4D76C447C4ABBFC45BBCE3F4C75CEF4255B64729204D225C47C8C3E3050D95AC7A1AEB5A6D90CB624B7A47E59E3ABF6AD95CFBC3BD5C2499C1F56421250E5D6610236FB8CB87650C1F49A9F948234310A80E42F86FD7D23DCCF90A73173AEAA42543EC789BD7E233DC48177A766C8B4C9FF2142A0BF59682B7B798572E112F545942E493B26A70967737DFBC82D1F2C6A461EE51 +78E5B46E8DD7DF7D +31862D27DB14EE9145BB54F1A010C9F86E282E5D154A71EE64A2DE170F31D86BD2133DFAEF4E75274F1C5DFD1C7D65E087C06840DA0FC055AC3C946117236657E468570D1D9620CB4E27EA063FBFED65DC1F9D87FDB7FF3C071C81E8ECA9CB5A0E373CE6DEEDEE8F4363697E1D5AF07F594EAAB557F4EDB723B23C03317A7629017F8283AA3829BDD11D59EDE0308CE8093B90CEF7F140717E041C45096045DA5D0E9B8B4F179B8E3B86A523B441E8D17E424960F5BAC2702C4245F6C5D5755992896CF4590982DBA694988A3A34F78250485985096081DFFD3DF3B561B4365F064335ECAE330D6CE7196F5E30E74D0DE591732D979E9F0244D5E1726F5D78AA +57C2A5E92EC6C384 +567DBE71A29DC686EC6B20AF16F682A54FA93ACE79C121103D21674274EABD58307814C51276581C1D5773B2AA13D89F2324307D85ADDEE29E1ADC2DB060DCA6D6A595387D1C377B81720A4905931CF702D7D633E965CEE7955DAB3CD2D9DEC9D82F84C1FB50F37177C25264A340E02193B4AC2F5275F4EA0FB429A5B3D7D290121D6D04B410C44F51C11BA47761A0E9D87C1144DA5A537E5D7AED5466244D0F9C1E25E5C2463C894EA55B908D2226F9FF072A3303619B6DFAF8D596DE00525C1BFB7B90BD8C76D920A00BBF878507F002669E1FDC22575ADFE3A4CBB6EB11A986A6821461913EE4F4B3F4EF974495DA6DCC74A5882496A5932EEDFEC6D2E006 +53C04C903C3CC3A5 +542E73578472DEBA4E45E4A80F822C8E24B4C1B7AF23D93A4C4C9A47AF551B73F9C46A6E959FCB98EE18A1C979B60FB3CF86F93037EF8BC3B8C7BE113D0E7C384DF8BF09D7513B7A60863AA61C7BFE7BC43D7D983CB5540CF32482919517032C4806F0594577E498EA624ACE57D198232996BE010E58322C2D26A936057971F75142099D2A4988E971294C72481FC523797F681821FE10D90EEA14E047B4F10F8ACD975C01CED883EE4ED06C6371B18857F93E1F10E935982F5BB17356BA28B5A4E5F4B4CB4F77976EDB7374FD2C93EDF7370F21E08BEC31AAD7F1178E766121CAE57D33601968AB7A85CE4A267423C5F249B78987AB2A3B52BC2DB47B3CE2 +00FDF2D41F60305BFD +00AFBE78D3B8D839EB90120B4DA3CE9DBE6C24366B48666571E9D19720C8D0429215C5172B803B281734D8017A554072D5AEE58775801BF34AE8EF71AD2149B0B80E10EB9531361D1C99D25DCAD09CA721012EA2BC550F8C7738993C46D46CA935CFC4E9F48D9F83C78F325F16FBC1D9877CC6BF2FDFF86723F9DFA3F45D348934442D530F33DCC6E1A566902AF8566E1E97091F0E626C459773C5AC37853AE57B1BD7C3405D46A9CDF08E3C76F84B4245D2FFABC94DCD961593CD5C12339188F3E723C3F6F87F9AD94584FEC35C0CC667AD79D2E8243DB3B11A9A6DE9FBCA264739D81C8A0A3DD1C35BD1129B1447C2899AAD236C53B3E95455E5B211D2B56BDE +415641A674296C20 +4772EECEE851F3EE8CF67E78411F791039476B3849AFCAD92DD7D5E757754024062C9601D8849732E02D9EBF83BDDC87DA0B13B2B08AC1788A83231ED37A6376768CCC5C464F659C70928B2A5E34E7BD91706CBD1E28D4596879BF037E2574ADDDBC70A14DC80C2D23930F3838E663BF62395089934B457A8973D2ACD5594989E0827DDC4734FE7270D2471267C1E43F4522509043CED5433D73B0683D06222D49F16AD11213881B4D0249863553F297EF3B1FB9FD5C82692876DDCF692D86D5665E627E18B19838AFA63CC95E44E3E49D1A2FC86F31BE57A0F02E9151BA141BFCBF508A2E02B1CB843B5AF5F3A3F3A0C0C170C7D37ABF242E44B3F93D136920 +263DF0F6987BC311 +0088A7FB48BD834CC4872EB59FDE7662B2385F7AB342E402429FE31A8B704B032121C54070A3A090BC186533EBB3FD344C247F45F1ABCA8EB135A23C12F57E6CE50FB60D01862F570B114A60DA504C63A21C7CDF0E15023F642F39178BD9B7A4C5D37D40FEDA1397A9F4B46EA772FBE5DA21450443E86736A5AA88B20E1DA3C12D0BF398C17D10B32551CD04397A0441E69123E3C069503940349A311DAEB6541862142CB7768C32E48095EA388C2D290973DC74EE6665770423ED14FB53EC992FDD81C7E096F19DD5B97F4AF20A4A9F3BA2A83351830FD8EC4B00CBE477889B4F84518609C1BBED4CEE0CD784AD53789D62E5E05BA9566D2598D929C95C7887FB +00CC3DCA62E3A33E90 +2E47F9D887321E63607CF95511A2953219A31EBEC359EBA26EFD3F4FAC644E7030D36CEA21DA47898DED96D18AACC27C25F542FF4A9A89F30130D78B0383F6AA2E8C104F7A7BF762BFBCD4E28CAEDD6610F56B5E09D644EC3F30CF842E8E9DF0522DA81508D18FEBA0C69354CCA72DADAF2E7EABE4D7919A7E8746EC67A623F9D7052A9C5021A6C03E1F2D4EDBAEE2CCFF2E424864DFDFC0F53EA80C3E03AEB26F311AD072105E1487497690C3E729426F234396375B81C2856E4CB5C2B1F6A436325776B3DA98CDF4C887CE6F4BF5145CFD53F05C549CBA8E7E5C0EAC54FF8B8901520AD25430202457CDE7DC133AA50352CDF8234804CA9C227A199AE2CF72 +0BEE0FAACEB52773 +53EBD9D4A10C7A66FAC2793C067A0F83D0659E4854493817C41BF26133E45DC40C7EC6DD9A378A8491CB1DDEEF708723FDAB86CBBB2744558BB45BDD72DD97419280CF308E8398CBED5E4AE631AE521E1A76F3E8C9CE8ACC46AFBBC71DABE65C8AA68960F6C6992D2785B59CDA0AFB9EA209D13B7BE1F10F7D61744B452DA69149839E578839FB3669ABEEA9867B53A598B92D24DD2C91166417FAEFF0C7EBBA3798F9F80646443E920A3754A1052F6CFCF21D28A1EAA80EEDB73835A4421D07CC02041CFAEF9380359591E788EAD40ED9533354C61B21CF9F0A9714C3C3CFCF20FBB9D6FDA2F00A10C0ACCFE16A5ADC39B9909BACF3683BF96773654883D303 +00A1A34C0151259F57 +429E47B51AF3E6B7E556A43D0A3FD42DD007631DBDBE3CC17F39E9DD3997BD81AFCC934D97F9A127E2F7B0FD90C9F46EEB403084B21E048E175F9F0EAAF10D27A67582B24452E56C35A01A97091AC8A4612AF24CF39EF3B0172E6E0D4FC4012745B96DF481D2E2AAF6460B70A84C519684038CA1C8B44F91E2AFF20BAB17434FC6CDA6E489ABC5A42AF2767F6D8ACCB462A47BB045AE7975276C78752BFA1AA327E6C788596CA16792E0CDA6A4E9D7B594E662CA907585508C71DBB322710210B0F8B4A5015408D27EC39D966A44A5814DE9F62C2457FAED55E880EA076E1F34A3F1CCE81EB3ACF68C52154398C8259F76B84F7276CCD640872C3C6788743EFA +00D114E792C580A36E +46EE4BB7BDAE8CBD2978A4A4FA7C03FD9A96E54DF864887442D4621AE3D4E0E47F37AEC9D89DFDB98B4EBDD742A09A48EBA2C316E6F338BD6C88901172D3E074E4F74F31A28CD4A387ECD7EB206C249B0647F6B27249FF213F1D6B3E774A74F69D597F5C0A75D9B65A96E071CC87A2A082FE632560B70B47907DEE1D75C7164C786FE4D7719398AB0E79FD2258A6DD5429FAD1E887E56DF670AA3CDE6110952B11099A4616EBE9AC2079B02F3FF1C2AAD153A6B8104C65029BF45FB44773BA8160E68BC76FCF21A1F7102889170DBA841ABCB641C95CCA13C3E3BAAA46C24E4A784EFB2884BA68FD476589BFF3BA3E1B173BBB3521EAB477013D63B1193EAF73 +00FA91AB8115B08275 +20D6E4F2A0CD8C74CE580017C8B825D03844328C0D2A97AF74B0F7160CCFA194DDB3329CE5C89AFD839F63C74BA9B8E2EC9D3836D0B557C892617324371B312EDECD6FC16341A5CFDF63BAA8B8F049D591EFFC70A827E907D86E6C8079F0B61ADFFF6CAF909AD878ED564421FDBA4CDDF7CEEA48CEFE411DBA88FCE8C6581AD86F61033D279D5D1F481CF26619DC9DA7E7E2DD60BA5B357799BEFADF15EC75B05346E437C240BC922FB96EE29526099AA001961320734A1A1EB7D108DFA1D65637E608872758B76BB9329BE67317AE32E771BCDD70A5C43EA2A9A8104DD724E88192A4C231F2588AA2D3CFC273215376BCB2F44E38EF138BA701C2DA73C812D4 +00CF4B51C297D89379 +0084412EC87B14640A613A83B9837709B427442F13333029E6D92534365FF6F88E6AC9D330663D7AB96A1BF09A06701DF125DBDB8F9D6F832835F34CCB4CE5D518D5FC244A90C5C549C9E15E4A676440E1954D74B4B2745C970A1B6CC2A558BA51E62FA140CB0B783485D04A58539167D5EC54CC7523C4932374A792206ED5CFE1ECD3CB79F44F366170644FCE1BFB7B2DBF58B84CEE4DD99BD142061629AE789AF8B79847E5F2F278BEABE8F7B2D5BCEB111B56BE68472126BA9E0D17662561DC19A6705142F26AB24DE2432773E147254B7524102267C8A959CE81F534E831CD8A5B84A7E474853412AD2897BDBEB932380D8AE1270EC7E2FB31EE3FD4C9D652 +00B7A7A96D3C82BE54 +4B527278114960D7C70359B0A5EFEDD02672FF5CB2C9DE345F40F68635B7D82EB43AA4C275B615C9402E4CAF4D218AC8457614DE965F3C26FB81A8B02161FC5380D70E39ADEF14A7805DE02880793A0AC2C3AEE40C58D221DB228C11882D5C03FBADF82B91646034272EEDFF062B5DE18FA68ECD2C33F25807A4CDEA6E8397E23E19BD1F1570362E68D9EF600127426E850768CD48DC9D8389E494919490D9E5D937639F56402C21E937A9DB5060C904C362DA360DC828896E51324D4966F979D12F25F5A37151A2A4B57074C87E00C28BD32F93525510CC3B57BB963918651135B4F9BC77AFC36CEEC3B8D8F3DE0B09D8A464DD3A8A38346C4F6FE8EEF74219 +3CBF3A34FC8300E0 +00B04CE658C3DA18D2BE4087B7A886BE7A6DF7193E754D2CE076148FA1A185018CAD81C6A6C9814A5A5DA3665570D6452F1B9521B1AB5097636FA05DB42AD8A3AD58CBC32A2B5998F2DBE40C6A4AA3DD758A15F9A706102081142BF3AF3B4556347CD8F3EBAA7C206FEAD275B2CDB62CEB0DD5F7545F6970D5737630E3C8A920998F1D1025C9122D40EC120541E23E13688905FADD85E2A0CE73183331117E4E4474EC425D27C62E33CD9FF6DD2ED3B49A8D9C2CD734E110ECBE3DBB562C0BB2C1221FEDBBA986DE094E161B7AD69C69905AFCEFA8F9F309213B14583A4BBEA175447BEB84133D1EA5DFCA9F9C03757746111EA13CB9B3259B9095CB485592F6B1 +395F33A6FE2E0799 +0090B6C9A38BB8C576908D1983111704C00626FED54AFBD84671054698F2F7952580EE084903497ED89B922668715AFCC8F61942A1A287CD4BFB78286A79B96A13D618FE28C8A3CD11E583F46687C39EE3FF854D00B90F217BBEDB592F25435F67CD2ECF0BE5A488E7EC134F2BE655D2E0970151DE58C1F0EFB3F2AAD0E05AAB10E10D70DD7F6A61C74509371731A0491E9729F20C4795FBF54A242A24B8D047F3CB8D90B9E6A0D004DB88A768C1D9C969E08565BDDD4E6D29B70E266DBCF5F6EDBF8E3D950EA1B110CAFC6CFDA79AAAC5C4936C41A2C7FD22E603A6643857FE6F804B8854DFBF98096AB3791F0BEB4B358E7F1FA37E37B86CE6A1E702DC5F9111 +00BE8597F14A024902 +1CE840F10D8578729BC1507A833ED36E86C505C1561FE56CFAAC1250045724268AE8B70FF3BB49248D9D8ACF1041FD7F5680177289B48020D16A38AC69BA78D26054AD4BD22C6E7AC21DB00FF2E91F1B49F8F022E7FD62F4F7AADC011525A18332424F37D90860887DC2729519D3D2A962805B98809D8BCBFB12E60328CA24CF069B4D8717A3FE465DC319D9BE6FB93A1B4DEA2E336E32C21D27D347202953BBDC74BA3BA08A8A91138DCB0F97EDB9866F141B76E84FE781B19FAE203C24EAACD5F6B7EC5DACD126D714D8A49A7C530269389FA8C83151A2F5818E7870FBAF35C3F6CC31AFAF8A062C4FAB6449D8C3C4B22B1481D1981617FA5EE17AE390FE56 +00AAB58CCC11F3D7B7 +797AC39B0F1B318EC03ABBDA2B24F92CBC0C4948ECA63E202851A7E96937FB1A0AC09C35DB5588C697053A7DC7C431EBE3A8A5FEA14AC1A3E291004B7350C2C6FF85B0B2E5018ADAD94AFE3FEC8658F22DDC51DCDF73AAE1092ABEF1536EE0625890C11ADFD48D350BAE1F95F705DB2643E784D632C603C3EBB5BC5F0F6D38D0F66F0F1CE1A5E0BA8AC94B6D0899E04CD4F498FD2FBD5D2B191C802B1024C5557FBFE8A1CCC87E4C6B6DF8A8243E2F4CDA6A5092D91649C0619D2150318DB09F55622E2799CEDF6CC0CA6CE9FFE0DFE7BEB1E47425116E7655AEE255B3ED75A4A96C109D76853688B7D36311A24BFABDE62B7641234504641985D7550314188E +00B1E5C528E3150A7A +00A6145DBE9F88914C8CAA1E49FF2B0B02C8A2F4F3D8B6A0E48A2A5EC9C6DE9DCF7E61AFBD0127632AD2FA9296A94804A6BB3F9B30F8C528B7AAB107754F357E017AF8518BC3889260E71139930DA6FD79EBA7FEEF3DF3320C31AF7B4DBC023B86CD79166120FD34D152F65F4033932E40A8E9D854210A154A551BC1FDD1B71FCAC7B5403FC29A51455D88D6FB0BDCA62444D0C52DF71400D1465A4FD3863A29A3AB853E7EE8B8D87E0970FEE88FB22A85D1B83D50F036ED6FA4EF53FBA91B8BABB94C19F2F8E8D49FD2B4409B427C8B856559BCE56A8C78CECCC9D97E26BDDF4516280DDDCB4E1CC79942CDC9FC336D4F6211025B2571972AE7F28FE063DACC1B +009EE7C919343837A9 +00A566F5D0AE53834855C4B1E6D4B46AAA3C591DD407C10B92CAE646D96ADCE27D7997E5CDDA348862326AC408DED65715EAA72E877A71B87890DA0860E6A1AAACB1086A36908DD4F00B9F3DD3EF777110B21CFAD3F54191CE5876026DEE11E4BA196633A337E86C00413951A2A682E97CC45E68F8CACEE76E17E7AB32145081FB502B2E4618400AA85D002A78174669912F679FAD3EEE7CD406A85EE09B14EB6C5321A55FC613DA7BC0B16623089C84B0154F53366410AC4500D610B39F20F642E80BF612F25E877096758284FD4E13C4A58A0931776810F05F203CD60C334702AD21D20469DF409E4A0EDED54C0CC21E1392A7B03655FA03625140E91A862575 +00DAE3025E757D432A +4C8FDE18A09668E17D9FEC4B74DEE31C1DED79ED606938FA12D572B0F6CFD53E03F5C0202CB96D7DFD949D34C959749C4970C2445A313167DD39DC14C4C0F02A97F6E25B03911876DDAE2680EF995082CAFDFAA9C8D92CD1B98E744AE807E89F25082DCFBC5B0B1A525CC5C0585F47AB5EE72E1E442549C7CF020097A0EE756686F1B8D1E88BCDB4729B74A01FFC686CFE3998753F6D874E2C9C79C498C3F4F4ED2B74A72D63D75C643C40DF2826F9465786FDC949A9552F7A162437C6C3CDB94CCBB94DEA709B424E985936337811A89F149ACB9200443F6DFA29E8F2E01DE50893972D238727808640F01E7421005DEC62AD8635B2BCC2391EFA048C950920 +7D800AC9BA2ACFE8 +5CCA26CBAF395F4B0535910B04893DF7DEE481C24AD589F39244857C7A4D3CE8BB09B37D6A5B93725FD167C5BD72ACFF404441AFE5203856E8FE8E55E54A3339EF25C41D5168FFBF0AA8FC6903327F2F26D37B0A499C67ABF82B401EBF7960A1DD3CF5D9991AEFDA9C4DA1E641A4B6DB77110B83117CF8310E522515DCE5AA120CB1B6BCB72013DF89DB54C6735C3722FF109C1C921FA7CE17198971CE06EB6F1014A4C187FB2907C7FBC3A9CE61393EB1BAA0C5DC855BE8CD820FE039757C98F282AFF592847D807AE969DC0478D67473B7828BF9EB59B9CC71E1CDD4652E3409440339A47B9F4E40CB7E0396F34063428CF8127B796641E2B74CC459DCF1EF +00DD29F87D7227A4BD +586153DCBCC20CFBA56319E867837069DA55B0CB918B3F1F9CE2E152EF07E12E7BB4400D4E74B507969B41D1AE4C13104FAA470C0B60BAACE91B67EF63612776AFC1D9A42B1562B675BE3129784A6B3D2E9E9AFF36AB7EF6ABCEF86309230429FA7548A2A956B5DE19F943B2FCE5FD99688ACC223225F7464EC3284585226F01F4A7A2D357D11F7840B8572EF65EDF10E2562644D88A257E4D36058E492943A7B8C731FB01463F839FB6DA6B1A83A383E253501EA1246B54F25A9B4954A91B6B43D52CF0E9223A5B998DCD0C54DF31AAD8A73DBA896FEBA1B6F3BA0B7F9E429DF9557A7E013728D721203AFEFBB8953255F5AB7D299E7D1C9B5EA5255464B47E +00E8E1C308AB64A058 +0080593E917112EA1DF7A20C1C57548ABFE1F51FEF8F5778B6B96433E691E1EE8E750B3124BE65085F8F7168DC05DA4F6C864DFD56415F08B01483A7A8A3B85AD3AF923CF6FD0BCAB5880068977465377D1EAB42CF67590E6700B34B5343E34BB8A68F168D7EA060FE34E6D27A5107C863E7D85534CBE8C072E553F3D933C2D21BB7F3172445805E4A1D3DCA5B450A5863152116BEAEDBBA70037B531471409030F9AE6C8434D153EC31FAC3F49D6BCE77C17D2924204077929BD97A0A84A962C215747596F3AE514C54DA3D951E2EC2A5A535AAF5AAAAC80AB5E87547B692609E279410E76FB91CB1CF8310F7EB8F90F06E279D810236DD99C679BFDA7584F0D2 +00A5CAFE20669183DB +4D50DF3E05BD47D40066D1A991AFEB3D1CFDB86D5518A1B93B49BA3BEE34960F4FC024C73027925E3C38F483E6BB2335EB534EE607A0C26EB2A0C08F79B726ED0FE386ED20855002D1A3C9F4D72068216ED3E025AD8E644EBA77F38DE120DBEB30A064C071BC9619B1634D79A246F447D248B038D8183AD544056205DDD76E8B4E50EAC174C0E7BE2F97A82987719138F57D561CF838AEF426E978BF4DD23377C4EBB2EDCDAD84E3F5033209CD608751EF8B4A5F6E4E5C1CB400E644C20C9CFE2DDA1A7B3A6F45A0F28CBF2968B5E931995A24BDB4D6498D49845A43D57C79ADBF305C33250C5F30AAB7E45A2190F54921434BF343717FD0B1091AE87E9ED1F4 +00998A05B1AC8BEB34 +00B6A5B0D23A9F6CCB76ACA2ECFC618C8C319E0405570142F2E5222F5A8A0464B70F8C58CBA65EC86EB5685B5BD0AA128377436BBD79F72B9C7377CDDAFED868700A9AAD82F7D1ECCDEC8B44085E67738AD0DB6C2B580B69181D464E82F6D109AF2EB04E1080E1C0FBE6C538475552984B8D4B99CDE4E157DFCEA12DA81A1933D97204AF73E31DB44BFE2F6BDFBB3DFDF7A5E84AE07FEDECA8000E283A69D29EA9A20595AE446E443940FE03934E5A6D99AF911EA4DF09AAA33A726D2B444EFA864045FD9D084102224278ECCDF092710FD6DE81FCB860D86816E08059206521537342E7276231969B3A8D942F17B8F19D79675EDB44B1082A69C65F573F29573C +008C456ABC19D1CBA7 +009437BAC19C37189D23FFE631567F68789BF9B5C89593C94CCCDC6A742C3F92F14367988CE328E3132B6F48CD26EF2231EC468560F8F6DA8B05E4FDEAC1A4937BAC608A43AD3ADD7D7A0A67848960386EB25F11B120B06820653491C3AEEFA365DF94EEB3947604EDF4038B6CAC85B01754D6B359850F05B6158D29293B42AD5A88AE85706A7C8E36E54A43A6E3C877165C92B2AF6A03C0F3F9017ABD4DB1FB724DB68FA9DD0365014AFBE0322D61BE779C4074EAA48A538CECBC6A1AE3E76EDD1AFD78E6633398B1A88DB84B49CB0F46B8AF22410B31246FF315AA038CF402B9D2107F00A2E162FB27F67825C10CCA887A8378C19150600581A679262A944355 diff --git a/xcode/lib-gpu-generate/publickey.txt b/xcode/lib-gpu-generate/publickey.txt @@ -1,2 +1,2 @@ -00C1D957F4DFBBA9A0D24C766F54CBC6C46E868CF9317E398690014454218004A68EB8CB69D6C819B60CA2097094BE7B52CAA39B1F7267940151C94208953E9AE244E58F4689994D09D780EF426718495597BBE5B8FAFDDFCA158FD3594B03C0B8FAF8F9D69E2AA598AABC36FDACE52EEEDD05ACFF0A368B68302EC41F264D3C49CB0ED468D10683EA43D80A19AC113753B4561323338A6F9C7945BE88D1864CCDFC6A13B0F398F83EBE7F60D32426FC1B80076ED89B2166215DA50543A64D21E6881B18F1C2E7C54D7DDF034293D96C755B47AAE86A4A4291DBB506F25BABCDA02BCA79FF31BF4F1142D43D08BC94CF6C6572CA2818221A326070D93F1F54CCA3 +00BB5175E55C2F1BBAE52B0C1225F43385FF54B3BFEA88B42B21044328815B8742E303C843ABE76D147861AE92D563592EFD748BF2E5BE4D76793FB32FCF6B38F755D408D114C9DF89B3FAA77EDF0C9358AC3BC23C90CDAA8337927A3530DCF2AD6EFC023C96A7932F8A7935B9B3F5C84668B41FB39059A1B723A40D59A7B1BD03F56933D641409F2A49E614BBAA9F2573ED24899840585B73329A01071793332BA92A0C9033D7004B45FD01C3A850125FA2E4A40818F8E233B7B7595ABAB04B84AE88E4F7B516359EAB7C285F399A3EFF467113DDBDB17981F2F4F2DE405BA18863046570C1621AD9446CE8A3884893CEF50933CB60053B6862E2443CC8554121 010001 diff --git a/xcode/lib-gpu-verify.xcodeproj/project.pbxproj b/xcode/lib-gpu-verify.xcodeproj/project.pbxproj @@ -7,6 +7,8 @@ objects = { /* Begin PBXBuildFile section */ + 6A7914CF2B0CF320001EDCC1 /* gmp.c in Sources */ = {isa = PBXBuildFile; fileRef = 6A7914CB2B0CF320001EDCC1 /* gmp.c */; }; + 6A7914D02B0CF320001EDCC1 /* montgomery.c in Sources */ = {isa = PBXBuildFile; fileRef = 6A7914CD2B0CF320001EDCC1 /* montgomery.c */; }; 6A8A795F2A89672700116D7D /* verify.cl in Sources */ = {isa = PBXBuildFile; fileRef = 6A8A795E2A89672700116D7D /* verify.cl */; }; 6AA38E5B2B0A97FC00E85243 /* main.c in Sources */ = {isa = PBXBuildFile; fileRef = 6AA38E5A2B0A97FC00E85243 /* main.c */; }; 6AD85E072AF71AD900662919 /* big-int-test.c in Sources */ = {isa = PBXBuildFile; fileRef = 6AF7487D2ADADF4500D58E08 /* big-int-test.c */; }; @@ -40,6 +42,10 @@ /* Begin PBXFileReference section */ 466E0F5F0C932E1A00ED01DB /* lib-gpu-verify */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = "lib-gpu-verify"; sourceTree = BUILT_PRODUCTS_DIR; }; + 6A7914CB2B0CF320001EDCC1 /* gmp.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; name = gmp.c; path = ../source/gmp.c; sourceTree = "<group>"; }; + 6A7914CC2B0CF320001EDCC1 /* montgomery.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; name = montgomery.h; path = ../source/montgomery.h; sourceTree = "<group>"; }; + 6A7914CD2B0CF320001EDCC1 /* montgomery.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; name = montgomery.c; path = ../source/montgomery.c; sourceTree = "<group>"; }; + 6A7914CE2B0CF320001EDCC1 /* gmp.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; name = gmp.h; path = ../source/gmp.h; sourceTree = "<group>"; }; 6A8A795E2A89672700116D7D /* verify.cl */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.opencl; path = verify.cl; sourceTree = "<group>"; }; 6AA38E582B0A97FC00E85243 /* lib-gpu-generate */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = "lib-gpu-generate"; sourceTree = BUILT_PRODUCTS_DIR; }; 6AA38E5A2B0A97FC00E85243 /* main.c */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.c.c; path = main.c; sourceTree = "<group>"; }; @@ -101,6 +107,8 @@ 6AF7487B2ADADF4500D58E08 /* big-int-test.h */, 6AF748842ADADFAD00D58E08 /* opencl-test.h */, 6AD85E0A2AFA510C00662919 /* openssl-test.h */, + 6A7914CE2B0CF320001EDCC1 /* gmp.h */, + 6A7914CC2B0CF320001EDCC1 /* montgomery.h */, ); name = Headers; sourceTree = "<group>"; @@ -124,6 +132,8 @@ 6AF7487D2ADADF4500D58E08 /* big-int-test.c */, 6AF7487F2ADADF4500D58E08 /* rsa-test.c */, 6AF748852ADADFAD00D58E08 /* opencl-test.c */, + 6A7914CB2B0CF320001EDCC1 /* gmp.c */, + 6A7914CD2B0CF320001EDCC1 /* montgomery.c */, ); name = Sources; sourceTree = "<group>"; @@ -217,7 +227,9 @@ 6AD85E072AF71AD900662919 /* big-int-test.c in Sources */, 6AF7487A2ADADEBD00D58E08 /* lib-gpu-verify.c in Sources */, 6A8A795F2A89672700116D7D /* verify.cl in Sources */, + 6A7914CF2B0CF320001EDCC1 /* gmp.c in Sources */, 6AF748832ADADF4500D58E08 /* rsa-test.c in Sources */, + 6A7914D02B0CF320001EDCC1 /* montgomery.c in Sources */, 6AF748862ADADFAD00D58E08 /* opencl-test.c in Sources */, ); runOnlyForDeploymentPostprocessing = 0; diff --git a/xcode/lib-gpu-verify.xcodeproj/project.xcworkspace/xcuserdata/cedriczwahlen.xcuserdatad/UserInterfaceState.xcuserstate b/xcode/lib-gpu-verify.xcodeproj/project.xcworkspace/xcuserdata/cedriczwahlen.xcuserdatad/UserInterfaceState.xcuserstate Binary files differ. diff --git a/xcode/lib-gpu-verify.xcodeproj/xcshareddata/xcschemes/lib-gpu-generate.xcscheme b/xcode/lib-gpu-verify.xcodeproj/xcshareddata/xcschemes/lib-gpu-generate.xcscheme @@ -53,7 +53,7 @@ </BuildableProductRunnable> <CommandLineArguments> <CommandLineArgument - argument = "16" + argument = "32" isEnabled = "YES"> </CommandLineArgument> </CommandLineArguments> diff --git a/xcode/lib-gpu-verify.xcodeproj/xcuserdata/cedriczwahlen.xcuserdatad/xcdebugger/Breakpoints_v2.xcbkptlist b/xcode/lib-gpu-verify.xcodeproj/xcuserdata/cedriczwahlen.xcuserdatad/xcdebugger/Breakpoints_v2.xcbkptlist @@ -1199,10 +1199,87 @@ filePath = "../source/montgomery.c" startingColumnNumber = "9223372036854775807" endingColumnNumber = "9223372036854775807" - startingLineNumber = "58" - endingLineNumber = "58" + startingLineNumber = "96" + endingLineNumber = "96" landmarkName = "mont_prepare(base, exponent, modulus)" landmarkType = "9"> + <Locations> + <Location + uuid = "985780EE-603E-4B6C-BF80-1BB11F65F6BA - 554a9b6fb8ae965a" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_modexp" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "80" + endingLineNumber = "80" + offsetFromSymbolStart = "4"> + </Location> + <Location + uuid = "985780EE-603E-4B6C-BF80-1BB11F65F6BA - 4382d64135f5dab5" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_prepare" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "78" + endingLineNumber = "78" + offsetFromSymbolStart = "613"> + </Location> + <Location + uuid = "985780EE-603E-4B6C-BF80-1BB11F65F6BA - 4382d64135f5db8c" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_prepare" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "85" + endingLineNumber = "85" + offsetFromSymbolStart = "671"> + </Location> + <Location + uuid = "985780EE-603E-4B6C-BF80-1BB11F65F6BA - 4382d64135f5d884" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_prepare" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "93" + endingLineNumber = "93" + offsetFromSymbolStart = "665"> + </Location> + <Location + uuid = "985780EE-603E-4B6C-BF80-1BB11F65F6BA - 4382d64135f5d8a5" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_prepare" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "94" + endingLineNumber = "94" + offsetFromSymbolStart = "677"> + </Location> + </Locations> </BreakpointContent> </BreakpointProxy> <BreakpointProxy @@ -1215,8 +1292,8 @@ filePath = "../source/montgomery.c" startingColumnNumber = "9223372036854775807" endingColumnNumber = "9223372036854775807" - startingLineNumber = "38" - endingLineNumber = "38" + startingLineNumber = "49" + endingLineNumber = "49" landmarkName = "mont_prepare(base, exponent, modulus)" landmarkType = "9"> <Locations> @@ -2094,8 +2171,246 @@ endingLineNumber = "615" offsetFromSymbolStart = "25"> </Location> + <Location + uuid = "9303A078-7EDA-414A-8825-97250FC649BC - b0b9078e770cf6a3" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "rsa_tests" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/rsa-test.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "611" + endingLineNumber = "611" + offsetFromSymbolStart = "52"> + </Location> </Locations> </BreakpointContent> </BreakpointProxy> + <BreakpointProxy + BreakpointExtensionID = "Xcode.Breakpoint.FileBreakpoint"> + <BreakpointContent + uuid = "97D703B2-D63F-4EC4-A612-0E2231CC441C" + shouldBeEnabled = "No" + ignoreCount = "0" + continueAfterRunningActions = "No" + filePath = "../source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "74" + endingLineNumber = "74" + landmarkName = "mont_prepare(base, exponent, modulus)" + landmarkType = "9"> + <Locations> + <Location + uuid = "97D703B2-D63F-4EC4-A612-0E2231CC441C - 4382d64135f42485" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_prepare" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "62" + endingLineNumber = "62" + offsetFromSymbolStart = "448"> + </Location> + <Location + uuid = "97D703B2-D63F-4EC4-A612-0E2231CC441C - 4382d64135f424a6" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_prepare" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "63" + endingLineNumber = "63" + offsetFromSymbolStart = "460"> + </Location> + </Locations> + </BreakpointContent> + </BreakpointProxy> + <BreakpointProxy + BreakpointExtensionID = "Xcode.Breakpoint.FileBreakpoint"> + <BreakpointContent + uuid = "3FBF3C5D-8FAE-42D2-AF5F-6D59CDF9A7D5" + shouldBeEnabled = "No" + ignoreCount = "0" + continueAfterRunningActions = "No" + filePath = "../source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "68" + endingLineNumber = "68" + landmarkName = "mont_prepare(base, exponent, modulus)" + landmarkType = "9"> + <Locations> + <Location + uuid = "3FBF3C5D-8FAE-42D2-AF5F-6D59CDF9A7D5 - 4382d64135f427e0" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_prepare" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "57" + endingLineNumber = "57" + offsetFromSymbolStart = "382"> + </Location> + <Location + uuid = "3FBF3C5D-8FAE-42D2-AF5F-6D59CDF9A7D5 - 4382d64135f427e0" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_prepare" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "57" + endingLineNumber = "57" + offsetFromSymbolStart = "387"> + </Location> + </Locations> + </BreakpointContent> + </BreakpointProxy> + <BreakpointProxy + BreakpointExtensionID = "Xcode.Breakpoint.FileBreakpoint"> + <BreakpointContent + uuid = "DB8454C8-42E4-4F53-A753-783E113BA61B" + shouldBeEnabled = "No" + ignoreCount = "0" + continueAfterRunningActions = "No" + filePath = "../source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "183" + endingLineNumber = "183" + landmarkName = "mont_product(ret, a, b, r, r_1, n, ni)" + landmarkType = "9"> + <Locations> + <Location + uuid = "DB8454C8-42E4-4F53-A753-783E113BA61B - 4382d80c0bf89e2c" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_product" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "135" + endingLineNumber = "135" + offsetFromSymbolStart = "288"> + </Location> + <Location + uuid = "DB8454C8-42E4-4F53-A753-783E113BA61B - 4382d80c0bf897b3" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_product" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "186" + endingLineNumber = "186" + offsetFromSymbolStart = "288"> + </Location> + </Locations> + </BreakpointContent> + </BreakpointProxy> + <BreakpointProxy + BreakpointExtensionID = "Xcode.Breakpoint.FileBreakpoint"> + <BreakpointContent + uuid = "5237754D-0F3B-47A1-B768-D4F7FD47830D" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + filePath = "../source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "123" + endingLineNumber = "123" + landmarkName = "mont_modexp(ret, a, e, M, n, ni, r, r_1)" + landmarkType = "9"> + <Locations> + <Location + uuid = "5237754D-0F3B-47A1-B768-D4F7FD47830D - 554a9b6fb8ae9009" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_modexp" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "125" + endingLineNumber = "125" + offsetFromSymbolStart = "231"> + </Location> + <Location + uuid = "5237754D-0F3B-47A1-B768-D4F7FD47830D - 554a9b6fb8ae9009" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_modexp" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "125" + endingLineNumber = "125" + offsetFromSymbolStart = "250"> + </Location> + <Location + uuid = "5237754D-0F3B-47A1-B768-D4F7FD47830D - 554a9b6fb8ae938d" + shouldBeEnabled = "Yes" + ignoreCount = "0" + continueAfterRunningActions = "No" + symbolName = "mont_modexp" + moduleName = "lib-gpu-verify" + usesParentBreakpointCondition = "Yes" + urlString = "file:///Users/cedriczwahlen/libgpuverify/source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "121" + endingLineNumber = "121" + offsetFromSymbolStart = "250"> + </Location> + </Locations> + </BreakpointContent> + </BreakpointProxy> + <BreakpointProxy + BreakpointExtensionID = "Xcode.Breakpoint.FileBreakpoint"> + <BreakpointContent + uuid = "611D5F4B-F0CE-4BB5-B380-7D1B42A74A10" + shouldBeEnabled = "No" + ignoreCount = "0" + continueAfterRunningActions = "No" + filePath = "../source/montgomery.c" + startingColumnNumber = "9223372036854775807" + endingColumnNumber = "9223372036854775807" + startingLineNumber = "92" + endingLineNumber = "92" + landmarkName = "mont_prepare(base, exponent, modulus)" + landmarkType = "9"> + </BreakpointContent> + </BreakpointProxy> </Breakpoints> </Bucket>