commit c44d1c803e6264930b07fb2d3cc79e4b3ce63e20
parent 7971a7a87e61d64445b0bc2fefc19e4c058ccad3
Author: Christian Grothoff <grothoff@gnunet.org>
Date: Thu, 7 Dec 2023 00:05:16 +0900
-restore to state with mp_alloc
Diffstat:
| M | source/montgomery.cl | | | 151 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------- |
1 file changed, 116 insertions(+), 35 deletions(-)
diff --git a/source/montgomery.cl b/source/montgomery.cl
@@ -201,6 +201,8 @@
#define assert(x){if((x)==0){printf((char __constant *)"assert reached\n");}}
+#define NULL ((void*)0)
+
typedef unsigned MINI_GMP_LIMB_TYPE mp_limb_t;
typedef long mp_size_t;
typedef unsigned long mp_bitcnt_t;
@@ -210,10 +212,14 @@ 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[256];
+ //mp_limb_t *_mp_d; /* Pointer to the limbs. */
+
+ mp_limb_t _mp_d[256];
} __mpz_struct;
@@ -262,6 +268,9 @@ void
mpz_init2 (mpz_t r, mp_bitcnt_t bits);
void
mpz_init_set_ui (mpz_t r, unsigned long int x);
+void
+gmp_die (const char *msg);
+
mp_size_t mpn_normalized_size (mp_srcptr xp, mp_size_t n);
void
@@ -371,6 +380,8 @@ int
mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index);
mp_bitcnt_t
mpn_limb_size_in_base_2 (mp_limb_t u);
+size_t
+mpz_sizeinbase (const mpz_t u, int base);
int
mpz_sgn (const mpz_t u);
mp_bitcnt_t
@@ -401,7 +412,8 @@ mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v);
unsigned
mpn_base_power_of_two_p (unsigned b);
-
+void
+mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b);
int isspace_gpu(unsigned char c);
int strlen_c(__global char *c);
mp_size_t mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
@@ -428,9 +440,12 @@ mpz_set_lg (unsigned long *r, __global unsigned long *x);
#define mpn_invert_limb(x) mpn_invert_3by2 ((x), 0)
+#define MPZ_REALLOC(z,n) (z)->_mp_d
+
void
mpz_init (mpz_t r)
{
+ r->_mp_alloc = 0;
r->_mp_size = 0;
for (int i = 0; i < 256; i++) {
r->_mp_d[i] = 0;
@@ -455,7 +470,7 @@ mpz_set (mpz_t r, const mpz_t x)
mp_ptr rp;
n = GMP_ABS (x->_mp_size);
- rp = r->_mp_d;
+ rp = MPZ_REALLOC (r, n);
mpn_copyi (rp, x->_mp_d, n);
r->_mp_size = x->_mp_size;
@@ -465,8 +480,14 @@ mpz_set (mpz_t r, const mpz_t x)
void
mpz_set_lg (unsigned long *r, __global unsigned long *x)
{
- for (int i = 0; i < 256; i++)
+ r[0] = x[0];
+ r[1] = x[1];
+
+ for (int i = 2; i < 256; i++) {
+
r[i] = x[i];
+
+ }
}
@@ -476,14 +497,14 @@ mpz_set_ui (mpz_t r, unsigned long int x)
if (x > 0)
{
r->_mp_size = 1;
- r->_mp_d[0] = x;
+ 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;
- r->_mp_d[r->_mp_size - 1] = x;
+ MPZ_REALLOC (r, r->_mp_size)[r->_mp_size - 1] = x;
}
}
}
@@ -514,7 +535,7 @@ mpz_set_si (mpz_t r, signed long int x)
else
{
r->_mp_size = -1;
- r->_mp_d[0] = GMP_NEG_CAST (unsigned long int, x);
+ MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x);
}
}
@@ -541,6 +562,7 @@ mpz_init2 (mpz_t r, mp_bitcnt_t bits)
bits -= (bits != 0); /* Round down, except if 0 */
rn = 1 + bits / GMP_LIMB_BITS;
+ r->_mp_alloc = rn;
r->_mp_size = 0;
}
@@ -552,13 +574,13 @@ mpz_init_set_ui (mpz_t r, unsigned long int x)
}
-#define gmp_die(msg) \
-do \
-{ \
- printf (msg); \
- assert(0); \
-} while (0)
-
+void
+gmp_die (const char *msg)
+{
+ //fprintf (stderr, "%s\n", msg);
+ //abort();
+}
+
mp_size_t mpn_normalized_size (mp_srcptr xp, mp_size_t n)
{
while (n > 0 && xp[n-1] == 0)
@@ -803,8 +825,9 @@ mpz_div_qr (mpz_t q, mpz_t r,
ds = d->_mp_size;
if (ds == 0) {
- gmp_die("mpz_div_qr: Divide by zero.");
+
}
+ //gmp_die("mpz_div_qr: Divide by zero.");
if (ns == 0)
{
@@ -940,6 +963,9 @@ mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v)
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_t temp;
mpz_init(temp);
@@ -1372,9 +1398,9 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
en = GMP_ABS (e->_mp_size);
mn = GMP_ABS (m->_mp_size);
- if (mn == 0) {
- gmp_die ("mpz_powm: Zero modulo.");
- }
+ if (mn == 0) {}
+ //gmp_die ("mpz_powm: Zero modulo.");
+
if (en == 0)
{
mpz_set_ui (r, mpz_cmpabs_ui (m, 1));
@@ -1400,9 +1426,8 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
if (e->_mp_size < 0)
{
- if (!mpz_invert (base, b, m)) {
- gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
- }
+ if (!mpz_invert (base, b, m)) {}
+ //gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
}
else
{
@@ -1421,7 +1446,7 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
m. */
if (b->_mp_size < 0)
{
- mp_ptr bp = base->_mp_d;
+ mp_ptr bp = MPZ_REALLOC (base, mn);
gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
bn = mn;
}
@@ -1481,13 +1506,13 @@ mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
if (cmp > 0)
{
- rp = r->_mp_d;
+ 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 = r->_mp_d;
+ 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);
}
@@ -1542,7 +1567,7 @@ mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
MP_SIZE_T_SWAP (an, bn);
}
- rp = r->_mp_d;
+ rp = MPZ_REALLOC (r, an + 1);
cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
rp[an] = cy;
@@ -1710,7 +1735,7 @@ mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
shift = bits % GMP_LIMB_BITS;
rn = un + limbs + (shift > 0);
- rp = r->_mp_d;
+ rp = MPZ_REALLOC (r, rn);
if (shift > 0)
{
mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
@@ -1765,7 +1790,7 @@ mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
qn = 0;
else
{
- qp = q->_mp_d;
+ qp = MPZ_REALLOC (q, qn);
if (bit_index != 0)
{
@@ -1864,13 +1889,22 @@ mpn_limb_size_in_base_2 (mp_limb_t u)
}
size_t
-mpz_sizeinbase2 (const mpz_t u)
+mpz_sizeinbase (const mpz_t u, int base)
{
- mp_size_t un;
+ mp_size_t un, tn;
mp_srcptr up;
+ //mp_ptr tp;
+ mpz_t tp;
mp_bitcnt_t bits;
+ struct gmp_div_inverse bi;
+ size_t ndigits;
+ mpz_init(tp);
+
+ assert (base >= 2);
+assert (base <= 62);
+
un = GMP_ABS (u->_mp_size);
if (un == 0)
return 1;
@@ -1878,7 +1912,39 @@ mpz_sizeinbase2 (const mpz_t u)
up = u->_mp_d;
bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]);
- return bits;
+ 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->_mp_d, up, un);
+ mpn_div_qr_1_invert (&bi, base);
+
+ tn = un;
+ ndigits = 0;
+ do
+ {
+ ndigits++;
+ mpn_div_qr_1_preinv (tp->_mp_d, tp->_mp_d, tn, &bi);
+ tn -= (tp->_mp_d[tn-1] == 0);
+ }
+ while (tn > 0);
+
+ //gmp_free_limbs (tp, un);
+ return ndigits;
}
int
@@ -1998,7 +2064,7 @@ mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
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 = d->_mp_d;
+ dp = MPZ_REALLOC (d, limb_index + 1);
dp[limb_index] = bit;
for (i = dn; i < limb_index; i++)
@@ -2014,7 +2080,7 @@ mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit);
if (cy > 0)
{
- dp = d->_mp_d;
+ dp = MPZ_REALLOC (d, dn + 1);
dp[dn++] = cy;
}
}
@@ -2258,6 +2324,21 @@ mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v)
}
+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;
+}
+
// Montgomery multiplication
@@ -2318,7 +2399,7 @@ void mont_prepare(mpz_t b, mpz_t e, mpz_t m,
mpz_init_set_si(one,1);
mpz_init_set_si(oo,0);
- //unsigned long len = mpz_sizeinbase2(m);
+ //unsigned long len = mpz_sizeinbase(m,2);
unsigned long len = 2049;
@@ -2370,7 +2451,7 @@ void mont_modexp(mpz_t ret,
mpz_init_set(aa, M);
mpz_init_set(xx, a);
- int k = (int)mpz_sizeinbase2(e);
+ int k = (int)mpz_sizeinbase(e,2);
for (int i = k - 1; i >= 0; i--) {
@@ -2430,7 +2511,7 @@ void mont_product(mpz_t ret,
mpz_add(ab, ab, mn);
- unsigned long sz = mpz_sizeinbase2(r) - 1;
+ unsigned long sz = mpz_sizeinbase(r,2) - 1;
mpz_tdiv_q_2exp(u, ab, sz); // this is essentially a bit shift, instead of a division
if (mpz_cmp(u, n) >= 0)
