github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpn/generic/gcd_1.c (about)

     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  }