contrib/gmp/mpn/generic/gcdext_lehmer.c

   1 /* mpn_gcdext -- Extended Greatest Common Divisor.
   2
   3 Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2005, 2008, 2009 Free Software
   4 Foundation, Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   8 The GNU MP Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MP Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  20
  21 #include "gmp.h"
  22 #include "gmp-impl.h"
  23 #include "longlong.h"
  24
  25 /* Temporary storage: 3*(n+1) for u. n+1 for the matrix-vector
  26    multiplications (if hgcd2 succeeds). If hgcd fails, n+1 limbs are
  27    needed for the division, with most n for the quotient, and n+1 for
  28    the product q u0. In all, 4n + 3. */
  29
  30 mp_size_t
  31 mpn_gcdext_lehmer_n (mp_ptr gp, mp_ptr up, mp_size_t *usize,
  32                      mp_ptr ap, mp_ptr bp, mp_size_t n,
  33                      mp_ptr tp)
  34 {
  35   mp_size_t ualloc = n + 1;
  36
  37   /* Keeps track of the second row of the reduction matrix
  38    *
  39    *   M = (v0, v1 ; u0, u1)
  40    *
  41    * which correspond to the first column of the inverse
  42    *
  43    *   M^{-1} = (u1, -v1; -u0, v0)
  44    */
  45
  46   mp_size_t un;
  47   mp_ptr u0;
  48   mp_ptr u1;
  49   mp_ptr u2;
  50
  51   MPN_ZERO (tp, 3*ualloc);
  52   u0 = tp; tp += ualloc;
  53   u1 = tp; tp += ualloc;
  54   u2 = tp; tp += ualloc;
  55
  56   u1[0] = 1; un = 1;
  57
  58   /* FIXME: Handle n == 2 differently, after the loop? */
  59   while (n >= 2)
  60     {
  61       struct hgcd_matrix1 M;
  62       mp_limb_t ah, al, bh, bl;
  63       mp_limb_t mask;
  64
  65       mask = ap[n-1] | bp[n-1];
  66       ASSERT (mask > 0);
  67
  68       if (mask & GMP_NUMB_HIGHBIT)
  69         {
  70           ah = ap[n-1]; al = ap[n-2];
  71           bh = bp[n-1]; bl = bp[n-2];
  72         }
  73       else if (n == 2)
  74         {
  75           /* We use the full inputs without truncation, so we can
  76              safely shift left. */
  77           int shift;
  78
  79           count_leading_zeros (shift, mask);
  80           ah = MPN_EXTRACT_NUMB (shift, ap[1], ap[0]);
  81           al = ap[0] << shift;
  82           bh = MPN_EXTRACT_NUMB (shift, bp[1], bp[0]);
  83           bl = bp[0] << shift;
  84         }
  85       else
  86         {
  87           int shift;
  88
  89           count_leading_zeros (shift, mask);
  90           ah = MPN_EXTRACT_NUMB (shift, ap[n-1], ap[n-2]);
  91           al = MPN_EXTRACT_NUMB (shift, ap[n-2], ap[n-3]);
  92           bh = MPN_EXTRACT_NUMB (shift, bp[n-1], bp[n-2]);
  93           bl = MPN_EXTRACT_NUMB (shift, bp[n-2], bp[n-3]);
  94         }
  95
  96       /* Try an mpn_nhgcd2 step */
  97       if (mpn_hgcd2 (ah, al, bh, bl, &M))
  98         {
  99           n = mpn_hgcd_mul_matrix1_inverse_vector (&M, tp, ap, bp, n);
 100           MP_PTR_SWAP (ap, tp);
 101           un = mpn_hgcd_mul_matrix1_vector(&M, u2, u0, u1, un);
 102           MP_PTR_SWAP (u0, u2);
 103         }
 104       else
 105         {
 106           /* mpn_hgcd2 has failed. Then either one of a or b is very
 107              small, or the difference is very small. Perform one
 108              subtraction followed by one division. */
 109           mp_size_t gn;
 110           mp_size_t updated_un = un;
 111
 112           /* Temporary storage n for the quotient and ualloc for the
 113              new cofactor. */
 114           n = mpn_gcdext_subdiv_step (gp, &gn, up, usize, ap, bp, n,
 115                                       u0, u1, &updated_un, tp, u2);
 116           if (n == 0)
 117             return gn;
 118
 119           un = updated_un;
 120         }
 121     }
 122   ASSERT_ALWAYS (ap[0] > 0);
 123   ASSERT_ALWAYS (bp[0] > 0);
 124
 125   if (ap[0] == bp[0])
 126     {
 127       int c;
 128
 129       /* Which cofactor to return now? Candidates are +u1 and -u0,
 130          depending on which of a and b was most recently reduced,
 131          which we don't keep track of. So compare and get the smallest
 132          one. */
 133
 134       gp[0] = ap[0];
 135
 136       MPN_CMP (c, u0, u1, un);
 137       ASSERT (c != 0 || (un == 1 && u0[0] == 1 && u1[0] == 1));
 138       if (c < 0)
 139         {
 140           MPN_NORMALIZE (u0, un);
 141           MPN_COPY (up, u0, un);
 142           *usize = -un;
 143         }
 144       else
 145         {
 146           MPN_NORMALIZE_NOT_ZERO (u1, un);
 147           MPN_COPY (up, u1, un);
 148           *usize = un;
 149         }
 150       return 1;
 151     }
 152   else
 153     {
 154       mp_limb_t uh, vh;
 155       mp_limb_signed_t u;
 156       mp_limb_signed_t v;
 157       int negate;
 158
 159       gp[0] = mpn_gcdext_1 (&u, &v, ap[0], bp[0]);
 160
 161       /* Set up = u u1 - v u0. Keep track of size, un grows by one or
 162          two limbs. */
 163
 164       if (u == 0)
 165         {
 166           ASSERT (v == 1);
 167           MPN_NORMALIZE (u0, un);
 168           MPN_COPY (up, u0, un);
 169           *usize = -un;
 170           return 1;
 171         }
 172       else if (v == 0)
 173         {
 174           ASSERT (u == 1);
 175           MPN_NORMALIZE (u1, un);
 176           MPN_COPY (up, u1, un);
 177           *usize = un;
 178           return 1;
 179         }
 180       else if (u > 0)
 181         {
 182           negate = 0;
 183           ASSERT (v < 0);
 184           v = -v;
 185         }
 186       else
 187         {
 188           negate = 1;
 189           ASSERT (v > 0);
 190           u = -u;
 191         }
 192
 193       uh = mpn_mul_1 (up, u1, un, u);
 194       vh = mpn_addmul_1 (up, u0, un, v);
 195
 196       if ( (uh | vh) > 0)
 197         {
 198           uh += vh;
 199           up[un++] = uh;
 200           if (uh < vh)
 201             up[un++] = 1;
 202         }
 203
 204       MPN_NORMALIZE_NOT_ZERO (up, un);
 205
 206       *usize = negate ? -un : un;
 207       return 1;
 208     }
 209 }