1 /* mpn_gcdext -- Extended Greatest Common Divisor. 2 3 Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2005, 2008, 2009 Free Software 4 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 the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MP Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along 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 25 /* Temporary storage: 3*(n+1) for u. n+1 for the matrix-vector 26 multiplications (if hgcd2 succeeds). If hgcd fails, n+1 limbs are 27 needed for the division, with most n for the quotient, and n+1 for 28 the product q u0. In all, 4n + 3. */ 29 30 mp_size_t 31 mpn_gcdext_lehmer_n (mp_ptr gp, mp_ptr up, mp_size_t *usize, 32 mp_ptr ap, mp_ptr bp, mp_size_t n, 33 mp_ptr tp) 34 { 35 mp_size_t ualloc = n + 1; 36 37 /* Keeps track of the second row of the reduction matrix 38 * 39 * M = (v0, v1 ; u0, u1) 40 * 41 * which correspond to the first column of the inverse 42 * 43 * M^{-1} = (u1, -v1; -u0, v0) 44 */ 45 46 mp_size_t un; 47 mp_ptr u0; 48 mp_ptr u1; 49 mp_ptr u2; 50 51 MPN_ZERO (tp, 3*ualloc); 52 u0 = tp; tp += ualloc; 53 u1 = tp; tp += ualloc; 54 u2 = tp; tp += ualloc; 55 56 u1[0] = 1; un = 1; 57 58 /* FIXME: Handle n == 2 differently, after the loop? */ 59 while (n >= 2) 60 { 61 struct hgcd_matrix1 M; 62 mp_limb_t ah, al, bh, bl; 63 mp_limb_t mask; 64 65 mask = ap[n-1] | bp[n-1]; 66 ASSERT (mask > 0); 67 68 if (mask & GMP_NUMB_HIGHBIT) 69 { 70 ah = ap[n-1]; al = ap[n-2]; 71 bh = bp[n-1]; bl = bp[n-2]; 72 } 73 else if (n == 2) 74 { 75 /* We use the full inputs without truncation, so we can 76 safely shift left. */ 77 int shift; 78 79 count_leading_zeros (shift, mask); 80 ah = MPN_EXTRACT_NUMB (shift, ap[1], ap[0]); 81 al = ap[0] << shift; 82 bh = MPN_EXTRACT_NUMB (shift, bp[1], bp[0]); 83 bl = bp[0] << shift; 84 } 85 else 86 { 87 int shift; 88 89 count_leading_zeros (shift, mask); 90 ah = MPN_EXTRACT_NUMB (shift, ap[n-1], ap[n-2]); 91 al = MPN_EXTRACT_NUMB (shift, ap[n-2], ap[n-3]); 92 bh = MPN_EXTRACT_NUMB (shift, bp[n-1], bp[n-2]); 93 bl = MPN_EXTRACT_NUMB (shift, bp[n-2], bp[n-3]); 94 } 95 96 /* Try an mpn_nhgcd2 step */ 97 if (mpn_hgcd2 (ah, al, bh, bl, &M)) 98 { 99 n = mpn_hgcd_mul_matrix1_inverse_vector (&M, tp, ap, bp, n); 100 MP_PTR_SWAP (ap, tp); 101 un = mpn_hgcd_mul_matrix1_vector(&M, u2, u0, u1, un); 102 MP_PTR_SWAP (u0, u2); 103 } 104 else 105 { 106 /* mpn_hgcd2 has failed. Then either one of a or b is very 107 small, or the difference is very small. Perform one 108 subtraction followed by one division. */ 109 mp_size_t gn; 110 mp_size_t updated_un = un; 111 112 /* Temporary storage n for the quotient and ualloc for the 113 new cofactor. */ 114 n = mpn_gcdext_subdiv_step (gp, &gn, up, usize, ap, bp, n, 115 u0, u1, &updated_un, tp, u2); 116 if (n == 0) 117 return gn; 118 119 un = updated_un; 120 } 121 } 122 ASSERT_ALWAYS (ap[0] > 0); 123 ASSERT_ALWAYS (bp[0] > 0); 124 125 if (ap[0] == bp[0]) 126 { 127 int c; 128 129 /* Which cofactor to return now? Candidates are +u1 and -u0, 130 depending on which of a and b was most recently reduced, 131 which we don't keep track of. So compare and get the smallest 132 one. */ 133 134 gp[0] = ap[0]; 135 136 MPN_CMP (c, u0, u1, un); 137 ASSERT (c != 0 || (un == 1 && u0[0] == 1 && u1[0] == 1)); 138 if (c < 0) 139 { 140 MPN_NORMALIZE (u0, un); 141 MPN_COPY (up, u0, un); 142 *usize = -un; 143 } 144 else 145 { 146 MPN_NORMALIZE_NOT_ZERO (u1, un); 147 MPN_COPY (up, u1, un); 148 *usize = un; 149 } 150 return 1; 151 } 152 else 153 { 154 mp_limb_t uh, vh; 155 mp_limb_signed_t u; 156 mp_limb_signed_t v; 157 int negate; 158 159 gp[0] = mpn_gcdext_1 (&u, &v, ap[0], bp[0]); 160 161 /* Set up = u u1 - v u0. Keep track of size, un grows by one or 162 two limbs. */ 163 164 if (u == 0) 165 { 166 ASSERT (v == 1); 167 MPN_NORMALIZE (u0, un); 168 MPN_COPY (up, u0, un); 169 *usize = -un; 170 return 1; 171 } 172 else if (v == 0) 173 { 174 ASSERT (u == 1); 175 MPN_NORMALIZE (u1, un); 176 MPN_COPY (up, u1, un); 177 *usize = un; 178 return 1; 179 } 180 else if (u > 0) 181 { 182 negate = 0; 183 ASSERT (v < 0); 184 v = -v; 185 } 186 else 187 { 188 negate = 1; 189 ASSERT (v > 0); 190 u = -u; 191 } 192 193 uh = mpn_mul_1 (up, u1, un, u); 194 vh = mpn_addmul_1 (up, u0, un, v); 195 196 if ( (uh | vh) > 0) 197 { 198 uh += vh; 199 up[un++] = uh; 200 if (uh < vh) 201 up[un++] = 1; 202 } 203 204 MPN_NORMALIZE_NOT_ZERO (up, un); 205 206 *usize = negate ? -un : un; 207 return 1; 208 } 209 } 210