Line data Source code
1 : /* Copyright (C) 2021 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_mat
19 :
20 : GEN
21 228567 : zero_F3v(long m)
22 : {
23 228567 : long l = nbits2nlong(2*m);
24 228567 : GEN v = const_vecsmall(l+1, 0);
25 228566 : v[1] = m;
26 228566 : return v;
27 : }
28 :
29 : GEN
30 85223 : zero_F3m_copy(long m, long n)
31 : {
32 : long i;
33 85223 : GEN M = cgetg(n+1, t_MAT);
34 239651 : for (i = 1; i <= n; i++)
35 154428 : gel(M,i)= zero_F3v(m);
36 85223 : return M;
37 : }
38 : #define TRITS_IN_LONG (BITS_IN_LONG>>1)
39 : #define TRITS_MASK (ULONG_MAX/3UL)
40 : #define TWOPOTTRITS_IN_LONG (TWOPOTBITS_IN_LONG-1)
41 :
42 : ulong
43 5600160 : F3v_coeff(GEN x,long v)
44 : {
45 5600160 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
46 5600160 : long r = (v-1)&(TRITS_IN_LONG-1);
47 5600160 : ulong u=(ulong)x[2+pos];
48 5600160 : return (u>>(2*r))&3UL;
49 : }
50 :
51 : void
52 431069 : F3v_clear(GEN x, long v)
53 : {
54 431069 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
55 431069 : long r = (v-1)&(TRITS_IN_LONG-1);
56 431069 : ulong *u=(ulong*)&x[2+pos];
57 431069 : *u&=~(3UL<<(2*r));
58 431069 : }
59 :
60 : void
61 1070902 : F3v_set(GEN x, long v, ulong n)
62 : {
63 1070902 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
64 1070902 : long r = (v-1)&(TRITS_IN_LONG-1);
65 1070902 : ulong *u=(ulong*)&x[2+pos];
66 1070902 : *u&=~(3UL<<(2*r));
67 1070902 : *u|=(n<<(2*r));
68 1070902 : }
69 :
70 : INLINE void
71 1111109 : F3v_setneg(GEN x, long v)
72 : {
73 1111109 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
74 1111109 : long r = (v-1)&(TRITS_IN_LONG-1);
75 1111109 : ulong *u=(ulong*)&x[2+pos];
76 1111109 : if ((*u>>(2*r))&3UL)
77 360275 : *u^=(3UL<<(2*r));
78 1111109 : }
79 :
80 : INLINE void
81 1111109 : F3m_setneg(GEN x, long a, long b) { F3v_setneg(gel(x,b), a); }
82 :
83 : static ulong
84 2929110 : bitswap(ulong a)
85 : {
86 2929110 : const ulong m = TRITS_MASK;
87 2929110 : return ((a&m)<<1)|((a>>1)&m);
88 : }
89 :
90 : static ulong
91 545553 : F3_add(ulong a, ulong b)
92 : {
93 545553 : ulong c = a^b^bitswap(a&b);
94 545553 : return c&~bitswap(c);
95 : }
96 :
97 : static ulong
98 612668 : F3_sub(ulong a, ulong b)
99 : {
100 612668 : ulong bi = bitswap(b);
101 612668 : ulong c = a^bi^bitswap(a&bi);
102 612668 : return c&~bitswap(c);
103 : }
104 :
105 : /* Allow lg(y)<lg(x) */
106 : static void
107 285293 : F3v_add_inplace(GEN x, GEN y)
108 : {
109 285293 : long n = lg(y);
110 : long i;
111 830846 : for (i = 2; i < n; i++)
112 545553 : x[i] = F3_add(x[i], y[i]);
113 285293 : }
114 :
115 : /* Allow lg(y)<lg(x) */
116 : static void
117 345504 : F3v_sub_inplace(GEN x, GEN y)
118 : {
119 345504 : long n = lg(y);
120 : long i;
121 958172 : for (i = 2; i < n; i++)
122 612668 : x[i] = F3_sub(x[i], y[i]);
123 345504 : }
124 :
125 : static ulong
126 0 : F3_dotsquare(ulong x)
127 0 : { return hammingu(x)%3UL; }
128 :
129 : ulong
130 0 : F3v_dotsquare(GEN x)
131 : {
132 0 : long i, lx = lg(x);
133 : ulong c;
134 0 : if (lx <= 2) return 0;
135 0 : c = F3_dotsquare(uel(x,2));
136 0 : for (i=3; i<lx; i++) c = Fl_add(c, F3_dotsquare(uel(x,i)), 3);
137 0 : return c;
138 : }
139 :
140 : static ulong
141 0 : F3_dotproduct(ulong x, ulong y)
142 : {
143 0 : ulong x1 = 3*(x&TRITS_MASK);
144 0 : ulong x2 = 3*((x>>1)&TRITS_MASK);
145 0 : ulong z =(x1&y) | bitswap(x2&y);
146 0 : return (hammingu(z&TRITS_MASK)+2*hammingu(z&(TRITS_MASK<<1)))%3;
147 : }
148 :
149 : ulong
150 0 : F3v_dotproduct(GEN x, GEN y)
151 : {
152 0 : long i, lx = lg(x);
153 : ulong c;
154 0 : if (lx <= 2) return 0;
155 0 : c = F3_dotproduct(uel(x,2), uel(y,2));
156 0 : for (i=3; i<lx; i++) c = Fl_add(c, F3_dotproduct(uel(x,i), uel(y,i)), 3);
157 0 : return c;
158 : }
159 :
160 : GEN
161 0 : Flv_to_F3v(GEN x)
162 : {
163 0 : long l = lg(x)-1;
164 0 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
165 : long i,j,k;
166 0 : z[1] = l;
167 0 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
168 : {
169 0 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
170 0 : z[k] |= (uel(x,i)%3)<<j;
171 : }
172 0 : return z;
173 : }
174 :
175 : GEN
176 0 : Flm_to_F3m(GEN x) { pari_APPLY_same(Flv_to_F3v(gel(x,i))) }
177 :
178 : GEN
179 741363 : ZV_to_F3v(GEN x)
180 : {
181 741363 : long l = lg(x)-1;
182 741363 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
183 : long i,j,k;
184 741365 : z[1] = l;
185 11987860 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
186 : {
187 11246496 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
188 11246496 : z[k] |= umodiu(gel(x,i),3)<<j;
189 : }
190 741364 : return z;
191 : }
192 :
193 : GEN
194 900324 : ZM_to_F3m(GEN x) { pari_APPLY_same(ZV_to_F3v(gel(x,i))) }
195 :
196 : GEN
197 1939 : RgV_to_F3v(GEN x)
198 : {
199 1939 : long l = lg(x)-1;
200 1939 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
201 : long i,j,k;
202 1939 : z[1] = l;
203 44905 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
204 : {
205 42966 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
206 42966 : z[k] |= Rg_to_Fl(gel(x,i),3)<<j;
207 : }
208 1939 : return z;
209 : }
210 :
211 : GEN
212 2324 : RgM_to_F3m(GEN x) { pari_APPLY_same(RgV_to_F3v(gel(x,i))) }
213 :
214 : GEN
215 0 : F3v_to_Flv(GEN x)
216 : {
217 0 : long l = x[1]+1, i, j, k;
218 0 : GEN z = cgetg(l, t_VECSMALL);
219 0 : for (i=2,k=1; i<lg(x); i++)
220 0 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
221 0 : z[k] = (uel(x,i)>>j)&3UL;
222 0 : return z;
223 : }
224 : GEN
225 227621 : F3c_to_ZC(GEN x)
226 : {
227 227621 : long l = x[1]+1, i, j, k;
228 227621 : GEN z = cgetg(l, t_COL);
229 457023 : for (i=2,k=1; i<lg(x); i++)
230 1448743 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
231 1219341 : switch((uel(x,i)>>j)&3UL)
232 : {
233 776243 : case 0: gel(z,k) = gen_0; break;
234 334313 : case 1: gel(z,k) = gen_1; break;
235 108785 : default:gel(z,k) = gen_2; break;
236 : }
237 227621 : return z;
238 : }
239 : GEN
240 945 : F3c_to_mod(GEN x)
241 : {
242 945 : long l = x[1]+1, i, j, k;
243 945 : GEN z = cgetg(l, t_COL), N = utoipos(3);
244 945 : GEN _0 = mkintmod(gen_0, N);
245 945 : GEN _1 = mkintmod(gen_1, N);
246 945 : GEN _2 = mkintmod(gen_2, N);
247 2978 : for (i=2,k=1; i<lg(x); i++)
248 38398 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
249 36365 : switch((uel(x,i)>>j)&3UL)
250 : {
251 34006 : case 0: gel(z,k) = _0; break;
252 1232 : case 1: gel(z,k) = _1; break;
253 1127 : default: gel(z,k)= _2; break;
254 : }
255 945 : return z;
256 : }
257 :
258 : GEN
259 238951 : F3m_to_ZM(GEN x) { pari_APPLY_same(F3c_to_ZC(gel(x,i))) }
260 : GEN
261 1176 : F3m_to_mod(GEN x) { pari_APPLY_same(F3c_to_mod(gel(x,i))) }
262 : GEN
263 0 : F3m_to_Flm(GEN x) { pari_APPLY_same(F3v_to_Flv(gel(x,i))) }
264 :
265 : /* in place, destroy x */
266 : GEN
267 159039 : F3m_ker_sp(GEN x, long deplin)
268 : {
269 : GEN y, c, d;
270 : long i, j, k, r, m, n;
271 :
272 159039 : n = lg(x)-1;
273 159039 : m = mael(x,1,1); r=0;
274 :
275 159039 : d = cgetg(n+1, t_VECSMALL);
276 159039 : c = const_F2v(m);
277 744536 : for (k=1; k<=n; k++)
278 : {
279 659313 : GEN xk = gel(x,k);
280 4659787 : for (j=1; j<=m; j++)
281 4431542 : if (F2v_coeff(c,j) && F3v_coeff(xk,j)) break;
282 659314 : if (j>m)
283 : {
284 228245 : if (deplin) {
285 73817 : GEN v = zero_F3v(n);
286 340550 : for (i=1; i<k; i++) F3v_set(v, i, F3v_coeff(xk, d[i]));
287 73817 : F3v_set(v, k, 1); return v;
288 : }
289 154428 : r++; d[k] = 0;
290 : }
291 : else
292 : {
293 431069 : ulong xkj = F3v_coeff(xk,j);
294 431069 : F3v_clear(xk, j);
295 431069 : F2v_clear(c,j); d[k] = j;
296 2077493 : for (i=k+1; i<=n; i++)
297 : {
298 1646426 : GEN xi = gel(x,i);
299 1646426 : ulong u = F3v_coeff(xi,j);
300 1646425 : if (u)
301 : {
302 630034 : if (u==xkj) F3v_sub_inplace(xi, xk);
303 285111 : else F3v_add_inplace(xi, xk);
304 : }
305 : }
306 431067 : F3v_set(xk, j, 2);
307 431067 : if (xkj==1)
308 1407220 : for (i=k+1; i<=n; i++) F3m_setneg(x,j,i);
309 : }
310 : }
311 85223 : if (deplin) return NULL;
312 85223 : y = zero_F3m_copy(n,r);
313 239651 : for (j=k=1; j<=r; j++,k++)
314 : {
315 154428 : GEN C = gel(y,j);
316 223936 : while (d[k]) k++;
317 529858 : for (i=1; i<k; i++)
318 375430 : if (d[i]) F3v_set(C,i,F3m_coeff(x,d[i],k));
319 154428 : F3v_set(C, k, 1);
320 : }
321 85223 : return y;
322 : }
323 :
324 : GEN
325 0 : F3m_ker(GEN x) { return F3m_ker_sp(F3m_copy(x), 0); }
326 :
327 : long
328 0 : F3m_rank(GEN x)
329 : {
330 0 : pari_sp av = avma;
331 0 : return gc_long (av, lg(x) - lg(F3m_ker(x)));
332 : }
333 :
334 : INLINE GEN
335 322 : F3m_F3c_mul_i(GEN x, GEN y, long lx, long l)
336 : {
337 : long j;
338 322 : GEN z = zero_F3v(l);
339 :
340 2114 : for (j=1; j<lx; j++)
341 : {
342 1792 : ulong c = F3v_coeff(y,j);
343 1792 : if (!c) continue;
344 763 : if (c==1)
345 182 : F3v_add_inplace(z,gel(x,j));
346 : else
347 581 : F3v_sub_inplace(z,gel(x,j));
348 : }
349 322 : return z;
350 : }
351 :
352 : GEN
353 154 : F3m_mul(GEN x, GEN y)
354 : {
355 154 : long i,j,l,lx=lg(x), ly=lg(y);
356 : GEN z;
357 154 : if (ly==1) return cgetg(1,t_MAT);
358 154 : z = cgetg(ly,t_MAT);
359 154 : if (lx==1)
360 : {
361 0 : for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
362 0 : return z;
363 : }
364 154 : l = coeff(x,1,1);
365 476 : for (j=1; j<ly; j++) gel(z,j) = F3m_F3c_mul_i(x, gel(y,j), lx, l);
366 154 : return z;
367 : }
368 :
369 : GEN
370 0 : F3m_row(GEN x, long j)
371 : {
372 0 : long i, l = lg(x);
373 0 : GEN V = zero_F3v(l-1);
374 0 : for(i = 1; i < l; i++) F3v_set(V, i, F3m_coeff(x,j,i));
375 0 : return V;
376 : }
377 :
378 : GEN
379 0 : F3m_transpose(GEN x)
380 : {
381 : long i, l;
382 : GEN y;
383 0 : if (lg(x) == 1) return cgetg(1,t_MAT);
384 0 : l = coeff(x,1,1) + 1; y = cgetg(l, t_MAT);
385 0 : for (i = 1; i < l; i++) gel(y,i) = F3m_row(x,i);
386 0 : return y;
387 : }
|
