contrib/libgmp/mpz/gcd.c

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