1/* hgcd.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-2005, 2008, 2011, 2012 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 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 18or 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 24or both in parallel, as here. 25 26The GNU MP Library is distributed in the hope that it will be useful, but 27WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 28or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 29for more details. 30 31You should have received copies of the GNU General Public License and the 32GNU Lesser General Public License along with the GNU MP Library. If not, 33see https://www.gnu.org/licenses/. */ 34 35#include "gmp-impl.h" 36#include "longlong.h" 37 38 39/* Size analysis for hgcd: 40 41 For the recursive calls, we have n1 <= ceil(n / 2). Then the 42 storage need is determined by the storage for the recursive call 43 computing M1, and hgcd_matrix_adjust and hgcd_matrix_mul calls that use M1 44 (after this, the storage needed for M1 can be recycled). 45 46 Let S(r) denote the required storage. For M1 we need 4 * (ceil(n1/2) + 1) 47 = 4 * (ceil(n/4) + 1), for the hgcd_matrix_adjust call, we need n + 2, 48 and for the hgcd_matrix_mul, we may need 3 ceil(n/2) + 8. In total, 49 4 * ceil(n/4) + 3 ceil(n/2) + 12 <= 10 ceil(n/4) + 12. 50 51 For the recursive call, we need S(n1) = S(ceil(n/2)). 52 53 S(n) <= 10*ceil(n/4) + 12 + S(ceil(n/2)) 54 <= 10*(ceil(n/4) + ... + ceil(n/2^(1+k))) + 12k + S(ceil(n/2^k)) 55 <= 10*(2 ceil(n/4) + k) + 12k + S(ceil(n/2^k)) 56 <= 20 ceil(n/4) + 22k + S(ceil(n/2^k)) 57*/ 58 59mp_size_t 60mpn_hgcd_itch (mp_size_t n) 61{ 62 unsigned k; 63 int count; 64 mp_size_t nscaled; 65 66 if (BELOW_THRESHOLD (n, HGCD_THRESHOLD)) 67 return n; 68 69 /* Get the recursion depth. */ 70 nscaled = (n - 1) / (HGCD_THRESHOLD - 1); 71 count_leading_zeros (count, nscaled); 72 k = GMP_LIMB_BITS - count; 73 74 return 20 * ((n+3) / 4) + 22 * k + HGCD_THRESHOLD; 75} 76 77/* Reduces a,b until |a-b| fits in n/2 + 1 limbs. Constructs matrix M 78 with elements of size at most (n+1)/2 - 1. Returns new size of a, 79 b, or zero if no reduction is possible. */ 80 81mp_size_t 82mpn_hgcd (mp_ptr ap, mp_ptr bp, mp_size_t n, 83 struct hgcd_matrix *M, mp_ptr tp) 84{ 85 mp_size_t s = n/2 + 1; 86 87 mp_size_t nn; 88 int success = 0; 89 90 if (n <= s) 91 /* Happens when n <= 2, a fairly uninteresting case but exercised 92 by the random inputs of the testsuite. */ 93 return 0; 94 95 ASSERT ((ap[n-1] | bp[n-1]) > 0); 96 97 ASSERT ((n+1)/2 - 1 < M->alloc); 98 99 if (ABOVE_THRESHOLD (n, HGCD_THRESHOLD)) 100 { 101 mp_size_t n2 = (3*n)/4 + 1; 102 mp_size_t p = n/2; 103 104 nn = mpn_hgcd_reduce (M, ap, bp, n, p, tp); 105 if (nn) 106 { 107 n = nn; 108 success = 1; 109 } 110 111 /* NOTE: It appears this loop never runs more than once (at 112 least when not recursing to hgcd_appr). */ 113 while (n > n2) 114 { 115 /* Needs n + 1 storage */ 116 nn = mpn_hgcd_step (n, ap, bp, s, M, tp); 117 if (!nn) 118 return success ? n : 0; 119 120 n = nn; 121 success = 1; 122 } 123 124 if (n > s + 2) 125 { 126 struct hgcd_matrix M1; 127 mp_size_t scratch; 128 129 p = 2*s - n + 1; 130 scratch = MPN_HGCD_MATRIX_INIT_ITCH (n-p); 131 132 mpn_hgcd_matrix_init(&M1, n - p, tp); 133 134 /* FIXME: Should use hgcd_reduce, but that may require more 135 scratch space, which requires review. */ 136 137 nn = mpn_hgcd (ap + p, bp + p, n - p, &M1, tp + scratch); 138 if (nn > 0) 139 { 140 /* We always have max(M) > 2^{-(GMP_NUMB_BITS + 1)} max(M1) */ 141 ASSERT (M->n + 2 >= M1.n); 142 143 /* Furthermore, assume M ends with a quotient (1, q; 0, 1), 144 then either q or q + 1 is a correct quotient, and M1 will 145 start with either (1, 0; 1, 1) or (2, 1; 1, 1). This 146 rules out the case that the size of M * M1 is much 147 smaller than the expected M->n + M1->n. */ 148 149 ASSERT (M->n + M1.n < M->alloc); 150 151 /* Needs 2 (p + M->n) <= 2 (2*s - n2 + 1 + n2 - s - 1) 152 = 2*s <= 2*(floor(n/2) + 1) <= n + 2. */ 153 n = mpn_hgcd_matrix_adjust (&M1, p + nn, ap, bp, p, tp + scratch); 154 155 /* We need a bound for of M->n + M1.n. Let n be the original 156 input size. Then 157 158 ceil(n/2) - 1 >= size of product >= M.n + M1.n - 2 159 160 and it follows that 161 162 M.n + M1.n <= ceil(n/2) + 1 163 164 Then 3*(M.n + M1.n) + 5 <= 3 * ceil(n/2) + 8 is the 165 amount of needed scratch space. */ 166 mpn_hgcd_matrix_mul (M, &M1, tp + scratch); 167 success = 1; 168 } 169 } 170 } 171 172 for (;;) 173 { 174 /* Needs s+3 < n */ 175 nn = mpn_hgcd_step (n, ap, bp, s, M, tp); 176 if (!nn) 177 return success ? n : 0; 178 179 n = nn; 180 success = 1; 181 } 182} 183