source/libs/gmp/gmp-src/mpn/generic/gcdext_1.c

   1 /* mpn_gcdext -- Extended Greatest Common Divisor.
   2
   3 Copyright 1996, 1998, 2000-2005, 2008, 2009 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of either:
   9
  10   * the GNU Lesser General Public License as published by the Free
  11     Software Foundation; either version 3 of the License, or (at your
  12     option) any later version.
  13
  14 or
  15
  16   * the GNU General Public License as published by the Free Software
  17     Foundation; either version 2 of the License, or (at your option) any
  18     later version.
  19
  20 or both in parallel, as here.
  21
  22 The GNU MP Library is distributed in the hope that it will be useful, but
  23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  24 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  25 for more details.
  26
  27 You should have received copies of the GNU General Public License and the
  28 GNU Lesser General Public License along with the GNU MP Library.  If not,
  29 see https://www.gnu.org/licenses/.  */
  30
  31 #include "gmp.h"
  32 #include "gmp-impl.h"
  33 #include "longlong.h"
  34
  35 #ifndef GCDEXT_1_USE_BINARY
  36 #define GCDEXT_1_USE_BINARY 0
  37 #endif
  38
  39 #ifndef GCDEXT_1_BINARY_METHOD
  40 #define GCDEXT_1_BINARY_METHOD 2
  41 #endif
  42
  43 #ifndef USE_ZEROTAB
  44 #define USE_ZEROTAB 1
  45 #endif
  46
  47 #if GCDEXT_1_USE_BINARY
  48
  49 #if USE_ZEROTAB
  50 static unsigned char zerotab[0x40] = {
  51   6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
  52   4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
  53   5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
  54   4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
  55 };
  56 #endif
  57
  58 mp_limb_t
  59 mpn_gcdext_1 (mp_limb_signed_t *sp, mp_limb_signed_t *tp,
  60               mp_limb_t u, mp_limb_t v)
  61 {
  62   /* Maintain
  63
  64      U = t1 u + t0 v
  65      V = s1 u + s0 v
  66
  67      where U, V are the inputs (without any shared power of two),
  68      and the matrix has determinant ± 2^{shift}.
  69   */
  70   mp_limb_t s0 = 1;
  71   mp_limb_t t0 = 0;
  72   mp_limb_t s1 = 0;
  73   mp_limb_t t1 = 1;
  74   mp_limb_t ug;
  75   mp_limb_t vg;
  76   mp_limb_t ugh;
  77   mp_limb_t vgh;
  78   unsigned zero_bits;
  79   unsigned shift;
  80   unsigned i;
  81 #if GCDEXT_1_BINARY_METHOD == 2
  82   mp_limb_t det_sign;
  83 #endif
  84
  85   ASSERT (u > 0);
  86   ASSERT (v > 0);
  87
  88   count_trailing_zeros (zero_bits, u | v);
  89   u >>= zero_bits;
  90   v >>= zero_bits;
  91
  92   if ((u & 1) == 0)
  93     {
  94       count_trailing_zeros (shift, u);
  95       u >>= shift;
  96       t1 <<= shift;
  97     }
  98   else if ((v & 1) == 0)
  99     {
 100       count_trailing_zeros (shift, v);
 101       v >>= shift;
 102       s0 <<= shift;
 103     }
 104   else
 105     shift = 0;
 106
 107 #if GCDEXT_1_BINARY_METHOD == 1
 108   while (u != v)
 109     {
 110       unsigned count;
 111       if (u > v)
 112         {
 113           u -= v;
 114 #if USE_ZEROTAB
 115           count = zerotab [u & 0x3f];
 116           u >>= count;
 117           if (UNLIKELY (count == 6))
 118             {
 119               unsigned c;
 120               do
 121                 {
 122                   c = zerotab[u & 0x3f];
 123                   u >>= c;
 124                   count += c;
 125                 }
 126               while (c == 6);
 127             }
 128 #else
 129           count_trailing_zeros (count, u);
 130           u >>= count;
 131 #endif
 132           t0 += t1; t1 <<= count;
 133           s0 += s1; s1 <<= count;
 134         }
 135       else
 136         {
 137           v -= u;
 138 #if USE_ZEROTAB
 139           count = zerotab [v & 0x3f];
 140           v >>= count;
 141           if (UNLIKELY (count == 6))
 142             {
 143               unsigned c;
 144               do
 145                 {
 146                   c = zerotab[v & 0x3f];
 147                   v >>= c;
 148                   count += c;
 149                 }
 150               while (c == 6);
 151             }
 152 #else
 153           count_trailing_zeros (count, v);
 154           v >>= count;
 155 #endif
 156           t1 += t0; t0 <<= count;
 157           s1 += s0; s0 <<= count;
 158         }
 159       shift += count;
 160     }
 161 #else
 162 # if GCDEXT_1_BINARY_METHOD == 2
 163   u >>= 1;
 164   v >>= 1;
 165
 166   det_sign = 0;
 167
 168   while (u != v)
 169     {
 170       unsigned count;
 171       mp_limb_t d =  u - v;
 172       mp_limb_t vgtu = LIMB_HIGHBIT_TO_MASK (d);
 173       mp_limb_t sx;
 174       mp_limb_t tx;
 175
 176       /* When v <= u (vgtu == 0), the updates are:
 177
 178            (u; v)   <-- ( (u - v) >> count; v)    (det = +(1<<count) for corr. M factor)
 179            (t1, t0) <-- (t1 << count, t0 + t1)
 180
 181          and when v > 0, the updates are
 182
 183            (u; v)   <-- ( (v - u) >> count; u)    (det = -(1<<count))
 184            (t1, t0) <-- (t0 << count, t0 + t1)
 185
 186          and similarly for s1, s0
 187       */
 188
 189       /* v <-- min (u, v) */
 190       v += (vgtu & d);
 191
 192       /* u <-- |u - v| */
 193       u = (d ^ vgtu) - vgtu;
 194
 195       /* Number of trailing zeros is the same no matter if we look at
 196        * d or u, but using d gives more parallelism. */
 197 #if USE_ZEROTAB
 198       count = zerotab[d & 0x3f];
 199       if (UNLIKELY (count == 6))
 200         {
 201           unsigned c = 6;
 202           do
 203             {
 204               d >>= c;
 205               c = zerotab[d & 0x3f];
 206               count += c;
 207             }
 208           while (c == 6);
 209         }
 210 #else
 211       count_trailing_zeros (count, d);
 212 #endif
 213       det_sign ^= vgtu;
 214
 215       tx = vgtu & (t0 - t1);
 216       sx = vgtu & (s0 - s1);
 217       t0 += t1;
 218       s0 += s1;
 219       t1 += tx;
 220       s1 += sx;
 221
 222       count++;
 223       u >>= count;
 224       t1 <<= count;
 225       s1 <<= count;
 226       shift += count;
 227     }
 228   u = (u << 1) + 1;
 229 # else /* GCDEXT_1_BINARY_METHOD == 2 */
 230 #  error Unknown GCDEXT_1_BINARY_METHOD
 231 # endif
 232 #endif
 233
 234   /* Now u = v = g = gcd (u,v). Compute U/g and V/g */
 235   ug = t0 + t1;
 236   vg = s0 + s1;
 237
 238   ugh = ug/2 + (ug & 1);
 239   vgh = vg/2 + (vg & 1);
 240
 241   /* Now ±2^{shift} g = s0 U - t0 V. Get rid of the power of two, using
 242      s0 U - t0 V = (s0 + V/g) U - (t0 + U/g) V. */
 243   for (i = 0; i < shift; i++)
 244     {
 245       mp_limb_t mask = - ( (s0 | t0) & 1);
 246
 247       s0 /= 2;
 248       t0 /= 2;
 249       s0 += mask & vgh;
 250       t0 += mask & ugh;
 251     }
 252   /* FIXME: Try simplifying this condition. */
 253   if ( (s0 > 1 && 2*s0 >= vg) || (t0 > 1 && 2*t0 >= ug) )
 254     {
 255       s0 -= vg;
 256       t0 -= ug;
 257     }
 258 #if GCDEXT_1_BINARY_METHOD == 2
 259   /* Conditional negation. */
 260   s0 = (s0 ^ det_sign) - det_sign;
 261   t0 = (t0 ^ det_sign) - det_sign;
 262 #endif
 263   *sp = s0;
 264   *tp = -t0;
 265
 266   return u << zero_bits;
 267 }
 268
 269 #else /* !GCDEXT_1_USE_BINARY */
 270
 271
 272 /* FIXME: Takes two single-word limbs. It could be extended to a
 273  * function that accepts a bignum for the first input, and only
 274  * returns the first co-factor. */
 275
 276 mp_limb_t
 277 mpn_gcdext_1 (mp_limb_signed_t *up, mp_limb_signed_t *vp,
 278               mp_limb_t a, mp_limb_t b)
 279 {
 280   /* Maintain
 281
 282      a =  u0 A + v0 B
 283      b =  u1 A + v1 B
 284
 285      where A, B are the original inputs.
 286   */
 287   mp_limb_signed_t u0 = 1;
 288   mp_limb_signed_t v0 = 0;
 289   mp_limb_signed_t u1 = 0;
 290   mp_limb_signed_t v1 = 1;
 291
 292   ASSERT (a > 0);
 293   ASSERT (b > 0);
 294
 295   if (a < b)
 296     goto divide_by_b;
 297
 298   for (;;)
 299     {
 300       mp_limb_t q;
 301
 302       q = a / b;
 303       a -= q * b;
 304
 305       if (a == 0)
 306         {
 307           *up = u1;
 308           *vp = v1;
 309           return b;
 310         }
 311       u0 -= q * u1;
 312       v0 -= q * v1;
 313
 314     divide_by_b:
 315       q = b / a;
 316       b -= q * a;
 317
 318       if (b == 0)
 319         {
 320           *up = u0;
 321           *vp = v0;
 322           return a;
 323         }
 324       u1 -= q * u0;
 325       v1 -= q * v0;
 326     }
 327 }
 328 #endif /* !GCDEXT_1_USE_BINARY */