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