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 223030 : zero_F3v(long m)
22 : {
23 223030 : long l = nbits2nlong(2*m);
24 223030 : GEN v = const_vecsmall(l+1, 0);
25 223030 : v[1] = m;
26 223030 : return v;
27 : }
28 :
29 : GEN
30 81373 : zero_F3m_copy(long m, long n)
31 : {
32 : long i;
33 81373 : GEN M = cgetg(n+1, t_MAT);
34 231515 : for (i = 1; i <= n; i++)
35 150142 : gel(M,i)= zero_F3v(m);
36 81373 : 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 4060188 : F3v_coeff(GEN x,long v)
44 : {
45 4060188 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
46 4060188 : long r = (v-1)&(TRITS_IN_LONG-1);
47 4060188 : ulong u=(ulong)x[2+pos];
48 4060188 : return (u>>(2*r))&3UL;
49 : }
50 :
51 : void
52 417044 : F3v_clear(GEN x, long v)
53 : {
54 417044 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
55 417044 : long r = (v-1)&(TRITS_IN_LONG-1);
56 417044 : ulong *u=(ulong*)&x[2+pos];
57 417044 : *u&=~(3UL<<(2*r));
58 417044 : }
59 :
60 : void
61 1034658 : F3v_set(GEN x, long v, ulong n)
62 : {
63 1034658 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
64 1034658 : long r = (v-1)&(TRITS_IN_LONG-1);
65 1034658 : ulong *u=(ulong*)&x[2+pos];
66 1034658 : *u&=~(3UL<<(2*r));
67 1034658 : *u|=(n<<(2*r));
68 1034658 : }
69 :
70 : INLINE void
71 1067119 : F3v_setneg(GEN x, long v)
72 : {
73 1067119 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
74 1067119 : long r = (v-1)&(TRITS_IN_LONG-1);
75 1067119 : ulong *u=(ulong*)&x[2+pos];
76 1067119 : if ((*u>>(2*r))&3UL)
77 352793 : *u^=(3UL<<(2*r));
78 1067119 : }
79 :
80 : INLINE void
81 1067119 : F3m_setneg(GEN x, long a, long b) { F3v_setneg(gel(x,b), a); }
82 :
83 : static ulong
84 2415990 : bitswap(ulong a)
85 : {
86 2415990 : const ulong m = TRITS_MASK;
87 2415990 : return ((a&m)<<1)|((a>>1)&m);
88 : }
89 :
90 : static ulong
91 448797 : F3_add(ulong a, ulong b)
92 : {
93 448797 : ulong c = a^b^bitswap(a&b);
94 448797 : return c&~bitswap(c);
95 : }
96 :
97 : static ulong
98 506132 : F3_sub(ulong a, ulong b)
99 : {
100 506132 : ulong bi = bitswap(b);
101 506132 : ulong c = a^bi^bitswap(a&bi);
102 506132 : return c&~bitswap(c);
103 : }
104 :
105 : /* Allow lg(y)<lg(x) */
106 : static void
107 278215 : F3v_add_inplace(GEN x, GEN y)
108 : {
109 278215 : long n = lg(y);
110 : long i;
111 727012 : for (i = 2; i < n; i++)
112 448797 : x[i] = F3_add(x[i], y[i]);
113 278215 : }
114 :
115 : /* Allow lg(y)<lg(x) */
116 : static void
117 337154 : F3v_sub_inplace(GEN x, GEN y)
118 : {
119 337154 : long n = lg(y);
120 : long i;
121 843286 : for (i = 2; i < n; i++)
122 506132 : x[i] = F3_sub(x[i], y[i]);
123 337154 : }
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 721307 : ZV_to_F3v(GEN x)
145 : {
146 721307 : long l = lg(x)-1;
147 721307 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
148 : long i,j,k;
149 721309 : z[1] = l;
150 7431970 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
151 : {
152 6710661 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
153 6710661 : z[k] |= umodiu(gel(x,i),3)<<j;
154 : }
155 721309 : return z;
156 : }
157 :
158 : GEN
159 875533 : ZM_to_F3m(GEN x) { pari_APPLY_same(ZV_to_F3v(gel(x,i))) }
160 :
161 : GEN
162 546 : RgV_to_F3v(GEN x)
163 : {
164 546 : long l = lg(x)-1;
165 546 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
166 : long i,j,k;
167 546 : z[1] = l;
168 34958 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
169 : {
170 34412 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
171 34412 : z[k] |= Rg_to_Fl(gel(x,i),3)<<j;
172 : }
173 546 : return z;
174 : }
175 :
176 : GEN
177 581 : 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 222526 : F3c_to_ZC(GEN x)
191 : {
192 222526 : long l = x[1]+1, i, j, k;
193 222526 : GEN z = cgetg(l, t_COL);
194 446347 : for (i=2,k=1; i<lg(x); i++)
195 1387197 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
196 1163376 : switch((uel(x,i)>>j)&3UL)
197 : {
198 727740 : case 0: gel(z,k) = gen_0; break;
199 328598 : case 1: gel(z,k) = gen_1; break;
200 107038 : default:gel(z,k) = gen_2; break;
201 : }
202 222526 : return z;
203 : }
204 : GEN
205 504 : F3c_to_mod(GEN x)
206 : {
207 504 : long l = x[1]+1, i, j, k;
208 504 : GEN z = cgetg(l, t_COL), N = utoipos(3);
209 504 : GEN _0 = mkintmod(gen_0, N);
210 504 : GEN _1 = mkintmod(gen_1, N);
211 504 : GEN _2 = mkintmod(gen_2, N);
212 2096 : for (i=2,k=1; i<lg(x); i++)
213 34968 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
214 33376 : switch((uel(x,i)>>j)&3UL)
215 : {
216 32228 : case 0: gel(z,k) = _0; break;
217 819 : case 1: gel(z,k) = _1; break;
218 329 : default: gel(z,k)= _2; break;
219 : }
220 504 : return z;
221 : }
222 :
223 : GEN
224 231032 : F3m_to_ZM(GEN x) { pari_APPLY_same(F3c_to_ZC(gel(x,i))) }
225 : GEN
226 525 : 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 154233 : F3m_ker_sp(GEN x, long deplin)
233 : {
234 : GEN y, c, d;
235 : long i, j, k, r, m, n;
236 :
237 154233 : n = lg(x)-1;
238 154233 : m = mael(x,1,1); r=0;
239 :
240 154233 : d = cgetg(n+1, t_VECSMALL);
241 154231 : c = const_F2v(m);
242 721417 : for (k=1; k<=n; k++)
243 : {
244 640044 : GEN xk = gel(x,k);
245 3161876 : for (j=1; j<=m; j++)
246 2938870 : if (F2v_coeff(c,j) && F3m_coeff(x,j,k)) break;
247 640050 : if (j>m)
248 : {
249 223002 : if (deplin) {
250 72860 : GEN v = zero_F3v(n);
251 337715 : for (i=1; i<k; i++) F3v_set(v, i, F3v_coeff(xk, d[i]));
252 72860 : F3v_set(v, k, 1); return v;
253 : }
254 150142 : r++; d[k] = 0;
255 : }
256 : else
257 : {
258 417048 : ulong xkj = F3v_coeff(xk,j);
259 417044 : F3v_clear(xk, j);
260 417044 : F2v_clear(c,j); d[k] = j;
261 1984090 : for (i=k+1; i<=n; i++)
262 : {
263 1567046 : GEN xi = gel(x,i);
264 1567046 : ulong u = F3v_coeff(xi,j);
265 1567046 : if (u)
266 : {
267 615320 : if (u==xkj) F3v_sub_inplace(xi, xk);
268 278201 : else F3v_add_inplace(xi, xk);
269 : }
270 : }
271 417044 : F3v_set(xk, j, 2);
272 417044 : if (xkj==1)
273 1354163 : for (i=k+1; i<=n; i++) F3m_setneg(x,j,i);
274 : }
275 : }
276 81373 : if (deplin) return NULL;
277 81373 : y = zero_F3m_copy(n,r);
278 231515 : for (j=k=1; j<=r; j++,k++)
279 : {
280 150142 : GEN C = gel(y,j);
281 216948 : while (d[k]) k++;
282 495968 : for (i=1; i<k; i++)
283 345826 : if (d[i]) F3v_set(C,i,F3m_coeff(x,d[i],k));
284 150142 : F3v_set(C, k, 1);
285 : }
286 81373 : 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 28 : F3m_F3c_mul_i(GEN x, GEN y, long lx, long l)
294 : {
295 : long j;
296 28 : GEN z = zero_F3v(l);
297 :
298 84 : for (j=1; j<lx; j++)
299 : {
300 56 : ulong c = F3v_coeff(y,j);
301 56 : if (!c) continue;
302 49 : if (c==1)
303 14 : F3v_add_inplace(z,gel(x,j));
304 : else
305 35 : F3v_sub_inplace(z,gel(x,j));
306 : }
307 28 : return z;
308 : }
309 :
310 : GEN
311 14 : F3m_mul(GEN x, GEN y)
312 : {
313 14 : long i,j,l,lx=lg(x), ly=lg(y);
314 : GEN z;
315 14 : if (ly==1) return cgetg(1,t_MAT);
316 14 : z = cgetg(ly,t_MAT);
317 14 : if (lx==1)
318 : {
319 0 : for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
320 0 : return z;
321 : }
322 14 : l = coeff(x,1,1);
323 42 : for (j=1; j<ly; j++) gel(z,j) = F3m_F3c_mul_i(x, gel(y,j), lx, l);
324 14 : 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 : }
|
