[Gmp-commit] /var/hg/gmp: New function mpn_sec_minvert.
mercurial at gmplib.org
mercurial at gmplib.org
Sun Dec 29 22:51:54 UTC 2013
details: /var/hg/gmp/rev/e0870b8bf13d
changeset: 16123:e0870b8bf13d
user: Niels M?ller <nisse at lysator.liu.se>
date: Sun Dec 29 23:49:58 2013 +0100
description:
New function mpn_sec_minvert.
diffstat:
ChangeLog | 13 +++
configure.ac | 1 +
gmp-h.in | 6 +
mpn/generic/sec_minvert.c | 195 ++++++++++++++++++++++++++++++++++++++++++++++
tests/mpn/Makefile.am | 2 +-
tests/mpn/t-minvert.c | 171 ++++++++++++++++++++++++++++++++++++++++
6 files changed, 387 insertions(+), 1 deletions(-)
diffs (truncated from 433 to 300 lines):
diff -r 4561887f5963 -r e0870b8bf13d ChangeLog
--- a/ChangeLog Sat Dec 28 19:38:42 2013 +0100
+++ b/ChangeLog Sun Dec 29 23:49:58 2013 +0100
@@ -1,3 +1,16 @@
+2013-12-29 Niels Möller <nisse at lysator.liu.se>
+
+ * tests/mpn/Makefile.am (check_PROGRAMS): Added t-minvert.
+ * tests/mpn/t-minvert.c: New file.
+
+ * configure.ac (gmp_mpn_functions): Added sec_minvert.
+ * gmp-h.in (mpn_sec_minvert, mpn_sec_minvert_itch): New
+ declarations.
+ * 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 <tege at gmplib.org>
* mpn/generic/sec_powm.c: Fix an ASSERT.
diff -r 4561887f5963 -r e0870b8bf13d configure.ac
--- a/configure.ac Sat Dec 28 19:38:42 2013 +0100
+++ b/configure.ac Sun Dec 29 23:49:58 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 -r e0870b8bf13d gmp-h.in
--- a/gmp-h.in Sat Dec 28 19:38:42 2013 +0100
+++ b/gmp-h.in Sun Dec 29 23:49:58 2013 +0100
@@ -1653,6 +1653,12 @@
#define mpn_sec_div_r_itch __MPN(sec_div_r_itch)
__GMP_DECLSPEC mp_size_t mpn_sec_div_r_itch (mp_size_t, 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);
+#define mpn_sec_minvert_itch __MPN(sec_minvert_itch)
+__GMP_DECLSPEC mp_size_t mpn_sec_minvert_itch (mp_size_t);
+
/**************** mpz inlines ****************/
diff -r 4561887f5963 -r e0870b8bf13d mpn/generic/sec_minvert.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpn/generic/sec_minvert.c Sun Dec 29 23:49:58 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 -r e0870b8bf13d 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 23:49:58 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 -r e0870b8bf13d tests/mpn/t-minvert.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/mpn/t-minvert.c Sun Dec 29 23:49:58 2013 +0100
@@ -0,0 +1,171 @@
+/* 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);
+}
+
More information about the gmp-commit
mailing list