1 /* mpn_gcdext -- Extended Greatest Common Divisor. 2 3 Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2005, 2008 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 if (ap[0] == 0) 123 { 124 gp[0] = bp[0]; 125 126 MPN_NORMALIZE_NOT_ZERO (u0, un); 127 MPN_COPY (up, u0, un); 128 129 *usize = -un; 130 return 1; 131 } 132 else if (bp[0] == 0) 133 { 134 gp[0] = ap[0]; 135 136 MPN_NORMALIZE_NOT_ZERO (u1, un); 137 MPN_COPY (up, u1, un); 138 139 *usize = un; 140 return 1; 141 } 142 else 143 { 144 mp_limb_t uh, vh; 145 mp_limb_t u; 146 mp_limb_t v; 147 148 gp[0] = mpn_gcdext_1 (&u, &v, ap[0], bp[0]); 149 150 /* Set up = u u1 + v u0. Keep track of size, un grows by one or 151 two limbs. */ 152 uh = mpn_mul_1 (up, u1, un, u); 153 vh = mpn_addmul_1 (up, u0, un, v); 154 155 if ( (uh | vh) > 0) 156 { 157 uh += vh; 158 up[un++] = uh; 159 if (uh < vh) 160 up[un++] = 1; 161 } 162 163 MPN_NORMALIZE_NOT_ZERO (up, un); 164 165 *usize = un; 166 return 1; 167 } 168 } 169