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 227195 : zero_F3v(long m)
22 : {
23 227195 : long l = nbits2nlong(2*m);
24 227195 : GEN v = const_vecsmall(l+1, 0);
25 227195 : v[1] = m;
26 227195 : return v;
27 : }
28 :
29 : GEN
30 83496 : zero_F3m_copy(long m, long n)
31 : {
32 : long i;
33 83496 : GEN M = cgetg(n+1, t_MAT);
34 237358 : for (i = 1; i <= n; i++)
35 153862 : gel(M,i)= zero_F3v(m);
36 83496 : 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 5518242 : F3v_coeff(GEN x,long v)
44 : {
45 5518242 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
46 5518242 : long r = (v-1)&(TRITS_IN_LONG-1);
47 5518242 : ulong u=(ulong)x[2+pos];
48 5518242 : return (u>>(2*r))&3UL;
49 : }
50 :
51 : void
52 423307 : F3v_clear(GEN x, long v)
53 : {
54 423307 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
55 423307 : long r = (v-1)&(TRITS_IN_LONG-1);
56 423307 : ulong *u=(ulong*)&x[2+pos];
57 423307 : *u&=~(3UL<<(2*r));
58 423307 : }
59 :
60 : void
61 1059917 : F3v_set(GEN x, long v, ulong n)
62 : {
63 1059917 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
64 1059917 : long r = (v-1)&(TRITS_IN_LONG-1);
65 1059917 : ulong *u=(ulong*)&x[2+pos];
66 1059917 : *u&=~(3UL<<(2*r));
67 1059917 : *u|=(n<<(2*r));
68 1059917 : }
69 :
70 : INLINE void
71 1095668 : F3v_setneg(GEN x, long v)
72 : {
73 1095668 : long pos = (v-1)>>TWOPOTTRITS_IN_LONG;
74 1095668 : long r = (v-1)&(TRITS_IN_LONG-1);
75 1095668 : ulong *u=(ulong*)&x[2+pos];
76 1095668 : if ((*u>>(2*r))&3UL)
77 356505 : *u^=(3UL<<(2*r));
78 1095668 : }
79 :
80 : INLINE void
81 1095668 : F3m_setneg(GEN x, long a, long b) { F3v_setneg(gel(x,b), a); }
82 :
83 : static ulong
84 3065229 : bitswap(ulong a)
85 : {
86 3065229 : const ulong m = TRITS_MASK;
87 3065229 : return ((a&m)<<1)|((a>>1)&m);
88 : }
89 :
90 : static ulong
91 572155 : F3_add(ulong a, ulong b)
92 : {
93 572155 : ulong c = a^b^bitswap(a&b);
94 572155 : return c&~bitswap(c);
95 : }
96 :
97 : static ulong
98 640307 : F3_sub(ulong a, ulong b)
99 : {
100 640307 : ulong bi = bitswap(b);
101 640307 : ulong c = a^bi^bitswap(a&bi);
102 640307 : return c&~bitswap(c);
103 : }
104 :
105 : /* Allow lg(y)<lg(x) */
106 : static void
107 282476 : F3v_add_inplace(GEN x, GEN y)
108 : {
109 282476 : long n = lg(y);
110 : long i;
111 854631 : for (i = 2; i < n; i++)
112 572155 : x[i] = F3_add(x[i], y[i]);
113 282476 : }
114 :
115 : /* Allow lg(y)<lg(x) */
116 : static void
117 341616 : F3v_sub_inplace(GEN x, GEN y)
118 : {
119 341616 : long n = lg(y);
120 : long i;
121 981923 : for (i = 2; i < n; i++)
122 640307 : x[i] = F3_sub(x[i], y[i]);
123 341616 : }
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 731031 : ZV_to_F3v(GEN x)
145 : {
146 731031 : long l = lg(x)-1;
147 731031 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
148 : long i,j,k;
149 731028 : z[1] = l;
150 11629686 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
151 : {
152 10898656 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
153 10898656 : z[k] |= umodiu(gel(x,i),3)<<j;
154 : }
155 731030 : return z;
156 : }
157 :
158 : GEN
159 887460 : ZM_to_F3m(GEN x) { pari_APPLY_same(ZV_to_F3v(gel(x,i))) }
160 :
161 : GEN
162 1939 : RgV_to_F3v(GEN x)
163 : {
164 1939 : long l = lg(x)-1;
165 1939 : GEN z = cgetg(nbits2lg(2*l), t_VECSMALL);
166 : long i,j,k;
167 1939 : z[1] = l;
168 44905 : for(i=1,k=1,j=BITS_IN_LONG; i<=l; i++,j+=2)
169 : {
170 42966 : if (j==BITS_IN_LONG) { j=0; z[++k]=0; }
171 42966 : z[k] |= Rg_to_Fl(gel(x,i),3)<<j;
172 : }
173 1939 : return z;
174 : }
175 :
176 : GEN
177 2324 : 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 226251 : F3c_to_ZC(GEN x)
191 : {
192 226251 : long l = x[1]+1, i, j, k;
193 226251 : GEN z = cgetg(l, t_COL);
194 454223 : for (i=2,k=1; i<lg(x); i++)
195 1437972 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
196 1210000 : switch((uel(x,i)>>j)&3UL)
197 : {
198 770125 : case 0: gel(z,k) = gen_0; break;
199 332071 : case 1: gel(z,k) = gen_1; break;
200 107804 : default:gel(z,k) = gen_2; break;
201 : }
202 226251 : return z;
203 : }
204 : GEN
205 945 : F3c_to_mod(GEN x)
206 : {
207 945 : long l = x[1]+1, i, j, k;
208 945 : GEN z = cgetg(l, t_COL), N = utoipos(3);
209 945 : GEN _0 = mkintmod(gen_0, N);
210 945 : GEN _1 = mkintmod(gen_1, N);
211 945 : GEN _2 = mkintmod(gen_2, N);
212 2978 : for (i=2,k=1; i<lg(x); i++)
213 38398 : for (j=0; j<BITS_IN_LONG && k<l; j+=2,k++)
214 36365 : switch((uel(x,i)>>j)&3UL)
215 : {
216 34006 : case 0: gel(z,k) = _0; break;
217 1232 : case 1: gel(z,k) = _1; break;
218 1127 : default: gel(z,k)= _2; break;
219 : }
220 945 : return z;
221 : }
222 :
223 : GEN
224 236658 : F3m_to_ZM(GEN x) { pari_APPLY_same(F3c_to_ZC(gel(x,i))) }
225 : GEN
226 1176 : 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 156508 : F3m_ker_sp(GEN x, long deplin)
233 : {
234 : GEN y, c, d;
235 : long i, j, k, r, m, n;
236 :
237 156508 : n = lg(x)-1;
238 156508 : m = mael(x,1,1); r=0;
239 :
240 156508 : d = cgetg(n+1, t_VECSMALL);
241 156508 : c = const_F2v(m);
242 733675 : for (k=1; k<=n; k++)
243 : {
244 650179 : GEN xk = gel(x,k);
245 4590086 : for (j=1; j<=m; j++)
246 4363213 : if (F2v_coeff(c,j) && F3m_coeff(x,j,k)) break;
247 650180 : if (j>m)
248 : {
249 226874 : if (deplin) {
250 73012 : GEN v = zero_F3v(n);
251 337502 : for (i=1; i<k; i++) F3v_set(v, i, F3v_coeff(xk, d[i]));
252 73012 : F3v_set(v, k, 1); return v;
253 : }
254 153862 : r++; d[k] = 0;
255 : }
256 : else
257 : {
258 423306 : ulong xkj = F3v_coeff(xk,j);
259 423307 : F3v_clear(xk, j);
260 423307 : F2v_clear(c,j); d[k] = j;
261 2045313 : for (i=k+1; i<=n; i++)
262 : {
263 1622006 : GEN xi = gel(x,i);
264 1622006 : ulong u = F3v_coeff(xi,j);
265 1622007 : if (u)
266 : {
267 623329 : if (u==xkj) F3v_sub_inplace(xi, xk);
268 282294 : else F3v_add_inplace(xi, xk);
269 : }
270 : }
271 423307 : F3v_set(xk, j, 2);
272 423307 : if (xkj==1)
273 1385797 : for (i=k+1; i<=n; i++) F3m_setneg(x,j,i);
274 : }
275 : }
276 83496 : if (deplin) return NULL;
277 83496 : y = zero_F3m_copy(n,r);
278 237358 : for (j=k=1; j<=r; j++,k++)
279 : {
280 153862 : GEN C = gel(y,j);
281 223542 : while (d[k]) k++;
282 528177 : for (i=1; i<k; i++)
283 374315 : if (d[i]) F3v_set(C,i,F3m_coeff(x,d[i],k));
284 153862 : F3v_set(C, k, 1);
285 : }
286 83496 : 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 322 : F3m_F3c_mul_i(GEN x, GEN y, long lx, long l)
294 : {
295 : long j;
296 322 : GEN z = zero_F3v(l);
297 :
298 2114 : for (j=1; j<lx; j++)
299 : {
300 1792 : ulong c = F3v_coeff(y,j);
301 1792 : if (!c) continue;
302 763 : if (c==1)
303 182 : F3v_add_inplace(z,gel(x,j));
304 : else
305 581 : F3v_sub_inplace(z,gel(x,j));
306 : }
307 322 : return z;
308 : }
309 :
310 : GEN
311 154 : F3m_mul(GEN x, GEN y)
312 : {
313 154 : long i,j,l,lx=lg(x), ly=lg(y);
314 : GEN z;
315 154 : if (ly==1) return cgetg(1,t_MAT);
316 154 : z = cgetg(ly,t_MAT);
317 154 : if (lx==1)
318 : {
319 0 : for (i=1; i<ly; i++) gel(z,i) = mkvecsmall(0);
320 0 : return z;
321 : }
322 154 : l = coeff(x,1,1);
323 476 : for (j=1; j<ly; j++) gel(z,j) = F3m_F3c_mul_i(x, gel(y,j), lx, l);
324 154 : 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 : }
|
