ghc/runtime/gmp/mpz_gcd.c

   1 /* mpz_gcd -- Calculate the greatest common divisior of two integers.
   2
   3 Copyright (C) 1991 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 General Public License as published by
   9 the Free Software Foundation; either version 2, or (at your option)
  10 any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with the GNU MP Library; see the file COPYING.  If not, write to
  19 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  20
  21 #include "gmp.h"
  22 #include "gmp-impl.h"
  23 #include "longlong.h"
  24
  25 #ifndef BERKELEY_MP
  26 void
  27 #ifdef __STDC__
  28 mpz_gcd (MP_INT *w, const MP_INT *u, const MP_INT *v)
  29 #else
  30 mpz_gcd (w, u, v)
  31      MP_INT *w;
  32      const MP_INT *u;
  33      const MP_INT *v;
  34 #endif
  35 #else /* BERKELEY_MP */
  36 void
  37 #ifdef __STDC__
  38 gcd (const MP_INT *u, const MP_INT *v, MP_INT *w)
  39 #else
  40 gcd (u, v, w)
  41      const MP_INT *u;
  42      const MP_INT *v;
  43      MP_INT *w;
  44 #endif
  45 #endif /* BERKELEY_MP */
  46 {
  47   mp_size usize, vsize, wsize;
  48   mp_ptr up_in, vp_in;
  49   mp_ptr up, vp;
  50   mp_ptr wp;
  51   mp_size i;
  52   mp_limb d;
  53   int bcnt;
  54   mp_size w_bcnt;
  55   mp_limb cy_digit;
  56
  57   usize = ABS (u->size);
  58   vsize = ABS (v->size);
  59
  60   /* GCD(0,v) == v.  */
  61   if (usize == 0)
  62     {
  63       if (w->alloc < vsize)
  64         _mpz_realloc (w, vsize);
  65
  66       w->size = vsize;
  67       MPN_COPY (w->d, v->d, vsize);
  68       return;
  69     }
  70
  71   /* GCD(0,u) == u.  */
  72   if (vsize == 0)
  73     {
  74       if (w->alloc < usize)
  75         _mpz_realloc (w, usize);
  76
  77       w->size = usize;
  78       MPN_COPY (w->d, u->d, usize);
  79       return;
  80     }
  81
  82   /* Make U odd by shifting it down as many bit positions as there
  83      are zero bits.  Put the result in temporary space.  */
  84   up = (mp_ptr) alloca (usize * BYTES_PER_MP_LIMB);
  85   up_in = u->d;
  86   for (i = 0; (d = up_in[i]) == 0; i++)
  87     ;
  88   count_leading_zeros (bcnt, d & -d);
  89   bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
  90   usize = mpn_rshift (up, up_in + i, usize - i, bcnt);
  91
  92   bcnt += i * BITS_PER_MP_LIMB;
  93   w_bcnt = bcnt;
  94
  95   /* Make V odd by shifting it down as many bit positions as there
  96      are zero bits.  Put the result in temporary space.  */
  97   vp = (mp_ptr) alloca (vsize * BYTES_PER_MP_LIMB);
  98   vp_in = v->d;
  99   for (i = 0; (d = vp_in[i]) == 0; i++)
 100     ;
 101   count_leading_zeros (bcnt, d & -d);
 102   bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
 103   vsize = mpn_rshift (vp, vp_in + i, vsize - i, bcnt);
 104
 105   /* W_BCNT is set to the minimum of the number of zero bits in U and V.
 106      Thus it represents the number of common 2 factors.  */
 107   bcnt += i * BITS_PER_MP_LIMB;
 108   if (bcnt < w_bcnt)
 109     w_bcnt = bcnt;
 110
 111   for (;;)
 112     {
 113       int cmp;
 114
 115       cmp = usize - vsize != 0 ? usize - vsize : mpn_cmp (up, vp, usize);
 116
 117       /* If U and V have become equal, we have found the GCD.  */
 118       if (cmp == 0)
 119         break;
 120
 121       if (cmp > 0)
 122         {
 123           /* Replace U by (U - V) >> cnt, with cnt being the least value
 124              making U odd again.  */
 125
 126           usize += mpn_sub (up, up, usize, vp, vsize);
 127           for (i = 0; (d = up[i]) == 0; i++)
 128             ;
 129           count_leading_zeros (bcnt, d & -d);
 130           bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
 131           usize = mpn_rshift (up, up + i, usize - i, bcnt);
 132         }
 133       else
 134         {
 135           /* Replace V by (V - U) >> cnt, with cnt being the least value
 136              making V odd again.  */
 137
 138           vsize += mpn_sub (vp, vp, vsize, up, usize);
 139           for (i = 0; (d = vp[i]) == 0; i++)
 140             ;
 141           count_leading_zeros (bcnt, d & -d);
 142           bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
 143           vsize = mpn_rshift (vp, vp + i, vsize - i, bcnt);
 144         }
 145     }
 146
 147   /* GCD(U_IN, V_IN) now is U * 2**W_BCNT.  */
 148
 149   wsize = usize + w_bcnt / BITS_PER_MP_LIMB + 1;
 150   if (w->alloc < wsize)
 151     _mpz_realloc (w, wsize);
 152
 153   wp = w->d;
 154
 155   MPN_ZERO (wp, w_bcnt / BITS_PER_MP_LIMB);
 156
 157   cy_digit = mpn_lshift (wp + w_bcnt / BITS_PER_MP_LIMB, up, usize,
 158                           w_bcnt % BITS_PER_MP_LIMB);
 159   wsize = usize + w_bcnt / BITS_PER_MP_LIMB;
 160   if (cy_digit != 0)
 161     {
 162       wp[wsize] = cy_digit;
 163       wsize++;
 164     }
 165
 166   w->size = wsize;
 167
 168   alloca (0);
 169 }