[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