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