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

   1 /* mpn_gcd_1 -- mpn and limb greatest common divisor.
   2
   3 Copyright 1994, 1996, 2000, 2001, 2009, 2012 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 GCD_1_METHOD
  36 #define GCD_1_METHOD 2
  37 #endif
  38
  39 #define USE_ZEROTAB 0
  40
  41 #if USE_ZEROTAB
  42 #define MAXSHIFT 4
  43 #define MASK ((1 << MAXSHIFT) - 1)
  44 static const unsigned char zerotab[1 << MAXSHIFT] =
  45 {
  46 #if MAXSHIFT > 4
  47   5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
  48 #endif
  49   4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
  50 };
  51 #endif
  52
  53 /* Does not work for U == 0 or V == 0.  It would be tough to make it work for
  54    V == 0 since gcd(x,0) = x, and U does not generally fit in an mp_limb_t.
  55
  56    The threshold for doing u%v when size==1 will vary by CPU according to
  57    the speed of a division and the code generated for the main loop.  Any
  58    tuning for this is left to a CPU specific implementation.  */
  59
  60 mp_limb_t
  61 mpn_gcd_1 (mp_srcptr up, mp_size_t size, mp_limb_t vlimb)
  62 {
  63   mp_limb_t      ulimb;
  64   unsigned long  zero_bits, u_low_zero_bits;
  65
  66   ASSERT (size >= 1);
  67   ASSERT (vlimb != 0);
  68   ASSERT_MPN_NONZERO_P (up, size);
  69
  70   ulimb = up[0];
  71
  72   /* Need vlimb odd for modexact, want it odd to get common zeros. */
  73   count_trailing_zeros (zero_bits, vlimb);
  74   vlimb >>= zero_bits;
  75
  76   if (size > 1)
  77     {
  78       /* Must get common zeros before the mod reduction.  If ulimb==0 then
  79          vlimb already gives the common zeros.  */
  80       if (ulimb != 0)
  81         {
  82           count_trailing_zeros (u_low_zero_bits, ulimb);
  83           zero_bits = MIN (zero_bits, u_low_zero_bits);
  84         }
  85
  86       ulimb = MPN_MOD_OR_MODEXACT_1_ODD (up, size, vlimb);
  87       if (ulimb == 0)
  88         goto done;
  89
  90       goto strip_u_maybe;
  91     }
  92
  93   /* size==1, so up[0]!=0 */
  94   count_trailing_zeros (u_low_zero_bits, ulimb);
  95   ulimb >>= u_low_zero_bits;
  96   zero_bits = MIN (zero_bits, u_low_zero_bits);
  97
  98   /* make u bigger */
  99   if (vlimb > ulimb)
 100     MP_LIMB_T_SWAP (ulimb, vlimb);
 101
 102   /* if u is much bigger than v, reduce using a division rather than
 103      chipping away at it bit-by-bit */
 104   if ((ulimb >> 16) > vlimb)
 105     {
 106       ulimb %= vlimb;
 107       if (ulimb == 0)
 108         goto done;
 109       goto strip_u_maybe;
 110     }
 111
 112   ASSERT (ulimb & 1);
 113   ASSERT (vlimb & 1);
 114
 115 #if GCD_1_METHOD == 1
 116   while (ulimb != vlimb)
 117     {
 118       ASSERT (ulimb & 1);
 119       ASSERT (vlimb & 1);
 120
 121       if (ulimb > vlimb)
 122         {
 123           ulimb -= vlimb;
 124           do
 125             {
 126               ulimb >>= 1;
 127               ASSERT (ulimb != 0);
 128             strip_u_maybe:
 129               ;
 130             }
 131           while ((ulimb & 1) == 0);
 132         }
 133       else /*  vlimb > ulimb.  */
 134         {
 135           vlimb -= ulimb;
 136           do
 137             {
 138               vlimb >>= 1;
 139               ASSERT (vlimb != 0);
 140             }
 141           while ((vlimb & 1) == 0);
 142         }
 143     }
 144 #else
 145 # if GCD_1_METHOD  == 2
 146
 147   ulimb >>= 1;
 148   vlimb >>= 1;
 149
 150   while (ulimb != vlimb)
 151     {
 152       int c;
 153       mp_limb_t t;
 154       mp_limb_t vgtu;
 155
 156       t = ulimb - vlimb;
 157       vgtu = LIMB_HIGHBIT_TO_MASK (t);
 158
 159       /* v <-- min (u, v) */
 160       vlimb += (vgtu & t);
 161
 162       /* u <-- |u - v| */
 163       ulimb = (t ^ vgtu) - vgtu;
 164
 165 #if USE_ZEROTAB
 166       /* Number of trailing zeros is the same no matter if we look at
 167        * t or ulimb, but using t gives more parallelism. */
 168       c = zerotab[t & MASK];
 169
 170       while (UNLIKELY (c == MAXSHIFT))
 171         {
 172           ulimb >>= MAXSHIFT;
 173           if (0)
 174           strip_u_maybe:
 175             vlimb >>= 1;
 176
 177           c = zerotab[ulimb & MASK];
 178         }
 179 #else
 180       if (0)
 181         {
 182         strip_u_maybe:
 183           vlimb >>= 1;
 184           t = ulimb;
 185         }
 186       count_trailing_zeros (c, t);
 187 #endif
 188       ulimb >>= (c + 1);
 189     }
 190
 191   vlimb = (vlimb << 1) | 1;
 192 # else
 193 #  error Unknown GCD_1_METHOD
 194 # endif
 195 #endif
 196
 197  done:
 198   return vlimb << zero_bits;
 199 }