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 7Copyright 2003, 2004, 2005, 2008, 2009 Free Software Foundation, Inc. 8 9This file is part of the GNU MP Library. 10 11The GNU MP Library is free software; you can redistribute it and/or modify 12it under the terms of the GNU Lesser General Public License as published by 13the Free Software Foundation; either version 3 of the License, or (at your 14option) any later version. 15 16The GNU MP Library is distributed in the hope that it will be useful, but 17WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 18or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 19License for more details. 20 21You should have received a copy of the GNU Lesser General Public License 22along 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. */ 38mp_size_t 39mpn_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