Torbjorn Granlund <[email protected]> writes: > * Finalise and commit mpn_sec_minvert.
Here's a new version, including tests. Seems to work. I'll try to get this committed fairly soon. It would be good to also add some support in speed, so performance can be compared easily to gcdext and powm. Regards, /Niels
diff -r 4561887f5963 ChangeLog --- a/ChangeLog Sat Dec 28 19:38:42 2013 +0100 +++ b/ChangeLog Sun Dec 29 16:17:28 2013 +0100 @@ -1,3 +1,14 @@ +2013-12-29 Niels Möller <[email protected]> + + * tests/mpn/Makefile.am (check_PROGRAMS): Added t-minvert. + * tests/mpn/t-minvert.c: New file. + + * configure.ac (gmp_mpn_functions): Added sec_minvert. + * mpn/generic/sec_minvert.c (mpn_sec_minvert) + (mpn_sec_minvert_itch): New functions. + (mpn_sec_add_1, mpn_cnd_neg, mpn_cnd_swap, mpn_sec_eq_ui): New + helper functions. + 2013-12-28 Torbjorn Granlund <[email protected]> * mpn/generic/sec_powm.c: Fix an ASSERT. diff -r 4561887f5963 configure.ac --- a/configure.ac Sat Dec 28 19:38:42 2013 +0100 +++ b/configure.ac Sun Dec 29 16:17:28 2013 +0100 @@ -2832,6 +2832,7 @@ bdiv_q bdiv_qr broot brootinv bsqrt bsqrtinv \ divexact bdiv_dbm1c redc_1 redc_2 redc_n powm powlo sec_powm \ sec_mul sec_sqr sec_div_qr sec_div_r sec_pi1_div_qr sec_pi1_div_r \ + sec_minvert \ trialdiv remove \ and_n andn_n nand_n ior_n iorn_n nior_n xor_n xnor_n \ copyi copyd zero sec_tabselect \ diff -r 4561887f5963 gmp-impl.h --- a/gmp-impl.h Sat Dec 28 19:38:42 2013 +0100 +++ b/gmp-impl.h Sun Dec 29 16:17:28 2013 +0100 @@ -1557,6 +1557,11 @@ #define mpn_sec_pi1_div_r __MPN(sec_pi1_div_r) __GMP_DECLSPEC void mpn_sec_pi1_div_r (mp_ptr, mp_size_t, mp_srcptr, mp_size_t, mp_limb_t, mp_ptr); +#define mpn_sec_minvert_itch __MPN(sec_minvert_itch) +__GMP_DECLSPEC mp_size_t mpn_sec_minvert_itch (mp_size_t); +#define mpn_sec_minvert __MPN(sec_minvert) +__GMP_DECLSPEC int mpn_sec_minvert (mp_ptr, mp_ptr ap, mp_srcptr mp, + mp_size_t, mp_bitcnt_t, mp_ptr); /* Override mpn_addlsh1_n, mpn_addlsh2_n, mpn_sublsh1_n, etc with mpn_addlsh_n, etc when !HAVE_NATIVE the former but HAVE_NATIVE_ the latter. We then lie diff -r 4561887f5963 mpn/generic/sec_minvert.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpn/generic/sec_minvert.c Sun Dec 29 16:17:28 2013 +0100 @@ -0,0 +1,195 @@ +/* mpn_sec_minvert + + Contributed to the GNU project by Niels Möller + +Copyright 2013 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 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. + +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 Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +static mp_size_t +mpn_sec_add_1_itch (mp_size_t n) +{ + return n; +} + +static mp_limb_t +mpn_sec_add_1 (mp_limb_t *rp, mp_limb_t *ap, mp_size_t n, mp_limb_t b, + mp_ptr scratch) +{ + scratch[0] = b; + MPN_ZERO (scratch + 1, n-1); + return mpn_add_n (rp, ap, scratch, n); +} + +mp_size_t +mpn_cnd_neg_itch (mp_size_t n) +{ + return n; +} + +/* FIXME: Ought to return carry */ +static void +mpn_cnd_neg (int cnd, mp_limb_t *rp, const mp_limb_t *ap, mp_size_t n, + mp_ptr scratch) +{ + mpn_lshift (scratch, ap, n, 1); + mpn_cnd_sub_n (cnd, rp, ap, scratch, n); +} + +static void +mpn_cnd_swap (int cnd, mp_limb_t *ap, mp_limb_t *bp, mp_size_t n) +{ + mp_limb_t mask = - (mp_limb_t) (cnd != 0); + mp_size_t i; + for (i = 0; i < n; i++) + { + mp_limb_t a, b, t; + a = ap[i]; + b = bp[i]; + t = (a ^ b) & mask; + ap[i] = a ^ t; + bp[i] = b ^ t; + } +} + +static int +mpn_sec_eq_ui (mp_srcptr ap, mp_size_t n, mp_limb_t b) +{ + mp_limb_t d; + ASSERT (n > 0); + + d = ap[0] ^ b; + + while (--n > 0) + d |= ap[n]; + + return d == 0; +} + +mp_size_t +mpn_sec_minvert_itch (mp_size_t n) +{ + return 4*n; +} + +/* Compute V <-- A^{-1} (mod M), in data-independent time. M must be + odd. Returns 1 on success, and 0 on failure (i.e., if gcd (A, m) != + 1). Inputs and outputs of size n, and no overlap allowed. The {ap, + n} area is destroyed. For arbitrary inputs, bit_size should be + 2*n*GMP_NUMB_BITS, but if A or M are known to be smaller, e.g., if + M = 2^521 - 1 and A < M, bit_size can be any bound on the sum of + the bit sizes of A and M. */ +int +mpn_sec_minvert (mp_ptr vp, mp_ptr ap, mp_srcptr mp, + mp_size_t n, mp_bitcnt_t bit_size, + mp_ptr scratch) +{ + ASSERT (n > 0); + ASSERT (bit_size > 0); + ASSERT (mp[0] & 1); + ASSERT (! MPN_OVERLAP_P (ap, n, vp, n)); +#define bp (scratch + n) +#define up (scratch + 2*n) +#define m1hp (scratch + 3*n) + + /* Maintain + + a = u * orig_a (mod m) + b = v * orig_a (mod m) + + and b odd at all times. Initially, + + a = a_orig, u = 1 + b = m, v = 0 + */ + + + up[0] = 1; + mpn_zero (up+1, n - 1); + mpn_copyi (bp, mp, n); + mpn_zero (vp, n); + + ASSERT_CARRY (mpn_rshift (m1hp, mp, n, 1)); + ASSERT_NOCARRY (mpn_sec_add_1 (m1hp, m1hp, n, 1, scratch)); + + while (bit_size-- > 0) + { + mp_limb_t odd, swap, cy; + + /* Always maintain b odd. The logic of the iteration is as + follows. For a, b: + + odd = a & 1 + a -= odd * b + if (underflow from a-b) + { + b += a, assigns old a + a = B^n-a + } + + a /= 2 + + For u, v: + + if (underflow from a - b) + swap u, v + u -= odd * v + if (underflow from u - v) + u += m + + u /= 2 + if (a one bit was shifted out) + u += (m+1)/2 + + As long as a > 0, the quantity + + (bitsize of a) + (bitsize of b) + + is reduced by at least one bit per iteration, hence after + (bit_size of orig_a) + (bit_size of m) - 1 iterations we + surely have a = 0. Then b = gcd(orig_a, m) and if b = 1 then + also v = orig_a^{-1} (mod m) + */ + + ASSERT (bp[0] & 1); + odd = ap[0] & 1; + + swap = mpn_cnd_sub_n (odd, ap, ap, bp, n); + mpn_cnd_add_n (swap, bp, bp, ap, n); + mpn_cnd_neg (swap, ap, ap, n, scratch); + + mpn_cnd_swap (swap, up, vp, n); + cy = mpn_cnd_sub_n (odd, up, up, vp, n); + cy -= mpn_cnd_add_n (cy, up, up, mp, n); + ASSERT (cy == 0); + + cy = mpn_rshift (ap, ap, n, 1); + ASSERT (cy == 0); + cy = mpn_rshift (up, up, n, 1); + cy = mpn_cnd_add_n (cy, up, up, m1hp, n); + ASSERT (cy == 0); + } + /* Should be all zeros, but check only extreme limbs */ + ASSERT ( (ap[0] | ap[n-1]) == 0); + /* Check if indeed gcd == 1. */ + return mpn_sec_eq_ui (bp, n, 1); +#undef bp +#undef up +#undef m1hp +} diff -r 4561887f5963 tests/mpn/Makefile.am --- a/tests/mpn/Makefile.am Sat Dec 28 19:38:42 2013 +0100 +++ b/tests/mpn/Makefile.am Sun Dec 29 16:17:28 2013 +0100 @@ -28,7 +28,7 @@ t-toom2-sqr t-toom3-sqr t-toom4-sqr t-toom6-sqr t-toom8-sqr \ t-mul t-mullo t-mulmod_bnm1 t-sqrmod_bnm1 t-mulmid \ t-hgcd t-hgcd_appr t-matrix22 t-invert t-div t-bdiv \ - t-broot t-brootinv + t-broot t-brootinv t-minvert EXTRA_DIST = toom-shared.h toom-sqr-shared.h diff -r 4561887f5963 tests/mpn/t-minvert.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/mpn/t-minvert.c Sun Dec 29 16:17:28 2013 +0100 @@ -0,0 +1,164 @@ +/* Copyright 2013 Free Software Foundation, Inc. + +This file is part of the GNU MP Library test suite. + +The GNU MP Library test suite is free software; you can redistribute it +and/or modify it under the terms of the GNU General Public License as +published by the Free Software Foundation; either version 3 of the License, +or (at your option) any later version. + +The GNU MP Library test suite 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 a copy of the GNU General Public License along with +the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ + +#include <stdio.h> +#include <stdlib.h> /* for strtol */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" +#include "tests/tests.h" + +#define MAX_SIZE 50 + +#define COUNT 200 + +static void +mpz_to_mpn (mp_ptr ap, mp_size_t an, const mpz_t b) +{ + mp_size_t bn = mpz_size (b); + ASSERT_ALWAYS (bn <= an); + MPN_COPY_INCR (ap, mpz_limbs_read (b), bn); + MPN_ZERO (ap + bn, an - bn); +} + +int +mpz_eq_mpn (mp_ptr ap, mp_size_t an, const mpz_t b) +{ + mp_size_t bn = mpz_size (b); + + return (bn >= 0 && bn <= an + && mpn_cmp (ap, mpz_limbs_read (b), bn) == 0 + && mpn_zero_p (ap + bn, an - bn)); +} + +int +main (int argc, char **argv) +{ + gmp_randstate_ptr rands; + long count = COUNT; + mp_ptr mp; + mp_ptr ap; + mp_ptr vp; + mp_ptr tp; + mp_ptr scratch; + mpz_t m, a, r, g; + int test; + tests_start (); + rands = RANDS; + + TMP_DECL; + + TMP_MARK; + mpz_init (m); + mpz_init (a); + mpz_init (r); + mpz_init (g); + + if (argc > 1) + { + char *end; + count = strtol (argv[1], &end, 0); + if (*end || count <= 0) + { + fprintf (stderr, "Invalid test count: %s.\n", argv[1]); + return 1; + } + } + + mp = TMP_ALLOC_LIMBS (MAX_SIZE); + ap = TMP_ALLOC_LIMBS (MAX_SIZE); + vp = TMP_ALLOC_LIMBS (MAX_SIZE); + tp = TMP_ALLOC_LIMBS (MAX_SIZE); + scratch = TMP_ALLOC_LIMBS (mpn_sec_minvert_itch (MAX_SIZE)); + + for (test = 0; test < count; test++) + { + mp_bitcnt_t bits; + int rres, tres; + mp_size_t n; + + bits = urandom () % (GMP_NUMB_BITS * MAX_SIZE) + 1; + + if (test & 1) + mpz_rrandomb (m, rands, bits); + else + mpz_urandomb (m, rands, bits); + if (test & 2) + mpz_rrandomb (a, rands, bits); + else + mpz_urandomb (a, rands, bits); + + mpz_setbit (m, 0); + if (test & 4) + { + /* Ensure it really is invertible */ + if (mpz_sgn (a) == 0) + mpz_set_ui (a, 1); + else + for (;;) + { + mpz_gcd (g, a, m); + if (mpz_cmp_ui (g, 1) == 0) + break; + mpz_remove (a, a, g); + } + } + + rres = mpz_invert (r, a, m); + if ( (test & 4) && !rres) + { + gmp_fprintf (stderr, "test %d: Not invertible!\n" + "m = %Zd\n" + "a = %Zd\n", test, m, a); + abort (); + } + ASSERT_ALWAYS (! (test & 4) || rres); + + n = (bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; + ASSERT_ALWAYS (n <= MAX_SIZE); + + mpz_to_mpn (ap, n, a); + mpz_to_mpn (mp, n, m); + tres = mpn_sec_minvert (tp, + ap, mp, n, 2*bits, scratch); + + if (rres != tres || (rres == 1 && !mpz_eq_mpn (tp, n, r))) + { + gmp_fprintf (stderr, "Test %d failed.\n" + "m = %Zd\n" + "a = %Zd\n", test, m, a); + fprintf (stderr, "ref ret: %d\n" + "got ret: %d\n", rres, tres); + if (rres) + gmp_fprintf (stderr, "ref: %Zd\n", r); + if (tres) + gmp_fprintf (stderr, "got: %Nd\n", tp, n); + abort (); + } + } + + TMP_FREE; + + mpz_clear (m); + mpz_clear (a); + mpz_clear (r); + mpz_clear (g); + + tests_end (); + return 0; +}
-- Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26. Internet email is subject to wholesale government surveillance.
_______________________________________________ gmp-devel mailing list [email protected] https://gmplib.org/mailman/listinfo/gmp-devel
