1 /* hgcd_appr.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 2011, 2012 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 "gmp-impl.h" 36 #include "longlong.h" 37 38 /* Identical to mpn_hgcd_itch. FIXME: Do we really need to add 39 HGCD_THRESHOLD at the end? */ 40 mp_size_t 41 mpn_hgcd_appr_itch (mp_size_t n) 42 { 43 if (BELOW_THRESHOLD (n, HGCD_APPR_THRESHOLD)) 44 return n; 45 else 46 { 47 unsigned k; 48 int count; 49 mp_size_t nscaled; 50 51 /* Get the recursion depth. */ 52 nscaled = (n - 1) / (HGCD_APPR_THRESHOLD - 1); 53 count_leading_zeros (count, nscaled); 54 k = GMP_LIMB_BITS - count; 55 56 return 20 * ((n+3) / 4) + 22 * k + HGCD_THRESHOLD; 57 } 58 } 59 60 /* Destroys inputs. */ 61 int 62 mpn_hgcd_appr (mp_ptr ap, mp_ptr bp, mp_size_t n, 63 struct hgcd_matrix *M, mp_ptr tp) 64 { 65 mp_size_t s; 66 int success = 0; 67 68 ASSERT (n > 0); 69 70 ASSERT ((ap[n-1] | bp[n-1]) != 0); 71 72 if (n <= 2) 73 /* Implies s = n. A fairly uninteresting case but exercised by the 74 random inputs of the testsuite. */ 75 return 0; 76 77 ASSERT ((n+1)/2 - 1 < M->alloc); 78 79 /* We aim for reduction of to GMP_NUMB_BITS * s bits. But each time 80 we discard some of the least significant limbs, we must keep one 81 additional bit to account for the truncation error. We maintain 82 the GMP_NUMB_BITS * s - extra_bits as the current target size. */ 83 84 s = n/2 + 1; 85 if (BELOW_THRESHOLD (n, HGCD_APPR_THRESHOLD)) 86 { 87 unsigned extra_bits = 0; 88 89 while (n > 2) 90 { 91 mp_size_t nn; 92 93 ASSERT (n > s); 94 ASSERT (n <= 2*s); 95 96 nn = mpn_hgcd_step (n, ap, bp, s, M, tp); 97 if (!nn) 98 break; 99 100 n = nn; 101 success = 1; 102 103 /* We can truncate and discard the lower p bits whenever nbits <= 104 2*sbits - p. To account for the truncation error, we must 105 adjust 106 107 sbits <-- sbits + 1 - p, 108 109 rather than just sbits <-- sbits - p. This adjustment makes 110 the produced matrix slightly smaller than it could be. */ 111 112 if (GMP_NUMB_BITS * (n + 1) + 2 * extra_bits <= 2*GMP_NUMB_BITS * s) 113 { 114 mp_size_t p = (GMP_NUMB_BITS * (2*s - n) - 2*extra_bits) / GMP_NUMB_BITS; 115 116 if (extra_bits == 0) 117 { 118 /* We cross a limb boundary and bump s. We can't do that 119 if the result is that it makes makes min(U, V) 120 smaller than 2^{GMP_NUMB_BITS} s. */ 121 if (s + 1 == n 122 || mpn_zero_p (ap + s + 1, n - s - 1) 123 || mpn_zero_p (bp + s + 1, n - s - 1)) 124 continue; 125 126 extra_bits = GMP_NUMB_BITS - 1; 127 s++; 128 } 129 else 130 { 131 extra_bits--; 132 } 133 134 /* Drop the p least significant limbs */ 135 ap += p; bp += p; n -= p; s -= p; 136 } 137 } 138 139 ASSERT (s > 0); 140 141 if (extra_bits > 0) 142 { 143 /* We can get here only of we have dropped at least one of the least 144 significant bits, so we can decrement ap and bp. We can then shift 145 left extra bits using mpn_rshift. */ 146 /* NOTE: In the unlikely case that n is large, it would be preferable 147 to do an initial subdiv step to reduce the size before shifting, 148 but that would mean duplicating mpn_gcd_subdiv_step with a bit 149 count rather than a limb count. */ 150 ap--; bp--; 151 ap[0] = mpn_rshift (ap+1, ap+1, n, GMP_NUMB_BITS - extra_bits); 152 bp[0] = mpn_rshift (bp+1, bp+1, n, GMP_NUMB_BITS - extra_bits); 153 n += (ap[n] | bp[n]) > 0; 154 155 ASSERT (success); 156 157 while (n > 2) 158 { 159 mp_size_t nn; 160 161 ASSERT (n > s); 162 ASSERT (n <= 2*s); 163 164 nn = mpn_hgcd_step (n, ap, bp, s, M, tp); 165 166 if (!nn) 167 return 1; 168 169 n = nn; 170 } 171 } 172 173 if (n == 2) 174 { 175 struct hgcd_matrix1 M1; 176 ASSERT (s == 1); 177 178 if (mpn_hgcd2 (ap[1], ap[0], bp[1], bp[0], &M1)) 179 { 180 /* Multiply M <- M * M1 */ 181 mpn_hgcd_matrix_mul_1 (M, &M1, tp); 182 success = 1; 183 } 184 } 185 return success; 186 } 187 else 188 { 189 mp_size_t n2 = (3*n)/4 + 1; 190 mp_size_t p = n/2; 191 mp_size_t nn; 192 193 nn = mpn_hgcd_reduce (M, ap, bp, n, p, tp); 194 if (nn) 195 { 196 n = nn; 197 /* FIXME: Discard some of the low limbs immediately? */ 198 success = 1; 199 } 200 201 while (n > n2) 202 { 203 mp_size_t nn; 204 205 /* Needs n + 1 storage */ 206 nn = mpn_hgcd_step (n, ap, bp, s, M, tp); 207 if (!nn) 208 return success; 209 210 n = nn; 211 success = 1; 212 } 213 if (n > s + 2) 214 { 215 struct hgcd_matrix M1; 216 mp_size_t scratch; 217 218 p = 2*s - n + 1; 219 scratch = MPN_HGCD_MATRIX_INIT_ITCH (n-p); 220 221 mpn_hgcd_matrix_init(&M1, n - p, tp); 222 if (mpn_hgcd_appr (ap + p, bp + p, n - p, &M1, tp + scratch)) 223 { 224 /* We always have max(M) > 2^{-(GMP_NUMB_BITS + 1)} max(M1) */ 225 ASSERT (M->n + 2 >= M1.n); 226 227 /* Furthermore, assume M ends with a quotient (1, q; 0, 1), 228 then either q or q + 1 is a correct quotient, and M1 will 229 start with either (1, 0; 1, 1) or (2, 1; 1, 1). This 230 rules out the case that the size of M * M1 is much 231 smaller than the expected M->n + M1->n. */ 232 233 ASSERT (M->n + M1.n < M->alloc); 234 235 /* We need a bound for of M->n + M1.n. Let n be the original 236 input size. Then 237 238 ceil(n/2) - 1 >= size of product >= M.n + M1.n - 2 239 240 and it follows that 241 242 M.n + M1.n <= ceil(n/2) + 1 243 244 Then 3*(M.n + M1.n) + 5 <= 3 * ceil(n/2) + 8 is the 245 amount of needed scratch space. */ 246 mpn_hgcd_matrix_mul (M, &M1, tp + scratch); 247 return 1; 248 } 249 } 250 251 for(;;) 252 { 253 mp_size_t nn; 254 255 ASSERT (n > s); 256 ASSERT (n <= 2*s); 257 258 nn = mpn_hgcd_step (n, ap, bp, s, M, tp); 259 260 if (!nn) 261 return success; 262 263 n = nn; 264 success = 1; 265 } 266 } 267 } 268