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