Home | History | Annotate | Line # | Download | only in mpz 
      1 /* mpz/gcd.c:   Calculate the greatest common divisor of two integers.
      2 
      3 Copyright 1991, 1993, 1994, 1996, 2000-2002, 2005, 2010, 2015, 2016
      4 Free Software 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-impl.h"
     33 #include "longlong.h"
     34 
     35 
     36 void
     37 mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v)
     38 {
     39   unsigned long int g_zero_bits, u_zero_bits, v_zero_bits;
     40   mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs;
     41   mp_ptr tp;
     42   mp_ptr up;
     43   mp_size_t usize;
     44   mp_ptr vp;
     45   mp_size_t vsize;
     46   mp_size_t gsize;
     47   TMP_DECL;
     48 
     49   up = PTR(u);
     50   usize = ABSIZ (u);
     51   vp = PTR(v);
     52   vsize = ABSIZ (v);
     53   /* GCD(0, V) == V.  */
     54   if (usize == 0)
     55     {
     56       SIZ (g) = vsize;
     57       if (g == v)
     58 	return;
     59       tp = MPZ_NEWALLOC (g, vsize);
     60       MPN_COPY (tp, vp, vsize);
     61       return;
     62     }
     63 
     64   /* GCD(U, 0) == U.  */
     65   if (vsize == 0)
     66     {
     67       SIZ (g) = usize;
     68       if (g == u)
     69 	return;
     70       tp = MPZ_NEWALLOC (g, usize);
     71       MPN_COPY (tp, up, usize);
     72       return;
     73     }
     74 
     75   if (usize == 1)
     76     {
     77       SIZ (g) = 1;
     78       MPZ_NEWALLOC (g, 1)[0] = mpn_gcd_1 (vp, vsize, up[0]);
     79       return;
     80     }
     81 
     82   if (vsize == 1)
     83     {
     84       SIZ(g) = 1;
     85       MPZ_NEWALLOC (g, 1)[0] = mpn_gcd_1 (up, usize, vp[0]);
     86       return;
     87     }
     88 
     89   TMP_MARK;
     90 
     91   /*  Eliminate low zero bits from U and V and move to temporary storage.  */
     92   tp = up;
     93   while (*tp == 0)
     94     tp++;
     95   u_zero_limbs = tp - up;
     96   usize -= u_zero_limbs;
     97   count_trailing_zeros (u_zero_bits, *tp);
     98   up = TMP_ALLOC_LIMBS (usize);
     99   if (u_zero_bits != 0)
    100     {
    101       mpn_rshift (up, tp, usize, u_zero_bits);
    102       usize -= up[usize - 1] == 0;
    103     }
    104   else
    105     MPN_COPY (up, tp, usize);
    106 
    107   tp = vp;
    108   while (*tp == 0)
    109     tp++;
    110   v_zero_limbs = tp - vp;
    111   vsize -= v_zero_limbs;
    112   count_trailing_zeros (v_zero_bits, *tp);
    113   vp = TMP_ALLOC_LIMBS (vsize);
    114   if (v_zero_bits != 0)
    115     {
    116       mpn_rshift (vp, tp, vsize, v_zero_bits);
    117       vsize -= vp[vsize - 1] == 0;
    118     }
    119   else
    120     MPN_COPY (vp, tp, vsize);
    121 
    122   if (u_zero_limbs > v_zero_limbs)
    123     {
    124       g_zero_limbs = v_zero_limbs;
    125       g_zero_bits = v_zero_bits;
    126     }
    127   else
    128     {
    129       g_zero_limbs = u_zero_limbs;
    130       if (u_zero_limbs < v_zero_limbs)
    131 	g_zero_bits = u_zero_bits;
    132       else  /*  Equal.  */
    133 	g_zero_bits = MIN (u_zero_bits, v_zero_bits);
    134     }
    135 
    136   /*  Call mpn_gcd.  The 2nd argument must not have more bits than the 1st.  */
    137   vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
    138     ? mpn_gcd (vp, vp, vsize, up, usize)
    139     : mpn_gcd (vp, up, usize, vp, vsize);
    140 
    141   /*  Here G <-- V << (g_zero_limbs*GMP_LIMB_BITS + g_zero_bits).  */
    142   gsize = vsize + g_zero_limbs;
    143   if (g_zero_bits != 0)
    144     {
    145       mp_limb_t cy_limb;
    146       gsize += (vp[vsize - 1] >> (GMP_NUMB_BITS - g_zero_bits)) != 0;
    147       tp = MPZ_NEWALLOC (g, gsize);
    148       MPN_ZERO (tp, g_zero_limbs);
    149 
    150       tp = tp + g_zero_limbs;
    151       cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
    152       if (cy_limb != 0)
    153 	tp[vsize] = cy_limb;
    154     }
    155   else
    156     {
    157       tp = MPZ_NEWALLOC (g, gsize);
    158       MPN_ZERO (tp, g_zero_limbs);
    159       MPN_COPY (tp + g_zero_limbs, vp, vsize);
    160     }
    161 
    162   SIZ (g) = gsize;
    163   TMP_FREE;
    164 }
    165