CoCalc -- FpV.c

Testing latest pari + WASM + node.js... and it works?! Wow.
pari / 2021-07-01 / pari-2.14.0.alpha / src / basemath / FpV.c
²⁸⁴⁷⁹ views
License: GPL3
ubuntu2004
1
/* Copyright (C) 2000  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
/********************************************************************/
21
/**                                                                **/
22
/**                           REDUCTION                            **/
23
/**                                                                **/
24
/********************************************************************/
25
/* z in Z^n, return lift(Col(z) * Mod(1,p)) */
26
GEN
27
FpC_red(GEN x, GEN p)
28
{ pari_APPLY_type(t_COL, modii(gel(x,i), p)) }
29

30
/* z in Z^n, return lift(Vec(z) * Mod(1,p)) */
31
GEN
32
FpV_red(GEN x, GEN p)
33
{ pari_APPLY_type(t_VEC, modii(gel(x,i), p)) }
34

35
GEN
36
FpC_center(GEN x, GEN p, GEN pov2)
37
{ pari_APPLY_type(t_COL, Fp_center(gel(x,i), p, pov2)) }
38

39
/* assume 0 <= u < p and ps2 = p>>1 */
40
INLINE void
41
Fp_center_inplace(GEN u, GEN p, GEN ps2)
42
{ if (abscmpii(u,ps2) > 0) subiiz(u,p,u); }
43

44
void
45
FpC_center_inplace(GEN z, GEN p, GEN ps2)
46
{
47
  long i,l = lg(z);
48
  for (i=1; i<l; i++)
49
    Fp_center_inplace(gel(z,i), p, ps2);
50
}
51

52
GEN
53
Flv_center(GEN x, ulong p, ulong ps2)
54
{ pari_APPLY_ulong(Fl_center(uel(x,i),p,ps2)) }
55

56
/* z in Mat m,n(Z), return lift(z * Mod(1,p)) */
57
GEN
58
FpM_red(GEN x, GEN p)
59
{ pari_APPLY_same(FpC_red(gel(x,i), p)) }
60

61
GEN
62
FpM_center(GEN x, GEN p, GEN pov2)
63
{ pari_APPLY_same(FpC_center(gel(x,i), p,pov2)) }
64

65
void
66
FpM_center_inplace(GEN z, GEN p, GEN pov2)
67
{
68
  long i, l = lg(z);
69
  for (i=1; i<l; i++) FpC_center_inplace(gel(z,i), p, pov2);
70
}
71
GEN
72
Flm_center(GEN x, ulong p, ulong ps2)
73
{ pari_APPLY_same(Flv_center(gel(x,i), p,ps2)) }
74

75
GEN
76
random_FpV(long d, GEN p)
77
{
78
  long i;
79
  GEN y = cgetg(d+1,t_VEC);
80
  for (i=1; i<=d; i++) gel(y,i) = randomi(p);
81
  return y;
82
}
83

84
GEN
85
random_FpC(long d, GEN p)
86
{
87
  long i;
88
  GEN y = cgetg(d+1,t_COL);
89
  for (i=1; i<=d; i++) gel(y,i) = randomi(p);
90
  return y;
91
}
92

93
GEN
94
random_Flv(long n, ulong p)
95
{
96
  GEN y = cgetg(n+1, t_VECSMALL);
97
  long i;
98
  for (i=1; i<=n; i++) uel(y,i) = random_Fl(p);
99
  return y;
100
}
101

102
/********************************************************************/
103
/**                                                                **/
104
/**                           ADD, SUB                             **/
105
/**                                                                **/
106
/********************************************************************/
107
GEN
108
FpC_add(GEN x, GEN y, GEN p)
109
{ pari_APPLY_type(t_COL, Fp_add(gel(x,i), gel(y,i), p)) }
110

111
GEN
112
FpV_add(GEN x, GEN y, GEN p)
113
{ pari_APPLY_type(t_VEC, Fp_add(gel(x,i), gel(y,i), p)) }
114

115
GEN
116
FpM_add(GEN x, GEN y, GEN p)
117
{ pari_APPLY_same(FpC_add(gel(x,i), gel(y,i), p)) }
118

119
GEN
120
Flv_add(GEN x, GEN y, ulong p)
121
{
122
  if (p==2)
123
    pari_APPLY_ulong(x[i]^y[i])
124
  else
125
    pari_APPLY_ulong(Fl_add(x[i], y[i], p))
126
}
127

128
void
129
Flv_add_inplace(GEN x, GEN y, ulong p)
130
{
131
  long i, l = lg(x);
132
  if (p==2)
133
    for (i = 1; i < l; i++) x[i] ^= y[i];
134
  else
135
    for (i = 1; i < l; i++) x[i] = Fl_add(x[i], y[i], p);
136
}
137

138
ulong
139
Flv_sum(GEN x, ulong p)
140
{
141
  long i, l = lg(x);
142
  ulong s = 0;
143
  if (p==2)
144
    for (i = 1; i < l; i++) s ^= x[i];
145
  else
146
    for (i = 1; i < l; i++) s = Fl_add(s, x[i], p);
147
  return s;
148
}
149

150
GEN
151
FpC_sub(GEN x, GEN y, GEN p)
152
{ pari_APPLY_type(t_COL, Fp_sub(gel(x,i), gel(y,i), p)) }
153

154
GEN
155
FpV_sub(GEN x, GEN y, GEN p)
156
{ pari_APPLY_type(t_VEC, Fp_sub(gel(x,i), gel(y,i), p)) }
157

158
GEN
159
FpM_sub(GEN x, GEN y, GEN p)
160
{ pari_APPLY_same(FpC_sub(gel(x,i), gel(y,i), p)) }
161

162
GEN
163
Flv_sub(GEN x, GEN y, ulong p)
164
{ pari_APPLY_ulong(Fl_sub(uel(x,i), uel(y,i), p)) }
165

166
void
167
Flv_sub_inplace(GEN x, GEN y, ulong p)
168
{
169
  long i, l = lg(x);
170
  for (i = 1; i < l; i++) x[i] = Fl_sub(x[i], y[i], p);
171
}
172

173
GEN
174
Flm_Fl_add(GEN x, ulong y, ulong p)
175
{
176
  long l = lg(x), i, j;
177
  GEN z = cgetg(l,t_MAT);
178

179
  if (l==1) return z;
180
  if (l != lgcols(x)) pari_err_OP( "+", x, utoi(y));
181
  for (i=1; i<l; i++)
182
  {
183
    GEN zi = cgetg(l,t_VECSMALL), xi = gel(x,i);
184
    gel(z,i) = zi;
185
    for (j=1; j<l; j++) zi[j] = xi[j];
186
    zi[i] = Fl_add(zi[i], y, p);
187
  }
188
  return z;
189
}
190
GEN
191
Flm_Fl_sub(GEN x, ulong y, ulong p) { return Flm_Fl_add(x, Fl_neg(y,p), p); }
192

193
GEN
194
Flm_add(GEN x, GEN y, ulong p)
195
{ pari_APPLY_same(Flv_add(gel(x,i), gel(y,i), p)) }
196

197
GEN
198
Flm_sub(GEN x, GEN y, ulong p)
199
{ pari_APPLY_same(Flv_sub(gel(x,i), gel(y,i), p)) }
200

201
/********************************************************************/
202
/**                                                                **/
203
/**                           MULTIPLICATION                       **/
204
/**                                                                **/
205
/********************************************************************/
206
GEN
207
FpC_Fp_mul(GEN x, GEN y, GEN p)
208
{ pari_APPLY_type(t_COL, Fp_mul(gel(x,i), y, p)) }
209

210
GEN
211
Flv_Fl_mul(GEN x, ulong y, ulong p)
212
{ pari_APPLY_ulong(Fl_mul(x[i], y, p)) }
213

214
GEN
215
Flv_Fl_div(GEN x, ulong y, ulong p)
216
{
217
  return Flv_Fl_mul(x, Fl_inv(y, p), p);
218
}
219
void
220
Flv_Fl_div_inplace(GEN x, ulong y, ulong p)
221
{
222
  Flv_Fl_mul_inplace(x, Fl_inv(y, p), p);
223
}
224
GEN
225
FpM_Fp_mul(GEN X, GEN c, GEN p) {
226
  long i, j, h, l = lg(X);
227
  GEN A = cgetg(l, t_MAT);
228
  if (l == 1) return A;
229
  h = lgcols(X);
230
  for (j=1; j<l; j++)
231
  {
232
    GEN a = cgetg(h, t_COL), x = gel(X, j);
233
    for (i = 1; i < h; i++) gel(a,i) = Fp_mul(gel(x,i), c, p);
234
    gel(A,j) = a;
235
  }
236
  return A;
237
}
238

239
/* x *= y */
240
void
241
Flv_Fl_mul_part_inplace(GEN x, ulong y, ulong p, long l)
242
{
243
  long i;
244
  if (HIGHWORD(y | p))
245
    for(i=1; i<=l; i++) x[i] = Fl_mul(x[i], y, p);
246
  else
247
    for(i=1; i<=l; i++) x[i] = (x[i] * y) % p;
248
}
249
void
250
Flv_Fl_mul_inplace(GEN x, ulong y, ulong p)
251
{ Flv_Fl_mul_part_inplace(x, y, p, lg(x)-1); }
252

253
/* set y *= x */
254
void
255
Flm_Fl_mul_inplace(GEN y, ulong x, ulong p)
256
{
257
  long i, j, m, l = lg(y);
258
  if (l == 1) return;
259
  m = lgcols(y);
260
  if (HIGHWORD(x | p))
261
    for(j=1; j<l; j++)
262
      for(i=1; i<m; i++) ucoeff(y,i,j) = Fl_mul(ucoeff(y,i,j), x, p);
263
  else
264
    for(j=1; j<l; j++)
265
      for(i=1; i<m; i++) ucoeff(y,i,j) = (ucoeff(y,i,j) * x) % p;
266
}
267

268
/* return x * y */
269
static GEN
270
Flm_Fl_mul_pre_i(GEN y, ulong x, ulong p, ulong pi)
271
{
272
  long i, j, m, l = lg(y);
273
  GEN z = cgetg(l, t_MAT);
274
  if (l == 1) return z;
275
  m = lgcols(y);
276
  for(j=1; j<l; j++) {
277
    GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
278
    for(i=1; i<m; i++) uel(c,i) = Fl_mul_pre(ucoeff(y,i,j), x, p, pi);
279
  }
280
  return z;
281
}
282

283
/* return x * y */
284
static GEN
285
Flm_Fl_mul_OK(GEN y, ulong x, ulong p)
286
{
287
  long i, j, m, l = lg(y);
288
  GEN z = cgetg(l, t_MAT);
289
  if (l == 1) return z;
290
  m = lgcols(y);
291
  for(j=1; j<l; j++) {
292
    GEN c = cgetg(m, t_VECSMALL); gel(z,j) = c;
293
    for(i=1; i<m; i++) uel(c,i) = (ucoeff(y,i,j) * x) % p;
294
  }
295
  return z;
296
}
297

298
/* return x * y */
299
GEN
300
Flm_Fl_mul_pre(GEN y, ulong x, ulong p, ulong pi)
301
{
302
  if (HIGHWORD(p))
303
    return Flm_Fl_mul_pre_i(y, x, p, pi);
304
  else
305
    return Flm_Fl_mul_OK(y, x, p);
306
}
307

308
/* return x * y */
309
GEN
310
Flm_Fl_mul(GEN y, ulong x, ulong p)
311
{
312
  if (HIGHWORD(p))
313
    return Flm_Fl_mul_pre_i(y, x, p, get_Fl_red(p));
314
  else
315
    return Flm_Fl_mul_OK(y, x, p);
316
}
317

318
GEN
319
Flv_neg(GEN x, ulong p)
320
{ pari_APPLY_ulong(Fl_neg(uel(x,i), p)) }
321

322
void
323
Flv_neg_inplace(GEN v, ulong p)
324
{
325
  long i;
326
  for (i = 1; i < lg(v); ++i)
327
    v[i] = Fl_neg(v[i], p);
328
}
329

330
GEN
331
Flm_neg(GEN x, ulong p)
332
{ pari_APPLY_same(Flv_neg(gel(x,i), p)) }
333

334
/* x[i,]*y. Assume lx > 1 and 0 < i < lgcols(x) */
335
static GEN
336
ZMrow_ZC_mul_i(GEN x, GEN y, long lx, long i)
337
{
338
  GEN c = mulii(gcoeff(x,i,1), gel(y,1));
339
  long k;
340
  for (k = 2; k < lx; k++)
341
  {
342
    GEN t = mulii(gcoeff(x,i,k), gel(y,k));
343
    if (signe(t)) c = addii(c, t);
344
  }
345
  return c;
346
}
347

348
static long
349
zmrow_zc_mul(GEN x, GEN y, long lx, long i)
350
{
351
  long k;
352
  long c = coeff(x,i,1) * y[1];
353
  for (k = 2; k < lx; k++)
354
    c += coeff(x,i,k) * y[k];
355
  return c;
356
}
357

358
GEN
359
zm_zc_mul(GEN x, GEN y)
360
{
361
  long lx = lg(x), l, i;
362
  GEN z;
363
  if (lx == 1) return cgetg(1, t_VECSMALL);
364
  l = lg(gel(x,1));
365
  z = cgetg(l,t_VECSMALL);
366
  for (i=1; i<l; i++) z[i] = zmrow_zc_mul(x, y, lx, i);
367
  return z;
368
}
369

370
GEN
371
zm_mul(GEN x, GEN y)
372
{
373
  long i,j,lx=lg(x), ly=lg(y);
374
  GEN z;
375
  if (ly==1) return cgetg(1,t_MAT);
376
  z = cgetg(ly,t_MAT);
377
  if (lx==1)
378
  {
379
    for (i=1; i<ly; i++) gel(z,i) = cgetg(1,t_VECSMALL);
380
    return z;
381
  }
382
  for (j=1; j<ly; j++)
383
    gel(z,j) = zm_zc_mul(x, gel(y,j));
384
  return z;
385
}
386

387
static ulong
388
Flmrow_Flc_mul_SMALL(GEN x, GEN y, ulong p, long lx, long i)
389
{
390
  ulong c = ucoeff(x,i,1) * uel(y,1);
391
  long k;
392
  for (k = 2; k < lx; k++) {
393
    c += ucoeff(x,i,k) * uel(y,k);
394
    if (c & HIGHBIT) c %= p;
395
  }
396
  return c % p;
397
}
398

399
static ulong
400
Flmrow_Flc_mul_i(GEN x, GEN y, ulong p, ulong pi, long lx, long i)
401
{
402
  ulong l0, l1, v1, h0, h1;
403
  long k = 1;
404
  LOCAL_OVERFLOW;
405
  LOCAL_HIREMAINDER;
406
  l1 = mulll(ucoeff(x,i,k), uel(y,k)); h1 = hiremainder; v1 = 0;
407
  while (++k < lx) {
408
    l0 = mulll(ucoeff(x,i,k), uel(y,k)); h0 = hiremainder;
409
    l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
410
  }
411
  if (v1 == 0) return remll_pre(h1, l1, p, pi);
412
  else return remlll_pre(v1, h1, l1, p, pi);
413
}
414

415
static GEN
416
Flm_Flc_mul_i_2(GEN x, GEN y, long lx, long l)
417
{
418
  long i,j;
419
  GEN z = NULL;
420

421
  for (j=1; j<lx; j++)
422
  {
423
    if (!y[j]) continue;
424
    if (!z) z = Flv_copy(gel(x,j));
425
    else for (i = 1; i < l; i++) z[i] ^= coeff(x,i,j);
426
  }
427
  if (!z) z = zero_zv(l-1);
428
  return z;
429
}
430

431
static GEN
432
FpM_FpC_mul_i(GEN x, GEN y, long lx, long l, GEN p)
433
{
434
  GEN z = cgetg(l,t_COL);
435
  long i;
436
  for (i = 1; i < l; i++)
437
  {
438
    pari_sp av = avma;
439
    GEN c = ZMrow_ZC_mul_i(x, y, lx, i);
440
    gel(z,i) = gerepileuptoint(av, modii(c,p));
441
  }
442
  return z;
443
}
444

445
static void
446
__Flm_Flc_mul_i_SMALL(GEN z, GEN x, GEN y, long lx, long l, ulong p)
447
{
448
  long i;
449
  for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_SMALL(x, y, p, lx, i);
450
}
451
static GEN
452
Flm_Flc_mul_i_SMALL(GEN x, GEN y, long lx, long l, ulong p)
453
{
454
  GEN z = cgetg(l,t_VECSMALL);
455
  __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
456
  return z;
457
}
458

459
static void
460
__Flm_Flc_mul_i(GEN z, GEN x, GEN y, long lx, long l, ulong p, ulong pi)
461
{
462
  long i;
463
  for (i = 1; i < l; i++) z[i] = Flmrow_Flc_mul_i(x, y, p, pi, lx, i);
464
}
465
static GEN
466
Flm_Flc_mul_i(GEN x, GEN y, long lx, long l, ulong p, ulong pi)
467
{
468
  GEN z = cgetg(l,t_VECSMALL);
469
  __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
470
  return z;
471
}
472

473
GEN
474
FpM_mul(GEN x, GEN y, GEN p)
475
{
476
  long lx=lg(x), ly=lg(y);
477
  GEN z;
478
  pari_sp av = avma;
479
  if (ly==1) return cgetg(1,t_MAT);
480
  if (lx==1) return zeromat(0, ly-1);
481
  if (lgefint(p) == 3)
482
  {
483
    ulong pp = uel(p,2);
484
    if (pp == 2)
485
    {
486
      x = ZM_to_F2m(x);
487
      y = ZM_to_F2m(y);
488
      z = F2m_to_ZM(F2m_mul(x,y));
489
    }
490
    else
491
    {
492
      x = ZM_to_Flm(x, pp);
493
      y = ZM_to_Flm(y, pp);
494
      z = Flm_to_ZM(Flm_mul(x,y, pp));
495
    }
496
  } else
497
    z = FpM_red(ZM_mul(x, y), p);
498
  return gerepileupto(av, z);
499
}
500

501
static GEN
502
Flm_mul_classical(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
503
{
504
  long j;
505
  GEN z = cgetg(ly, t_MAT);
506
  if (SMALL_ULONG(p))
507
    for (j=1; j<ly; j++)
508
      gel(z,j) = Flm_Flc_mul_i_SMALL(x, gel(y,j), lx, l, p);
509
  else
510
    for (j=1; j<ly; j++)
511
      gel(z,j) = Flm_Flc_mul_i(x, gel(y,j), lx, l, p, pi);
512
  return z;
513
}
514

515
/*
516
  Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
517
  as an (m x n)-matrix, padding the input with zeroes as necessary.
518
*/
519
static void
520
add_slices_ip(long m, long n,
521
           GEN A, long ma, long da, long na, long ea,
522
           GEN B, long mb, long db, long nb, long eb,
523
           GEN M, long dy, long dx, ulong p)
524
{
525
  long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
526
  GEN C;
527

528
  for (j = 1; j <= min_e; j++) {
529
    C = gel(M, j + dx) + dy;
530
    for (i = 1; i <= min_d; i++)
531
      uel(C, i) = Fl_add(ucoeff(A, ma + i, na + j),
532
                        ucoeff(B, mb + i, nb + j), p);
533
    for (; i <= da; i++)
534
      uel(C, i) = ucoeff(A, ma + i, na + j);
535
    for (; i <= db; i++)
536
      uel(C, i) = ucoeff(B, mb + i, nb + j);
537
    for (; i <= m; i++)
538
      uel(C, i) = 0;
539
  }
540
  for (; j <= ea; j++) {
541
    C = gel(M, j + dx) + dy;
542
    for (i = 1; i <= da; i++)
543
      uel(C, i) = ucoeff(A, ma + i, na + j);
544
    for (; i <= m; i++)
545
      uel(C, i) = 0;
546
  }
547
  for (; j <= eb; j++) {
548
    C = gel(M, j + dx) + dy;
549
    for (i = 1; i <= db; i++)
550
      uel(C, i) = ucoeff(B, mb + i, nb + j);
551
    for (; i <= m; i++)
552
      uel(C, i) = 0;
553
  }
554
  for (; j <= n; j++) {
555
    C = gel(M, j + dx) + dy;
556
    for (i = 1; i <= m; i++)
557
      uel(C, i) = 0;
558
  }
559
}
560

561
static GEN
562
add_slices(long m, long n,
563
           GEN A, long ma, long da, long na, long ea,
564
           GEN B, long mb, long db, long nb, long eb, ulong p)
565
{
566
  GEN M;
567
  long j;
568
  M = cgetg(n + 1, t_MAT);
569
  for (j = 1; j <= n; j++)
570
    gel(M, j) = cgetg(m + 1, t_VECSMALL);
571
  add_slices_ip(m, n, A, ma, da, na, ea, B, mb, db, nb, eb, M, 0, 0, p);
572
  return M;
573
}
574

575
/*
576
  Return A[ma+1..ma+da, na+1..na+ea] - B[mb+1..mb+db, nb+1..nb+eb]
577
  as an (m x n)-matrix, padding the input with zeroes as necessary.
578
*/
579
static GEN
580
subtract_slices(long m, long n,
581
                GEN A, long ma, long da, long na, long ea,
582
                GEN B, long mb, long db, long nb, long eb, ulong p)
583
{
584
  long min_d = minss(da, db), min_e = minss(ea, eb), i, j;
585
  GEN M = cgetg(n + 1, t_MAT), C;
586

587
  for (j = 1; j <= min_e; j++) {
588
    gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
589
    for (i = 1; i <= min_d; i++)
590
      uel(C, i) = Fl_sub(ucoeff(A, ma + i, na + j),
591
                        ucoeff(B, mb + i, nb + j), p);
592
    for (; i <= da; i++)
593
      uel(C, i) = ucoeff(A, ma + i, na + j);
594
    for (; i <= db; i++)
595
      uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
596
    for (; i <= m; i++)
597
      uel(C, i) = 0;
598
  }
599
  for (; j <= ea; j++) {
600
    gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
601
    for (i = 1; i <= da; i++)
602
      uel(C, i) = ucoeff(A, ma + i, na + j);
603
    for (; i <= m; i++)
604
      uel(C, i) = 0;
605
  }
606
  for (; j <= eb; j++) {
607
    gel(M, j) = C = cgetg(m + 1, t_VECSMALL);
608
    for (i = 1; i <= db; i++)
609
      uel(C, i) = Fl_neg(ucoeff(B, mb + i, nb + j), p);
610
    for (; i <= m; i++)
611
      uel(C, i) = 0;
612
  }
613
  for (; j <= n; j++)
614
    gel(M, j) = zero_Flv(m);
615
  return M;
616
}
617

618
static GEN Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi);
619

620
/* Strassen-Winograd matrix product A (m x n) * B (n x p) */
621
static GEN
622
Flm_mul_sw(GEN A, GEN B, long m, long n, long p, ulong l, ulong li)
623
{
624
  pari_sp av;
625
  long m1 = (m + 1)/2, m2 = m/2,
626
    n1 = (n + 1)/2, n2 = n/2,
627
    p1 = (p + 1)/2, p2 = p/2;
628
  GEN A11, A12, A22, B11, B21, B22,
629
    S1, S2, S3, S4, T1, T2, T3, T4,
630
    M1, M2, M3, M4, M5, M6, M7,
631
    V1, V2, V3, C;
632
  long j;
633
  C = cgetg(p + 1, t_MAT);
634
  for (j = 1; j <= p; j++)
635
    gel(C, j) = cgetg(m + 1, t_VECSMALL);
636
  av = avma;
637
  T2 = subtract_slices(n1, p2, B, 0, n1, p1, p2, B, n1, n2, p1, p2, l);
638
  S1 = subtract_slices(m2, n1, A, m1, m2, 0, n1, A, 0, m2, 0, n1, l);
639
  M2 = Flm_mul_i(S1, T2, m2 + 1, n1 + 1, p2 + 1, l, li);
640
  if (gc_needed(av, 1))
641
    gerepileall(av, 2, &T2, &M2);  /* destroy S1 */
642
  T3 = subtract_slices(n1, p1, T2, 0, n1, 0, p2, B, 0, n1, 0, p1, l);
643
  if (gc_needed(av, 1))
644
    gerepileall(av, 2, &M2, &T3);  /* destroy T2 */
645
  S2 = add_slices(m2, n1, A, m1, m2, 0, n1, A, m1, m2, n1, n2, l);
646
  T1 = subtract_slices(n1, p1, B, 0, n1, p1, p2, B, 0, n1, 0, p2, l);
647
  M3 = Flm_mul_i(S2, T1, m2 + 1, n1 + 1, p2 + 1, l, li);
648
  if (gc_needed(av, 1))
649
    gerepileall(av, 4, &M2, &T3, &S2, &M3);  /* destroy T1 */
650
  S3 = subtract_slices(m1, n1, S2, 0, m2, 0, n1, A, 0, m1, 0, n1, l);
651
  if (gc_needed(av, 1))
652
    gerepileall(av, 4, &M2, &T3, &M3, &S3);  /* destroy S2 */
653
  A11 = matslice(A, 1, m1, 1, n1);
654
  B11 = matslice(B, 1, n1, 1, p1);
655
  M1 = Flm_mul_i(A11, B11, m1 + 1, n1 + 1, p1 + 1, l, li);
656
  if (gc_needed(av, 1))
657
    gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy A11, B11 */
658
  A12 = matslice(A, 1, m1, n1 + 1, n);
659
  B21 = matslice(B, n1 + 1, n, 1, p1);
660
  M4 = Flm_mul_i(A12, B21, m1 + 1, n2 + 1, p1 + 1, l, li);
661
  if (gc_needed(av, 1))
662
    gerepileall(av, 6, &M2, &T3, &M3, &S3, &M1, &M4);  /* destroy A12, B21 */
663
  add_slices_ip(m1, p1, M1, 0, m1, 0, p1, M4, 0, m1, 0, p1, C, 0, 0, l);
664
  if (gc_needed(av, 1))
665
    gerepileall(av, 5, &M2, &T3, &M3, &S3, &M1);  /* destroy M4 */
666
  M5 = Flm_mul_i(S3, T3, m1 + 1, n1 + 1, p1 + 1, l, li);
667
  S4 = subtract_slices(m1, n2, A, 0, m1, n1, n2, S3, 0, m1, 0, n2, l);
668
  if (gc_needed(av, 1))
669
    gerepileall(av, 6, &M2, &T3, &M3, &M1, &M5, &S4);  /* destroy S3 */
670
  T4 = add_slices(n2, p1, B, n1, n2, 0, p1, T3, 0, n2, 0, p1, l);
671
  if (gc_needed(av, 1))
672
    gerepileall(av, 6, &M2, &M3, &M1, &M5, &S4, &T4);  /* destroy T3 */
673
  V1 = subtract_slices(m1, p1, M1, 0, m1, 0, p1, M5, 0, m1, 0, p1, l);
674
  if (gc_needed(av, 1))
675
    gerepileall(av, 5, &M2, &M3, &S4, &T4, &V1);  /* destroy M1, M5 */
676
  B22 = matslice(B, n1 + 1, n, p1 + 1, p);
677
  M6 = Flm_mul_i(S4, B22, m1 + 1, n2 + 1, p2 + 1, l, li);
678
  if (gc_needed(av, 1))
679
    gerepileall(av, 5, &M2, &M3, &T4, &V1, &M6);  /* destroy S4, B22 */
680
  A22 = matslice(A, m1 + 1, m, n1 + 1, n);
681
  M7 = Flm_mul_i(A22, T4, m2 + 1, n2 + 1, p1 + 1, l, li);
682
  if (gc_needed(av, 1))
683
    gerepileall(av, 5, &M2, &M3, &V1, &M6, &M7);  /* destroy A22, T4 */
684
  V3 = add_slices(m1, p2, V1, 0, m1, 0, p2, M3, 0, m2, 0, p2, l);
685
  add_slices_ip(m1, p2, V3, 0, m1, 0, p2, M6, 0, m1, 0, p2, C, 0, p1, l);
686
  if (gc_needed(av, 1))
687
    gerepileall(av, 4, &M2, &M3, &V1, &M7);  /* destroy V3, M6 */
688
  V2 = add_slices(m2, p1, V1, 0, m2, 0, p1, M2, 0, m2, 0, p2, l);
689
  if (gc_needed(av, 1))
690
    gerepileall(av, 3, &M3, &M7, &V2);  /* destroy V1, M2 */
691
  add_slices_ip(m2, p1, V2, 0, m2, 0, p1, M7, 0, m2, 0, p1, C, m1, 0, l);
692
  add_slices_ip(m2, p2, V2, 0, m2, 0, p2, M3, 0, m2, 0, p2, C, m1, p1, l);
693
  set_avma(av); return C;
694
}
695

696
/* Strassen-Winograd used for dim >= ZM_sw_bound */
697
static GEN
698
Flm_mul_i(GEN x, GEN y, long l, long lx, long ly, ulong p, ulong pi)
699
{
700
  ulong e = expu(p);
701
#ifdef LONG_IS_64BIT
702
  long ZM_sw_bound = e <= 29 ? 140: e <=62 ? 40: 70;
703
#else
704
  long ZM_sw_bound = e <= 12 ? 230: e <=14 ? 170 : e <=17 ? 110: 120;
705
#endif
706
  if (l <= ZM_sw_bound || lx <= ZM_sw_bound || ly <= ZM_sw_bound)
707
    return Flm_mul_classical(x, y, l, lx, ly, p, pi);
708
  else
709
    return Flm_mul_sw(x, y, l - 1, lx - 1, ly - 1, p, pi);
710
}
711

712
GEN
713
Flm_mul_pre(GEN x, GEN y, ulong p, ulong pi)
714
{
715
  long lx=lg(x), ly=lg(y);
716
  if (ly==1) return cgetg(1,t_MAT);
717
  if (lx==1) return zero_Flm(0, ly-1);
718
  return Flm_mul_i(x, y, lgcols(x), lx, ly, p, pi);
719
}
720

721
GEN
722
Flm_mul(GEN x, GEN y, ulong p)
723
{
724
  long lx=lg(x), ly=lg(y);
725
  if (ly==1) return cgetg(1,t_MAT);
726
  if (lx==1) return zero_Flm(0, ly-1);
727
  return Flm_mul_i(x, y, lgcols(x), lx, ly, p, get_Fl_red(p));
728
}
729

730
struct _Flm
731
{
732
  ulong p;
733
  long n;
734
};
735

736
static GEN
737
_Flm_mul(void *E , GEN x, GEN y)
738
{ return Flm_mul(x,y,((struct _Flm*)E)->p); }
739
static GEN
740
_Flm_sqr(void *E, GEN x)
741
{ return Flm_mul(x,x,((struct _Flm*)E)->p); }
742
static GEN
743
_Flm_one(void *E)
744
{ return matid_Flm(((struct _Flm*)E)->n); }
745
GEN
746
Flm_powu(GEN x, ulong n, ulong p)
747
{
748
  struct _Flm d;
749
  if (!n) return matid(lg(x)-1);
750
  d.p = p;
751
  return gen_powu(x, n, (void*)&d, &_Flm_sqr, &_Flm_mul);
752
}
753
GEN
754
Flm_powers(GEN x, ulong n, ulong p)
755
{
756
  struct _Flm d;
757
  d.p = p;
758
  d.n = lg(x)-1;
759
  return gen_powers(x, n, 1, (void*)&d,
760
                    &_Flm_sqr, &_Flm_mul, &_Flm_one);
761
}
762

763
static GEN
764
_FpM_mul(void *p , GEN x, GEN y)
765
{ return FpM_mul(x,y,(GEN)p); }
766
static GEN
767
_FpM_sqr(void *p, GEN x)
768
{ return FpM_mul(x,x,(GEN)p); }
769
GEN
770
FpM_powu(GEN x, ulong n, GEN p)
771
{
772
  if (!n) return matid(lg(x)-1);
773
  if (lgefint(p) == 3)
774
  {
775
    pari_sp av = avma;
776
    ulong pp = uel(p,2);
777
    GEN z;
778
    if (pp == 2)
779
      z = F2m_to_ZM(F2m_powu(ZM_to_F2m(x),n));
780
    else
781
      z = Flm_to_ZM(Flm_powu(ZM_to_Flm(x, pp), n, pp));
782
    return gerepileupto(av, z);
783
  }
784
  return gen_powu(x, n, (void*)p, &_FpM_sqr, &_FpM_mul);
785
}
786

787
/*Multiple a column vector by a line vector to make a matrix*/
788
GEN
789
Flc_Flv_mul(GEN x, GEN y, ulong p)
790
{
791
  long i,j, lx=lg(x), ly=lg(y);
792
  GEN z;
793
  if (ly==1) return cgetg(1,t_MAT);
794
  z = cgetg(ly, t_MAT);
795
  for (j=1; j < ly; j++)
796
  {
797
    GEN zj = cgetg(lx,t_VECSMALL);
798
    for (i=1; i<lx; i++)
799
      uel(zj,i) = Fl_mul(uel(x,i), uel(y,j), p);
800
    gel(z,j) = zj;
801
  }
802
  return z;
803
}
804

805
/*Multiple a column vector by a line vector to make a matrix*/
806
GEN
807
FpC_FpV_mul(GEN x, GEN y, GEN p)
808
{
809
  long i,j, lx=lg(x), ly=lg(y);
810
  GEN z;
811
  if (ly==1) return cgetg(1,t_MAT);
812
  z = cgetg(ly,t_MAT);
813
  for (j=1; j < ly; j++)
814
  {
815
    GEN zj = cgetg(lx,t_COL);
816
    for (i=1; i<lx; i++) gel(zj,i) = Fp_mul(gel(x,i),gel(y,j), p);
817
    gel(z, j) = zj;
818
  }
819
  return z;
820
}
821

822
/* Multiply a line vector by a column and return a scalar (t_INT) */
823
GEN
824
FpV_dotproduct(GEN x, GEN y, GEN p)
825
{
826
  long i, lx = lg(x);
827
  pari_sp av;
828
  GEN c;
829
  if (lx == 1) return gen_0;
830
  av = avma; c = mulii(gel(x,1),gel(y,1));
831
  for (i=2; i<lx; i++) c = addii(c, mulii(gel(x,i),gel(y,i)));
832
  return gerepileuptoint(av, modii(c,p));
833
}
834
GEN
835
FpV_dotsquare(GEN x, GEN p)
836
{
837
  long i, lx = lg(x);
838
  pari_sp av;
839
  GEN c;
840
  if (lx == 1) return gen_0;
841
  av = avma; c = sqri(gel(x,1));
842
  for (i=2; i<lx; i++) c = addii(c, sqri(gel(x,i)));
843
  return gerepileuptoint(av, modii(c,p));
844
}
845

846
INLINE ulong
847
Flv_dotproductspec_SMALL(GEN x, GEN y, ulong p, long lx)
848
{
849
  ulong c = uel(x,0) * uel(y,0);
850
  long k;
851
  for (k = 1; k < lx; k++) {
852
    c += uel(x,k) * uel(y,k);
853
    if (c & HIGHBIT) c %= p;
854
  }
855
  return c % p;
856
}
857

858
INLINE ulong
859
Flv_dotproductspec_i(GEN x, GEN y, ulong p, ulong pi, long lx)
860
{
861
  ulong l0, l1, v1, h0, h1;
862
  long i = 0;
863
  LOCAL_OVERFLOW;
864
  LOCAL_HIREMAINDER;
865
  l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
866
  while (++i < lx) {
867
    l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
868
    l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
869
  }
870
  if (v1 == 0) return remll_pre(h1, l1, p, pi);
871
  else return remlll_pre(v1, h1, l1, p, pi);
872
}
873

874
ulong
875
Flv_dotproduct(GEN x, GEN y, ulong p)
876
{
877
  long lx = lg(x)-1;
878
  if (lx == 0) return 0;
879
  if (SMALL_ULONG(p))
880
    return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
881
  else
882
    return Flv_dotproductspec_i(x+1, y+1, p, get_Fl_red(p), lx);
883
}
884

885
ulong
886
Flv_dotproduct_pre(GEN x, GEN y, ulong p, ulong pi)
887
{
888
  long lx = lg(x)-1;
889
  if (lx == 0) return 0;
890
  if (SMALL_ULONG(p))
891
    return Flv_dotproductspec_SMALL(x+1, y+1, p, lx);
892
  else
893
    return Flv_dotproductspec_i(x+1, y+1, p, pi, lx);
894
}
895

896
ulong
897
Flx_dotproduct(GEN x, GEN y, ulong p)
898
{
899
  long lx = minss(lgpol(x), lgpol(y));
900
  if (lx == 0) return 0;
901
  if (SMALL_ULONG(p))
902
    return Flv_dotproductspec_SMALL(x+2, y+2, p, lx);
903
  else
904
    return Flv_dotproductspec_i(x+2, y+2, p, get_Fl_red(p), lx);
905
}
906

907
GEN
908
FpM_FpC_mul(GEN x, GEN y, GEN p)
909
{
910
  long lx = lg(x);
911
  return lx==1? cgetg(1,t_COL): FpM_FpC_mul_i(x, y, lx, lgcols(x), p);
912
}
913
GEN
914
Flm_Flc_mul(GEN x, GEN y, ulong p)
915
{
916
  long l, lx = lg(x);
917
  if (lx==1) return cgetg(1,t_VECSMALL);
918
  l = lgcols(x);
919
  if (p==2)
920
    return Flm_Flc_mul_i_2(x, y, lx, l);
921
  else if (SMALL_ULONG(p))
922
    return Flm_Flc_mul_i_SMALL(x, y, lx, l, p);
923
  else
924
    return Flm_Flc_mul_i(x, y, lx, l, p, get_Fl_red(p));
925
}
926

927
GEN
928
Flm_Flc_mul_pre(GEN x, GEN y, ulong p, ulong pi)
929
{
930
  long l, lx = lg(x);
931
  GEN z;
932
  if (lx==1) return cgetg(1,t_VECSMALL);
933
  l = lgcols(x);
934
  z = cgetg(l, t_VECSMALL);
935
  if (SMALL_ULONG(p))
936
    __Flm_Flc_mul_i_SMALL(z, x, y, lx, l, p);
937
  else
938
    __Flm_Flc_mul_i(z, x, y, lx, l, p, pi);
939
  return z;
940
}
941

942
GEN
943
Flm_Flc_mul_pre_Flx(GEN x, GEN y, ulong p, ulong pi, long sv)
944
{
945
  long l, lx = lg(x);
946
  GEN z;
947
  if (lx==1) return pol0_Flx(sv);
948
  l = lgcols(x);
949
  z = cgetg(l + 1, t_VECSMALL); z[1] = sv;
950
  if (SMALL_ULONG(p))
951
    __Flm_Flc_mul_i_SMALL(z + 1, x, y, lx, l, p);
952
  else
953
    __Flm_Flc_mul_i(z + 1, x, y, lx, l, p, pi);
954
  return Flx_renormalize(z, l + 1);
955
}
956

957
/* RgV_to_RgX(FpM_FpC_mul(x,y,p), v), p != NULL, memory clean */
958
GEN
959
FpM_FpC_mul_FpX(GEN x, GEN y, GEN p, long v)
960
{
961
  long i, l, lx = lg(x);
962
  GEN z;
963
  if (lx==1) return pol_0(v);
964
  l = lgcols(x);
965
  z = new_chunk(l+1);
966
  for (i=l-1; i; i--)
967
  {
968
    pari_sp av = avma;
969
    GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
970
    p1 = modii(p1, p);
971
    if (signe(p1))
972
    {
973
      if (i != l-1) stackdummy((pari_sp)(z + l+1), (pari_sp)(z + i+2));
974
      gel(z,i+1) = gerepileuptoint(av, p1);
975
      break;
976
    }
977
    set_avma(av);
978
  }
979
  if (!i) { set_avma((pari_sp)(z + l+1)); return pol_0(v); }
980
  z[0] = evaltyp(t_POL) | evallg(i+2);
981
  z[1] = evalsigne(1) | evalvarn(v);
982
  for (; i; i--)
983
  {
984
    pari_sp av = avma;
985
    GEN p1 = ZMrow_ZC_mul_i(x,y,lx,i);
986
    gel(z,i+1) = gerepileuptoint(av, modii(p1,p));
987
  }
988
  return z;
989
}
990

991
/********************************************************************/
992
/**                                                                **/
993
/**                           TRANSPOSITION                        **/
994
/**                                                                **/
995
/********************************************************************/
996

997
/* == zm_transpose */
998
GEN
999
Flm_transpose(GEN x)
1000
{
1001
  long i, dx, lx = lg(x);
1002
  GEN y;
1003
  if (lx == 1) return cgetg(1,t_MAT);
1004
  dx = lgcols(x); y = cgetg(dx,t_MAT);
1005
  for (i=1; i<dx; i++) gel(y,i) = Flm_row(x,i);
1006
  return y;
1007
}
1008

1009
/********************************************************************/
1010
/**                                                                **/
1011
/**                           SCALAR MATRICES                      **/
1012
/**                                                                **/
1013
/********************************************************************/
1014

1015
GEN
1016
gen_matid(long n, void *E, const struct bb_field *S)
1017
{
1018
  GEN y = cgetg(n+1,t_MAT), _0, _1;
1019
  long i;
1020
  if (n < 0) pari_err_DOMAIN("gen_matid", "dimension","<",gen_0,stoi(n));
1021
  _0 = S->s(E,0);
1022
  _1 = S->s(E,1);
1023
  for (i=1; i<=n; i++)
1024
  {
1025
    GEN z = const_col(n, _0); gel(z,i) = _1;
1026
    gel(y, i) = z;
1027
  }
1028
  return y;
1029
}
1030

1031
GEN
1032
matid_F2xqM(long n, GEN T)
1033
{
1034
  void *E;
1035
  const struct bb_field *S = get_F2xq_field(&E, T);
1036
  return gen_matid(n, E, S);
1037
}
1038
GEN
1039
matid_FlxqM(long n, GEN T, ulong p)
1040
{
1041
  void *E;
1042
  const struct bb_field *S = get_Flxq_field(&E, T, p);
1043
  return gen_matid(n, E, S);
1044
}
1045

1046
GEN
1047
matid_Flm(long n)
1048
{
1049
  GEN y = cgetg(n+1,t_MAT);
1050
  long i;
1051
  if (n < 0) pari_err_DOMAIN("matid_Flm", "dimension","<",gen_0,stoi(n));
1052
  for (i=1; i<=n; i++) { gel(y,i) = zero_zv(n); ucoeff(y, i,i) = 1; }
1053
  return y;
1054
}
1055

1056
GEN
1057
scalar_Flm(long s, long n)
1058
{
1059
  long i;
1060
  GEN y = cgetg(n+1,t_MAT);
1061
  for (i=1; i<=n; i++) { gel(y,i) = zero_Flv(n); coeff(y, i,i) = s; }
1062
  return y;
1063
}
1064

1065
/********************************************************************/
1066
/**                                                                **/
1067
/**                           CONVERSIONS                          **/
1068
/**                                                                **/
1069
/********************************************************************/
1070
GEN
1071
ZV_to_Flv(GEN x, ulong p)
1072
{ pari_APPLY_ulong(umodiu(gel(x,i), p)) }
1073

1074
GEN
1075
ZM_to_Flm(GEN x, ulong p)
1076
{ pari_APPLY_same(ZV_to_Flv(gel(x,i), p)) }
1077

1078
GEN
1079
ZMV_to_FlmV(GEN x, ulong m)
1080
{ pari_APPLY_type(t_VEC,ZM_to_Flm(gel(x,i), m)) }
1081

1082
/*                          TO INTMOD                        */
1083
static GEN
1084
to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
1085
static GEN
1086
Fl_to_intmod(ulong x, GEN p) { retmkintmod(utoi(x), p); }
1087

1088
GEN
1089
Fp_to_mod(GEN z, GEN p)
1090
{
1091
  retmkintmod(modii(z, p), icopy(p));
1092
}
1093

1094
GEN
1095
FpX_to_mod_raw(GEN z, GEN p)
1096
{
1097
  long i,l = lg(z);
1098
  GEN x;
1099
  if (l == 2)
1100
  {
1101
    x = cgetg(3,t_POL); x[1] = z[1];
1102
    gel(x,2) = mkintmod(gen_0,p); return x;
1103
  }
1104
  x = cgetg(l,t_POL);
1105
  for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1106
  x[1] = z[1]; return normalizepol_lg(x,l);
1107
}
1108

1109
/* z in Z[X], return z * Mod(1,p), normalized*/
1110
GEN
1111
FpX_to_mod(GEN z, GEN p)
1112
{
1113
  long i,l = lg(z);
1114
  GEN x;
1115
  if (l == 2)
1116
  {
1117
    x = cgetg(3,t_POL); x[1] = z[1];
1118
    gel(x,2) = mkintmod(gen_0,icopy(p)); return x;
1119
  }
1120
  x = cgetg(l,t_POL);
1121
  if (l >2) p = icopy(p);
1122
  for (i=2; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1123
  x[1] = z[1]; return normalizepol_lg(x,l);
1124
}
1125

1126
GEN
1127
FpXC_to_mod(GEN z, GEN p)
1128
{
1129
  long i,l = lg(z);
1130
  GEN x = cgetg(l, t_COL);
1131
  p = icopy(p);
1132
  for (i=1; i<l; i++)
1133
    gel(x,i) = FpX_to_mod_raw(gel(z,i), p);
1134
  return x;
1135
}
1136

1137
static GEN
1138
FpXC_to_mod_raw(GEN x, GEN p)
1139
{ pari_APPLY_type(t_COL, FpX_to_mod_raw(gel(x,i), p)) }
1140

1141
GEN
1142
FpXM_to_mod(GEN z, GEN p)
1143
{
1144
  long i,l = lg(z);
1145
  GEN x = cgetg(l, t_MAT);
1146
  p = icopy(p);
1147
  for (i=1; i<l; i++)
1148
    gel(x,i) = FpXC_to_mod_raw(gel(z,i), p);
1149
  return x;
1150
}
1151

1152
/* z in Z^n, return z * Mod(1,p), normalized*/
1153
GEN
1154
FpV_to_mod(GEN z, GEN p)
1155
{
1156
  long i,l = lg(z);
1157
  GEN x = cgetg(l, t_VEC);
1158
  if (l == 1) return x;
1159
  p = icopy(p);
1160
  for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1161
  return x;
1162
}
1163
/* z in Z^n, return z * Mod(1,p), normalized*/
1164
GEN
1165
FpC_to_mod(GEN z, GEN p)
1166
{
1167
  long i, l = lg(z);
1168
  GEN x = cgetg(l, t_COL);
1169
  if (l == 1) return x;
1170
  p = icopy(p);
1171
  for (i=1; i<l; i++) gel(x,i) = to_intmod(gel(z,i), p);
1172
  return x;
1173
}
1174
/* z in Mat m,n(Z), return z * Mod(1,p), normalized*/
1175
GEN
1176
FpM_to_mod(GEN z, GEN p)
1177
{
1178
  long i, j, m, l = lg(z);
1179
  GEN  x = cgetg(l,t_MAT), y, zi;
1180
  if (l == 1) return x;
1181
  m = lgcols(z);
1182
  p = icopy(p);
1183
  for (i=1; i<l; i++)
1184
  {
1185
    gel(x,i) = cgetg(m,t_COL);
1186
    y = gel(x,i); zi= gel(z,i);
1187
    for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1188
  }
1189
  return x;
1190
}
1191
GEN
1192
Flc_to_mod(GEN z, ulong pp)
1193
{
1194
  long i, l = lg(z);
1195
  GEN p, x = cgetg(l, t_COL);
1196
  if (l == 1) return x;
1197
  p = utoipos(pp);
1198
  for (i=1; i<l; i++) gel(x,i) = Fl_to_intmod(z[i], p);
1199
  return x;
1200
}
1201
GEN
1202
Flm_to_mod(GEN z, ulong pp)
1203
{
1204
  long i, j, m, l = lg(z);
1205
  GEN p, x = cgetg(l,t_MAT), y, zi;
1206
  if (l == 1) return x;
1207
  m = lgcols(z);
1208
  p = utoipos(pp);
1209
  for (i=1; i<l; i++)
1210
  {
1211
    gel(x,i) = cgetg(m,t_COL);
1212
    y = gel(x,i); zi= gel(z,i);
1213
    for (j=1; j<m; j++) gel(y,j) = Fl_to_intmod(zi[j], p);
1214
  }
1215
  return x;
1216
}
1217

1218
GEN
1219
FpVV_to_mod(GEN z, GEN p)
1220
{
1221
  long i, j, m, l = lg(z);
1222
  GEN  x = cgetg(l,t_VEC), y, zi;
1223
  if (l == 1) return x;
1224
  m = lgcols(z);
1225
  p = icopy(p);
1226
  for (i=1; i<l; i++)
1227
  {
1228
    gel(x,i) = cgetg(m,t_VEC);
1229
    y = gel(x,i); zi= gel(z,i);
1230
    for (j=1; j<m; j++) gel(y,j) = to_intmod(gel(zi,j), p);
1231
  }
1232
  return x;
1233
}
1234

1235
/* z in Z^n, return z * Mod(1,p), normalized*/
1236
GEN
1237
FpXQC_to_mod(GEN z, GEN T, GEN p)
1238
{
1239
  long i,l = lg(z);
1240
  GEN x = cgetg(l, t_COL);
1241
  if (l == 1) return x;
1242
  p = icopy(p);
1243
  T = FpX_to_mod_raw(T, p);
1244
  for (i=1; i<l; i++)
1245
    gel(x,i) = mkpolmod(FpX_to_mod_raw(gel(z,i), p), T);
1246
  return x;
1247
}
1248

1249
static GEN
1250
Fq_to_mod_raw(GEN x, GEN T, GEN p)
1251
{
1252
  GEN z = typ(x)==t_INT ? mkintmod(modii(x, p), p): FpX_to_mod_raw(x, p);
1253
  return mkpolmod(z, T);
1254
}
1255

1256
/* z in Z^n, return z * Mod(1,p), normalized*/
1257
GEN
1258
FqC_to_mod(GEN z, GEN T, GEN p)
1259
{
1260
  GEN x;
1261
  long i,l = lg(z);
1262
  if (!T) return FpC_to_mod(z, p);
1263
  x = cgetg(l, t_COL);
1264
  if (l == 1) return x;
1265
  p = icopy(p);
1266
  T = FpX_to_mod_raw(T, p);
1267
  for (i=1; i<l; i++)
1268
    gel(x,i) = Fq_to_mod_raw(gel(z, i), T, p);
1269
  return x;
1270
}
1271

1272
/* z in Z^n, return z * Mod(1,p), normalized*/
1273
static GEN
1274
FqC_to_mod_raw(GEN x, GEN T, GEN p)
1275
{ pari_APPLY_type(t_COL, Fq_to_mod_raw(gel(x,i), T, p)) }
1276

1277
/* z in Z^n, return z * Mod(1,p), normalized*/
1278
GEN
1279
FqM_to_mod(GEN z, GEN T, GEN p)
1280
{
1281
  GEN x;
1282
  long i,l = lg(z);
1283
  if (!T) return FpM_to_mod(z, p);
1284
  x = cgetg(l, t_MAT);
1285
  if (l == 1) return x;
1286
  p = icopy(p);
1287
  T = FpX_to_mod_raw(T, p);
1288
  for (i=1; i<l; i++)
1289
    gel(x,i) = FqC_to_mod_raw(gel(z, i), T, p);
1290
  return x;
1291
}
1292

1293
/********************************************************************/
1294
/*                                                                  */
1295
/*                     Blackbox linear algebra                      */
1296
/*                                                                  */
1297
/********************************************************************/
1298

1299
/* A sparse column (zCs) v is a t_COL with two components C and E which are
1300
 * t_VECSMALL of the same length, representing sum_i E[i]*e_{C[i]}, where
1301
 * (e_j) is the canonical basis.
1302
 * A sparse matrix (zMs) is a t_VEC of zCs */
1303

1304
/* FpCs and FpMs are identical but E[i] is interpreted as a _signed_ C long
1305
 * integer representing an element of Fp. This is important since p can be
1306
 * large and p+E[i] would not fit in a C long.  */
1307

1308
/* RgCs and RgMs are similar, except that the type of the component is
1309
 * unspecified. Functions handling RgCs/RgMs must be independent of the type
1310
 * of E. */
1311

1312
/* Most functions take an argument nbrow which is the number of lines of the
1313
 * column/matrix, which cannot be derived from the data. */
1314

1315
GEN
1316
zCs_to_ZC(GEN R, long nbrow)
1317
{
1318
  GEN C = gel(R,1), E = gel(R,2), c = zerocol(nbrow);
1319
  long j, l = lg(C);
1320
  for (j = 1; j < l; ++j) gel(c, C[j]) = stoi(E[j]);
1321
  return c;
1322
}
1323

1324
GEN
1325
zMs_to_ZM(GEN x, long nbrow)
1326
{ pari_APPLY_same(zCs_to_ZC(gel(x, i), nbrow)) }
1327

1328
/* Solve equation f(X) = B (mod p) where B is a FpV, and f is an endomorphism.
1329
 * Return either a solution as a t_COL, or a kernel vector as a t_VEC. */
1330
GEN
1331
gen_FpM_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p)
1332
{
1333
  pari_sp ltop = avma;
1334
  long col = 0, n = lg(B)-1, m = 2*n+1;
1335
  if (ZV_equal0(B)) return zerocol(n);
1336
  while (++col <= n)
1337
  {
1338
    pari_sp btop = avma, av;
1339
    long i, lQ;
1340
    GEN V, Q, M, W = B;
1341
    GEN b = cgetg(m+2, t_POL);
1342
    b[1] = evalsigne(1)|evalvarn(0);
1343
    gel(b, 2) = gel(W, col);
1344
    for (i = 3; i<m+2; i++)
1345
      gel(b, i) = cgeti(lgefint(p));
1346
    av = avma;
1347
    for (i = 3; i<m+2; i++)
1348
    {
1349
      W = f(E, W);
1350
      affii(gel(W, col),gel(b, i));
1351
      if (gc_needed(av,1))
1352
      {
1353
        if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: first loop, %ld",i);
1354
        W = gerepileupto(av, W);
1355
      }
1356
    }
1357
    b = FpX_renormalize(b, m+2);
1358
    if (lgpol(b)==0) {set_avma(btop); continue; }
1359
    M = FpX_halfgcd(b, pol_xn(m, 0), p);
1360
    Q = FpX_neg(FpX_normalize(gcoeff(M, 2, 1),p),p);
1361
    W = B; lQ =lg(Q);
1362
    if (DEBUGLEVEL) err_printf("Wiedemann: deg. minpoly: %ld\n",lQ-3);
1363
    V = FpC_Fp_mul(W, gel(Q, lQ-2), p);
1364
    av = avma;
1365
    for (i = lQ-3; i > 1; i--)
1366
    {
1367
      W = f(E, W);
1368
      V = ZC_lincomb(gen_1, gel(Q,i), V, W);
1369
      if (gc_needed(av,1))
1370
      {
1371
        if (DEBUGMEM>1) pari_warn(warnmem,"Wiedemann: second loop, %ld",i);
1372
        gerepileall(av, 2, &V, &W);
1373
      }
1374
    }
1375
    V = FpC_red(V, p);
1376
    W = FpC_sub(f(E,V), B, p);
1377
    if (ZV_equal0(W)) return gerepilecopy(ltop, V);
1378
    av = avma;
1379
    for (i = 1; i <= n; ++i)
1380
    {
1381
      V = W;
1382
      W = f(E, V);
1383
      if (ZV_equal0(W))
1384
        return gerepilecopy(ltop, shallowtrans(V));
1385
      gerepileall(av, 2, &V, &W);
1386
    }
1387
    set_avma(btop);
1388
  }
1389
  return NULL;
1390
}
1391

1392
GEN
1393
zMs_ZC_mul(GEN M, GEN B)
1394
{
1395
  long i, j;
1396
  long n = lg(B)-1;
1397
  GEN V = zerocol(n);
1398
  for (i = 1; i <= n; ++i)
1399
    if (signe(gel(B, i)))
1400
    {
1401
      GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1402
      long l = lg(C);
1403
      for (j = 1; j < l; ++j)
1404
      {
1405
        long k = C[j];
1406
        switch(E[j])
1407
        {
1408
        case 1:
1409
          gel(V, k) = gel(V,k)==gen_0 ? gel(B,i) : addii(gel(V, k), gel(B,i));
1410
          break;
1411
        case -1:
1412
          gel(V, k) = gel(V,k)==gen_0 ? negi(gel(B,i)) : subii(gel(V, k), gel(B,i));
1413
          break;
1414
        default:
1415
          gel(V, k) = gel(V,k)==gen_0 ? mulis(gel(B, i), E[j]) : addii(gel(V, k), mulis(gel(B, i), E[j]));
1416
          break;
1417
        }
1418
      }
1419
    }
1420
  return V;
1421
}
1422

1423
GEN
1424
FpMs_FpC_mul(GEN M, GEN B, GEN p) { return FpC_red(zMs_ZC_mul(M, B), p); }
1425

1426
GEN
1427
ZV_zMs_mul(GEN B, GEN M)
1428
{
1429
  long i, j;
1430
  long m = lg(M)-1;
1431
  GEN V = cgetg(m+1,t_VEC);
1432
  for (i = 1; i <= m; ++i)
1433
  {
1434
    GEN R = gel(M, i), C = gel(R, 1), E = gel(R, 2);
1435
    long l = lg(C);
1436
    GEN z;
1437
    if (l == 1)
1438
    {
1439
      gel(V,i) = gen_0;
1440
      continue;
1441
    }
1442
    z = mulis(gel(B, C[1]), E[1]);
1443
    for (j = 2; j < l; ++j)
1444
    {
1445
      long k = C[j];
1446
      switch(E[j])
1447
      {
1448
      case 1:
1449
        z = addii(z, gel(B,k));
1450
        break;
1451
      case -1:
1452
        z = subii(z, gel(B,k));
1453
        break;
1454
      default:
1455
        z = addii(z, mulis(gel(B,k), E[j]));
1456
        break;
1457
      }
1458
    }
1459
    gel(V,i) = z;
1460
  }
1461
  return V;
1462
}
1463

1464
GEN
1465
FpV_FpMs_mul(GEN B, GEN M, GEN p) { return FpV_red(ZV_zMs_mul(B, M), p); }
1466

1467
GEN
1468
ZlM_gauss(GEN a, GEN b, ulong p, long e, GEN C)
1469
{
1470
  pari_sp av = avma, av2;
1471
  GEN xi, xb, pi = gen_1, P;
1472
  long i;
1473
  if (!C) {
1474
    C = Flm_inv(ZM_to_Flm(a, p), p);
1475
    if (!C) pari_err_INV("ZlM_gauss", a);
1476
  }
1477
  P = utoipos(p);
1478
  av2 = avma;
1479
  xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1480
  xb = Flm_to_ZM(xi);
1481
  for (i = 2; i <= e; i++)
1482
  {
1483
    pi = muliu(pi, p); /* = p^(i-1) */
1484
    b = ZM_Z_divexact(ZM_sub(b, ZM_nm_mul(a, xi)), P);
1485
    if (gc_needed(av,2))
1486
    {
1487
      if(DEBUGMEM>1) pari_warn(warnmem,"ZlM_gauss. i=%ld/%ld",i,e);
1488
      gerepileall(av2,3, &pi,&b,&xb);
1489
    }
1490
    xi = Flm_mul(C, ZM_to_Flm(b, p), p);
1491
    xb = ZM_add(xb, nm_Z_mul(xi, pi));
1492
  }
1493
  return gerepileupto(av, xb);
1494
}
1495

1496
struct wrapper_modp_s {
1497
  GEN (*f)(void*, GEN);
1498
  void *E;
1499
  GEN p;
1500
};
1501

1502
/* compute f(x) mod p */
1503
static GEN
1504
wrap_relcomb_modp(void *E, GEN x)
1505
{
1506
  struct wrapper_modp_s *W = (struct wrapper_modp_s*)E;
1507
  return FpC_red(W->f(W->E, x), W->p);
1508
}
1509
static GEN
1510
wrap_relcomb(void*E, GEN x) { return zMs_ZC_mul((GEN)E, x); }
1511

1512
static GEN
1513
wrap_relker(void*E, GEN x) { return ZV_zMs_mul(x, (GEN)E); }
1514

1515
/* Solve f(X) = B (mod p^e); blackbox version of ZlM_gauss */
1516
GEN
1517
gen_ZpM_Dixon_Wiedemann(void *E, GEN (*f)(void*, GEN), GEN B, GEN p, long e)
1518
{
1519
  struct wrapper_modp_s W;
1520
  pari_sp av = avma;
1521
  GEN xb, xi, pi = gen_1;
1522
  long i;
1523
  W.E = E;
1524
  W.f = f;
1525
  W.p = p;
1526
  xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p); /* f^(-1) B */
1527
  if (!xi || e == 1 || typ(xi) == t_VEC) return xi;
1528
  xb = xi;
1529
  for (i = 2; i <= e; i++)
1530
  {
1531
    pi = mulii(pi, p); /* = p^(i-1) */
1532
    B = ZC_Z_divexact(ZC_sub(B, f(E, xi)), p);
1533
    if (gc_needed(av,2))
1534
    {
1535
      if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Dixon_Wiedemann. i=%ld",i);
1536
      gerepileall(av,3, &pi,&B,&xb);
1537
    }
1538
    xi = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, FpC_red(B, p), p);
1539
    if (!xi) return NULL;
1540
    if (typ(xi) == t_VEC) return gerepileupto(av, xi);
1541
    xb = ZC_add(xb, ZC_Z_mul(xi, pi));
1542
  }
1543
  return gerepileupto(av, xb);
1544
}
1545

1546
static GEN
1547
vecprow(GEN A, GEN prow)
1548
{
1549
  return mkvec2(vecsmallpermute(prow,gel(A,1)), gel(A,2));
1550
}
1551

1552
/* Solve the equation MX = A. Return either a solution as a t_COL,
1553
 * or the index of a column which is linearly dependent from the others as a
1554
 * t_VECSMALL with a single component. */
1555
GEN
1556
ZpMs_ZpCs_solve(GEN M, GEN A, long nbrow, GEN p, long e)
1557
{
1558
  pari_sp av = avma;
1559
  GEN pcol, prow;
1560
  long nbi=lg(M)-1, lR;
1561
  long i, n;
1562
  GEN Mp, Ap, Rp;
1563
  pari_timer ti;
1564
  if (DEBUGLEVEL) timer_start(&ti);
1565
  RgMs_structelim(M, nbrow, gel(A, 1), &pcol, &prow);
1566
  if (!pcol) return gc_NULL(av);
1567
  if (DEBUGLEVEL)
1568
    timer_printf(&ti,"structured elimination (%ld -> %ld)",nbi,lg(pcol)-1);
1569
  n = lg(pcol)-1;
1570
  Mp = cgetg(n+1, t_MAT);
1571
  for(i=1; i<=n; i++)
1572
    gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1573
  Ap = zCs_to_ZC(vecprow(A, prow), n);
1574
  if (DEBUGLEVEL) timer_start(&ti);
1575
  Rp = gen_ZpM_Dixon_Wiedemann((void*)Mp,wrap_relcomb, Ap, p, e);
1576
  if (DEBUGLEVEL) timer_printf(&ti,"linear algebra");
1577
  if (!Rp) return gc_NULL(av);
1578
  lR = lg(Rp)-1;
1579
  if (typ(Rp) == t_COL)
1580
  {
1581
    GEN R = zerocol(nbi+1);
1582
    for (i=1; i<=lR; i++)
1583
       gel(R,pcol[i]) = gel(Rp,i);
1584
    return gerepilecopy(av, R);
1585
  }
1586
  for (i = 1; i <= lR; ++i)
1587
    if (signe(gel(Rp, i)))
1588
      return gerepileuptoleaf(av, mkvecsmall(pcol[i]));
1589
  return NULL; /* LCOV_EXCL_LINE */
1590
}
1591

1592
GEN
1593
FpMs_FpCs_solve(GEN M, GEN A, long nbrow, GEN p)
1594
{
1595
  return ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1596
}
1597

1598
GEN
1599
FpMs_FpCs_solve_safe(GEN M, GEN A, long nbrow, GEN p)
1600
{
1601
  GEN res;
1602
  pari_CATCH(e_INV)
1603
  {
1604
    GEN E = pari_err_last();
1605
    GEN x = err_get_compo(E,2);
1606
    if (typ(x) != t_INTMOD) pari_err(0,E);
1607
    if (DEBUGLEVEL)
1608
      pari_warn(warner,"FpMs_FpCs_solve_safe, impossible inverse %Ps", x);
1609
    res = NULL;
1610
  } pari_TRY {
1611
    res = ZpMs_ZpCs_solve(M, A, nbrow, p, 1);
1612
  } pari_ENDCATCH
1613
  return res;
1614
}
1615

1616
static GEN
1617
FpMs_structelim_back(GEN M, GEN V, GEN prow, GEN p)
1618
{
1619
  long i, j, oldf = 0, f = 0;
1620
  long lrow = lg(prow), lM = lg(M);
1621
  GEN W = const_vecsmall(lM-1,1);
1622
  GEN R = cgetg(lrow, t_VEC);
1623
  for (i=1; i<lrow; i++)
1624
    gel(R,i) = prow[i]==0 ? NULL: gel(V,prow[i]);
1625
  do
1626
  {
1627
    oldf = f;
1628
    for(i=1; i<lM; i++)
1629
      if (W[i])
1630
      {
1631
        GEN C = gel(M,i), F = gel(C,1), E = gel(C,2);
1632
        long c=0, cj=0, lF = lg(F);
1633
        for(j=1; j<lF; j++)
1634
          if (!gel(R,F[j])) { c++; cj=j; }
1635
        if (c>=2) continue;
1636
        if (c==1)
1637
        {
1638
          pari_sp av = avma;
1639
          GEN s = gen_0;
1640
          for(j=1; j<lF; j++)
1641
            if (j!=cj) s = Fp_add(s, mulis(gel(R,F[j]), E[j]), p);
1642
          /* s /= -E[cj] mod p */
1643
          s = E[cj] < 0? Fp_divu(s, -E[cj], p): Fp_divu(Fp_neg(s,p), E[cj], p);
1644
          gel(R,F[cj]) = gerepileuptoint(av, s);
1645
          f++;
1646
        }
1647
        W[i]=0;
1648
      }
1649
  } while(oldf!=f);
1650
  for (i=1; i<lrow; i++)
1651
    if (!gel(R,i)) gel(R,i) = gen_0;
1652
  return R;
1653
}
1654

1655
/* Return a linear form Y such that YM is essentially 0 */
1656
GEN
1657
FpMs_leftkernel_elt_col(GEN M, long nbcol, long nbrow, GEN p)
1658
{
1659
  pari_sp av = avma, av2;
1660
  GEN pcol, prow;
1661
  long i, n;
1662
  GEN Mp, B, MB, R, Rp;
1663
  pari_timer ti;
1664
  struct wrapper_modp_s W;
1665
  if (DEBUGLEVEL) timer_start(&ti);
1666
  RgMs_structelim_col(M, nbcol, nbrow, cgetg(1,t_VECSMALL), &pcol, &prow);
1667
  if (!pcol) return gc_NULL(av);
1668
  if (DEBUGLEVEL)
1669
    timer_printf(&ti,"structured elimination (%ld -> %ld)",nbcol,lg(pcol)-1);
1670
  n = lg(pcol)-1;
1671
  Mp = cgetg(n+1, t_MAT);
1672
  for(i=1; i<=n; i++)
1673
    gel(Mp, i) = vecprow(gel(M,pcol[i]), prow);
1674
  W.E = (void*) Mp;
1675
  W.f = wrap_relker;
1676
  W.p = p;
1677
  av2 = avma;
1678
  for(;;)
1679
  {
1680
    set_avma(av2);
1681
    B = random_FpV(n, p);
1682
    MB = FpV_FpMs_mul(B, Mp, p);
1683
    if (DEBUGLEVEL) timer_start(&ti);
1684
    pari_CATCH(e_INV)
1685
    {
1686
      GEN E = pari_err_last();
1687
      GEN x = err_get_compo(E,2);
1688
      if (typ(x) != t_INTMOD) pari_err(0,E);
1689
      if (DEBUGLEVEL)
1690
        pari_warn(warner,"FpMs_leftkernel_elt, impossible inverse %Ps", x);
1691
      Rp = NULL;
1692
    } pari_TRY {
1693
      Rp = gen_FpM_Wiedemann((void*)&W, wrap_relcomb_modp, MB, p);
1694
    } pari_ENDCATCH
1695
    if (!Rp) continue;
1696
    if (typ(Rp)==t_COL)
1697
    {
1698
      Rp = FpC_sub(Rp,B,p);
1699
      if (ZV_equal0(Rp)) continue;
1700
    }
1701
    R = FpMs_structelim_back(M, Rp, prow, p);
1702
    if (DEBUGLEVEL) timer_printf(&ti,"Wiedemann left kernel");
1703
    return gerepilecopy(av, R);
1704
  }
1705
}
1706

1707
GEN
1708
FpMs_leftkernel_elt(GEN M, long nbrow, GEN p)
1709
{
1710
  return FpMs_leftkernel_elt_col(M, lg(M)-1, nbrow, p);
1711
}
1712

1713
Product

Resources

Company
