1 /* gcd_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-2005, 2008, 2010, 2011 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 either: 13 14 * the GNU Lesser General Public License as published by the Free 15 Software Foundation; either version 3 of the License, or (at your 16 option) any later version. 17 18 or 19 20 * the GNU General Public License as published by the Free Software 21 Foundation; either version 2 of the License, or (at your option) any 22 later version. 23 24 or both in parallel, as here. 25 26 The GNU MP Library is distributed in the hope that it will be useful, but 27 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 28 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 29 for more details. 30 31 You should have received copies of the GNU General Public License and the 32 GNU Lesser General Public License along with the GNU MP Library. If not, 33 see https://www.gnu.org/licenses/. */ 34 35 #include <stdlib.h> /* for NULL */ 36 37 #include "gmp-impl.h" 38 #include "longlong.h" 39 40 /* Used when mpn_hgcd or mpn_hgcd2 has failed. Then either one of a or 41 b is small, or the difference is small. Perform one subtraction 42 followed by one division. The normal case is to compute the reduced 43 a and b, and return the new size. 44 45 If s == 0 (used for gcd and gcdext), returns zero if the gcd is 46 found. 47 48 If s > 0, don't reduce to size <= s, and return zero if no 49 reduction is possible (if either a, b or |a-b| is of size <= s). */ 50 51 /* The hook function is called as 52 53 hook(ctx, gp, gn, qp, qn, d) 54 55 in the following cases: 56 57 + If A = B at the start, G is the gcd, Q is NULL, d = -1. 58 59 + If one input is zero at the start, G is the gcd, Q is NULL, 60 d = 0 if A = G and d = 1 if B = G. 61 62 Otherwise, if d = 0 we have just subtracted a multiple of A from B, 63 and if d = 1 we have subtracted a multiple of B from A. 64 65 + If A = B after subtraction, G is the gcd, Q is NULL. 66 67 + If we get a zero remainder after division, G is the gcd, Q is the 68 quotient. 69 70 + Otherwise, G is NULL, Q is the quotient (often 1). 71 72 */ 73 74 mp_size_t 75 mpn_gcd_subdiv_step (mp_ptr ap, mp_ptr bp, mp_size_t n, mp_size_t s, 76 gcd_subdiv_step_hook *hook, void *ctx, 77 mp_ptr tp) 78 { 79 static const mp_limb_t one = CNST_LIMB(1); 80 mp_size_t an, bn, qn; 81 82 int swapped; 83 84 ASSERT (n > 0); 85 ASSERT (ap[n-1] > 0 || bp[n-1] > 0); 86 87 an = bn = n; 88 MPN_NORMALIZE (ap, an); 89 MPN_NORMALIZE (bp, bn); 90 91 swapped = 0; 92 93 /* Arrange so that a < b, subtract b -= a, and maintain 94 normalization. */ 95 if (an == bn) 96 { 97 int c; 98 MPN_CMP (c, ap, bp, an); 99 if (UNLIKELY (c == 0)) 100 { 101 /* For gcdext, return the smallest of the two cofactors, so 102 pass d = -1. */ 103 if (s == 0) 104 hook (ctx, ap, an, NULL, 0, -1); 105 return 0; 106 } 107 else if (c > 0) 108 { 109 MP_PTR_SWAP (ap, bp); 110 swapped ^= 1; 111 } 112 } 113 else 114 { 115 if (an > bn) 116 { 117 MPN_PTR_SWAP (ap, an, bp, bn); 118 swapped ^= 1; 119 } 120 } 121 if (an <= s) 122 { 123 if (s == 0) 124 hook (ctx, bp, bn, NULL, 0, swapped ^ 1); 125 return 0; 126 } 127 128 ASSERT_NOCARRY (mpn_sub (bp, bp, bn, ap, an)); 129 MPN_NORMALIZE (bp, bn); 130 ASSERT (bn > 0); 131 132 if (bn <= s) 133 { 134 /* Undo subtraction. */ 135 mp_limb_t cy = mpn_add (bp, ap, an, bp, bn); 136 if (cy > 0) 137 bp[an] = cy; 138 return 0; 139 } 140 141 /* Arrange so that a < b */ 142 if (an == bn) 143 { 144 int c; 145 MPN_CMP (c, ap, bp, an); 146 if (UNLIKELY (c == 0)) 147 { 148 if (s > 0) 149 /* Just record subtraction and return */ 150 hook (ctx, NULL, 0, &one, 1, swapped); 151 else 152 /* Found gcd. */ 153 hook (ctx, bp, bn, NULL, 0, swapped); 154 return 0; 155 } 156 157 hook (ctx, NULL, 0, &one, 1, swapped); 158 159 if (c > 0) 160 { 161 MP_PTR_SWAP (ap, bp); 162 swapped ^= 1; 163 } 164 } 165 else 166 { 167 hook (ctx, NULL, 0, &one, 1, swapped); 168 169 if (an > bn) 170 { 171 MPN_PTR_SWAP (ap, an, bp, bn); 172 swapped ^= 1; 173 } 174 } 175 176 mpn_tdiv_qr (tp, bp, 0, bp, bn, ap, an); 177 qn = bn - an + 1; 178 bn = an; 179 MPN_NORMALIZE (bp, bn); 180 181 if (UNLIKELY (bn <= s)) 182 { 183 if (s == 0) 184 { 185 hook (ctx, ap, an, tp, qn, swapped); 186 return 0; 187 } 188 189 /* Quotient is one too large, so decrement it and add back A. */ 190 if (bn > 0) 191 { 192 mp_limb_t cy = mpn_add (bp, ap, an, bp, bn); 193 if (cy) 194 bp[an++] = cy; 195 } 196 else 197 MPN_COPY (bp, ap, an); 198 199 MPN_DECR_U (tp, qn, 1); 200 } 201 202 hook (ctx, NULL, 0, tp, qn, swapped); 203 return an; 204 } 205