1 /* gcdext_subdiv_step.c. 2 3 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY 4 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST 5 GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE. 6 7 Copyright 2003, 2004, 2005, 2008, 2009 Free Software Foundation, Inc. 8 9 This file is part of the GNU MP Library. 10 11 The GNU MP Library is free software; you can redistribute it and/or modify 12 it under the terms of the GNU Lesser General Public License as published by 13 the Free Software Foundation; either version 3 of the License, or (at your 14 option) any later version. 15 16 The GNU MP Library is distributed in the hope that it will be useful, but 17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 19 License for more details. 20 21 You should have received a copy of the GNU Lesser General Public License 22 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 23 24 #include "gmp.h" 25 #include "gmp-impl.h" 26 #include "longlong.h" 27 28 /* Used when mpn_hgcd or mpn_hgcd2 has failed. Then either one of a or 29 b is small, or the difference is small. Perform one subtraction 30 followed by one division. If the gcd is found, stores it in gp and 31 *gn, and returns zero. Otherwise, compute the reduced a and b, 32 return the new size, and cofactors. */ 33 34 /* Temporary storage: Needs n limbs for the quotient, at qp. tp must 35 point to an area large enough for the resulting cofactor, plus one 36 limb extra. All in all, 2N + 1 if N is a bound for both inputs and 37 outputs. */ 38 mp_size_t 39 mpn_gcdext_subdiv_step (mp_ptr gp, mp_size_t *gn, mp_ptr up, mp_size_t *usizep, 40 mp_ptr ap, mp_ptr bp, mp_size_t n, 41 mp_ptr u0, mp_ptr u1, mp_size_t *unp, 42 mp_ptr qp, mp_ptr tp) 43 { 44 mp_size_t an, bn, un; 45 mp_size_t qn; 46 mp_size_t u0n; 47 48 int swapped; 49 50 an = bn = n; 51 52 ASSERT (an > 0); 53 ASSERT (ap[an-1] > 0 || bp[an-1] > 0); 54 55 MPN_NORMALIZE (ap, an); 56 MPN_NORMALIZE (bp, bn); 57 58 un = *unp; 59 60 swapped = 0; 61 62 if (UNLIKELY (an == 0)) 63 { 64 return_b: 65 MPN_COPY (gp, bp, bn); 66 *gn = bn; 67 68 MPN_NORMALIZE (u0, un); 69 MPN_COPY (up, u0, un); 70 71 *usizep = swapped ? un : -un; 72 73 return 0; 74 } 75 else if (UNLIKELY (bn == 0)) 76 { 77 MPN_COPY (gp, ap, an); 78 *gn = an; 79 80 MPN_NORMALIZE (u1, un); 81 MPN_COPY (up, u1, un); 82 83 *usizep = swapped ? -un : un; 84 85 return 0; 86 } 87 88 /* Arrange so that a > b, subtract an -= bn, and maintain 89 normalization. */ 90 if (an < bn) 91 { 92 MPN_PTR_SWAP (ap, an, bp, bn); 93 MP_PTR_SWAP (u0, u1); 94 swapped ^= 1; 95 } 96 else if (an == bn) 97 { 98 int c; 99 MPN_CMP (c, ap, bp, an); 100 if (UNLIKELY (c == 0)) 101 { 102 MPN_COPY (gp, ap, an); 103 *gn = an; 104 105 /* Must return the smallest cofactor, +u1 or -u0 */ 106 MPN_CMP (c, u0, u1, un); 107 ASSERT (c != 0 || (un == 1 && u0[0] == 1 && u1[0] == 1)); 108 109 if (c < 0) 110 { 111 MPN_NORMALIZE (u0, un); 112 MPN_COPY (up, u0, un); 113 swapped ^= 1; 114 } 115 else 116 { 117 MPN_NORMALIZE_NOT_ZERO (u1, un); 118 MPN_COPY (up, u1, un); 119 } 120 121 *usizep = swapped ? -un : un; 122 return 0; 123 } 124 else if (c < 0) 125 { 126 MP_PTR_SWAP (ap, bp); 127 MP_PTR_SWAP (u0, u1); 128 swapped ^= 1; 129 } 130 } 131 /* Reduce a -= b, u1 += u0 */ 132 ASSERT_NOCARRY (mpn_sub (ap, ap, an, bp, bn)); 133 MPN_NORMALIZE (ap, an); 134 ASSERT (an > 0); 135 136 u1[un] = mpn_add_n (u1, u1, u0, un); 137 un += (u1[un] > 0); 138 139 /* Arrange so that a > b, and divide a = q b + r */ 140 if (an < bn) 141 { 142 MPN_PTR_SWAP (ap, an, bp, bn); 143 MP_PTR_SWAP (u0, u1); 144 swapped ^= 1; 145 } 146 else if (an == bn) 147 { 148 int c; 149 MPN_CMP (c, ap, bp, an); 150 if (UNLIKELY (c == 0)) 151 goto return_b; 152 else if (c < 0) 153 { 154 MP_PTR_SWAP (ap, bp); 155 MP_PTR_SWAP (u0, u1); 156 swapped ^= 1; 157 } 158 } 159 160 /* Reduce a -= q b, u1 += q u0 */ 161 qn = an - bn + 1; 162 mpn_tdiv_qr (qp, ap, 0, ap, an, bp, bn); 163 164 if (mpn_zero_p (ap, bn)) 165 goto return_b; 166 167 n = bn; 168 169 /* Update u1 += q u0 */ 170 u0n = un; 171 MPN_NORMALIZE (u0, u0n); 172 173 if (u0n > 0) 174 { 175 qn -= (qp[qn - 1] == 0); 176 177 if (qn > u0n) 178 mpn_mul (tp, qp, qn, u0, u0n); 179 else 180 mpn_mul (tp, u0, u0n, qp, qn); 181 182 if (qn + u0n > un) 183 { 184 mp_size_t u1n = un; 185 un = qn + u0n; 186 un -= (tp[un-1] == 0); 187 u1[un] = mpn_add (u1, tp, un, u1, u1n); 188 } 189 else 190 { 191 u1[un] = mpn_add (u1, u1, un, tp, qn + u0n); 192 } 193 194 un += (u1[un] > 0); 195 } 196 197 *unp = un; 198 return n; 199 } 200