1 /* hgcd_matrix.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, 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 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 "gmp.h" 25 #include "gmp-impl.h" 26 #include "longlong.h" 27 28 /* For input of size n, matrix elements are of size at most ceil(n/2) 29 - 1, but we need two limbs extra. */ 30 void 31 mpn_hgcd_matrix_init (struct hgcd_matrix *M, mp_size_t n, mp_ptr p) 32 { 33 mp_size_t s = (n+1)/2 + 1; 34 M->alloc = s; 35 M->n = 1; 36 MPN_ZERO (p, 4 * s); 37 M->p[0][0] = p; 38 M->p[0][1] = p + s; 39 M->p[1][0] = p + 2 * s; 40 M->p[1][1] = p + 3 * s; 41 42 M->p[0][0][0] = M->p[1][1][0] = 1; 43 } 44 45 /* Update column COL, adding in Q * column (1-COL). Temporary storage: 46 * qn + n <= M->alloc, where n is the size of the largest element in 47 * column 1 - COL. */ 48 void 49 mpn_hgcd_matrix_update_q (struct hgcd_matrix *M, mp_srcptr qp, mp_size_t qn, 50 unsigned col, mp_ptr tp) 51 { 52 ASSERT (col < 2); 53 54 if (qn == 1) 55 { 56 mp_limb_t q = qp[0]; 57 mp_limb_t c0, c1; 58 59 c0 = mpn_addmul_1 (M->p[0][col], M->p[0][1-col], M->n, q); 60 c1 = mpn_addmul_1 (M->p[1][col], M->p[1][1-col], M->n, q); 61 62 M->p[0][col][M->n] = c0; 63 M->p[1][col][M->n] = c1; 64 65 M->n += (c0 | c1) != 0; 66 } 67 else 68 { 69 unsigned row; 70 71 /* Carries for the unlikely case that we get both high words 72 from the multiplication and carries from the addition. */ 73 mp_limb_t c[2]; 74 mp_size_t n; 75 76 /* The matrix will not necessarily grow in size by qn, so we 77 need normalization in order not to overflow M. */ 78 79 for (n = M->n; n + qn > M->n; n--) 80 { 81 ASSERT (n > 0); 82 if (M->p[0][1-col][n-1] > 0 || M->p[1][1-col][n-1] > 0) 83 break; 84 } 85 86 ASSERT (qn + n <= M->alloc); 87 88 for (row = 0; row < 2; row++) 89 { 90 if (qn <= n) 91 mpn_mul (tp, M->p[row][1-col], n, qp, qn); 92 else 93 mpn_mul (tp, qp, qn, M->p[row][1-col], n); 94 95 ASSERT (n + qn >= M->n); 96 c[row] = mpn_add (M->p[row][col], tp, n + qn, M->p[row][col], M->n); 97 } 98 99 n += qn; 100 101 if (c[0] | c[1]) 102 { 103 M->p[0][col][n] = c[0]; 104 M->p[1][col][n] = c[1]; 105 n++; 106 } 107 else 108 { 109 n -= (M->p[0][col][n-1] | M->p[1][col][n-1]) == 0; 110 ASSERT (n >= M->n); 111 } 112 M->n = n; 113 } 114 115 ASSERT (M->n < M->alloc); 116 } 117 118 /* Multiply M by M1 from the right. Since the M1 elements fit in 119 GMP_NUMB_BITS - 1 bits, M grows by at most one limb. Needs 120 temporary space M->n */ 121 void 122 mpn_hgcd_matrix_mul_1 (struct hgcd_matrix *M, const struct hgcd_matrix1 *M1, 123 mp_ptr tp) 124 { 125 mp_size_t n0, n1; 126 127 /* Could avoid copy by some swapping of pointers. */ 128 MPN_COPY (tp, M->p[0][0], M->n); 129 n0 = mpn_hgcd_mul_matrix1_vector (M1, M->p[0][0], tp, M->p[0][1], M->n); 130 MPN_COPY (tp, M->p[1][0], M->n); 131 n1 = mpn_hgcd_mul_matrix1_vector (M1, M->p[1][0], tp, M->p[1][1], M->n); 132 133 /* Depends on zero initialization */ 134 M->n = MAX(n0, n1); 135 ASSERT (M->n < M->alloc); 136 } 137 138 /* Multiply M by M1 from the right. Needs 3*(M->n + M1->n) + 5 limbs 139 of temporary storage (see mpn_matrix22_mul_itch). */ 140 void 141 mpn_hgcd_matrix_mul (struct hgcd_matrix *M, const struct hgcd_matrix *M1, 142 mp_ptr tp) 143 { 144 mp_size_t n; 145 146 /* About the new size of M:s elements. Since M1's diagonal elements 147 are > 0, no element can decrease. The new elements are of size 148 M->n + M1->n, one limb more or less. The computation of the 149 matrix product produces elements of size M->n + M1->n + 1. But 150 the true size, after normalization, may be three limbs smaller. 151 152 The reason that the product has normalized size >= M->n + M1->n - 153 2 is subtle. It depends on the fact that M and M1 can be factored 154 as products of (1,1; 0,1) and (1,0; 1,1), and that we can't have 155 M ending with a large power and M1 starting with a large power of 156 the same matrix. */ 157 158 /* FIXME: Strassen multiplication gives only a small speedup. In FFT 159 multiplication range, this function could be sped up quite a lot 160 using invariance. */ 161 ASSERT (M->n + M1->n < M->alloc); 162 163 ASSERT ((M->p[0][0][M->n-1] | M->p[0][1][M->n-1] 164 | M->p[1][0][M->n-1] | M->p[1][1][M->n-1]) > 0); 165 166 ASSERT ((M1->p[0][0][M1->n-1] | M1->p[0][1][M1->n-1] 167 | M1->p[1][0][M1->n-1] | M1->p[1][1][M1->n-1]) > 0); 168 169 mpn_matrix22_mul (M->p[0][0], M->p[0][1], 170 M->p[1][0], M->p[1][1], M->n, 171 M1->p[0][0], M1->p[0][1], 172 M1->p[1][0], M1->p[1][1], M1->n, tp); 173 174 /* Index of last potentially non-zero limb, size is one greater. */ 175 n = M->n + M1->n; 176 177 n -= ((M->p[0][0][n] | M->p[0][1][n] | M->p[1][0][n] | M->p[1][1][n]) == 0); 178 n -= ((M->p[0][0][n] | M->p[0][1][n] | M->p[1][0][n] | M->p[1][1][n]) == 0); 179 n -= ((M->p[0][0][n] | M->p[0][1][n] | M->p[1][0][n] | M->p[1][1][n]) == 0); 180 181 ASSERT ((M->p[0][0][n] | M->p[0][1][n] | M->p[1][0][n] | M->p[1][1][n]) > 0); 182 183 M->n = n + 1; 184 } 185 186 /* Multiplies the least significant p limbs of (a;b) by M^-1. 187 Temporary space needed: 2 * (p + M->n)*/ 188 mp_size_t 189 mpn_hgcd_matrix_adjust (const struct hgcd_matrix *M, 190 mp_size_t n, mp_ptr ap, mp_ptr bp, 191 mp_size_t p, mp_ptr tp) 192 { 193 /* M^-1 (a;b) = (r11, -r01; -r10, r00) (a ; b) 194 = (r11 a - r01 b; - r10 a + r00 b */ 195 196 mp_ptr t0 = tp; 197 mp_ptr t1 = tp + p + M->n; 198 mp_limb_t ah, bh; 199 mp_limb_t cy; 200 201 ASSERT (p + M->n < n); 202 203 /* First compute the two values depending on a, before overwriting a */ 204 205 if (M->n >= p) 206 { 207 mpn_mul (t0, M->p[1][1], M->n, ap, p); 208 mpn_mul (t1, M->p[1][0], M->n, ap, p); 209 } 210 else 211 { 212 mpn_mul (t0, ap, p, M->p[1][1], M->n); 213 mpn_mul (t1, ap, p, M->p[1][0], M->n); 214 } 215 216 /* Update a */ 217 MPN_COPY (ap, t0, p); 218 ah = mpn_add (ap + p, ap + p, n - p, t0 + p, M->n); 219 220 if (M->n >= p) 221 mpn_mul (t0, M->p[0][1], M->n, bp, p); 222 else 223 mpn_mul (t0, bp, p, M->p[0][1], M->n); 224 225 cy = mpn_sub (ap, ap, n, t0, p + M->n); 226 ASSERT (cy <= ah); 227 ah -= cy; 228 229 /* Update b */ 230 if (M->n >= p) 231 mpn_mul (t0, M->p[0][0], M->n, bp, p); 232 else 233 mpn_mul (t0, bp, p, M->p[0][0], M->n); 234 235 MPN_COPY (bp, t0, p); 236 bh = mpn_add (bp + p, bp + p, n - p, t0 + p, M->n); 237 cy = mpn_sub (bp, bp, n, t1, p + M->n); 238 ASSERT (cy <= bh); 239 bh -= cy; 240 241 if (ah > 0 || bh > 0) 242 { 243 ap[n] = ah; 244 bp[n] = bh; 245 n++; 246 } 247 else 248 { 249 /* The subtraction can reduce the size by at most one limb. */ 250 if (ap[n-1] == 0 && bp[n-1] == 0) 251 n--; 252 } 253 ASSERT (ap[n-1] > 0 || bp[n-1] > 0); 254 return n; 255 } 256