github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpz/gcd.c (about)

     1  /* mpz/gcd.c:   Calculate the greatest common divisor of two integers.
     2  
     3  Copyright 1991, 1993, 1994, 1996, 2000-2002, 2005, 2010 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 either:
    10  
    11    * the GNU Lesser General Public License as published by the Free
    12      Software Foundation; either version 3 of the License, or (at your
    13      option) any later version.
    14  
    15  or
    16  
    17    * the GNU General Public License as published by the Free Software
    18      Foundation; either version 2 of the License, or (at your option) any
    19      later version.
    20  
    21  or both in parallel, as here.
    22  
    23  The GNU MP Library is distributed in the hope that it will be useful, but
    24  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    25  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    26  for more details.
    27  
    28  You should have received copies of the GNU General Public License and the
    29  GNU Lesser General Public License along with the GNU MP Library.  If not,
    30  see https://www.gnu.org/licenses/.  */
    31  
    32  #include "gmp.h"
    33  #include "gmp-impl.h"
    34  #include "longlong.h"
    35  
    36  
    37  void
    38  mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v)
    39  {
    40    unsigned long int g_zero_bits, u_zero_bits, v_zero_bits;
    41    mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs;
    42    mp_ptr tp;
    43    mp_ptr up;
    44    mp_size_t usize;
    45    mp_ptr vp;
    46    mp_size_t vsize;
    47    mp_size_t gsize;
    48    TMP_DECL;
    49  
    50    up = PTR(u);
    51    usize = ABSIZ (u);
    52    vp = PTR(v);
    53    vsize = ABSIZ (v);
    54    /* GCD(0, V) == V.  */
    55    if (usize == 0)
    56      {
    57        SIZ (g) = vsize;
    58        if (g == v)
    59  	return;
    60        MPZ_REALLOC (g, vsize);
    61        MPN_COPY (PTR (g), vp, vsize);
    62        return;
    63      }
    64  
    65    /* GCD(U, 0) == U.  */
    66    if (vsize == 0)
    67      {
    68        SIZ (g) = usize;
    69        if (g == u)
    70  	return;
    71        MPZ_REALLOC (g, usize);
    72        MPN_COPY (PTR (g), up, usize);
    73        return;
    74      }
    75  
    76    if (usize == 1)
    77      {
    78        SIZ (g) = 1;
    79        PTR (g)[0] = mpn_gcd_1 (vp, vsize, up[0]);
    80        return;
    81      }
    82  
    83    if (vsize == 1)
    84      {
    85        SIZ(g) = 1;
    86        PTR (g)[0] = mpn_gcd_1 (up, usize, vp[0]);
    87        return;
    88      }
    89  
    90    TMP_MARK;
    91  
    92    /*  Eliminate low zero bits from U and V and move to temporary storage.  */
    93    while (*up == 0)
    94      up++;
    95    u_zero_limbs = up - PTR(u);
    96    usize -= u_zero_limbs;
    97    count_trailing_zeros (u_zero_bits, *up);
    98    tp = up;
    99    up = TMP_ALLOC_LIMBS (usize);
   100    if (u_zero_bits != 0)
   101      {
   102        mpn_rshift (up, tp, usize, u_zero_bits);
   103        usize -= up[usize - 1] == 0;
   104      }
   105    else
   106      MPN_COPY (up, tp, usize);
   107  
   108    while (*vp == 0)
   109      vp++;
   110    v_zero_limbs = vp - PTR (v);
   111    vsize -= v_zero_limbs;
   112    count_trailing_zeros (v_zero_bits, *vp);
   113    tp = vp;
   114    vp = TMP_ALLOC_LIMBS (vsize);
   115    if (v_zero_bits != 0)
   116      {
   117        mpn_rshift (vp, tp, vsize, v_zero_bits);
   118        vsize -= vp[vsize - 1] == 0;
   119      }
   120    else
   121      MPN_COPY (vp, tp, vsize);
   122  
   123    if (u_zero_limbs > v_zero_limbs)
   124      {
   125        g_zero_limbs = v_zero_limbs;
   126        g_zero_bits = v_zero_bits;
   127      }
   128    else if (u_zero_limbs < v_zero_limbs)
   129      {
   130        g_zero_limbs = u_zero_limbs;
   131        g_zero_bits = u_zero_bits;
   132      }
   133    else  /*  Equal.  */
   134      {
   135        g_zero_limbs = u_zero_limbs;
   136        g_zero_bits = MIN (u_zero_bits, v_zero_bits);
   137      }
   138  
   139    /*  Call mpn_gcd.  The 2nd argument must not have more bits than the 1st.  */
   140    vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
   141      ? mpn_gcd (vp, vp, vsize, up, usize)
   142      : mpn_gcd (vp, up, usize, vp, vsize);
   143  
   144    /*  Here G <-- V << (g_zero_limbs*GMP_LIMB_BITS + g_zero_bits).  */
   145    gsize = vsize + g_zero_limbs;
   146    if (g_zero_bits != 0)
   147      {
   148        mp_limb_t cy_limb;
   149        gsize += (vp[vsize - 1] >> (GMP_NUMB_BITS - g_zero_bits)) != 0;
   150        MPZ_REALLOC (g, gsize);
   151        MPN_ZERO (PTR (g), g_zero_limbs);
   152  
   153        tp = PTR(g) + g_zero_limbs;
   154        cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
   155        if (cy_limb != 0)
   156  	tp[vsize] = cy_limb;
   157      }
   158    else
   159      {
   160        MPZ_REALLOC (g, gsize);
   161        MPN_ZERO (PTR (g), g_zero_limbs);
   162        MPN_COPY (PTR (g) + g_zero_limbs, vp, vsize);
   163      }
   164  
   165    SIZ (g) = gsize;
   166    TMP_FREE;
   167  }