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 194251 : zero_F3v(long m)
22 : {
23 194251 : long l = nbits2nlong(2*m);
24 194251 : GEN v = const_vecsmall(l+1, 0);
25 194251 : v[1] = m;
26 194251 : return v;
27 : }
28 :
29 : GEN
30 72338 : zero_F3m_copy(long m, long n)
31 : {
32 : long i;
33 72338 : GEN M = cgetg(n+1, t_MAT);
34 203191 : for (i = 1; i <= n; i++)
35 130853 : gel(M,i)= zero_F3v(m);
36 72338 : 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 4764208 : F3v_coeff(GEN x,long v)
44 : {
45 4764208 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
46 4764208 : long r = (v-1)&(TRITS_IN_LONG-1);
47 4764208 : ulong u=(ulong)x[2+pos];
48 4764208 : return (u>>(2*r))&3UL;
49 : }
50 :
51 : void
52 366382 : F3v_clear(GEN x, long v)
53 : {
54 366382 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
55 366382 : long r = (v-1)&(TRITS_IN_LONG-1);
56 366382 : ulong *u=(ulong*)&x[2+pos];
57 366382 : *u&=~(3UL<<(2*r));
58 366382 : }
59 :
60 : void
61 909622 : F3v_set(GEN x, long v, ulong n)
62 : {
63 909622 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
64 909622 : long r = (v-1)&(TRITS_IN_LONG-1);
65 909622 : ulong *u=(ulong*)&x[2+pos];
66 909622 : *u&=~(3UL<<(2*r));
67 909622 : *u|=(n<<(2*r));
68 909622 : }
69 :
70 : INLINE void
71 945678 : F3v_setneg(GEN x, long v)
72 : {
73 945678 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
74 945678 : long r = (v-1)&(TRITS_IN_LONG-1);
75 945678 : ulong *u=(ulong*)&x[2+pos];
76 945678 : if ((*u>>(2*r))&3UL)
77 304921 : *u^=(3UL<<(2*r));
78 945678 : }
79 :
80 : INLINE void
81 945678 : F3m_setneg(GEN x, long a, long b) { F3v_setneg(gel(x,b), a); }
82 :
83 : static ulong
84 2526085 : bitswap(ulong a)
85 : {
86 2526085 : const ulong m = TRITS_MASK;
87 2526085 : return ((a&m)<<1)|((a>>1)&m);
88 : }
89 :
90 : static ulong
91 470183 : F3_add(ulong a, ulong b)
92 : {
93 470183 : ulong c = a^b^bitswap(a&b);
94 470183 : return c&~bitswap(c);
95 : }
96 :
97 : static ulong
98 528573 : F3_sub(ulong a, ulong b)
99 : {
100 528573 : ulong bi = bitswap(b);
101 528573 : ulong c = a^bi^bitswap(a&bi);
102 528573 : return c&~bitswap(c);
103 : }
104 :
105 : /* Allow lg(y)<lg(x) */
106 : static void
107 241155 : F3v_add_inplace(GEN x, GEN y)
108 : {
109 241155 : long n = lg(y);
110 : long i;
111 711338 : for (i = 2; i < n; i++)
112 470183 : x[i] = F3_add(x[i], y[i]);
113 241155 : }
114 :
115 : /* Allow lg(y)<lg(x) */
116 : static void
117 292865 : F3v_sub_inplace(GEN x, GEN y)
118 : {
119 292865 : long n = lg(y);
120 : long i;
121 821438 : for (i = 2; i < n; i++)
122 528573 : x[i] = F3_sub(x[i], y[i]);
123 292865 : }
124 :
125 : GEN
126 0 : Flv_to_F3v(GEN x)
127 : {
128 0 : long l = lg(x)-1;
129 0 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
130 : long i,j,k;
131 0 : z[1] = l;
132 0 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
133 : {
134 0 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
135 0 : z[k] |= (uel(x,i)%3)<<j;
136 : }
137 0 : return z;
138 : }
139 :
140 : GEN
141 0 : Flm_to_F3m(GEN x) { pari_APPLY_same(Flv_to_F3v(gel(x,i))) }
142 :
143 : GEN
144 630470 : ZV_to_F3v(GEN x)
145 : {
146 630470 : long l = lg(x)-1;
147 630470 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
148 : long i,j,k;
149 630470 : z[1] = l;
150 10229466 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
151 : {
152 9598996 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
153 9598996 : z[k] |= umodiu(gel(x,i),3)<<j;
154 : }
155 630470 : return z;
156 : }
157 :
158 : GEN
159 765864 : ZM_to_F3m(GEN x) { pari_APPLY_same(ZV_to_F3v(gel(x,i))) }
160 :
161 : GEN
162 1662 : RgV_to_F3v(GEN x)
163 : {
164 1662 : long l = lg(x)-1;
165 1662 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
166 : long i,j,k;
167 1662 : z[1] = l;
168 38490 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
169 : {
170 36828 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
171 36828 : z[k] |= Rg_to_Fl(gel(x,i),3)<<j;
172 : }
173 1662 : return z;
174 : }
175 :
176 : GEN
177 1992 : RgM_to_F3m(GEN x) { pari_APPLY_same(RgV_to_F3v(gel(x,i))) }
178 :
179 : GEN
180 0 : F3v_to_Flv(GEN x)
181 : {
182 0 : long l = x[1]+1, i, j, k;
183 0 : GEN z = cgetg(l, t_VECSMALL);
184 0 : for (i=2,k=1; i<lg(x); i++)
185 0 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
186 0 : z[k] = (uel(x,i)>>j)&3UL;
187 0 : return z;
188 : }
189 : GEN
190 193441 : F3c_to_ZC(GEN x)
191 : {
192 193441 : long l = x[1]+1, i, j, k;
193 193441 : GEN z = cgetg(l, t_COL);
194 388490 : for (i=2,k=1; i<lg(x); i++)
195 1227061 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
196 1032012 : switch((uel(x,i)>>j)&3UL)
197 : {
198 655624 : case 0: gel(z,k) = gen_0; break;
199 283905 : case 1: gel(z,k) = gen_1; break;
200 92483 : default:gel(z,k) = gen_2; break;
201 : }
202 193441 : return z;
203 : }
204 : GEN
205 810 : F3c_to_mod(GEN x)
206 : {
207 810 : long l = x[1]+1, i, j, k;
208 810 : GEN z = cgetg(l, t_COL), N = utoipos(3);
209 810 : GEN _0 = mkintmod(gen_0, N);
210 810 : GEN _1 = mkintmod(gen_1, N);
211 810 : GEN _2 = mkintmod(gen_2, N);
212 2572 : for (i=2,k=1; i<lg(x); i++)
213 32932 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
214 31170 : switch((uel(x,i)>>j)&3UL)
215 : {
216 29148 : case 0: gel(z,k) = _0; break;
217 1056 : case 1: gel(z,k) = _1; break;
218 966 : default: gel(z,k)= _2; break;
219 : }
220 810 : return z;
221 : }
222 :
223 : GEN
224 202591 : F3m_to_ZM(GEN x) { pari_APPLY_same(F3c_to_ZC(gel(x,i))) }
225 : GEN
226 1008 : F3m_to_mod(GEN x) { pari_APPLY_same(F3c_to_mod(gel(x,i))) }
227 : GEN
228 0 : F3m_to_Flm(GEN x) { pari_APPLY_same(F3v_to_Flv(gel(x,i))) }
229 :
230 : /* in place, destroy x */
231 : GEN
232 135460 : F3m_ker_sp(GEN x, long deplin)
233 : {
234 : GEN y, c, d;
235 : long i, j, k, r, m, n;
236 :
237 135460 : n = lg(x)-1;
238 135460 : m = mael(x,1,1); r=0;
239 :
240 135460 : d = cgetg(n+1, t_VECSMALL);
241 135460 : c = const_F2v(m);
242 632695 : for (k=1; k<=n; k++)
243 : {
244 560357 : GEN xk = gel(x,k);
245 3966282 : for (j=1; j<=m; j++)
246 3772307 : if (F2v_coeff(c,j) && F3m_coeff(x,j,k)) break;
247 560357 : if (j>m)
248 : {
249 193975 : if (deplin) {
250 63122 : GEN v = zero_F3v(n);
251 291153 : for (i=1; i<k; i++) F3v_set(v, i, F3v_coeff(xk, d[i]));
252 63122 : F3v_set(v, k, 1); return v;
253 : }
254 130853 : r++; d[k] = 0;
255 : }
256 : else
257 : {
258 366382 : ulong xkj = F3v_coeff(xk,j);
259 366382 : F3v_clear(xk, j);
260 366382 : F2v_clear(c,j); d[k] = j;
261 1765183 : for (i=k+1; i<=n; i++)
262 : {
263 1398801 : GEN xi = gel(x,i);
264 1398801 : ulong u = F3v_coeff(xi,j);
265 1398801 : if (u)
266 : {
267 533366 : if (u==xkj) F3v_sub_inplace(xi, xk);
268 240999 : else F3v_add_inplace(xi, xk);
269 : }
270 : }
271 366382 : F3v_set(xk, j, 2);
272 366382 : if (xkj==1)
273 1198017 : for (i=k+1; i<=n; i++) F3m_setneg(x,j,i);
274 : }
275 : }
276 72338 : if (deplin) return NULL;
277 72338 : y = zero_F3m_copy(n,r);
278 203191 : for (j=k=1; j<=r; j++,k++)
279 : {
280 130853 : GEN C = gel(y,j);
281 188711 : while (d[k]) k++;
282 446506 : for (i=1; i<k; i++)
283 315653 : if (d[i]) F3v_set(C,i,F3m_coeff(x,d[i],k));
284 130853 : F3v_set(C, k, 1);
285 : }
286 72338 : return y;
287 : }
288 :
289 : GEN
290 0 : F3m_ker(GEN x) { return F3m_ker_sp(F3m_copy(x), 0); }
291 :
292 : INLINE GEN
293 276 : F3m_F3c_mul_i(GEN x, GEN y, long lx, long l)
294 : {
295 : long j;
296 276 : GEN z = zero_F3v(l);
297 :
298 1812 : for (j=1; j<lx; j++)
299 : {
300 1536 : ulong c = F3v_coeff(y,j);
301 1536 : if (!c) continue;
302 654 : if (c==1)
303 156 : F3v_add_inplace(z,gel(x,j));
304 : else
305 498 : F3v_sub_inplace(z,gel(x,j));
306 : }
307 276 : return z;
308 : }
309 :
310 : GEN
311 132 : F3m_mul(GEN x, GEN y)
312 : {
313 132 : long i,j,l,lx=lg(x), ly=lg(y);
314 : GEN z;
315 132 : if (ly==1) return cgetg(1,t_MAT);
316 132 : z = cgetg(ly,t_MAT);
317 132 : if (lx==1)
318 : {
319 0 : for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
320 0 : return z;
321 : }
322 132 : l = coeff(x,1,1);
323 408 : for (j=1; j<ly; j++) gel(z,j) = F3m_F3c_mul_i(x, gel(y,j), lx, l);
324 132 : return z;
325 : }
326 :
327 : GEN
328 0 : F3m_row(GEN x, long j)
329 : {
330 0 : long i, l = lg(x);
331 0 : GEN V = zero_F3v(l-1);
332 0 : for(i = 1; i < l; i++) F3v_set(V, i, F3m_coeff(x,j,i));
333 0 : return V;
334 : }
335 :
336 : GEN
337 0 : F3m_transpose(GEN x)
338 : {
339 : long i, l;
340 : GEN y;
341 0 : if (lg(x) == 1) return cgetg(1,t_MAT);
342 0 : l = coeff(x,1,1) + 1; y = cgetg(l, t_MAT);
343 0 : for (i = 1; i < l; i++) gel(y,i) = F3m_row(x,i);
344 0 : return y;
345 : }
|
