X-Git-Url: https://oss.titaniummirror.com/gitweb?a=blobdiff_plain;f=gmp%2Fmpn%2Fgeneric%2Fgcdext_lehmer.c;fp=gmp%2Fmpn%2Fgeneric%2Fgcdext_lehmer.c;h=90c6412c002a3e87a19acfbf47376db1918ac716;hb=6fed43773c9b0ce596dca5686f37ac3fc0fa11c0;hp=0000000000000000000000000000000000000000;hpb=27b11d56b743098deb193d510b337ba22dc52e5c;p=msp430-gcc.git

diff --git a/gmp/mpn/generic/gcdext_lehmer.c b/gmp/mpn/generic/gcdext_lehmer.c
new file mode 100644
index 00000000..90c6412c
--- /dev/null
+++ b/gmp/mpn/generic/gcdext_lehmer.c
@@ -0,0 +1,168 @@
+/* mpn_gcdext -- Extended Greatest Common Divisor.
+
+Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2005, 2008 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 http://www.gnu.org/licenses/.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+/* Temporary storage: 3*(n+1) for u. n+1 for the matrix-vector
+   multiplications (if hgcd2 succeeds). If hgcd fails, n+1 limbs are
+   needed for the division, with most n for the quotient, and n+1 for
+   the product q u0. In all, 4n + 3. */
+
+mp_size_t
+mpn_gcdext_lehmer_n (mp_ptr gp, mp_ptr up, mp_size_t *usize,
+		     mp_ptr ap, mp_ptr bp, mp_size_t n,
+		     mp_ptr tp)
+{
+  mp_size_t ualloc = n + 1;
+
+  /* Keeps track of the second row of the reduction matrix
+   *
+   *   M = (v0, v1 ; u0, u1)
+   *
+   * which correspond to the first column of the inverse
+   *
+   *   M^{-1} = (u1, -v1; -u0, v0)
+   */
+
+  mp_size_t un;
+  mp_ptr u0;
+  mp_ptr u1;
+  mp_ptr u2;
+
+  MPN_ZERO (tp, 3*ualloc);
+  u0 = tp; tp += ualloc;
+  u1 = tp; tp += ualloc;
+  u2 = tp; tp += ualloc;
+
+  u1[0] = 1; un = 1;
+
+  /* FIXME: Handle n == 2 differently, after the loop? */
+  while (n >= 2)
+    {
+      struct hgcd_matrix1 M;
+      mp_limb_t ah, al, bh, bl;
+      mp_limb_t mask;
+
+      mask = ap[n-1] | bp[n-1];
+      ASSERT (mask > 0);
+
+      if (mask & GMP_NUMB_HIGHBIT)
+	{
+	  ah = ap[n-1]; al = ap[n-2];
+	  bh = bp[n-1]; bl = bp[n-2];
+	}
+      else if (n == 2)
+	{
+	  /* We use the full inputs without truncation, so we can
+	     safely shift left. */
+	  int shift;
+
+	  count_leading_zeros (shift, mask);
+	  ah = MPN_EXTRACT_NUMB (shift, ap[1], ap[0]);
+	  al = ap[0] << shift;
+	  bh = MPN_EXTRACT_NUMB (shift, bp[1], bp[0]);
+	  bl = bp[0] << shift;
+	}
+      else
+	{
+	  int shift;
+
+	  count_leading_zeros (shift, mask);
+	  ah = MPN_EXTRACT_NUMB (shift, ap[n-1], ap[n-2]);
+	  al = MPN_EXTRACT_NUMB (shift, ap[n-2], ap[n-3]);
+	  bh = MPN_EXTRACT_NUMB (shift, bp[n-1], bp[n-2]);
+	  bl = MPN_EXTRACT_NUMB (shift, bp[n-2], bp[n-3]);
+	}
+
+      /* Try an mpn_nhgcd2 step */
+      if (mpn_hgcd2 (ah, al, bh, bl, &M))
+	{
+	  n = mpn_hgcd_mul_matrix1_inverse_vector (&M, tp, ap, bp, n);
+	  MP_PTR_SWAP (ap, tp);
+	  un = mpn_hgcd_mul_matrix1_vector(&M, u2, u0, u1, un);
+	  MP_PTR_SWAP (u0, u2);
+	}
+      else
+	{
+	  /* mpn_hgcd2 has failed. Then either one of a or b is very
+	     small, or the difference is very small. Perform one
+	     subtraction followed by one division. */
+	  mp_size_t gn;
+	  mp_size_t updated_un = un;
+
+	  /* Temporary storage n for the quotient and ualloc for the
+	     new cofactor. */
+	  n = mpn_gcdext_subdiv_step (gp, &gn, up, usize, ap, bp, n,
+				      u0, u1, &updated_un, tp, u2);
+	  if (n == 0)
+	    return gn;
+
+	  un = updated_un;
+	}
+    }
+  if (ap[0] == 0)
+    {
+      gp[0] = bp[0];
+
+      MPN_NORMALIZE_NOT_ZERO (u0, un);
+      MPN_COPY (up, u0, un);
+
+      *usize = -un;
+      return 1;
+    }
+  else if (bp[0] == 0)
+    {
+      gp[0] = ap[0];
+
+      MPN_NORMALIZE_NOT_ZERO (u1, un);
+      MPN_COPY (up, u1, un);
+
+      *usize = un;
+      return 1;
+    }
+  else
+    {
+      mp_limb_t uh, vh;
+      mp_limb_t u;
+      mp_limb_t v;
+
+      gp[0] = mpn_gcdext_1 (&u, &v, ap[0], bp[0]);
+
+      /* Set up = u u1 + v u0. Keep track of size, un grows by one or
+	 two limbs. */
+      uh = mpn_mul_1 (up, u1, un, u);
+      vh = mpn_addmul_1 (up, u0, un, v);
+
+      if ( (uh | vh) > 0)
+	{
+	  uh += vh;
+	  up[un++] = uh;
+	  if (uh < vh)
+	    up[un++] = 1;
+	}
+
+      MPN_NORMALIZE_NOT_ZERO (up, un);
+
+      *usize = un;
+      return 1;
+    }
+}