rts/gmp/mpz/gcd.c

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