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 /*******************************************************************/
16 /* */
17 /* BASIC NF OPERATIONS */
18 /* */
19 /*******************************************************************/
20 #include "pari.h"
21 #include "paripriv.h"
22
23 /*******************************************************************/
24 /* */
25 /* OPERATIONS OVER NUMBER FIELD ELEMENTS. */
26 /* represented as column vectors over the integral basis */
27 /* */
28 /*******************************************************************/
29 static GEN
get_tab(GEN nf,long * N)30 get_tab(GEN nf, long *N)
31 {
32 GEN tab = (typ(nf) == t_MAT)? nf: gel(nf,9);
33 *N = nbrows(tab); return tab;
34 }
35
36 /* x != 0, y t_INT. Return x * y (not memory clean if x = 1) */
37 static GEN
_mulii(GEN x,GEN y)38 _mulii(GEN x, GEN y) {
39 return is_pm1(x)? (signe(x) < 0)? negi(y): y
40 : mulii(x, y);
41 }
42
43 GEN
tablemul_ei_ej(GEN M,long i,long j)44 tablemul_ei_ej(GEN M, long i, long j)
45 {
46 long N;
47 GEN tab = get_tab(M, &N);
48 tab += (i-1)*N; return gel(tab,j);
49 }
50
51 /* Outputs x.ei, where ei is the i-th elt of the algebra basis.
52 * x an RgV of correct length and arbitrary content (polynomials, scalars...).
53 * M is the multiplication table ei ej = sum_k M_k^(i,j) ek */
54 GEN
tablemul_ei(GEN M,GEN x,long i)55 tablemul_ei(GEN M, GEN x, long i)
56 {
57 long j, k, N;
58 GEN v, tab;
59
60 if (i==1) return gcopy(x);
61 tab = get_tab(M, &N);
62 if (typ(x) != t_COL) { v = zerocol(N); gel(v,i) = gcopy(x); return v; }
63 tab += (i-1)*N; v = cgetg(N+1,t_COL);
64 /* wi . x = [ sum_j tab[k,j] x[j] ]_k */
65 for (k=1; k<=N; k++)
66 {
67 pari_sp av = avma;
68 GEN s = gen_0;
69 for (j=1; j<=N; j++)
70 {
71 GEN c = gcoeff(tab,k,j);
72 if (!gequal0(c)) s = gadd(s, gmul(c, gel(x,j)));
73 }
74 gel(v,k) = gerepileupto(av,s);
75 }
76 return v;
77 }
78 /* as tablemul_ei, assume x a ZV of correct length */
79 GEN
zk_ei_mul(GEN nf,GEN x,long i)80 zk_ei_mul(GEN nf, GEN x, long i)
81 {
82 long j, k, N;
83 GEN v, tab;
84
85 if (i==1) return ZC_copy(x);
86 tab = get_tab(nf, &N); tab += (i-1)*N;
87 v = cgetg(N+1,t_COL);
88 for (k=1; k<=N; k++)
89 {
90 pari_sp av = avma;
91 GEN s = gen_0;
92 for (j=1; j<=N; j++)
93 {
94 GEN c = gcoeff(tab,k,j);
95 if (signe(c)) s = addii(s, _mulii(c, gel(x,j)));
96 }
97 gel(v,k) = gerepileuptoint(av, s);
98 }
99 return v;
100 }
101
102 /* table of multiplication by wi in R[w1,..., wN] */
103 GEN
ei_multable(GEN TAB,long i)104 ei_multable(GEN TAB, long i)
105 {
106 long k,N;
107 GEN m, tab = get_tab(TAB, &N);
108 tab += (i-1)*N;
109 m = cgetg(N+1,t_MAT);
110 for (k=1; k<=N; k++) gel(m,k) = gel(tab,k);
111 return m;
112 }
113
114 GEN
zk_multable(GEN nf,GEN x)115 zk_multable(GEN nf, GEN x)
116 {
117 long i, l = lg(x);
118 GEN mul = cgetg(l,t_MAT);
119 gel(mul,1) = x; /* assume w_1 = 1 */
120 for (i=2; i<l; i++) gel(mul,i) = zk_ei_mul(nf,x,i);
121 return mul;
122 }
123 GEN
multable(GEN M,GEN x)124 multable(GEN M, GEN x)
125 {
126 long i, N;
127 GEN mul;
128 if (typ(x) == t_MAT) return x;
129 M = get_tab(M, &N);
130 if (typ(x) != t_COL) return scalarmat(x, N);
131 mul = cgetg(N+1,t_MAT);
132 gel(mul,1) = x; /* assume w_1 = 1 */
133 for (i=2; i<=N; i++) gel(mul,i) = tablemul_ei(M,x,i);
134 return mul;
135 }
136
137 /* x integral in nf; table of multiplication by x in ZK = Z[w1,..., wN].
138 * Return a t_INT if x is scalar, and a ZM otherwise */
139 GEN
zk_scalar_or_multable(GEN nf,GEN x)140 zk_scalar_or_multable(GEN nf, GEN x)
141 {
142 long tx = typ(x);
143 if (tx == t_MAT || tx == t_INT) return x;
144 x = nf_to_scalar_or_basis(nf, x);
145 return (typ(x) == t_COL)? zk_multable(nf, x): x;
146 }
147
148 GEN
nftrace(GEN nf,GEN x)149 nftrace(GEN nf, GEN x)
150 {
151 pari_sp av = avma;
152 nf = checknf(nf);
153 x = nf_to_scalar_or_basis(nf, x);
154 x = (typ(x) == t_COL)? RgV_dotproduct(x, gel(nf_get_Tr(nf),1))
155 : gmulgs(x, nf_get_degree(nf));
156 return gerepileupto(av, x);
157 }
158 GEN
rnfelttrace(GEN rnf,GEN x)159 rnfelttrace(GEN rnf, GEN x)
160 {
161 pari_sp av = avma;
162 checkrnf(rnf);
163 x = rnfeltabstorel(rnf, x);
164 x = (typ(x) == t_POLMOD)? rnfeltdown(rnf, gtrace(x))
165 : gmulgs(x, rnf_get_degree(rnf));
166 return gerepileupto(av, x);
167 }
168
169 /* assume nf is a genuine nf, fa a famat */
170 static GEN
famat_norm(GEN nf,GEN fa)171 famat_norm(GEN nf, GEN fa)
172 {
173 pari_sp av = avma;
174 GEN g = gel(fa,1), e = gel(fa,2), N = gen_1;
175 long i, l = lg(g);
176 for (i = 1; i < l; i++)
177 N = gmul(N, powgi(nfnorm(nf, gel(g,i)), gel(e,i)));
178 return gerepileupto(av, N);
179 }
180 GEN
nfnorm(GEN nf,GEN x)181 nfnorm(GEN nf, GEN x)
182 {
183 pari_sp av = avma;
184 GEN c, den;
185 long n;
186 nf = checknf(nf);
187 n = nf_get_degree(nf);
188 if (typ(x) == t_MAT) return famat_norm(nf, x);
189 x = nf_to_scalar_or_basis(nf, x);
190 if (typ(x)!=t_COL)
191 return gerepileupto(av, gpowgs(x, n));
192 x = nf_to_scalar_or_alg(nf, Q_primitive_part(x, &c));
193 x = Q_remove_denom(x, &den);
194 x = ZX_resultant_all(nf_get_pol(nf), x, den, 0);
195 return gerepileupto(av, c ? gmul(x, gpowgs(c, n)): x);
196 }
197
198 static GEN
to_RgX(GEN P,long vx)199 to_RgX(GEN P, long vx)
200 {
201 return varn(P) == vx ? P: scalarpol_shallow(P, vx);
202 }
203
204 GEN
rnfeltnorm(GEN rnf,GEN x)205 rnfeltnorm(GEN rnf, GEN x)
206 {
207 pari_sp av = avma;
208 GEN nf, pol;
209 long v = rnf_get_varn(rnf);
210 checkrnf(rnf);
211 x = liftpol_shallow(rnfeltabstorel(rnf, x));
212 nf = rnf_get_nf(rnf); pol = rnf_get_pol(rnf);
213 x = (typ(x) == t_POL)
214 ? rnfeltdown(rnf, nfX_resultant(nf,pol,to_RgX(x,v)))
215 : gpowgs(x, rnf_get_degree(rnf));
216 return gerepileupto(av, x);
217 }
218
219 /* x + y in nf */
220 GEN
nfadd(GEN nf,GEN x,GEN y)221 nfadd(GEN nf, GEN x, GEN y)
222 {
223 pari_sp av = avma;
224 GEN z;
225
226 nf = checknf(nf);
227 x = nf_to_scalar_or_basis(nf, x);
228 y = nf_to_scalar_or_basis(nf, y);
229 if (typ(x) != t_COL)
230 { z = (typ(y) == t_COL)? RgC_Rg_add(y, x): gadd(x,y); }
231 else
232 { z = (typ(y) == t_COL)? RgC_add(x, y): RgC_Rg_add(x, y); }
233 return gerepileupto(av, z);
234 }
235 /* x - y in nf */
236 GEN
nfsub(GEN nf,GEN x,GEN y)237 nfsub(GEN nf, GEN x, GEN y)
238 {
239 pari_sp av = avma;
240 GEN z;
241
242 nf = checknf(nf);
243 x = nf_to_scalar_or_basis(nf, x);
244 y = nf_to_scalar_or_basis(nf, y);
245 if (typ(x) != t_COL)
246 { z = (typ(y) == t_COL)? Rg_RgC_sub(x,y): gsub(x,y); }
247 else
248 { z = (typ(y) == t_COL)? RgC_sub(x,y): RgC_Rg_sub(x,y); }
249 return gerepileupto(av, z);
250 }
251
252 /* product of ZC x,y in nf; ( sum_i x_i sum_j y_j m^{i,j}_k )_k */
253 static GEN
nfmuli_ZC(GEN nf,GEN x,GEN y)254 nfmuli_ZC(GEN nf, GEN x, GEN y)
255 {
256 long i, j, k, N;
257 GEN TAB = get_tab(nf, &N), v = cgetg(N+1,t_COL);
258
259 for (k = 1; k <= N; k++)
260 {
261 pari_sp av = avma;
262 GEN s, TABi = TAB;
263 if (k == 1)
264 s = mulii(gel(x,1),gel(y,1));
265 else
266 s = addii(mulii(gel(x,1),gel(y,k)),
267 mulii(gel(x,k),gel(y,1)));
268 for (i=2; i<=N; i++)
269 {
270 GEN t, xi = gel(x,i);
271 TABi += N;
272 if (!signe(xi)) continue;
273
274 t = NULL;
275 for (j=2; j<=N; j++)
276 {
277 GEN p1, c = gcoeff(TABi, k, j); /* m^{i,j}_k */
278 if (!signe(c)) continue;
279 p1 = _mulii(c, gel(y,j));
280 t = t? addii(t, p1): p1;
281 }
282 if (t) s = addii(s, mulii(xi, t));
283 }
284 gel(v,k) = gerepileuptoint(av,s);
285 }
286 return v;
287 }
288 static int
is_famat(GEN x)289 is_famat(GEN x) { return typ(x) == t_MAT && lg(x) == 3; }
290 /* product of x and y in nf */
291 GEN
nfmul(GEN nf,GEN x,GEN y)292 nfmul(GEN nf, GEN x, GEN y)
293 {
294 GEN z;
295 pari_sp av = avma;
296
297 if (x == y) return nfsqr(nf,x);
298
299 nf = checknf(nf);
300 if (is_famat(x) || is_famat(y)) return famat_mul(x, y);
301 x = nf_to_scalar_or_basis(nf, x);
302 y = nf_to_scalar_or_basis(nf, y);
303 if (typ(x) != t_COL)
304 {
305 if (isintzero(x)) return gen_0;
306 z = (typ(y) == t_COL)? RgC_Rg_mul(y, x): gmul(x,y); }
307 else
308 {
309 if (typ(y) != t_COL)
310 {
311 if (isintzero(y)) return gen_0;
312 z = RgC_Rg_mul(x, y);
313 }
314 else
315 {
316 GEN dx, dy;
317 x = Q_remove_denom(x, &dx);
318 y = Q_remove_denom(y, &dy);
319 z = nfmuli_ZC(nf,x,y);
320 dx = mul_denom(dx,dy);
321 if (dx) z = ZC_Z_div(z, dx);
322 }
323 }
324 return gerepileupto(av, z);
325 }
326 /* square of ZC x in nf */
327 static GEN
nfsqri_ZC(GEN nf,GEN x)328 nfsqri_ZC(GEN nf, GEN x)
329 {
330 long i, j, k, N;
331 GEN TAB = get_tab(nf, &N), v = cgetg(N+1,t_COL);
332
333 for (k = 1; k <= N; k++)
334 {
335 pari_sp av = avma;
336 GEN s, TABi = TAB;
337 if (k == 1)
338 s = sqri(gel(x,1));
339 else
340 s = shifti(mulii(gel(x,1),gel(x,k)), 1);
341 for (i=2; i<=N; i++)
342 {
343 GEN p1, c, t, xi = gel(x,i);
344 TABi += N;
345 if (!signe(xi)) continue;
346
347 c = gcoeff(TABi, k, i);
348 t = signe(c)? _mulii(c,xi): NULL;
349 for (j=i+1; j<=N; j++)
350 {
351 c = gcoeff(TABi, k, j);
352 if (!signe(c)) continue;
353 p1 = _mulii(c, shifti(gel(x,j),1));
354 t = t? addii(t, p1): p1;
355 }
356 if (t) s = addii(s, mulii(xi, t));
357 }
358 gel(v,k) = gerepileuptoint(av,s);
359 }
360 return v;
361 }
362 /* square of x in nf */
363 GEN
nfsqr(GEN nf,GEN x)364 nfsqr(GEN nf, GEN x)
365 {
366 pari_sp av = avma;
367 GEN z;
368
369 nf = checknf(nf);
370 if (is_famat(x)) return famat_sqr(x);
371 x = nf_to_scalar_or_basis(nf, x);
372 if (typ(x) != t_COL) z = gsqr(x);
373 else
374 {
375 GEN dx;
376 x = Q_remove_denom(x, &dx);
377 z = nfsqri_ZC(nf,x);
378 if (dx) z = RgC_Rg_div(z, sqri(dx));
379 }
380 return gerepileupto(av, z);
381 }
382
383 /* x a ZC, v a t_COL of ZC/Z */
384 GEN
zkC_multable_mul(GEN v,GEN x)385 zkC_multable_mul(GEN v, GEN x)
386 {
387 long i, l = lg(v);
388 GEN y = cgetg(l, t_COL);
389 for (i = 1; i < l; i++)
390 {
391 GEN c = gel(v,i);
392 if (typ(c)!=t_COL) {
393 if (!isintzero(c)) c = ZC_Z_mul(gel(x,1), c);
394 } else {
395 c = ZM_ZC_mul(x,c);
396 if (ZV_isscalar(c)) c = gel(c,1);
397 }
398 gel(y,i) = c;
399 }
400 return y;
401 }
402
403 GEN
nfC_multable_mul(GEN v,GEN x)404 nfC_multable_mul(GEN v, GEN x)
405 {
406 long i, l = lg(v);
407 GEN y = cgetg(l, t_COL);
408 for (i = 1; i < l; i++)
409 {
410 GEN c = gel(v,i);
411 if (typ(c)!=t_COL) {
412 if (!isintzero(c)) c = RgC_Rg_mul(gel(x,1), c);
413 } else {
414 c = RgM_RgC_mul(x,c);
415 if (QV_isscalar(c)) c = gel(c,1);
416 }
417 gel(y,i) = c;
418 }
419 return y;
420 }
421
422 GEN
nfC_nf_mul(GEN nf,GEN v,GEN x)423 nfC_nf_mul(GEN nf, GEN v, GEN x)
424 {
425 long tx;
426 GEN y;
427
428 x = nf_to_scalar_or_basis(nf, x);
429 tx = typ(x);
430 if (tx != t_COL)
431 {
432 long l, i;
433 if (tx == t_INT)
434 {
435 long s = signe(x);
436 if (!s) return zerocol(lg(v)-1);
437 if (is_pm1(x)) return s > 0? leafcopy(v): RgC_neg(v);
438 }
439 l = lg(v); y = cgetg(l, t_COL);
440 for (i=1; i < l; i++)
441 {
442 GEN c = gel(v,i);
443 if (typ(c) != t_COL) c = gmul(c, x); else c = RgC_Rg_mul(c, x);
444 gel(y,i) = c;
445 }
446 return y;
447 }
448 else
449 {
450 GEN dx;
451 x = zk_multable(nf, Q_remove_denom(x,&dx));
452 y = nfC_multable_mul(v, x);
453 return dx? RgC_Rg_div(y, dx): y;
454 }
455 }
456 static GEN
mulbytab(GEN M,GEN c)457 mulbytab(GEN M, GEN c)
458 { return typ(c) == t_COL? RgM_RgC_mul(M,c): RgC_Rg_mul(gel(M,1), c); }
459 GEN
tablemulvec(GEN M,GEN x,GEN v)460 tablemulvec(GEN M, GEN x, GEN v)
461 {
462 long l, i;
463 GEN y;
464
465 if (typ(x) == t_COL && RgV_isscalar(x))
466 {
467 x = gel(x,1);
468 return typ(v) == t_POL? RgX_Rg_mul(v,x): RgV_Rg_mul(v,x);
469 }
470 x = multable(M, x); /* multiplication table by x */
471 y = cgetg_copy(v, &l);
472 if (typ(v) == t_POL)
473 {
474 y[1] = v[1];
475 for (i=2; i < l; i++) gel(y,i) = mulbytab(x, gel(v,i));
476 y = normalizepol(y);
477 }
478 else
479 {
480 for (i=1; i < l; i++) gel(y,i) = mulbytab(x, gel(v,i));
481 }
482 return y;
483 }
484
485 GEN
zkmultable_capZ(GEN mx)486 zkmultable_capZ(GEN mx) { return Q_denom(zkmultable_inv(mx)); }
487 GEN
zkmultable_inv(GEN mx)488 zkmultable_inv(GEN mx) { return ZM_gauss(mx, col_ei(lg(mx)-1,1)); }
489 /* nf a true nf, x a ZC */
490 GEN
zk_inv(GEN nf,GEN x)491 zk_inv(GEN nf, GEN x) { return zkmultable_inv(zk_multable(nf,x)); }
492
493 /* inverse of x in nf */
494 GEN
nfinv(GEN nf,GEN x)495 nfinv(GEN nf, GEN x)
496 {
497 pari_sp av = avma;
498 GEN z;
499
500 nf = checknf(nf);
501 if (is_famat(x)) return famat_inv(x);
502 x = nf_to_scalar_or_basis(nf, x);
503 if (typ(x) == t_COL)
504 {
505 GEN d;
506 x = Q_remove_denom(x, &d);
507 z = zk_inv(nf, x);
508 if (d) z = RgC_Rg_mul(z, d);
509 }
510 else
511 z = ginv(x);
512 return gerepileupto(av, z);
513 }
514
515 /* quotient of x and y in nf */
516 GEN
nfdiv(GEN nf,GEN x,GEN y)517 nfdiv(GEN nf, GEN x, GEN y)
518 {
519 pari_sp av = avma;
520 GEN z;
521
522 nf = checknf(nf);
523 if (is_famat(x) || is_famat(y)) return famat_div(x,y);
524 y = nf_to_scalar_or_basis(nf, y);
525 if (typ(y) != t_COL)
526 {
527 x = nf_to_scalar_or_basis(nf, x);
528 z = (typ(x) == t_COL)? RgC_Rg_div(x, y): gdiv(x,y);
529 }
530 else
531 {
532 GEN d;
533 y = Q_remove_denom(y, &d);
534 z = nfmul(nf, x, zk_inv(nf,y));
535 if (d) z = typ(z) == t_COL? RgC_Rg_mul(z, d): gmul(z, d);
536 }
537 return gerepileupto(av, z);
538 }
539
540 /* product of INTEGERS (t_INT or ZC) x and y in nf */
541 GEN
nfmuli(GEN nf,GEN x,GEN y)542 nfmuli(GEN nf, GEN x, GEN y)
543 {
544 if (typ(x) == t_INT) return (typ(y) == t_COL)? ZC_Z_mul(y, x): mulii(x,y);
545 if (typ(y) == t_INT) return ZC_Z_mul(x, y);
546 return nfmuli_ZC(nf, x, y);
547 }
548 GEN
nfsqri(GEN nf,GEN x)549 nfsqri(GEN nf, GEN x)
550 { return (typ(x) == t_INT)? sqri(x): nfsqri_ZC(nf, x); }
551
552 /* both x and y are RgV */
553 GEN
tablemul(GEN TAB,GEN x,GEN y)554 tablemul(GEN TAB, GEN x, GEN y)
555 {
556 long i, j, k, N;
557 GEN s, v;
558 if (typ(x) != t_COL) return gmul(x, y);
559 if (typ(y) != t_COL) return gmul(y, x);
560 N = lg(x)-1;
561 v = cgetg(N+1,t_COL);
562 for (k=1; k<=N; k++)
563 {
564 pari_sp av = avma;
565 GEN TABi = TAB;
566 if (k == 1)
567 s = gmul(gel(x,1),gel(y,1));
568 else
569 s = gadd(gmul(gel(x,1),gel(y,k)),
570 gmul(gel(x,k),gel(y,1)));
571 for (i=2; i<=N; i++)
572 {
573 GEN t, xi = gel(x,i);
574 TABi += N;
575 if (gequal0(xi)) continue;
576
577 t = NULL;
578 for (j=2; j<=N; j++)
579 {
580 GEN p1, c = gcoeff(TABi, k, j); /* m^{i,j}_k */
581 if (gequal0(c)) continue;
582 p1 = gmul(c, gel(y,j));
583 t = t? gadd(t, p1): p1;
584 }
585 if (t) s = gadd(s, gmul(xi, t));
586 }
587 gel(v,k) = gerepileupto(av,s);
588 }
589 return v;
590 }
591 GEN
tablesqr(GEN TAB,GEN x)592 tablesqr(GEN TAB, GEN x)
593 {
594 long i, j, k, N;
595 GEN s, v;
596
597 if (typ(x) != t_COL) return gsqr(x);
598 N = lg(x)-1;
599 v = cgetg(N+1,t_COL);
600
601 for (k=1; k<=N; k++)
602 {
603 pari_sp av = avma;
604 GEN TABi = TAB;
605 if (k == 1)
606 s = gsqr(gel(x,1));
607 else
608 s = gmul2n(gmul(gel(x,1),gel(x,k)), 1);
609 for (i=2; i<=N; i++)
610 {
611 GEN p1, c, t, xi = gel(x,i);
612 TABi += N;
613 if (gequal0(xi)) continue;
614
615 c = gcoeff(TABi, k, i);
616 t = !gequal0(c)? gmul(c,xi): NULL;
617 for (j=i+1; j<=N; j++)
618 {
619 c = gcoeff(TABi, k, j);
620 if (gequal0(c)) continue;
621 p1 = gmul(gmul2n(c,1), gel(x,j));
622 t = t? gadd(t, p1): p1;
623 }
624 if (t) s = gadd(s, gmul(xi, t));
625 }
626 gel(v,k) = gerepileupto(av,s);
627 }
628 return v;
629 }
630
631 static GEN
_mul(void * data,GEN x,GEN y)632 _mul(void *data, GEN x, GEN y) { return nfmuli((GEN)data,x,y); }
633 static GEN
_sqr(void * data,GEN x)634 _sqr(void *data, GEN x) { return nfsqri((GEN)data,x); }
635
636 /* Compute z^n in nf, left-shift binary powering */
637 GEN
nfpow(GEN nf,GEN z,GEN n)638 nfpow(GEN nf, GEN z, GEN n)
639 {
640 pari_sp av = avma;
641 long s;
642 GEN x, cx;
643
644 if (typ(n)!=t_INT) pari_err_TYPE("nfpow",n);
645 nf = checknf(nf);
646 s = signe(n); if (!s) return gen_1;
647 if (is_famat(z)) return famat_pow(z, n);
648 x = nf_to_scalar_or_basis(nf, z);
649 if (typ(x) != t_COL) return powgi(x,n);
650 if (s < 0)
651 { /* simplified nfinv */
652 GEN d;
653 x = Q_remove_denom(x, &d);
654 x = zk_inv(nf, x);
655 x = primitive_part(x, &cx);
656 cx = mul_content(cx, d);
657 n = negi(n);
658 }
659 else
660 x = primitive_part(x, &cx);
661 x = gen_pow_i(x, n, (void*)nf, _sqr, _mul);
662 if (cx)
663 x = gerepileupto(av, gmul(x, powgi(cx, n)));
664 else
665 x = gerepilecopy(av, x);
666 return x;
667 }
668 /* Compute z^n in nf, left-shift binary powering */
669 GEN
nfpow_u(GEN nf,GEN z,ulong n)670 nfpow_u(GEN nf, GEN z, ulong n)
671 {
672 pari_sp av = avma;
673 GEN x, cx;
674
675 nf = checknf(nf);
676 if (!n) return gen_1;
677 x = nf_to_scalar_or_basis(nf, z);
678 if (typ(x) != t_COL) return gpowgs(x,n);
679 x = primitive_part(x, &cx);
680 x = gen_powu_i(x, n, (void*)nf, _sqr, _mul);
681 if (cx)
682 {
683 x = gmul(x, powgi(cx, utoipos(n)));
684 return gerepileupto(av,x);
685 }
686 return gerepilecopy(av, x);
687 }
688
689 static GEN
_nf_red(void * E,GEN x)690 _nf_red(void *E, GEN x) { (void)E; return x; }
691
692 static GEN
_nf_add(void * E,GEN x,GEN y)693 _nf_add(void *E, GEN x, GEN y) { return nfadd((GEN)E,x,y); }
694
695 static GEN
_nf_neg(void * E,GEN x)696 _nf_neg(void *E, GEN x) { (void)E; return gneg(x); }
697
698 static GEN
_nf_mul(void * E,GEN x,GEN y)699 _nf_mul(void *E, GEN x, GEN y) { return nfmul((GEN)E,x,y); }
700
701 static GEN
_nf_inv(void * E,GEN x)702 _nf_inv(void *E, GEN x) { return nfinv((GEN)E,x); }
703
704 static GEN
_nf_s(void * E,long x)705 _nf_s(void *E, long x) { (void)E; return stoi(x); }
706
707 static const struct bb_field nf_field={_nf_red,_nf_add,_nf_mul,_nf_neg,
708 _nf_inv,&gequal0,_nf_s };
709
get_nf_field(void ** E,GEN nf)710 const struct bb_field *get_nf_field(void **E, GEN nf)
711 { *E = (void*)nf; return &nf_field; }
712
713 GEN
nfM_det(GEN nf,GEN M)714 nfM_det(GEN nf, GEN M)
715 {
716 void *E;
717 const struct bb_field *S = get_nf_field(&E, nf);
718 return gen_det(M, E, S);
719 }
720 GEN
nfM_inv(GEN nf,GEN M)721 nfM_inv(GEN nf, GEN M)
722 {
723 void *E;
724 const struct bb_field *S = get_nf_field(&E, nf);
725 return gen_Gauss(M, matid(lg(M)-1), E, S);
726 }
727 GEN
nfM_mul(GEN nf,GEN A,GEN B)728 nfM_mul(GEN nf, GEN A, GEN B)
729 {
730 void *E;
731 const struct bb_field *S = get_nf_field(&E, nf);
732 return gen_matmul(A, B, E, S);
733 }
734 GEN
nfM_nfC_mul(GEN nf,GEN A,GEN B)735 nfM_nfC_mul(GEN nf, GEN A, GEN B)
736 {
737 void *E;
738 const struct bb_field *S = get_nf_field(&E, nf);
739 return gen_matcolmul(A, B, E, S);
740 }
741
742 /* valuation of integral x (ZV), with resp. to prime ideal pr */
743 long
ZC_nfvalrem(GEN x,GEN pr,GEN * newx)744 ZC_nfvalrem(GEN x, GEN pr, GEN *newx)
745 {
746 long i, v, l;
747 GEN r, y, p = pr_get_p(pr), mul = pr_get_tau(pr);
748
749 /* p inert */
750 if (typ(mul) == t_INT) return newx? ZV_pvalrem(x, p, newx):ZV_pval(x, p);
751 y = cgetg_copy(x, &l); /* will hold the new x */
752 x = leafcopy(x);
753 for(v=0;; v++)
754 {
755 for (i=1; i<l; i++)
756 { /* is (x.b)[i] divisible by p ? */
757 gel(y,i) = dvmdii(ZMrow_ZC_mul(mul,x,i),p,&r);
758 if (r != gen_0) { if (newx) *newx = x; return v; }
759 }
760 swap(x, y);
761 }
762 }
763 long
ZC_nfval(GEN x,GEN P)764 ZC_nfval(GEN x, GEN P)
765 { return ZC_nfvalrem(x, P, NULL); }
766
767 /* v_P(x) != 0, x a ZV. Simpler version of ZC_nfvalrem */
768 int
ZC_prdvd(GEN x,GEN P)769 ZC_prdvd(GEN x, GEN P)
770 {
771 pari_sp av = avma;
772 long i, l;
773 GEN p = pr_get_p(P), mul = pr_get_tau(P);
774 if (typ(mul) == t_INT) return ZV_Z_dvd(x, p);
775 l = lg(x);
776 for (i=1; i<l; i++)
777 if (!dvdii(ZMrow_ZC_mul(mul,x,i), p)) return gc_bool(av,0);
778 return gc_bool(av,1);
779 }
780
781 int
pr_equal(GEN P,GEN Q)782 pr_equal(GEN P, GEN Q)
783 {
784 GEN gQ, p = pr_get_p(P);
785 long e = pr_get_e(P), f = pr_get_f(P), n;
786 if (!equalii(p, pr_get_p(Q)) || e != pr_get_e(Q) || f != pr_get_f(Q))
787 return 0;
788 gQ = pr_get_gen(Q); n = lg(gQ)-1;
789 if (2*e*f > n) return 1; /* room for only one such pr */
790 return ZV_equal(pr_get_gen(P), gQ) || ZC_prdvd(gQ, P);
791 }
792
793 GEN
famat_nfvalrem(GEN nf,GEN x,GEN pr,GEN * py)794 famat_nfvalrem(GEN nf, GEN x, GEN pr, GEN *py)
795 {
796 pari_sp av = avma;
797 GEN P = gel(x,1), E = gel(x,2), V = gen_0, y = NULL;
798 long l = lg(P), simplify = 0, i;
799 if (py) { *py = gen_1; y = cgetg(l, t_COL); }
800
801 for (i = 1; i < l; i++)
802 {
803 GEN e = gel(E,i);
804 long v;
805 if (!signe(e))
806 {
807 if (py) gel(y,i) = gen_1;
808 simplify = 1; continue;
809 }
810 v = nfvalrem(nf, gel(P,i), pr, py? &gel(y,i): NULL);
811 if (v == LONG_MAX) { set_avma(av); if (py) *py = gen_0; return mkoo(); }
812 V = addmulii(V, stoi(v), e);
813 }
814 if (!py) V = gerepileuptoint(av, V);
815 else
816 {
817 y = mkmat2(y, gel(x,2));
818 if (simplify) y = famat_remove_trivial(y);
819 gerepileall(av, 2, &V, &y); *py = y;
820 }
821 return V;
822 }
823 long
nfval(GEN nf,GEN x,GEN pr)824 nfval(GEN nf, GEN x, GEN pr)
825 {
826 pari_sp av = avma;
827 long w, e;
828 GEN cx, p;
829
830 if (gequal0(x)) return LONG_MAX;
831 nf = checknf(nf);
832 checkprid(pr);
833 p = pr_get_p(pr);
834 e = pr_get_e(pr);
835 x = nf_to_scalar_or_basis(nf, x);
836 if (typ(x) != t_COL) return e*Q_pval(x,p);
837 x = Q_primitive_part(x, &cx);
838 w = ZC_nfval(x,pr);
839 if (cx) w += e*Q_pval(cx,p);
840 return gc_long(av,w);
841 }
842
843 /* want to write p^v = uniformizer^(e*v) * z^v, z coprime to pr */
844 /* z := tau^e / p^(e-1), algebraic integer coprime to pr; return z^v */
845 static GEN
powp(GEN nf,GEN pr,long v)846 powp(GEN nf, GEN pr, long v)
847 {
848 GEN b, z;
849 long e;
850 if (!v) return gen_1;
851 b = pr_get_tau(pr);
852 if (typ(b) == t_INT) return gen_1;
853 e = pr_get_e(pr);
854 z = gel(b,1);
855 if (e != 1) z = gdiv(nfpow_u(nf, z, e), powiu(pr_get_p(pr),e-1));
856 if (v < 0) { v = -v; z = nfinv(nf, z); }
857 if (v != 1) z = nfpow_u(nf, z, v);
858 return z;
859 }
860 long
nfvalrem(GEN nf,GEN x,GEN pr,GEN * py)861 nfvalrem(GEN nf, GEN x, GEN pr, GEN *py)
862 {
863 pari_sp av = avma;
864 long w, e;
865 GEN cx, p, t;
866
867 if (!py) return nfval(nf,x,pr);
868 if (gequal0(x)) { *py = gen_0; return LONG_MAX; }
869 nf = checknf(nf);
870 checkprid(pr);
871 p = pr_get_p(pr);
872 e = pr_get_e(pr);
873 x = nf_to_scalar_or_basis(nf, x);
874 if (typ(x) != t_COL) {
875 w = Q_pvalrem(x,p, py);
876 if (!w) { *py = gerepilecopy(av, x); return 0; }
877 *py = gerepileupto(av, gmul(powp(nf, pr, w), *py));
878 return e*w;
879 }
880 x = Q_primitive_part(x, &cx);
881 w = ZC_nfvalrem(x,pr, py);
882 if (cx)
883 {
884 long v = Q_pvalrem(cx,p, &t);
885 *py = nfmul(nf, *py, gmul(powp(nf,pr,v), t));
886 *py = gerepileupto(av, *py);
887 w += e*v;
888 }
889 else
890 *py = gerepilecopy(av, *py);
891 return w;
892 }
893 GEN
gpnfvalrem(GEN nf,GEN x,GEN pr,GEN * py)894 gpnfvalrem(GEN nf, GEN x, GEN pr, GEN *py)
895 {
896 long v;
897 if (is_famat(x)) return famat_nfvalrem(nf, x, pr, py);
898 v = nfvalrem(nf,x,pr,py);
899 return v == LONG_MAX? mkoo(): stoi(v);
900 }
901
902 /* true nf */
903 GEN
coltoalg(GEN nf,GEN x)904 coltoalg(GEN nf, GEN x)
905 {
906 return mkpolmod( nf_to_scalar_or_alg(nf, x), nf_get_pol(nf) );
907 }
908
909 GEN
basistoalg(GEN nf,GEN x)910 basistoalg(GEN nf, GEN x)
911 {
912 GEN T;
913
914 nf = checknf(nf);
915 switch(typ(x))
916 {
917 case t_COL: {
918 pari_sp av = avma;
919 return gerepilecopy(av, coltoalg(nf, x));
920 }
921 case t_POLMOD:
922 T = nf_get_pol(nf);
923 if (!RgX_equal_var(T,gel(x,1)))
924 pari_err_MODULUS("basistoalg", T,gel(x,1));
925 return gcopy(x);
926 case t_POL:
927 T = nf_get_pol(nf);
928 if (varn(T) != varn(x)) pari_err_VAR("basistoalg",x,T);
929 retmkpolmod(RgX_rem(x, T), ZX_copy(T));
930 case t_INT:
931 case t_FRAC:
932 T = nf_get_pol(nf);
933 retmkpolmod(gcopy(x), ZX_copy(T));
934 default:
935 pari_err_TYPE("basistoalg",x);
936 return NULL; /* LCOV_EXCL_LINE */
937 }
938 }
939
940 /* true nf, x a t_POL */
941 static GEN
pol_to_scalar_or_basis(GEN nf,GEN x)942 pol_to_scalar_or_basis(GEN nf, GEN x)
943 {
944 GEN T = nf_get_pol(nf);
945 long l = lg(x);
946 if (varn(x) != varn(T)) pari_err_VAR("nf_to_scalar_or_basis", x,T);
947 if (l >= lg(T)) { x = RgX_rem(x, T); l = lg(x); }
948 if (l == 2) return gen_0;
949 if (l == 3)
950 {
951 x = gel(x,2);
952 if (!is_rational_t(typ(x))) pari_err_TYPE("nf_to_scalar_or_basis",x);
953 return x;
954 }
955 return poltobasis(nf,x);
956 }
957 /* Assume nf is a genuine nf. */
958 GEN
nf_to_scalar_or_basis(GEN nf,GEN x)959 nf_to_scalar_or_basis(GEN nf, GEN x)
960 {
961 switch(typ(x))
962 {
963 case t_INT: case t_FRAC:
964 return x;
965 case t_POLMOD:
966 x = checknfelt_mod(nf,x,"nf_to_scalar_or_basis");
967 switch(typ(x))
968 {
969 case t_INT: case t_FRAC: return x;
970 case t_POL: return pol_to_scalar_or_basis(nf,x);
971 }
972 break;
973 case t_POL: return pol_to_scalar_or_basis(nf,x);
974 case t_COL:
975 if (lg(x)-1 != nf_get_degree(nf)) break;
976 return QV_isscalar(x)? gel(x,1): x;
977 }
978 pari_err_TYPE("nf_to_scalar_or_basis",x);
979 return NULL; /* LCOV_EXCL_LINE */
980 }
981 /* Let x be a polynomial with coefficients in Q or nf. Return the same
982 * polynomial with coefficients expressed as vectors (on the integral basis).
983 * No consistency checks, not memory-clean. */
984 GEN
RgX_to_nfX(GEN nf,GEN x)985 RgX_to_nfX(GEN nf, GEN x)
986 {
987 long i, l;
988 GEN y = cgetg_copy(x, &l); y[1] = x[1];
989 for (i=2; i<l; i++) gel(y,i) = nf_to_scalar_or_basis(nf, gel(x,i));
990 return y;
991 }
992
993 /* Assume nf is a genuine nf. */
994 GEN
nf_to_scalar_or_alg(GEN nf,GEN x)995 nf_to_scalar_or_alg(GEN nf, GEN x)
996 {
997 switch(typ(x))
998 {
999 case t_INT: case t_FRAC:
1000 return x;
1001 case t_POLMOD:
1002 x = checknfelt_mod(nf,x,"nf_to_scalar_or_alg");
1003 if (typ(x) != t_POL) return x;
1004 /* fall through */
1005 case t_POL:
1006 {
1007 GEN T = nf_get_pol(nf);
1008 long l = lg(x);
1009 if (varn(x) != varn(T)) pari_err_VAR("nf_to_scalar_or_alg", x,T);
1010 if (l >= lg(T)) { x = RgX_rem(x, T); l = lg(x); }
1011 if (l == 2) return gen_0;
1012 if (l == 3) return gel(x,2);
1013 return x;
1014 }
1015 case t_COL:
1016 {
1017 GEN dx;
1018 if (lg(x)-1 != nf_get_degree(nf)) break;
1019 if (QV_isscalar(x)) return gel(x,1);
1020 x = Q_remove_denom(x, &dx);
1021 x = RgV_RgC_mul(nf_get_zkprimpart(nf), x);
1022 dx = mul_denom(dx, nf_get_zkden(nf));
1023 return gdiv(x,dx);
1024 }
1025 }
1026 pari_err_TYPE("nf_to_scalar_or_alg",x);
1027 return NULL; /* LCOV_EXCL_LINE */
1028 }
1029
1030 /* gmul(A, RgX_to_RgC(x)), A t_MAT of compatible dimensions */
1031 GEN
RgM_RgX_mul(GEN A,GEN x)1032 RgM_RgX_mul(GEN A, GEN x)
1033 {
1034 long i, l = lg(x)-1;
1035 GEN z;
1036 if (l == 1) return zerocol(nbrows(A));
1037 z = gmul(gel(x,2), gel(A,1));
1038 for (i = 2; i < l; i++)
1039 if (!gequal0(gel(x,i+1))) z = gadd(z, gmul(gel(x,i+1), gel(A,i)));
1040 return z;
1041 }
1042 GEN
ZM_ZX_mul(GEN A,GEN x)1043 ZM_ZX_mul(GEN A, GEN x)
1044 {
1045 long i, l = lg(x)-1;
1046 GEN z;
1047 if (l == 1) return zerocol(nbrows(A));
1048 z = ZC_Z_mul(gel(A,1), gel(x,2));
1049 for (i = 2; i < l ; i++)
1050 if (signe(gel(x,i+1))) z = ZC_add(z, ZC_Z_mul(gel(A,i), gel(x,i+1)));
1051 return z;
1052 }
1053 /* x a t_POL, nf a genuine nf. No garbage collecting. No check. */
1054 GEN
poltobasis(GEN nf,GEN x)1055 poltobasis(GEN nf, GEN x)
1056 {
1057 GEN d, T = nf_get_pol(nf);
1058 if (varn(x) != varn(T)) pari_err_VAR( "poltobasis", x,T);
1059 if (degpol(x) >= degpol(T)) x = RgX_rem(x,T);
1060 x = Q_remove_denom(x, &d);
1061 if (!RgX_is_ZX(x)) pari_err_TYPE("poltobasis",x);
1062 x = ZM_ZX_mul(nf_get_invzk(nf), x);
1063 if (d) x = RgC_Rg_div(x, d);
1064 return x;
1065 }
1066
1067 GEN
algtobasis(GEN nf,GEN x)1068 algtobasis(GEN nf, GEN x)
1069 {
1070 pari_sp av;
1071
1072 nf = checknf(nf);
1073 switch(typ(x))
1074 {
1075 case t_POLMOD:
1076 if (!RgX_equal_var(nf_get_pol(nf),gel(x,1)))
1077 pari_err_MODULUS("algtobasis", nf_get_pol(nf),gel(x,1));
1078 x = gel(x,2);
1079 switch(typ(x))
1080 {
1081 case t_INT:
1082 case t_FRAC: return scalarcol(x, nf_get_degree(nf));
1083 case t_POL:
1084 av = avma;
1085 return gerepileupto(av,poltobasis(nf,x));
1086 }
1087 break;
1088
1089 case t_POL:
1090 av = avma;
1091 return gerepileupto(av,poltobasis(nf,x));
1092
1093 case t_COL:
1094 if (!RgV_is_QV(x)) pari_err_TYPE("nfalgtobasis",x);
1095 if (lg(x)-1 != nf_get_degree(nf)) pari_err_DIM("nfalgtobasis");
1096 return gcopy(x);
1097
1098 case t_INT:
1099 case t_FRAC: return scalarcol(x, nf_get_degree(nf));
1100 }
1101 pari_err_TYPE("algtobasis",x);
1102 return NULL; /* LCOV_EXCL_LINE */
1103 }
1104
1105 GEN
rnfbasistoalg(GEN rnf,GEN x)1106 rnfbasistoalg(GEN rnf,GEN x)
1107 {
1108 const char *f = "rnfbasistoalg";
1109 long lx, i;
1110 pari_sp av = avma;
1111 GEN z, nf, R, T;
1112
1113 checkrnf(rnf);
1114 nf = rnf_get_nf(rnf);
1115 T = nf_get_pol(nf);
1116 R = QXQX_to_mod_shallow(rnf_get_pol(rnf), T);
1117 switch(typ(x))
1118 {
1119 case t_COL:
1120 z = cgetg_copy(x, &lx);
1121 for (i=1; i<lx; i++)
1122 {
1123 GEN c = nf_to_scalar_or_alg(nf, gel(x,i));
1124 if (typ(c) == t_POL) c = mkpolmod(c,T);
1125 gel(z,i) = c;
1126 }
1127 z = RgV_RgC_mul(gel(rnf_get_zk(rnf),1), z);
1128 return gerepileupto(av, gmodulo(z,R));
1129
1130 case t_POLMOD:
1131 x = polmod_nffix(f, rnf, x, 0);
1132 if (typ(x) != t_POL) break;
1133 retmkpolmod(RgX_copy(x), RgX_copy(R));
1134 case t_POL:
1135 if (varn(x) == varn(T)) { RgX_check_QX(x,f); x = gmodulo(x,T); break; }
1136 if (varn(x) == varn(R))
1137 {
1138 x = RgX_nffix(f,nf_get_pol(nf),x,0);
1139 return gmodulo(x, R);
1140 }
1141 pari_err_VAR(f, x,R);
1142 }
1143 retmkpolmod(scalarpol(x, varn(R)), RgX_copy(R));
1144 }
1145
1146 GEN
matbasistoalg(GEN nf,GEN x)1147 matbasistoalg(GEN nf,GEN x)
1148 {
1149 long i, j, li, lx;
1150 GEN z = cgetg_copy(x, &lx);
1151
1152 if (lx == 1) return z;
1153 switch(typ(x))
1154 {
1155 case t_VEC: case t_COL:
1156 for (i=1; i<lx; i++) gel(z,i) = basistoalg(nf, gel(x,i));
1157 return z;
1158 case t_MAT: break;
1159 default: pari_err_TYPE("matbasistoalg",x);
1160 }
1161 li = lgcols(x);
1162 for (j=1; j<lx; j++)
1163 {
1164 GEN c = cgetg(li,t_COL), xj = gel(x,j);
1165 gel(z,j) = c;
1166 for (i=1; i<li; i++) gel(c,i) = basistoalg(nf, gel(xj,i));
1167 }
1168 return z;
1169 }
1170
1171 GEN
matalgtobasis(GEN nf,GEN x)1172 matalgtobasis(GEN nf,GEN x)
1173 {
1174 long i, j, li, lx;
1175 GEN z = cgetg_copy(x, &lx);
1176
1177 if (lx == 1) return z;
1178 switch(typ(x))
1179 {
1180 case t_VEC: case t_COL:
1181 for (i=1; i<lx; i++) gel(z,i) = algtobasis(nf, gel(x,i));
1182 return z;
1183 case t_MAT: break;
1184 default: pari_err_TYPE("matalgtobasis",x);
1185 }
1186 li = lgcols(x);
1187 for (j=1; j<lx; j++)
1188 {
1189 GEN c = cgetg(li,t_COL), xj = gel(x,j);
1190 gel(z,j) = c;
1191 for (i=1; i<li; i++) gel(c,i) = algtobasis(nf, gel(xj,i));
1192 }
1193 return z;
1194 }
1195 GEN
RgM_to_nfM(GEN nf,GEN x)1196 RgM_to_nfM(GEN nf,GEN x)
1197 {
1198 long i, j, li, lx;
1199 GEN z = cgetg_copy(x, &lx);
1200
1201 if (lx == 1) return z;
1202 li = lgcols(x);
1203 for (j=1; j<lx; j++)
1204 {
1205 GEN c = cgetg(li,t_COL), xj = gel(x,j);
1206 gel(z,j) = c;
1207 for (i=1; i<li; i++) gel(c,i) = nf_to_scalar_or_basis(nf, gel(xj,i));
1208 }
1209 return z;
1210 }
1211 GEN
RgC_to_nfC(GEN nf,GEN x)1212 RgC_to_nfC(GEN nf, GEN x)
1213 { pari_APPLY_type(t_COL, nf_to_scalar_or_basis(nf, gel(x,i))) }
1214
1215 /* x a t_POLMOD, supposedly in rnf = K[z]/(T), K = Q[y]/(Tnf) */
1216 GEN
polmod_nffix(const char * f,GEN rnf,GEN x,int lift)1217 polmod_nffix(const char *f, GEN rnf, GEN x, int lift)
1218 { return polmod_nffix2(f, rnf_get_nfpol(rnf), rnf_get_pol(rnf), x,lift); }
1219 GEN
polmod_nffix2(const char * f,GEN T,GEN R,GEN x,int lift)1220 polmod_nffix2(const char *f, GEN T, GEN R, GEN x, int lift)
1221 {
1222 if (RgX_equal_var(gel(x,1), R))
1223 {
1224 x = gel(x,2);
1225 if (typ(x) == t_POL && varn(x) == varn(R))
1226 {
1227 x = RgX_nffix(f, T, x, lift);
1228 switch(lg(x))
1229 {
1230 case 2: return gen_0;
1231 case 3: return gel(x,2);
1232 }
1233 return x;
1234 }
1235 }
1236 return Rg_nffix(f, T, x, lift);
1237 }
1238 GEN
rnfalgtobasis(GEN rnf,GEN x)1239 rnfalgtobasis(GEN rnf,GEN x)
1240 {
1241 const char *f = "rnfalgtobasis";
1242 pari_sp av = avma;
1243 GEN T, R;
1244
1245 checkrnf(rnf);
1246 R = rnf_get_pol(rnf);
1247 T = rnf_get_nfpol(rnf);
1248 switch(typ(x))
1249 {
1250 case t_COL:
1251 if (lg(x)-1 != rnf_get_degree(rnf)) pari_err_DIM(f);
1252 x = RgV_nffix(f, T, x, 0);
1253 return gerepilecopy(av, x);
1254
1255 case t_POLMOD:
1256 x = polmod_nffix(f, rnf, x, 0);
1257 if (typ(x) != t_POL) break;
1258 return gerepileupto(av, RgM_RgX_mul(rnf_get_invzk(rnf), x));
1259 case t_POL:
1260 if (varn(x) == varn(T))
1261 {
1262 RgX_check_QX(x,f);
1263 if (degpol(x) >= degpol(T)) x = RgX_rem(x,T);
1264 x = mkpolmod(x,T); break;
1265 }
1266 x = RgX_nffix(f, T, x, 0);
1267 if (degpol(x) >= degpol(R)) x = RgX_rem(x, R);
1268 return gerepileupto(av, RgM_RgX_mul(rnf_get_invzk(rnf), x));
1269 }
1270 return gerepileupto(av, scalarcol(x, rnf_get_degree(rnf)));
1271 }
1272
1273 /* Given a and b in nf, gives an algebraic integer y in nf such that a-b.y
1274 * is "small" */
1275 GEN
nfdiveuc(GEN nf,GEN a,GEN b)1276 nfdiveuc(GEN nf, GEN a, GEN b)
1277 {
1278 pari_sp av = avma;
1279 a = nfdiv(nf,a,b);
1280 return gerepileupto(av, ground(a));
1281 }
1282
1283 /* Given a and b in nf, gives a "small" algebraic integer r in nf
1284 * of the form a-b.y */
1285 GEN
nfmod(GEN nf,GEN a,GEN b)1286 nfmod(GEN nf, GEN a, GEN b)
1287 {
1288 pari_sp av = avma;
1289 GEN p1 = gneg_i(nfmul(nf,b,ground(nfdiv(nf,a,b))));
1290 return gerepileupto(av, nfadd(nf,a,p1));
1291 }
1292
1293 /* Given a and b in nf, gives a two-component vector [y,r] in nf such
1294 * that r=a-b.y is "small". */
1295 GEN
nfdivrem(GEN nf,GEN a,GEN b)1296 nfdivrem(GEN nf, GEN a, GEN b)
1297 {
1298 pari_sp av = avma;
1299 GEN p1,z, y = ground(nfdiv(nf,a,b));
1300
1301 p1 = gneg_i(nfmul(nf,b,y));
1302 z = cgetg(3,t_VEC);
1303 gel(z,1) = gcopy(y);
1304 gel(z,2) = nfadd(nf,a,p1); return gerepileupto(av, z);
1305 }
1306
1307 /*************************************************************************/
1308 /** **/
1309 /** LOGARITHMIC EMBEDDINGS **/
1310 /** **/
1311 /*************************************************************************/
1312
1313 static int
low_prec(GEN x)1314 low_prec(GEN x)
1315 {
1316 switch(typ(x))
1317 {
1318 case t_INT: return !signe(x);
1319 case t_REAL: return !signe(x) || realprec(x) <= DEFAULTPREC;
1320 default: return 0;
1321 }
1322 }
1323
1324 static GEN
cxlog_1(GEN nf)1325 cxlog_1(GEN nf) { return zerocol(lg(nf_get_roots(nf))-1); }
1326 static GEN
cxlog_m1(GEN nf,long prec)1327 cxlog_m1(GEN nf, long prec)
1328 {
1329 long i, l = lg(nf_get_roots(nf)), r1 = nf_get_r1(nf);
1330 GEN v = cgetg(l, t_COL), p, P;
1331 p = mppi(prec); P = mkcomplex(gen_0, p);
1332 for (i = 1; i <= r1; i++) gel(v,i) = P; /* IPi*/
1333 if (i < l) P = gmul2n(P,1);
1334 for ( ; i < l; i++) gel(v,i) = P; /* 2IPi */
1335 return v;
1336 }
1337 static GEN
famat_cxlog(GEN nf,GEN fa,long prec)1338 famat_cxlog(GEN nf, GEN fa, long prec)
1339 {
1340 GEN g, e, y = NULL;
1341 long i, l;
1342
1343 if (typ(fa) != t_MAT) pari_err_TYPE("famat_cxlog",fa);
1344 if (lg(fa) == 1) return cxlog_1(nf);
1345 g = gel(fa,1);
1346 e = gel(fa,2); l = lg(e);
1347 for (i = 1; i < l; i++)
1348 {
1349 GEN t, x = nf_to_scalar_or_basis(nf, gel(g,i));
1350 /* multiplicative arch would be better (save logs), but exponents overflow
1351 * [ could keep track of expo separately, but not worth it ] */
1352 t = nf_cxlog(nf,x,prec); if (!t) return NULL;
1353 if (gel(t,1) == gen_0) continue; /* positive rational */
1354 t = RgC_Rg_mul(t, gel(e,i));
1355 y = y? RgV_add(y,t): t;
1356 }
1357 return y ? y: cxlog_1(nf);
1358 }
1359 /* Archimedean components: [e_i Log( sigma_i(X) )], where X = primpart(x),
1360 * and e_i = 1 (resp 2.) for i <= R1 (resp. > R1) */
1361 GEN
nf_cxlog(GEN nf,GEN x,long prec)1362 nf_cxlog(GEN nf, GEN x, long prec)
1363 {
1364 long i, l, r1;
1365 GEN v;
1366 if (typ(x) == t_MAT) return famat_cxlog(nf,x,prec);
1367 x = nf_to_scalar_or_basis(nf,x);
1368 if (typ(x) != t_COL) return gsigne(x) > 0? cxlog_1(nf): cxlog_m1(nf, prec);
1369 x = RgM_RgC_mul(nf_get_M(nf), Q_primpart(x));
1370 l = lg(x); r1 = nf_get_r1(nf);
1371 for (i = 1; i <= r1; i++)
1372 if (low_prec(gel(x,i))) return NULL;
1373 for ( ; i < l; i++)
1374 if (low_prec(gnorm(gel(x,i)))) return NULL;
1375 v = cgetg(l,t_COL);
1376 for (i = 1; i <= r1; i++) gel(v,i) = glog(gel(x,i),prec);
1377 for ( ; i < l; i++) gel(v,i) = gmul2n(glog(gel(x,i),prec),1);
1378 return v;
1379 }
1380 GEN
nfV_cxlog(GEN nf,GEN x,long prec)1381 nfV_cxlog(GEN nf, GEN x, long prec)
1382 {
1383 long i, l;
1384 GEN v = cgetg_copy(x, &l);
1385 for (i = 1; i < l; i++)
1386 if (!(gel(v,i) = nf_cxlog(nf, gel(x,i), prec))) return NULL;
1387 return v;
1388 }
1389
1390 static GEN
scalar_logembed(GEN nf,GEN u,GEN * emb)1391 scalar_logembed(GEN nf, GEN u, GEN *emb)
1392 {
1393 GEN v, logu;
1394 long i, s = signe(u), RU = lg(nf_get_roots(nf))-1, R1 = nf_get_r1(nf);
1395
1396 if (!s) pari_err_DOMAIN("nflogembed","argument","=",gen_0,u);
1397 v = cgetg(RU+1, t_COL); logu = logr_abs(u);
1398 for (i = 1; i <= R1; i++) gel(v,i) = logu;
1399 if (i <= RU)
1400 {
1401 GEN logu2 = shiftr(logu,1);
1402 for ( ; i <= RU; i++) gel(v,i) = logu2;
1403 }
1404 if (emb) *emb = const_col(RU, u);
1405 return v;
1406 }
1407
1408 static GEN
famat_logembed(GEN nf,GEN x,GEN * emb,long prec)1409 famat_logembed(GEN nf,GEN x,GEN *emb,long prec)
1410 {
1411 GEN A, M, T, a, t, g = gel(x,1), e = gel(x,2);
1412 long i, l = lg(e);
1413
1414 if (l == 1) return scalar_logembed(nf, real_1(prec), emb);
1415 A = NULL; T = emb? cgetg(l, t_COL): NULL;
1416 if (emb) *emb = M = mkmat2(T, e);
1417 for (i = 1; i < l; i++)
1418 {
1419 a = nflogembed(nf, gel(g,i), &t, prec);
1420 if (!a) return NULL;
1421 a = RgC_Rg_mul(a, gel(e,i));
1422 A = A? RgC_add(A, a): a;
1423 if (emb) gel(T,i) = t;
1424 }
1425 return A;
1426 }
1427
1428 /* Get archimedean components: [e_i log( | sigma_i(x) | )], with e_i = 1
1429 * (resp 2.) for i <= R1 (resp. > R1) and set emb to the embeddings of x.
1430 * Return NULL if precision problem */
1431 GEN
nflogembed(GEN nf,GEN x,GEN * emb,long prec)1432 nflogembed(GEN nf, GEN x, GEN *emb, long prec)
1433 {
1434 long i, l, r1;
1435 GEN v, t;
1436
1437 if (typ(x) == t_MAT) return famat_logembed(nf,x,emb,prec);
1438 x = nf_to_scalar_or_basis(nf,x);
1439 if (typ(x) != t_COL) return scalar_logembed(nf, gtofp(x,prec), emb);
1440 x = RgM_RgC_mul(nf_get_M(nf), x);
1441 l = lg(x); r1 = nf_get_r1(nf); v = cgetg(l,t_COL);
1442 for (i = 1; i <= r1; i++)
1443 {
1444 t = gabs(gel(x,i),prec); if (low_prec(t)) return NULL;
1445 gel(v,i) = glog(t,prec);
1446 }
1447 for ( ; i < l; i++)
1448 {
1449 t = gnorm(gel(x,i)); if (low_prec(t)) return NULL;
1450 gel(v,i) = glog(t,prec);
1451 }
1452 if (emb) *emb = x;
1453 return v;
1454 }
1455
1456 /*************************************************************************/
1457 /** **/
1458 /** REAL EMBEDDINGS **/
1459 /** **/
1460 /*************************************************************************/
1461 static GEN
sarch_get_cyc(GEN sarch)1462 sarch_get_cyc(GEN sarch) { return gel(sarch,1); }
1463 static GEN
sarch_get_archp(GEN sarch)1464 sarch_get_archp(GEN sarch) { return gel(sarch,2); }
1465 static GEN
sarch_get_MI(GEN sarch)1466 sarch_get_MI(GEN sarch) { return gel(sarch,3); }
1467 static GEN
sarch_get_lambda(GEN sarch)1468 sarch_get_lambda(GEN sarch) { return gel(sarch,4); }
1469 static GEN
sarch_get_F(GEN sarch)1470 sarch_get_F(GEN sarch) { return gel(sarch,5); }
1471
1472 /* x not a scalar, true nf, return number of positive roots of char_x */
1473 static long
num_positive(GEN nf,GEN x)1474 num_positive(GEN nf, GEN x)
1475 {
1476 GEN T = nf_get_pol(nf), B, charx;
1477 long dnf, vnf, N;
1478 x = nf_to_scalar_or_alg(nf, x); /* not a scalar */
1479 charx = ZXQ_charpoly(x, T, 0);
1480 charx = ZX_radical(charx);
1481 N = degpol(T) / degpol(charx);
1482 /* real places are unramified ? */
1483 if (N == 1 || ZX_sturm(charx) * N == nf_get_r1(nf))
1484 return ZX_sturmpart(charx, mkvec2(gen_0,mkoo())) * N;
1485 /* painful case, multiply by random square until primitive */
1486 dnf = nf_get_degree(nf);
1487 vnf = varn(T);
1488 B = int2n(10);
1489 for(;;)
1490 {
1491 GEN y = RgXQ_sqr(random_FpX(dnf, vnf, B), T);
1492 y = RgXQ_mul(x, y, T);
1493 charx = ZXQ_charpoly(y, T, 0);
1494 if (ZX_is_squarefree(charx))
1495 return ZX_sturmpart(charx, mkvec2(gen_0,mkoo())) * N;
1496 }
1497 }
1498
1499 /* x a QC: return sigma_k(x) where 1 <= k <= r1+r2; correct but inefficient
1500 * if x in Q. M = nf_get_M(nf) */
1501 static GEN
nfembed_i(GEN M,GEN x,long k)1502 nfembed_i(GEN M, GEN x, long k)
1503 {
1504 long i, l = lg(M);
1505 GEN z = gel(x,1);
1506 for (i = 2; i < l; i++) z = gadd(z, gmul(gcoeff(M,k,i), gel(x,i)));
1507 return z;
1508 }
1509 GEN
nfembed(GEN nf,GEN x,long k)1510 nfembed(GEN nf, GEN x, long k)
1511 {
1512 pari_sp av = avma;
1513 nf = checknf(nf);
1514 x = nf_to_scalar_or_basis(nf,x);
1515 if (typ(x) != t_COL) return gerepilecopy(av, x);
1516 return gerepileupto(av, nfembed_i(nf_get_M(nf),x,k));
1517 }
1518
1519 /* x a ZC */
1520 static GEN
zk_embed(GEN M,GEN x,long k)1521 zk_embed(GEN M, GEN x, long k)
1522 {
1523 long i, l = lg(x);
1524 GEN z = gel(x,1); /* times M[k,1], which is 1 */
1525 for (i = 2; i < l; i++) z = mpadd(z, mpmul(gcoeff(M,k,i), gel(x,i)));
1526 return z;
1527 }
1528
1529 /* Given floating point approximation z of sigma_k(x), decide its sign
1530 * [0/+, 1/- and -1 for FAIL] */
1531 static long
eval_sign_embed(GEN z)1532 eval_sign_embed(GEN z)
1533 { /* dubious, fail */
1534 if (typ(z) == t_REAL && realprec(z) <= LOWDEFAULTPREC) return -1;
1535 return (signe(z) < 1)? 1: 0;
1536 }
1537 /* return v such that (-1)^v = sign(sigma_k(x)), x primitive ZC */
1538 static long
eval_sign(GEN M,GEN x,long k)1539 eval_sign(GEN M, GEN x, long k)
1540 { return eval_sign_embed( zk_embed(M, x, k) ); }
1541
1542 /* check that signs[i..#signs] == s; signs = NULL encodes "totally positive" */
1543 static int
oksigns(long l,GEN signs,long i,long s)1544 oksigns(long l, GEN signs, long i, long s)
1545 {
1546 if (!signs) return s == 0;
1547 for (; i < l; i++)
1548 if (signs[i] != s) return 0;
1549 return 1;
1550 }
1551 /* check that signs[i] = s and signs[i+1..#signs] = 1-s */
1552 static int
oksigns2(long l,GEN signs,long i,long s)1553 oksigns2(long l, GEN signs, long i, long s)
1554 {
1555 if (!signs) return s == 0 && i == l-1;
1556 return signs[i] == s && oksigns(l, signs, i+1, 1-s);
1557 }
1558
1559 /* true nf, x a ZC (primitive for efficiency), embx its embeddings or NULL */
1560 static int
nfchecksigns_i(GEN nf,GEN x,GEN embx,GEN signs,GEN archp)1561 nfchecksigns_i(GEN nf, GEN x, GEN embx, GEN signs, GEN archp)
1562 {
1563 long l = lg(archp), i;
1564 GEN M = nf_get_M(nf), sarch = NULL;
1565 long np = -1;
1566 for (i = 1; i < l; i++)
1567 {
1568 long s;
1569 if (embx)
1570 s = eval_sign_embed(gel(embx,i));
1571 else
1572 s = eval_sign(M, x, archp[i]);
1573 /* 0 / + or 1 / -; -1 for FAIL */
1574 if (s < 0) /* failure */
1575 {
1576 long ni, r1 = nf_get_r1(nf);
1577 GEN xi;
1578 if (np < 0)
1579 {
1580 np = num_positive(nf, x);
1581 if (np == 0) return oksigns(l, signs, i, 1);
1582 if (np == r1) return oksigns(l, signs, i, 0);
1583 sarch = nfarchstar(nf, NULL, identity_perm(r1));
1584 }
1585 xi = set_sign_mod_divisor(nf, vecsmall_ei(r1, archp[i]), gen_1, sarch);
1586 xi = Q_primpart(xi);
1587 ni = num_positive(nf, nfmuli(nf,x,xi));
1588 if (ni == 0) return oksigns2(l, signs, i, 0);
1589 if (ni == r1) return oksigns2(l, signs, i, 1);
1590 s = ni < np? 0: 1;
1591 }
1592 if (s != (signs? signs[i]: 0)) return 0;
1593 }
1594 return 1;
1595 }
1596 static void
pl_convert(GEN pl,GEN * psigns,GEN * parchp)1597 pl_convert(GEN pl, GEN *psigns, GEN *parchp)
1598 {
1599 long i, j, l = lg(pl);
1600 GEN signs = cgetg(l, t_VECSMALL);
1601 GEN archp = cgetg(l, t_VECSMALL);
1602 for (i = j = 1; i < l; i++)
1603 {
1604 if (!pl[i]) continue;
1605 archp[j] = i;
1606 signs[j] = (pl[i] < 0)? 1: 0;
1607 j++;
1608 }
1609 setlg(archp, j); *parchp = archp;
1610 setlg(signs, j); *psigns = signs;
1611 }
1612 /* pl : requested signs for real embeddings, 0 = no sign constraint */
1613 int
nfchecksigns(GEN nf,GEN x,GEN pl)1614 nfchecksigns(GEN nf, GEN x, GEN pl)
1615 {
1616 pari_sp av = avma;
1617 GEN signs, archp;
1618 nf = checknf(nf);
1619 x = nf_to_scalar_or_basis(nf,x);
1620 if (typ(x) != t_COL)
1621 {
1622 long i, l = lg(pl), s = gsigne(x);
1623 for (i = 1; i < l; i++)
1624 if (pl[i] && pl[i] != s) return gc_bool(av,0);
1625 return gc_bool(av,1);
1626 }
1627 pl_convert(pl, &signs, &archp);
1628 return gc_bool(av, nfchecksigns_i(nf, x, NULL, signs, archp));
1629 }
1630
1631 /* signs = NULL: totally positive, else sign[i] = 0 (+) or 1 (-) */
1632 static GEN
get_C(GEN lambda,long l,GEN signs)1633 get_C(GEN lambda, long l, GEN signs)
1634 {
1635 long i;
1636 GEN C, mlambda;
1637 if (!signs) return const_vec(l-1, lambda);
1638 C = cgetg(l, t_COL); mlambda = gneg(lambda);
1639 for (i = 1; i < l; i++) gel(C,i) = signs[i]? mlambda: lambda;
1640 return C;
1641 }
1642 /* signs = NULL: totally positive at archp */
1643 static GEN
nfsetsigns(GEN nf,GEN signs,GEN x,GEN sarch)1644 nfsetsigns(GEN nf, GEN signs, GEN x, GEN sarch)
1645 {
1646 long i, l = lg(sarch_get_archp(sarch));
1647 GEN ex;
1648 /* Is signature already correct ? */
1649 if (typ(x) != t_COL)
1650 {
1651 long s = gsigne(x);
1652 if (!s) i = 1;
1653 else if (!signs)
1654 i = (s < 0)? 1: l;
1655 else
1656 {
1657 s = s < 0? 1: 0;
1658 for (i = 1; i < l; i++)
1659 if (signs[i] != s) break;
1660 }
1661 ex = (i < l)? const_col(l-1, x): NULL;
1662 }
1663 else
1664 {
1665 pari_sp av = avma;
1666 GEN cex, M = nf_get_M(nf), archp = sarch_get_archp(sarch);
1667 GEN xp = Q_primitive_part(x,&cex);
1668 ex = cgetg(l,t_COL);
1669 for (i = 1; i < l; i++) gel(ex,i) = zk_embed(M,xp,archp[i]);
1670 if (nfchecksigns_i(nf, xp, ex, signs, archp)) { ex = NULL; set_avma(av); }
1671 else if (cex) ex = RgC_Rg_mul(ex, cex); /* put back content */
1672 }
1673 if (ex)
1674 { /* If no, fix it */
1675 GEN MI = sarch_get_MI(sarch), F = sarch_get_F(sarch);
1676 GEN lambda = sarch_get_lambda(sarch);
1677 GEN t = RgC_sub(get_C(lambda, l, signs), ex);
1678 long e;
1679 t = grndtoi(RgM_RgC_mul(MI,t), &e);
1680 if (lg(F) != 1) t = ZM_ZC_mul(F, t);
1681 x = typ(x) == t_COL? RgC_add(t, x): RgC_Rg_add(t, x);
1682 }
1683 return x;
1684 }
1685 /* - sarch = nfarchstar(nf, F);
1686 * - x encodes a vector of signs at arch.archp: either a t_VECSMALL
1687 * (vector of signs as {0,1}-vector), NULL (totally positive at archp),
1688 * or a nonzero number field element (replaced by its signature at archp);
1689 * - y is a nonzero number field element
1690 * Return z = y (mod F) with signs(y, archp) = signs(x) (a {0,1}-vector) */
1691 GEN
set_sign_mod_divisor(GEN nf,GEN x,GEN y,GEN sarch)1692 set_sign_mod_divisor(GEN nf, GEN x, GEN y, GEN sarch)
1693 {
1694 GEN archp = sarch_get_archp(sarch);
1695 if (lg(archp) == 1) return y;
1696 nf = checknf(nf);
1697 if (x && typ(x) != t_VECSMALL) x = nfsign_arch(nf, x, archp);
1698 y = nf_to_scalar_or_basis(nf,y);
1699 return nfsetsigns(nf, x, y, sarch);
1700 }
1701
1702 static GEN
setsigns_init(GEN nf,GEN archp,GEN F,GEN DATA)1703 setsigns_init(GEN nf, GEN archp, GEN F, GEN DATA)
1704 {
1705 GEN lambda, Mr = rowpermute(nf_get_M(nf), archp), MI = F? RgM_mul(Mr,F): Mr;
1706 lambda = gmul2n(matrixnorm(MI,DEFAULTPREC), -1);
1707 if (typ(lambda) != t_REAL) lambda = gmul(lambda, sstoQ(1001,1000));
1708 if (lg(archp) < lg(MI))
1709 {
1710 GEN perm = gel(indexrank(MI), 2);
1711 if (!F) F = matid(nf_get_degree(nf));
1712 MI = vecpermute(MI, perm);
1713 F = vecpermute(F, perm);
1714 }
1715 if (!F) F = cgetg(1,t_MAT);
1716 MI = RgM_inv(MI);
1717 return mkvec5(DATA, archp, MI, lambda, F);
1718 }
1719 /* F nonzero integral ideal in HNF (or NULL: Z_K), compute elements in 1+F
1720 * whose sign matrix at archp is identity; archp in 'indices' format */
1721 GEN
nfarchstar(GEN nf,GEN F,GEN archp)1722 nfarchstar(GEN nf, GEN F, GEN archp)
1723 {
1724 long nba = lg(archp) - 1;
1725 if (!nba) return mkvec2(cgetg(1,t_VEC), archp);
1726 if (F && equali1(gcoeff(F,1,1))) F = NULL;
1727 if (F) F = idealpseudored(F, nf_get_roundG(nf));
1728 return setsigns_init(nf, archp, F, const_vec(nba, gen_2));
1729 }
1730
1731 /*************************************************************************/
1732 /** **/
1733 /** IDEALCHINESE **/
1734 /** **/
1735 /*************************************************************************/
1736 static int
isprfact(GEN x)1737 isprfact(GEN x)
1738 {
1739 long i, l;
1740 GEN L, E;
1741 if (typ(x) != t_MAT || lg(x) != 3) return 0;
1742 L = gel(x,1); l = lg(L);
1743 E = gel(x,2);
1744 for(i=1; i<l; i++)
1745 {
1746 checkprid(gel(L,i));
1747 if (typ(gel(E,i)) != t_INT) return 0;
1748 }
1749 return 1;
1750 }
1751
1752 /* initialize projectors mod pr[i]^e[i] for idealchinese */
1753 static GEN
pr_init(GEN nf,GEN fa,GEN w,GEN dw)1754 pr_init(GEN nf, GEN fa, GEN w, GEN dw)
1755 {
1756 GEN U, E, F, L = gel(fa,1), E0 = gel(fa,2);
1757 long i, r = lg(L);
1758
1759 if (w && lg(w) != r) pari_err_TYPE("idealchinese", w);
1760 if (r == 1 && !dw) return cgetg(1,t_VEC);
1761 E = leafcopy(E0); /* do not destroy fa[2] */
1762 for (i = 1; i < r; i++)
1763 if (signe(gel(E,i)) < 0) gel(E,i) = gen_0;
1764 F = factorbackprime(nf, L, E);
1765 if (dw)
1766 {
1767 F = ZM_Z_mul(F, dw);
1768 for (i = 1; i < r; i++)
1769 {
1770 GEN pr = gel(L,i);
1771 long e = itos(gel(E0,i)), v = idealval(nf, dw, pr);
1772 if (e >= 0)
1773 gel(E,i) = addiu(gel(E,i), v);
1774 else if (v + e <= 0)
1775 F = idealmulpowprime(nf, F, pr, stoi(-v)); /* coprime to pr */
1776 else
1777 {
1778 F = idealmulpowprime(nf, F, pr, stoi(e));
1779 gel(E,i) = stoi(v + e);
1780 }
1781 }
1782 }
1783 U = cgetg(r, t_VEC);
1784 for (i = 1; i < r; i++)
1785 {
1786 GEN u;
1787 if (w && gequal0(gel(w,i))) u = gen_0; /* unused */
1788 else
1789 {
1790 GEN pr = gel(L,i), e = gel(E,i), t;
1791 t = idealdivpowprime(nf,F, pr, e);
1792 u = hnfmerge_get_1(t, idealpow(nf, pr, e));
1793 if (!u) pari_err_COPRIME("idealchinese", t,pr);
1794 }
1795 gel(U,i) = u;
1796 }
1797 F = idealpseudored(F, nf_get_roundG(nf));
1798 return mkvec2(F, U);
1799 }
1800
1801 static GEN
pl_normalize(GEN nf,GEN pl)1802 pl_normalize(GEN nf, GEN pl)
1803 {
1804 const char *fun = "idealchinese";
1805 if (lg(pl)-1 != nf_get_r1(nf)) pari_err_TYPE(fun,pl);
1806 switch(typ(pl))
1807 {
1808 case t_VEC: RgV_check_ZV(pl,fun); pl = ZV_to_zv(pl);
1809 /* fall through */
1810 case t_VECSMALL: break;
1811 default: pari_err_TYPE(fun,pl);
1812 }
1813 return pl;
1814 }
1815
1816 static int
is_chineseinit(GEN x)1817 is_chineseinit(GEN x)
1818 {
1819 GEN fa, pl;
1820 long l;
1821 if (typ(x) != t_VEC || lg(x)!=3) return 0;
1822 fa = gel(x,1);
1823 pl = gel(x,2);
1824 if (typ(fa) != t_VEC || typ(pl) != t_VEC) return 0;
1825 l = lg(fa);
1826 if (l != 1)
1827 {
1828 if (l != 3 || typ(gel(fa,1)) != t_MAT || typ(gel(fa,2)) != t_VEC)
1829 return 0;
1830 }
1831 l = lg(pl);
1832 if (l != 1)
1833 {
1834 if (l != 6 || typ(gel(pl,3)) != t_MAT || typ(gel(pl,1)) != t_VECSMALL
1835 || typ(gel(pl,2)) != t_VECSMALL)
1836 return 0;
1837 }
1838 return 1;
1839 }
1840
1841 /* nf a true 'nf' */
1842 static GEN
chineseinit_i(GEN nf,GEN fa,GEN w,GEN dw)1843 chineseinit_i(GEN nf, GEN fa, GEN w, GEN dw)
1844 {
1845 const char *fun = "idealchineseinit";
1846 GEN archp = NULL, pl = NULL;
1847 switch(typ(fa))
1848 {
1849 case t_VEC:
1850 if (is_chineseinit(fa))
1851 {
1852 if (dw) pari_err_DOMAIN(fun, "denom(y)", "!=", gen_1, w);
1853 return fa;
1854 }
1855 if (lg(fa) != 3) pari_err_TYPE(fun, fa);
1856 /* of the form [x,s] */
1857 pl = pl_normalize(nf, gel(fa,2));
1858 fa = gel(fa,1);
1859 archp = vecsmall01_to_indices(pl);
1860 /* keep pr_init, reset pl */
1861 if (is_chineseinit(fa)) { fa = gel(fa,1); break; }
1862 /* fall through */
1863 case t_MAT: /* factorization? */
1864 if (isprfact(fa)) { fa = pr_init(nf, fa, w, dw); break; }
1865 default: pari_err_TYPE(fun,fa);
1866 }
1867
1868 if (!pl) pl = cgetg(1,t_VEC);
1869 else
1870 {
1871 long r = lg(archp);
1872 if (r == 1) pl = cgetg(1, t_VEC);
1873 else
1874 {
1875 GEN F = (lg(fa) == 1)? NULL: gel(fa,1), signs = cgetg(r, t_VECSMALL);
1876 long i;
1877 for (i = 1; i < r; i++) signs[i] = (pl[archp[i]] < 0)? 1: 0;
1878 pl = setsigns_init(nf, archp, F, signs);
1879 }
1880 }
1881 return mkvec2(fa, pl);
1882 }
1883
1884 /* Given a prime ideal factorization x, possibly with 0 or negative exponents,
1885 * and a vector w of elements of nf, gives b such that
1886 * v_p(b-w_p)>=v_p(x) for all prime ideals p in the ideal factorization
1887 * and v_p(b)>=0 for all other p, using the standard proof given in GTM 138. */
1888 GEN
idealchinese(GEN nf,GEN x,GEN w)1889 idealchinese(GEN nf, GEN x, GEN w)
1890 {
1891 const char *fun = "idealchinese";
1892 pari_sp av = avma;
1893 GEN x1, x2, s, dw, F;
1894
1895 nf = checknf(nf);
1896 if (!w) return gerepilecopy(av, chineseinit_i(nf,x,NULL,NULL));
1897
1898 if (typ(w) != t_VEC) pari_err_TYPE(fun,w);
1899 w = Q_remove_denom(matalgtobasis(nf,w), &dw);
1900 if (!is_chineseinit(x)) x = chineseinit_i(nf,x,w,dw);
1901 /* x is a 'chineseinit' */
1902 x1 = gel(x,1); s = NULL;
1903 x2 = gel(x,2);
1904 if (lg(x1) == 1) F = NULL;
1905 else
1906 {
1907 GEN U = gel(x1,2);
1908 long i, r = lg(w);
1909 F = gel(x1,1);
1910 for (i=1; i<r; i++)
1911 if (!gequal0(gel(w,i)))
1912 {
1913 GEN t = nfmuli(nf, gel(U,i), gel(w,i));
1914 s = s? ZC_add(s,t): t;
1915 }
1916 if (s) s = ZC_reducemodmatrix(s, F);
1917 }
1918 if (lg(x2) != 1) s = nfsetsigns(nf, gel(x2,1), s? s: gen_0, x2);
1919 if (!s) { s = zerocol(nf_get_degree(nf)); dw = NULL; }
1920
1921 if (dw) s = RgC_Rg_div(s,dw);
1922 return gerepileupto(av, s);
1923 }
1924
1925 /*************************************************************************/
1926 /** **/
1927 /** (Z_K/I)^* **/
1928 /** **/
1929 /*************************************************************************/
1930 GEN
vecsmall01_to_indices(GEN v)1931 vecsmall01_to_indices(GEN v)
1932 {
1933 long i, k, l = lg(v);
1934 GEN p = new_chunk(l) + l;
1935 for (k=1, i=l-1; i; i--)
1936 if (v[i]) { *--p = i; k++; }
1937 *--p = evallg(k) | evaltyp(t_VECSMALL);
1938 set_avma((pari_sp)p); return p;
1939 }
1940 GEN
vec01_to_indices(GEN v)1941 vec01_to_indices(GEN v)
1942 {
1943 long i, k, l;
1944 GEN p;
1945
1946 switch (typ(v))
1947 {
1948 case t_VECSMALL: return v;
1949 case t_VEC: break;
1950 default: pari_err_TYPE("vec01_to_indices",v);
1951 }
1952 l = lg(v);
1953 p = new_chunk(l) + l;
1954 for (k=1, i=l-1; i; i--)
1955 if (signe(gel(v,i))) { *--p = i; k++; }
1956 *--p = evallg(k) | evaltyp(t_VECSMALL);
1957 set_avma((pari_sp)p); return p;
1958 }
1959 GEN
indices_to_vec01(GEN p,long r)1960 indices_to_vec01(GEN p, long r)
1961 {
1962 long i, l = lg(p);
1963 GEN v = zerovec(r);
1964 for (i = 1; i < l; i++) gel(v, p[i]) = gen_1;
1965 return v;
1966 }
1967
1968 /* return (column) vector of R1 signatures of x (0 or 1) */
1969 GEN
nfsign_arch(GEN nf,GEN x,GEN arch)1970 nfsign_arch(GEN nf, GEN x, GEN arch)
1971 {
1972 GEN sarch, M, V, archp = vec01_to_indices(arch);
1973 long i, s, np, n = lg(archp)-1;
1974 pari_sp av;
1975
1976 if (!n) return cgetg(1,t_VECSMALL);
1977 nf = checknf(nf);
1978 if (typ(x) == t_MAT)
1979 { /* factorisation */
1980 GEN g = gel(x,1), e = gel(x,2);
1981 long l = lg(g);
1982 V = zero_zv(n);
1983 for (i = 1; i < l; i++)
1984 if (mpodd(gel(e,i)))
1985 Flv_add_inplace(V, nfsign_arch(nf,gel(g,i),archp), 2);
1986 set_avma((pari_sp)V); return V;
1987 }
1988 av = avma; V = cgetg(n+1,t_VECSMALL);
1989 x = nf_to_scalar_or_basis(nf, x);
1990 switch(typ(x))
1991 {
1992 case t_INT:
1993 s = signe(x);
1994 if (!s) pari_err_DOMAIN("nfsign_arch","element","=",gen_0,x);
1995 set_avma(av); return const_vecsmall(n, (s < 0)? 1: 0);
1996 case t_FRAC:
1997 s = signe(gel(x,1));
1998 set_avma(av); return const_vecsmall(n, (s < 0)? 1: 0);
1999 }
2000 x = Q_primpart(x); M = nf_get_M(nf); sarch = NULL; np = -1;
2001 for (i = 1; i <= n; i++)
2002 {
2003 long s = eval_sign(M, x, archp[i]);
2004 if (s < 0) /* failure */
2005 {
2006 long ni, r1 = nf_get_r1(nf);
2007 GEN xi;
2008 if (np < 0)
2009 {
2010 np = num_positive(nf, x);
2011 if (np == 0) { set_avma(av); return const_vecsmall(n, 1); }
2012 if (np == r1){ set_avma(av); return const_vecsmall(n, 0); }
2013 sarch = nfarchstar(nf, NULL, identity_perm(r1));
2014 }
2015 xi = set_sign_mod_divisor(nf, vecsmall_ei(r1, archp[i]), gen_1, sarch);
2016 xi = Q_primpart(xi);
2017 ni = num_positive(nf, nfmuli(nf,x,xi));
2018 if (ni == 0) { set_avma(av); V = const_vecsmall(n, 1); V[i] = 0; return V; }
2019 if (ni == r1){ set_avma(av); V = const_vecsmall(n, 0); V[i] = 1; return V; }
2020 s = ni < np? 0: 1;
2021 }
2022 V[i] = s;
2023 }
2024 set_avma((pari_sp)V); return V;
2025 }
2026 static void
chk_ind(const char * s,long i,long r1)2027 chk_ind(const char *s, long i, long r1)
2028 {
2029 if (i <= 0) pari_err_DOMAIN(s, "index", "<=", gen_0, stoi(i));
2030 if (i > r1) pari_err_DOMAIN(s, "index", ">", utoi(r1), utoi(i));
2031 }
2032 static GEN
parse_embed(GEN ind,long r,const char * f)2033 parse_embed(GEN ind, long r, const char *f)
2034 {
2035 long l, i;
2036 if (!ind) return identity_perm(r);
2037 switch(typ(ind))
2038 {
2039 case t_INT: case t_VEC: case t_COL: ind = gtovecsmall(ind); break;
2040 case t_VECSMALL: break;
2041 default: pari_err_TYPE(f, ind);
2042 }
2043 l = lg(ind);
2044 for (i = 1; i < l; i++) chk_ind(f, ind[i], r);
2045 return ind;
2046 }
2047 GEN
nfeltsign(GEN nf,GEN x,GEN ind0)2048 nfeltsign(GEN nf, GEN x, GEN ind0)
2049 {
2050 pari_sp av = avma;
2051 long i, l;
2052 GEN v, ind;
2053 nf = checknf(nf);
2054 ind = parse_embed(ind0, nf_get_r1(nf), "nfeltsign");
2055 l = lg(ind);
2056 if (is_rational_t(typ(x)))
2057 { /* nfsign_arch would test this, but avoid converting t_VECSMALL -> t_VEC */
2058 GEN s;
2059 switch(gsigne(x))
2060 {
2061 case -1:s = gen_m1; break;
2062 case 1: s = gen_1; break;
2063 default: s = gen_0; break;
2064 }
2065 set_avma(av);
2066 return (ind0 && typ(ind0) == t_INT)? s: const_vec(l-1, s);
2067 }
2068 v = nfsign_arch(nf, x, ind);
2069 if (ind0 && typ(ind0) == t_INT) { set_avma(av); return v[1]? gen_m1: gen_1; }
2070 settyp(v, t_VEC);
2071 for (i = 1; i < l; i++) gel(v,i) = v[i]? gen_m1: gen_1;
2072 return gerepileupto(av, v);
2073 }
2074
2075 GEN
nfeltembed(GEN nf,GEN x,GEN ind0,long prec0)2076 nfeltembed(GEN nf, GEN x, GEN ind0, long prec0)
2077 {
2078 pari_sp av = avma;
2079 long i, e, l, r1, r2, prec, prec1;
2080 GEN v, ind, cx;
2081 nf = checknf(nf); nf_get_sign(nf,&r1,&r2);
2082 x = nf_to_scalar_or_basis(nf, x);
2083 ind = parse_embed(ind0, r1+r2, "nfeltembed");
2084 l = lg(ind);
2085 if (typ(x) != t_COL)
2086 {
2087 if (!(ind0 && typ(ind0) == t_INT)) x = const_vec(l-1, x);
2088 return gerepilecopy(av, x);
2089 }
2090 x = Q_primitive_part(x, &cx);
2091 prec1 = prec0; e = gexpo(x);
2092 if (e > 8) prec1 += nbits2extraprec(e);
2093 prec = prec1;
2094 if (nf_get_prec(nf) < prec) nf = nfnewprec_shallow(nf, prec);
2095 v = cgetg(l, t_VEC);
2096 for(;;)
2097 {
2098 GEN M = nf_get_M(nf);
2099 for (i = 1; i < l; i++)
2100 {
2101 GEN t = nfembed_i(M, x, ind[i]);
2102 long e = gexpo(t);
2103 if (gequal0(t) || precision(t) < prec0
2104 || (e < 0 && prec < prec1 + nbits2extraprec(-e)) ) break;
2105 if (cx) t = gmul(t, cx);
2106 gel(v,i) = t;
2107 }
2108 if (i == l) break;
2109 prec = precdbl(prec);
2110 if (DEBUGLEVEL>1) pari_warn(warnprec,"eltnfembed", prec);
2111 nf = nfnewprec_shallow(nf, prec);
2112 }
2113 if (ind0 && typ(ind0) == t_INT) v = gel(v,1);
2114 return gerepilecopy(av, v);
2115 }
2116
2117 /* number of distinct roots of sigma(f) */
2118 GEN
nfpolsturm(GEN nf,GEN f,GEN ind0)2119 nfpolsturm(GEN nf, GEN f, GEN ind0)
2120 {
2121 pari_sp av = avma;
2122 long d, l, r1, single;
2123 GEN ind, u, v, vr1, T, s, t;
2124
2125 nf = checknf(nf); T = nf_get_pol(nf); r1 = nf_get_r1(nf);
2126 ind = parse_embed(ind0, r1, "nfpolsturm");
2127 single = ind0 && typ(ind0) == t_INT;
2128 l = lg(ind);
2129
2130 if (gequal0(f)) pari_err_ROOTS0("nfpolsturm");
2131 if (typ(f) == t_POL && varn(f) != varn(T))
2132 {
2133 f = RgX_nffix("nfpolsturm", T, f,1);
2134 if (lg(f) == 3) f = NULL;
2135 }
2136 else
2137 {
2138 (void)Rg_nffix("nfpolsturm", T, f, 0);
2139 f = NULL;
2140 }
2141 if (!f) { set_avma(av); return single? gen_0: zerovec(l-1); }
2142 d = degpol(f);
2143 if (d == 1) { set_avma(av); return single? gen_1: const_vec(l-1,gen_1); }
2144
2145 vr1 = const_vecsmall(l-1, 1);
2146 u = Q_primpart(f); s = ZV_to_zv(nfeltsign(nf, gel(u,d+2), ind));
2147 v = RgX_deriv(u); t = odd(d)? leafcopy(s): zv_neg(s);
2148 for(;;)
2149 {
2150 GEN r = RgX_neg( Q_primpart(RgX_pseudorem(u, v)) ), sr;
2151 long i, dr = degpol(r);
2152 if (dr < 0) break;
2153 sr = ZV_to_zv(nfeltsign(nf, gel(r,dr+2), ind));
2154 for (i = 1; i < l; i++)
2155 if (sr[i] != s[i]) { s[i] = sr[i], vr1[i]--; }
2156 if (odd(dr)) sr = zv_neg(sr);
2157 for (i = 1; i < l; i++)
2158 if (sr[i] != t[i]) { t[i] = sr[i], vr1[i]++; }
2159 if (!dr) break;
2160 u = v; v = r;
2161 }
2162 if (single) { set_avma(av); return stoi(vr1[1]); }
2163 return gerepileupto(av, zv_to_ZV(vr1));
2164 }
2165
2166 /* return the vector of signs of x; the matrix of such if x is a vector
2167 * of nf elements */
2168 GEN
nfsign(GEN nf,GEN x)2169 nfsign(GEN nf, GEN x)
2170 {
2171 long i, l;
2172 GEN archp, S;
2173
2174 nf = checknf(nf);
2175 archp = identity_perm( nf_get_r1(nf) );
2176 if (typ(x) != t_VEC) return nfsign_arch(nf, x, archp);
2177 l = lg(x); S = cgetg(l, t_MAT);
2178 for (i=1; i<l; i++) gel(S,i) = nfsign_arch(nf, gel(x,i), archp);
2179 return S;
2180 }
2181
2182 /* x integral elt, A integral ideal in HNF; reduce x mod A */
2183 static GEN
zk_modHNF(GEN x,GEN A)2184 zk_modHNF(GEN x, GEN A)
2185 { return (typ(x) == t_COL)? ZC_hnfrem(x, A): modii(x, gcoeff(A,1,1)); }
2186
2187 /* given an element x in Z_K and an integral ideal y in HNF, coprime with x,
2188 outputs an element inverse of x modulo y */
2189 GEN
nfinvmodideal(GEN nf,GEN x,GEN y)2190 nfinvmodideal(GEN nf, GEN x, GEN y)
2191 {
2192 pari_sp av = avma;
2193 GEN a, yZ = gcoeff(y,1,1);
2194
2195 if (equali1(yZ)) return gen_0;
2196 x = nf_to_scalar_or_basis(nf, x);
2197 if (typ(x) == t_INT) return gerepileupto(av, Fp_inv(x, yZ));
2198
2199 a = hnfmerge_get_1(idealhnf_principal(nf,x), y);
2200 if (!a) pari_err_INV("nfinvmodideal", x);
2201 return gerepileupto(av, zk_modHNF(nfdiv(nf,a,x), y));
2202 }
2203
2204 static GEN
nfsqrmodideal(GEN nf,GEN x,GEN id)2205 nfsqrmodideal(GEN nf, GEN x, GEN id)
2206 { return zk_modHNF(nfsqri(nf,x), id); }
2207 static GEN
nfmulmodideal(GEN nf,GEN x,GEN y,GEN id)2208 nfmulmodideal(GEN nf, GEN x, GEN y, GEN id)
2209 { return x? zk_modHNF(nfmuli(nf,x,y), id): y; }
2210 /* assume x integral, k integer, A in HNF */
2211 GEN
nfpowmodideal(GEN nf,GEN x,GEN k,GEN A)2212 nfpowmodideal(GEN nf,GEN x,GEN k,GEN A)
2213 {
2214 long s = signe(k);
2215 pari_sp av;
2216 GEN y;
2217
2218 if (!s) return gen_1;
2219 av = avma;
2220 x = nf_to_scalar_or_basis(nf, x);
2221 if (typ(x) != t_COL) return Fp_pow(x, k, gcoeff(A,1,1));
2222 if (s < 0) { x = nfinvmodideal(nf, x,A); k = negi(k); }
2223 for(y = NULL;;)
2224 {
2225 if (mpodd(k)) y = nfmulmodideal(nf,y,x,A);
2226 k = shifti(k,-1); if (!signe(k)) break;
2227 x = nfsqrmodideal(nf,x,A);
2228 }
2229 return gerepileupto(av, y);
2230 }
2231
2232 /* a * g^n mod id */
2233 static GEN
nfmulpowmodideal(GEN nf,GEN a,GEN g,GEN n,GEN id)2234 nfmulpowmodideal(GEN nf, GEN a, GEN g, GEN n, GEN id)
2235 {
2236 return nfmulmodideal(nf, a, nfpowmodideal(nf,g,n,id), id);
2237 }
2238
2239 /* assume (num(g[i]), id) = 1 for all i. Return prod g[i]^e[i] mod id.
2240 * EX = multiple of exponent of (O_K/id)^* */
2241 GEN
famat_to_nf_modideal_coprime(GEN nf,GEN g,GEN e,GEN id,GEN EX)2242 famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN id, GEN EX)
2243 {
2244 GEN EXo2, plus = NULL, minus = NULL, idZ = gcoeff(id,1,1);
2245 long i, lx = lg(g);
2246
2247 if (equali1(idZ)) return gen_1; /* id = Z_K */
2248 EXo2 = (expi(EX) > 10)? shifti(EX,-1): NULL;
2249 for (i = 1; i < lx; i++)
2250 {
2251 GEN h, n = centermodii(gel(e,i), EX, EXo2);
2252 long sn = signe(n);
2253 if (!sn) continue;
2254
2255 h = nf_to_scalar_or_basis(nf, gel(g,i));
2256 switch(typ(h))
2257 {
2258 case t_INT: break;
2259 case t_FRAC:
2260 h = Fp_div(gel(h,1), gel(h,2), idZ); break;
2261 default:
2262 {
2263 GEN dh;
2264 h = Q_remove_denom(h, &dh);
2265 if (dh) h = FpC_Fp_mul(h, Fp_inv(dh,idZ), idZ);
2266 }
2267 }
2268 if (sn > 0)
2269 plus = nfmulpowmodideal(nf, plus, h, n, id);
2270 else /* sn < 0 */
2271 minus = nfmulpowmodideal(nf, minus, h, negi(n), id);
2272 }
2273 if (minus) plus = nfmulmodideal(nf, plus, nfinvmodideal(nf,minus,id), id);
2274 return plus? plus: gen_1;
2275 }
2276
2277 /* given 2 integral ideals x, y in HNF s.t x | y | x^2, compute (1+x)/(1+y) in
2278 * the form [[cyc],[gen], U], where U := ux^-1 as a pair [ZM, denom(U)] */
2279 static GEN
zidealij(GEN x,GEN y)2280 zidealij(GEN x, GEN y)
2281 {
2282 GEN U, G, cyc, xp = gcoeff(x,1,1), xi = hnf_invscale(x, xp);
2283 long j, N;
2284
2285 /* x^(-1) y = relations between the 1 + x_i (HNF) */
2286 cyc = ZM_snf_group(ZM_Z_divexact(ZM_mul(xi, y), xp), &U, &G);
2287 N = lg(cyc); G = ZM_mul(x,G); settyp(G, t_VEC); /* new generators */
2288 for (j=1; j<N; j++)
2289 {
2290 GEN c = gel(G,j);
2291 gel(c,1) = addiu(gel(c,1), 1); /* 1 + g_j */
2292 if (ZV_isscalar(c)) gel(G,j) = gel(c,1);
2293 }
2294 return mkvec4(cyc, G, ZM_mul(U,xi), xp);
2295 }
2296
2297 /* lg(x) > 1, x + 1; shallow */
2298 static GEN
ZC_add1(GEN x)2299 ZC_add1(GEN x)
2300 {
2301 long i, l = lg(x);
2302 GEN y = cgetg(l, t_COL);
2303 for (i = 2; i < l; i++) gel(y,i) = gel(x,i);
2304 gel(y,1) = addiu(gel(x,1), 1); return y;
2305 }
2306 /* lg(x) > 1, x - 1; shallow */
2307 static GEN
ZC_sub1(GEN x)2308 ZC_sub1(GEN x)
2309 {
2310 long i, l = lg(x);
2311 GEN y = cgetg(l, t_COL);
2312 for (i = 2; i < l; i++) gel(y,i) = gel(x,i);
2313 gel(y,1) = subiu(gel(x,1), 1); return y;
2314 }
2315
2316 /* x,y are t_INT or ZC */
2317 static GEN
zkadd(GEN x,GEN y)2318 zkadd(GEN x, GEN y)
2319 {
2320 long tx = typ(x);
2321 if (tx == typ(y))
2322 return tx == t_INT? addii(x,y): ZC_add(x,y);
2323 else
2324 return tx == t_INT? ZC_Z_add(y,x): ZC_Z_add(x,y);
2325 }
2326 /* x a t_INT or ZC, x+1; shallow */
2327 static GEN
zkadd1(GEN x)2328 zkadd1(GEN x)
2329 {
2330 long tx = typ(x);
2331 return tx == t_INT? addiu(x,1): ZC_add1(x);
2332 }
2333 /* x a t_INT or ZC, x-1; shallow */
2334 static GEN
zksub1(GEN x)2335 zksub1(GEN x)
2336 {
2337 long tx = typ(x);
2338 return tx == t_INT? subiu(x,1): ZC_sub1(x);
2339 }
2340 /* x,y are t_INT or ZC; x - y */
2341 static GEN
zksub(GEN x,GEN y)2342 zksub(GEN x, GEN y)
2343 {
2344 long tx = typ(x), ty = typ(y);
2345 if (tx == ty)
2346 return tx == t_INT? subii(x,y): ZC_sub(x,y);
2347 else
2348 return tx == t_INT? Z_ZC_sub(x,y): ZC_Z_sub(x,y);
2349 }
2350 /* x is t_INT or ZM (mult. map), y is t_INT or ZC; x * y */
2351 static GEN
zkmul(GEN x,GEN y)2352 zkmul(GEN x, GEN y)
2353 {
2354 long tx = typ(x), ty = typ(y);
2355 if (ty == t_INT)
2356 return tx == t_INT? mulii(x,y): ZC_Z_mul(gel(x,1),y);
2357 else
2358 return tx == t_INT? ZC_Z_mul(y,x): ZM_ZC_mul(x,y);
2359 }
2360
2361 /* (U,V) = 1 coprime ideals. Want z = x mod U, = y mod V; namely
2362 * z =vx + uy = v(x-y) + y, where u + v = 1, u in U, v in V.
2363 * zkc = [v, UV], v a t_INT or ZM (mult. by v map), UV a ZM (ideal in HNF);
2364 * shallow */
2365 GEN
zkchinese(GEN zkc,GEN x,GEN y)2366 zkchinese(GEN zkc, GEN x, GEN y)
2367 {
2368 GEN v = gel(zkc,1), UV = gel(zkc,2), z = zkadd(zkmul(v, zksub(x,y)), y);
2369 return zk_modHNF(z, UV);
2370 }
2371 /* special case z = x mod U, = 1 mod V; shallow */
2372 GEN
zkchinese1(GEN zkc,GEN x)2373 zkchinese1(GEN zkc, GEN x)
2374 {
2375 GEN v = gel(zkc,1), UV = gel(zkc,2), z = zkadd1(zkmul(v, zksub1(x)));
2376 return (typ(z) == t_INT)? z: ZC_hnfrem(z, UV);
2377 }
2378 static GEN
zkVchinese1(GEN zkc,GEN v)2379 zkVchinese1(GEN zkc, GEN v)
2380 {
2381 long i, ly;
2382 GEN y = cgetg_copy(v, &ly);
2383 for (i=1; i<ly; i++) gel(y,i) = zkchinese1(zkc, gel(v,i));
2384 return y;
2385 }
2386
2387 /* prepare to solve z = x (mod A), z = y mod (B) [zkchinese or zkchinese1] */
2388 GEN
zkchineseinit(GEN nf,GEN A,GEN B,GEN AB)2389 zkchineseinit(GEN nf, GEN A, GEN B, GEN AB)
2390 {
2391 GEN v;
2392 long e;
2393 nf = checknf(nf);
2394 v = idealaddtoone_raw(nf, A, B);
2395 if ((e = gexpo(v)) > 5)
2396 {
2397 GEN b = (typ(v) == t_COL)? v: scalarcol_shallow(v, nf_get_degree(nf));
2398 b= ZC_reducemodlll(b, AB);
2399 if (gexpo(b) < e) v = b;
2400 }
2401 return mkvec2(zk_scalar_or_multable(nf,v), AB);
2402 }
2403 /* prepare to solve z = x (mod A), z = 1 mod (B)
2404 * and then z = 1 (mod A), z = y mod (B) [zkchinese1 twice] */
2405 static GEN
zkchinese1init2(GEN nf,GEN A,GEN B,GEN AB)2406 zkchinese1init2(GEN nf, GEN A, GEN B, GEN AB)
2407 {
2408 GEN zkc = zkchineseinit(nf, A, B, AB);
2409 GEN mv = gel(zkc,1), mu;
2410 if (typ(mv) == t_INT) return mkvec2(zkc, mkvec2(subui(1,mv),AB));
2411 mu = RgM_Rg_add_shallow(ZM_neg(mv), gen_1);
2412 return mkvec2(mkvec2(mv,AB), mkvec2(mu,AB));
2413 }
2414
2415 static GEN
apply_U(GEN L,GEN a)2416 apply_U(GEN L, GEN a)
2417 {
2418 GEN e, U = gel(L,3), dU = gel(L,4);
2419 if (typ(a) == t_INT)
2420 e = ZC_Z_mul(gel(U,1), subiu(a, 1));
2421 else
2422 { /* t_COL */
2423 GEN t = shallowcopy(a);
2424 gel(t,1) = subiu(gel(t,1), 1); /* t = a - 1 */
2425 e = ZM_ZC_mul(U, t);
2426 }
2427 return gdiv(e, dU);
2428 }
2429
2430 /* true nf; vectors of [[cyc],[g],U.X^-1]. Assume k > 1. */
2431 static GEN
principal_units(GEN nf,GEN pr,long k,GEN prk)2432 principal_units(GEN nf, GEN pr, long k, GEN prk)
2433 {
2434 GEN list, prb;
2435 ulong mask = quadratic_prec_mask(k);
2436 long a = 1;
2437
2438 prb = pr_hnf(nf,pr);
2439 list = vectrunc_init(k);
2440 while (mask > 1)
2441 {
2442 GEN pra = prb;
2443 long b = a << 1;
2444
2445 if (mask & 1) b--;
2446 mask >>= 1;
2447 /* compute 1 + pr^a / 1 + pr^b, 2a <= b */
2448 prb = (b >= k)? prk: idealpows(nf,pr,b);
2449 vectrunc_append(list, zidealij(pra, prb));
2450 a = b;
2451 }
2452 return list;
2453 }
2454 /* a = 1 mod (pr) return log(a) on local-gens of 1+pr/1+pr^k */
2455 static GEN
log_prk1(GEN nf,GEN a,long nh,GEN L2,GEN prk)2456 log_prk1(GEN nf, GEN a, long nh, GEN L2, GEN prk)
2457 {
2458 GEN y = cgetg(nh+1, t_COL);
2459 long j, iy, c = lg(L2)-1;
2460 for (j = iy = 1; j <= c; j++)
2461 {
2462 GEN L = gel(L2,j), cyc = gel(L,1), gen = gel(L,2), E = apply_U(L,a);
2463 long i, nc = lg(cyc)-1;
2464 int last = (j == c);
2465 for (i = 1; i <= nc; i++, iy++)
2466 {
2467 GEN t, e = gel(E,i);
2468 if (typ(e) != t_INT) pari_err_COPRIME("zlog_prk1", a, prk);
2469 t = Fp_neg(e, gel(cyc,i));
2470 gel(y,iy) = negi(t);
2471 if (!last && signe(t)) a = nfmulpowmodideal(nf, a, gel(gen,i), t, prk);
2472 }
2473 }
2474 return y;
2475 }
2476 /* true nf */
2477 static GEN
principal_units_relations(GEN nf,GEN L2,GEN prk,long nh)2478 principal_units_relations(GEN nf, GEN L2, GEN prk, long nh)
2479 {
2480 GEN h = cgetg(nh+1,t_MAT);
2481 long ih, j, c = lg(L2)-1;
2482 for (j = ih = 1; j <= c; j++)
2483 {
2484 GEN L = gel(L2,j), F = gel(L,1), G = gel(L,2);
2485 long k, lG = lg(G);
2486 for (k = 1; k < lG; k++,ih++)
2487 { /* log(g^f) mod pr^e */
2488 GEN a = nfpowmodideal(nf,gel(G,k),gel(F,k),prk);
2489 gel(h,ih) = ZC_neg(log_prk1(nf, a, nh, L2, prk));
2490 gcoeff(h,ih,ih) = gel(F,k);
2491 }
2492 }
2493 return h;
2494 }
2495 /* true nf; k > 1; multiplicative group (1 + pr) / (1 + pr^k) */
2496 static GEN
idealprincipalunits_i(GEN nf,GEN pr,long k,GEN * pU)2497 idealprincipalunits_i(GEN nf, GEN pr, long k, GEN *pU)
2498 {
2499 GEN cyc, gen, L2, prk = idealpows(nf, pr, k);
2500
2501 L2 = principal_units(nf, pr, k, prk);
2502 if (k == 2)
2503 {
2504 GEN L = gel(L2,1);
2505 cyc = gel(L,1);
2506 gen = gel(L,2);
2507 if (pU) *pU = matid(lg(gen)-1);
2508 }
2509 else
2510 {
2511 long c = lg(L2), j;
2512 GEN EX, h, Ui, vg = cgetg(c, t_VEC);
2513 for (j = 1; j < c; j++) gel(vg, j) = gmael(L2,j,2);
2514 vg = shallowconcat1(vg);
2515 h = principal_units_relations(nf, L2, prk, lg(vg)-1);
2516 h = ZM_hnfall_i(h, NULL, 0);
2517 cyc = ZM_snf_group(h, pU, &Ui);
2518 c = lg(Ui); gen = cgetg(c, t_VEC); EX = cyc_get_expo(cyc);
2519 for (j = 1; j < c; j++)
2520 gel(gen,j) = famat_to_nf_modideal_coprime(nf, vg, gel(Ui,j), prk, EX);
2521 }
2522 return mkvec4(cyc, gen, prk, L2);
2523 }
2524 GEN
idealprincipalunits(GEN nf,GEN pr,long k)2525 idealprincipalunits(GEN nf, GEN pr, long k)
2526 {
2527 pari_sp av;
2528 GEN v;
2529 nf = checknf(nf);
2530 if (k == 1) { checkprid(pr); retmkvec3(gen_1,cgetg(1,t_VEC),cgetg(1,t_VEC)); }
2531 av = avma; v = idealprincipalunits_i(nf, pr, k, NULL);
2532 return gerepilecopy(av, mkvec3(powiu(pr_norm(pr), k-1), gel(v,1), gel(v,2)));
2533 }
2534
2535 /* true nf; given an ideal pr^k dividing an integral ideal x (in HNF form)
2536 * compute an 'sprk', the structure of G = (Z_K/pr^k)^* [ x = NULL for x=pr^k ]
2537 * Return a vector with at least 4 components [cyc],[gen],[HNF pr^k,pr,k],ff,
2538 * where
2539 * cyc : type of G as abelian group (SNF)
2540 * gen : generators of G, coprime to x
2541 * pr^k: in HNF
2542 * ff : data for log_g in (Z_K/pr)^*
2543 * Two extra components are present iff k > 1: L2, U
2544 * L2 : list of data structures to compute local DL in (Z_K/pr)^*,
2545 * and 1 + pr^a/ 1 + pr^b for various a < b <= min(2a, k)
2546 * U : base change matrices to convert a vector of local DL to DL wrt gen
2547 * If MOD is not NULL, initialize G / G^MOD instead */
2548 static GEN
sprkinit(GEN nf,GEN pr,long k,GEN x,GEN MOD)2549 sprkinit(GEN nf, GEN pr, long k, GEN x, GEN MOD)
2550 {
2551 GEN T, p, Ld, modpr, cyc, gen, g, g0, A, prk, U, L2, ord0 = NULL;
2552 long f = pr_get_f(pr);
2553
2554 if(DEBUGLEVEL>3) err_printf("treating pr^%ld, pr = %Ps\n",k,pr);
2555 modpr = nf_to_Fq_init(nf, &pr,&T,&p);
2556 if (MOD)
2557 {
2558 GEN A = subiu(powiu(p,f), 1), d = gcdii(A, MOD), fa = Z_factor(d);
2559 ord0 = mkvec2(A, fa); /* true order, factorization of order in G/G^MOD */
2560 Ld = gel(fa,1);
2561 if (lg(Ld) > 1 && equaliu(gel(Ld,1),2)) Ld = vecslice(Ld,2,lg(Ld)-1);
2562 }
2563 /* (Z_K / pr)^* */
2564 if (f == 1)
2565 {
2566 g0 = g = MOD? pgener_Fp_local(p, Ld): pgener_Fp(p);
2567 if (!ord0) ord0 = get_arith_ZZM(subiu(p,1));
2568 }
2569 else
2570 {
2571 g0 = g = MOD? gener_FpXQ_local(T, p, Ld): gener_FpXQ(T,p, &ord0);
2572 g = Fq_to_nf(g, modpr);
2573 if (typ(g) == t_POL) g = poltobasis(nf, g);
2574 }
2575 A = gel(ord0, 1); /* Norm(pr)-1 */
2576 /* If MOD != NULL, d = gcd(A, MOD): g^(A/d) has order d */
2577 if (k == 1)
2578 {
2579 cyc = mkvec(A);
2580 gen = mkvec(g);
2581 prk = pr_hnf(nf,pr);
2582 L2 = U = NULL;
2583 }
2584 else
2585 { /* local-gens of (1 + pr)/(1 + pr^k) = SNF-gens * U */
2586 GEN AB, B, u, v, w;
2587 long j, l;
2588 w = idealprincipalunits_i(nf, pr, k, &U);
2589 /* incorporate (Z_K/pr)^*, order A coprime to B = expo(1+pr/1+pr^k)*/
2590 cyc = leafcopy(gel(w,1)); B = cyc_get_expo(cyc); AB = mulii(A,B);
2591 gen = leafcopy(gel(w,2));
2592 prk = gel(w,3);
2593 g = nfpowmodideal(nf, g, B, prk);
2594 g0 = Fq_pow(g0, modii(B,A), T, p); /* update primitive root */
2595 L2 = mkvec3(A, g, gel(w,4));
2596 gel(cyc,1) = AB;
2597 gel(gen,1) = nfmulmodideal(nf, gel(gen,1), g, prk);
2598 u = mulii(Fp_inv(A,B), A);
2599 v = subui(1, u); l = lg(U);
2600 for (j = 1; j < l; j++) gcoeff(U,1,j) = Fp_mul(u, gcoeff(U,1,j), AB);
2601 U = mkvec2(Rg_col_ei(v, lg(gen)-1, 1), U);
2602 }
2603 /* local-gens of (Z_K/pr^k)^* = SNF-gens * U */
2604 if (x)
2605 {
2606 GEN uv = zkchineseinit(nf, idealmulpowprime(nf,x,pr,utoineg(k)), prk, x);
2607 gen = zkVchinese1(uv, gen);
2608 }
2609 return mkvecn(U? 6: 4, cyc, gen, prk, mkvec3(modpr,g0,ord0), L2, U);
2610 }
2611 static GEN
sprk_get_cyc(GEN s)2612 sprk_get_cyc(GEN s) { return gel(s,1); }
2613 static GEN
sprk_get_expo(GEN s)2614 sprk_get_expo(GEN s) { return cyc_get_expo(sprk_get_cyc(s)); }
2615 static GEN
sprk_get_gen(GEN s)2616 sprk_get_gen(GEN s) { return gel(s,2); }
2617 static GEN
sprk_get_prk(GEN s)2618 sprk_get_prk(GEN s) { return gel(s,3); }
2619 static GEN
sprk_get_ff(GEN s)2620 sprk_get_ff(GEN s) { return gel(s,4); }
2621 static GEN
sprk_get_pr(GEN s)2622 sprk_get_pr(GEN s) { GEN ff = gel(s,4); return modpr_get_pr(gel(ff,1)); }
2623 /* L2 to 1 + pr / 1 + pr^k */
2624 static GEN
sprk_get_L2(GEN s)2625 sprk_get_L2(GEN s) { return gmael(s,5,3); }
2626 /* lift to nf of primitive root of k(pr) */
2627 static GEN
sprk_get_gnf(GEN s)2628 sprk_get_gnf(GEN s) { return gmael(s,5,2); }
2629 static void
sprk_get_U2(GEN s,GEN * U1,GEN * U2)2630 sprk_get_U2(GEN s, GEN *U1, GEN *U2)
2631 { GEN v = gel(s,6); *U1 = gel(v,1); *U2 = gel(v,2); }
2632 static int
sprk_is_prime(GEN s)2633 sprk_is_prime(GEN s) { return lg(s) == 5; }
2634
2635 static GEN
famat_zlog_pr(GEN nf,GEN g,GEN e,GEN sprk,GEN mod)2636 famat_zlog_pr(GEN nf, GEN g, GEN e, GEN sprk, GEN mod)
2637 {
2638 GEN x, expo = sprk_get_expo(sprk);
2639 if (mod) expo = gcdii(expo,mod);
2640 x = famat_makecoprime(nf, g, e, sprk_get_pr(sprk), sprk_get_prk(sprk), expo);
2641 return log_prk(nf, x, sprk, mod);
2642 }
2643 /* famat_zlog_pr assuming (g,sprk.pr) = 1 */
2644 static GEN
famat_zlog_pr_coprime(GEN nf,GEN g,GEN e,GEN sprk,GEN MOD)2645 famat_zlog_pr_coprime(GEN nf, GEN g, GEN e, GEN sprk, GEN MOD)
2646 {
2647 GEN x = famat_to_nf_modideal_coprime(nf, g, e, sprk_get_prk(sprk),
2648 sprk_get_expo(sprk));
2649 return log_prk(nf, x, sprk, MOD);
2650 }
2651
2652 /* o t_INT, O = [ord,fa] format for multiple of o (for Fq_log);
2653 * return o in [ord,fa] format */
2654 static GEN
order_update(GEN o,GEN O)2655 order_update(GEN o, GEN O)
2656 {
2657 GEN p = gmael(O,2,1), z = o, P, E;
2658 long i, j, l = lg(p);
2659 P = cgetg(l, t_COL);
2660 E = cgetg(l, t_COL);
2661 for (i = j = 1; i < l; i++)
2662 {
2663 long v = Z_pvalrem(z, gel(p,i), &z);
2664 if (v)
2665 {
2666 gel(P,j) = gel(p,i);
2667 gel(E,j) = utoipos(v); j++;
2668 if (is_pm1(z)) break;
2669 }
2670 }
2671 setlg(P, j);
2672 setlg(E, j); return mkvec2(o, mkmat2(P,E));
2673 }
2674
2675 /* a in Z_K (t_COL or t_INT), pr prime ideal, sprk = sprkinit(nf,pr,k,x),
2676 * mod positive t_INT or NULL (meaning mod=0).
2677 * return log(a) modulo mod on SNF-generators of (Z_K/pr^k)^* */
2678 GEN
log_prk(GEN nf,GEN a,GEN sprk,GEN mod)2679 log_prk(GEN nf, GEN a, GEN sprk, GEN mod)
2680 {
2681 GEN e, prk, g, U1, U2, y, ff, O, o, oN, gN, N, T, p, modpr, pr, cyc;
2682
2683 if (typ(a) == t_MAT) return famat_zlog_pr(nf, gel(a,1), gel(a,2), sprk, mod);
2684 N = NULL;
2685 ff = sprk_get_ff(sprk);
2686 pr = gel(ff,1); /* modpr */
2687 g = gN = gel(ff,2);
2688 O = gel(ff,3); /* order of g = |Fq^*|, in [ord, fa] format */
2689 o = oN = gel(O,1); /* order as a t_INT */
2690 prk = sprk_get_prk(sprk);
2691 modpr = nf_to_Fq_init(nf, &pr, &T, &p);
2692 if (mod)
2693 {
2694 GEN d = gcdii(o,mod);
2695 if (!equalii(o, d))
2696 {
2697 N = diviiexact(o,d); /* > 1, coprime to p */
2698 a = nfpowmodideal(nf, a, N, prk);
2699 oN = d; /* order of g^N mod pr */
2700 }
2701 }
2702 if (equali1(oN))
2703 e = gen_0;
2704 else
2705 {
2706 if (N) { O = order_update(oN, O); gN = Fq_pow(g, N, T, p); }
2707 e = Fq_log(nf_to_Fq(nf,a,modpr), gN, O, T, p);
2708 }
2709 /* 0 <= e < oN is correct modulo oN */
2710 if (sprk_is_prime(sprk)) return mkcol(e); /* k = 1 */
2711
2712 sprk_get_U2(sprk, &U1,&U2);
2713 cyc = sprk_get_cyc(sprk);
2714 if (mod)
2715 {
2716 cyc = ZV_snf_gcd(cyc, mod);
2717 if (signe(remii(mod,p))) return vecmodii(ZC_Z_mul(U1,e), cyc);
2718 }
2719 if (signe(e))
2720 {
2721 GEN E = N? mulii(e, N): e;
2722 a = nfmulpowmodideal(nf, a, sprk_get_gnf(sprk), Fp_neg(E, o), prk);
2723 }
2724 /* a = 1 mod pr */
2725 y = log_prk1(nf, a, lg(U2)-1, sprk_get_L2(sprk), prk);
2726 if (N)
2727 { /* from DL(a^N) to DL(a) */
2728 GEN E = gel(sprk_get_cyc(sprk), 1), q = powiu(p, Z_pval(E, p));
2729 y = ZC_Z_mul(y, Fp_inv(N, q));
2730 }
2731 y = ZC_lincomb(gen_1, e, ZM_ZC_mul(U2,y), U1);
2732 return vecmodii(y, cyc);
2733 }
2734 GEN
log_prk_init(GEN nf,GEN pr,long k,GEN MOD)2735 log_prk_init(GEN nf, GEN pr, long k, GEN MOD)
2736 { return sprkinit(checknf(nf),pr,k,NULL,MOD);}
2737 GEN
veclog_prk(GEN nf,GEN v,GEN sprk)2738 veclog_prk(GEN nf, GEN v, GEN sprk)
2739 {
2740 long l = lg(v), i;
2741 GEN w = cgetg(l, t_MAT);
2742 for (i = 1; i < l; i++) gel(w,i) = log_prk(nf, gel(v,i), sprk, NULL);
2743 return w;
2744 }
2745
2746 static GEN
famat_zlog(GEN nf,GEN fa,GEN sgn,zlog_S * S)2747 famat_zlog(GEN nf, GEN fa, GEN sgn, zlog_S *S)
2748 {
2749 long i, l0, l = lg(S->U);
2750 GEN g = gel(fa,1), e = gel(fa,2), y = cgetg(l, t_COL);
2751 l0 = lg(S->sprk); /* = l (trivial arch. part), or l-1 */
2752 for (i=1; i < l0; i++) gel(y,i) = famat_zlog_pr(nf, g, e, gel(S->sprk,i), S->mod);
2753 if (l0 != l)
2754 {
2755 if (!sgn) sgn = nfsign_arch(nf, fa, S->archp);
2756 gel(y,l0) = Flc_to_ZC(sgn);
2757 }
2758 return y;
2759 }
2760
2761 /* assume that cyclic factors are normalized, in particular != [1] */
2762 static GEN
split_U(GEN U,GEN Sprk)2763 split_U(GEN U, GEN Sprk)
2764 {
2765 long t = 0, k, n, l = lg(Sprk);
2766 GEN vU = cgetg(l+1, t_VEC);
2767 for (k = 1; k < l; k++)
2768 {
2769 n = lg(sprk_get_cyc(gel(Sprk,k))) - 1; /* > 0 */
2770 gel(vU,k) = vecslice(U, t+1, t+n);
2771 t += n;
2772 }
2773 /* t+1 .. lg(U)-1 */
2774 n = lg(U) - t - 1; /* can be 0 */
2775 if (!n) setlg(vU,l); else gel(vU,l) = vecslice(U, t+1, t+n);
2776 return vU;
2777 }
2778
2779 static void
init_zlog_mod(zlog_S * S,GEN bid,GEN mod)2780 init_zlog_mod(zlog_S *S, GEN bid, GEN mod)
2781 {
2782 GEN fa2 = bid_get_fact2(bid);
2783 S->U = bid_get_U(bid);
2784 S->hU = lg(bid_get_cyc(bid))-1;
2785 S->archp = bid_get_archp(bid);
2786 S->sprk = bid_get_sprk(bid);
2787 S->bid = bid;
2788 S->mod = mod;
2789 S->P = gel(fa2,1);
2790 S->k = gel(fa2,2);
2791 S->no2 = lg(S->P) == lg(gel(bid_get_fact(bid),1));
2792 }
2793 void
init_zlog(zlog_S * S,GEN bid)2794 init_zlog(zlog_S *S, GEN bid)
2795 {
2796 return init_zlog_mod(S, bid, NULL);
2797 }
2798
2799 /* a a t_FRAC/t_INT, reduce mod bid */
2800 static GEN
Q_mod_bid(GEN bid,GEN a)2801 Q_mod_bid(GEN bid, GEN a)
2802 {
2803 GEN xZ = gcoeff(bid_get_ideal(bid),1,1);
2804 GEN b = Rg_to_Fp(a, xZ);
2805 if (gsigne(a) < 0) b = subii(b, xZ);
2806 return signe(b)? b: xZ;
2807 }
2808 /* Return decomposition of a on the CRT generators blocks attached to the
2809 * S->sprk and sarch; sgn = sign(a, S->arch), NULL if unknown */
2810 static GEN
zlog(GEN nf,GEN a,GEN sgn,zlog_S * S)2811 zlog(GEN nf, GEN a, GEN sgn, zlog_S *S)
2812 {
2813 long k, l;
2814 GEN y;
2815 a = nf_to_scalar_or_basis(nf, a);
2816 switch(typ(a))
2817 {
2818 case t_INT: break;
2819 case t_FRAC: a = Q_mod_bid(S->bid, a); break;
2820 default: /* case t_COL: */
2821 {
2822 GEN den;
2823 check_nfelt(a, &den);
2824 if (den)
2825 {
2826 a = Q_muli_to_int(a, den);
2827 a = mkmat2(mkcol2(a, den), mkcol2(gen_1, gen_m1));
2828 return famat_zlog(nf, a, sgn, S);
2829 }
2830 }
2831 }
2832 if (sgn)
2833 sgn = (lg(sgn) == 1)? NULL: leafcopy(sgn);
2834 else
2835 sgn = (lg(S->archp) == 1)? NULL: nfsign_arch(nf, a, S->archp);
2836 l = lg(S->sprk);
2837 y = cgetg(sgn? l+1: l, t_COL);
2838 for (k = 1; k < l; k++)
2839 {
2840 GEN sprk = gel(S->sprk,k);
2841 gel(y,k) = log_prk(nf, a, sprk, S->mod);
2842 }
2843 if (sgn) gel(y,l) = Flc_to_ZC(sgn);
2844 return y;
2845 }
2846
2847 /* true nf */
2848 GEN
pr_basis_perm(GEN nf,GEN pr)2849 pr_basis_perm(GEN nf, GEN pr)
2850 {
2851 long f = pr_get_f(pr);
2852 GEN perm;
2853 if (f == nf_get_degree(nf)) return identity_perm(f);
2854 perm = cgetg(f+1, t_VECSMALL);
2855 perm[1] = 1;
2856 if (f > 1)
2857 {
2858 GEN H = pr_hnf(nf,pr);
2859 long i, k;
2860 for (i = k = 2; k <= f; i++)
2861 if (!equali1(gcoeff(H,i,i))) perm[k++] = i;
2862 }
2863 return perm;
2864 }
2865
2866 /* \sum U[i]*y[i], U[i] ZM, y[i] ZC. We allow lg(y) > lg(U). */
2867 static GEN
ZMV_ZCV_mul(GEN U,GEN y)2868 ZMV_ZCV_mul(GEN U, GEN y)
2869 {
2870 long i, l = lg(U);
2871 GEN z = NULL;
2872 if (l == 1) return cgetg(1,t_COL);
2873 for (i = 1; i < l; i++)
2874 {
2875 GEN u = ZM_ZC_mul(gel(U,i), gel(y,i));
2876 z = z? ZC_add(z, u): u;
2877 }
2878 return z;
2879 }
2880 /* A * (U[1], ..., U[d] */
2881 static GEN
ZM_ZMV_mul(GEN A,GEN U)2882 ZM_ZMV_mul(GEN A, GEN U)
2883 {
2884 long i, l;
2885 GEN V = cgetg_copy(U,&l);
2886 for (i = 1; i < l; i++) gel(V,i) = ZM_mul(A,gel(U,i));
2887 return V;
2888 }
2889
2890 /* a = 1 mod pr, sprk mod pr^e, e >= 1 */
2891 static GEN
sprk_log_prk1_2(GEN nf,GEN a,GEN sprk)2892 sprk_log_prk1_2(GEN nf, GEN a, GEN sprk)
2893 {
2894 GEN U1, U2, y, L2 = sprk_get_L2(sprk);
2895 sprk_get_U2(sprk, &U1,&U2);
2896 y = ZM_ZC_mul(U2, log_prk1(nf, a, lg(U2)-1, L2, sprk_get_prk(sprk)));
2897 return vecmodii(y, sprk_get_cyc(sprk));
2898 }
2899 /* assume e >= 2 */
2900 static GEN
sprk_log_gen_pr2(GEN nf,GEN sprk,long e)2901 sprk_log_gen_pr2(GEN nf, GEN sprk, long e)
2902 {
2903 GEN M, G, pr = sprk_get_pr(sprk);
2904 long i, l;
2905 if (e == 2)
2906 {
2907 GEN L2 = sprk_get_L2(sprk), L = gel(L2,1);
2908 G = gel(L,2); l = lg(G);
2909 }
2910 else
2911 {
2912 GEN perm = pr_basis_perm(nf,pr), PI = nfpow_u(nf, pr_get_gen(pr), e-1);
2913 l = lg(perm);
2914 G = cgetg(l, t_VEC);
2915 if (typ(PI) == t_INT)
2916 { /* zk_ei_mul doesn't allow t_INT */
2917 long N = nf_get_degree(nf);
2918 gel(G,1) = addiu(PI,1);
2919 for (i = 2; i < l; i++)
2920 {
2921 GEN z = col_ei(N, 1);
2922 gel(G,i) = z; gel(z, perm[i]) = PI;
2923 }
2924 }
2925 else
2926 {
2927 gel(G,1) = nfadd(nf, gen_1, PI);
2928 for (i = 2; i < l; i++)
2929 gel(G,i) = nfadd(nf, gen_1, zk_ei_mul(nf, PI, perm[i]));
2930 }
2931 }
2932 M = cgetg(l, t_MAT);
2933 for (i = 1; i < l; i++) gel(M,i) = sprk_log_prk1_2(nf, gel(G,i), sprk);
2934 return M;
2935 }
2936 /* Log on bid.gen of generators of P_{1,I pr^{e-1}} / P_{1,I pr^e} (I,pr) = 1,
2937 * defined implicitly via CRT. 'ind' is the index of pr in modulus
2938 * factorization */
2939 GEN
log_gen_pr(zlog_S * S,long ind,GEN nf,long e)2940 log_gen_pr(zlog_S *S, long ind, GEN nf, long e)
2941 {
2942 GEN Uind = gel(S->U, ind);
2943 if (e == 1) retmkmat( gel(Uind,1) );
2944 return ZM_mul(Uind, sprk_log_gen_pr2(nf, gel(S->sprk,ind), e));
2945 }
2946 GEN
sprk_log_gen_pr(GEN nf,GEN sprk,long e)2947 sprk_log_gen_pr(GEN nf, GEN sprk, long e)
2948 {
2949 if (e == 1)
2950 {
2951 long n = lg(sprk_get_cyc(sprk))-1;
2952 retmkmat(col_ei(n, 1));
2953 }
2954 return sprk_log_gen_pr2(nf, sprk, e);
2955 }
2956 /* a = 1 mod pr */
2957 GEN
sprk_log_prk1(GEN nf,GEN a,GEN sprk)2958 sprk_log_prk1(GEN nf, GEN a, GEN sprk)
2959 {
2960 if (lg(sprk) == 5) return mkcol(gen_0); /* mod pr */
2961 return sprk_log_prk1_2(nf, a, sprk);
2962 }
2963 /* Log on bid.gen of generator of P_{1,f} / P_{1,f v[index]}
2964 * v = vector of r1 real places */
2965 GEN
log_gen_arch(zlog_S * S,long index)2966 log_gen_arch(zlog_S *S, long index)
2967 {
2968 GEN U = gel(S->U, lg(S->U)-1);
2969 return gel(U, index);
2970 }
2971
2972 /* compute bid.clgp: [h,cyc] or [h,cyc,gen] */
2973 static GEN
bid_grp(GEN nf,GEN U,GEN cyc,GEN g,GEN F,GEN sarch)2974 bid_grp(GEN nf, GEN U, GEN cyc, GEN g, GEN F, GEN sarch)
2975 {
2976 GEN G, h = ZV_prod(cyc);
2977 long c;
2978 if (!U) return mkvec2(h,cyc);
2979 c = lg(U);
2980 G = cgetg(c,t_VEC);
2981 if (c > 1)
2982 {
2983 GEN U0, Uoo, EX = cyc_get_expo(cyc); /* exponent of bid */
2984 long i, hU = nbrows(U), nba = lg(sarch_get_cyc(sarch))-1; /* #f_oo */
2985 if (!nba) { U0 = U; Uoo = NULL; }
2986 else if (nba == hU) { U0 = NULL; Uoo = U; }
2987 else
2988 {
2989 U0 = rowslice(U, 1, hU-nba);
2990 Uoo = rowslice(U, hU-nba+1, hU);
2991 }
2992 for (i = 1; i < c; i++)
2993 {
2994 GEN t = gen_1;
2995 if (U0) t = famat_to_nf_modideal_coprime(nf, g, gel(U0,i), F, EX);
2996 if (Uoo) t = set_sign_mod_divisor(nf, ZV_to_Flv(gel(Uoo,i),2), t, sarch);
2997 gel(G,i) = t;
2998 }
2999 }
3000 return mkvec3(h, cyc, G);
3001 }
3002
3003 /* remove prime ideals of norm 2 with exponent 1 from factorization */
3004 static GEN
famat_strip2(GEN fa)3005 famat_strip2(GEN fa)
3006 {
3007 GEN P = gel(fa,1), E = gel(fa,2), Q, F;
3008 long l = lg(P), i, j;
3009 Q = cgetg(l, t_COL);
3010 F = cgetg(l, t_COL);
3011 for (i = j = 1; i < l; i++)
3012 {
3013 GEN pr = gel(P,i), e = gel(E,i);
3014 if (!absequaliu(pr_get_p(pr), 2) || itou(e) != 1 || pr_get_f(pr) != 1)
3015 {
3016 gel(Q,j) = pr;
3017 gel(F,j) = e; j++;
3018 }
3019 }
3020 setlg(Q,j);
3021 setlg(F,j); return mkmat2(Q,F);
3022 }
3023
3024 /* Compute [[ideal,arch], [h,[cyc],[gen]], idealfact, [liste], U]
3025 flag may include nf_GEN | nf_INIT */
3026 static GEN
Idealstarmod_i(GEN nf,GEN ideal,long flag,GEN MOD)3027 Idealstarmod_i(GEN nf, GEN ideal, long flag, GEN MOD)
3028 {
3029 long i, k, nbp, R1;
3030 GEN y, cyc, U, u1 = NULL, fa, fa2, sprk, x, arch, archp, E, P, sarch, gen;
3031
3032 nf = checknf(nf);
3033 R1 = nf_get_r1(nf);
3034 if (typ(ideal) == t_VEC && lg(ideal) == 3)
3035 {
3036 arch = gel(ideal,2);
3037 ideal= gel(ideal,1);
3038 switch(typ(arch))
3039 {
3040 case t_VEC:
3041 if (lg(arch) != R1+1)
3042 pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
3043 archp = vec01_to_indices(arch);
3044 break;
3045 case t_VECSMALL:
3046 archp = arch;
3047 k = lg(archp)-1;
3048 if (k && archp[k] > R1)
3049 pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
3050 arch = indices_to_vec01(archp, R1);
3051 break;
3052 default:
3053 pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
3054 return NULL;/*LCOV_EXCL_LINE*/
3055 }
3056 }
3057 else
3058 {
3059 arch = zerovec(R1);
3060 archp = cgetg(1, t_VECSMALL);
3061 }
3062 if (MOD)
3063 {
3064 if (typ(MOD) != t_INT) pari_err_TYPE("bnrinit [incorrect cycmod]", MOD);
3065 if (mpodd(MOD) && lg(archp) != 1)
3066 MOD = shifti(MOD, 1); /* ensure elements of G^MOD are >> 0 */
3067 }
3068 if (is_nf_factor(ideal))
3069 {
3070 fa = ideal;
3071 x = factorbackprime(nf, gel(fa,1), gel(fa,2));
3072 }
3073 else
3074 {
3075 fa = idealfactor(nf, ideal);
3076 x = ideal;
3077 }
3078 if (typ(x) != t_MAT) x = idealhnf_shallow(nf, x);
3079 if (lg(x) == 1) pari_err_DOMAIN("Idealstar", "ideal","=",gen_0,x);
3080 if (typ(gcoeff(x,1,1)) != t_INT)
3081 pari_err_DOMAIN("Idealstar","denominator(ideal)", "!=",gen_1,x);
3082 sarch = nfarchstar(nf, x, archp);
3083 fa2 = famat_strip2(fa);
3084 P = gel(fa2,1);
3085 E = gel(fa2,2);
3086 nbp = lg(P)-1;
3087 sprk = cgetg(nbp+1,t_VEC);
3088 if (nbp)
3089 {
3090 GEN t = (lg(gel(fa,1))==2)? NULL: x; /* beware fa != fa2 */
3091 cyc = cgetg(nbp+2,t_VEC);
3092 gen = cgetg(nbp+1,t_VEC);
3093 for (i = 1; i <= nbp; i++)
3094 {
3095 GEN L = sprkinit(nf, gel(P,i), itou(gel(E,i)), t, MOD);
3096 gel(sprk,i) = L;
3097 gel(cyc,i) = sprk_get_cyc(L);
3098 /* true gens are congruent to those mod x AND positive at archp */
3099 gel(gen,i) = sprk_get_gen(L);
3100 }
3101 gel(cyc,i) = sarch_get_cyc(sarch);
3102 cyc = shallowconcat1(cyc);
3103 gen = shallowconcat1(gen);
3104 cyc = ZV_snf_group(cyc, &U, (flag & nf_GEN)? &u1: NULL);
3105 }
3106 else
3107 {
3108 cyc = sarch_get_cyc(sarch);
3109 gen = cgetg(1,t_VEC);
3110 U = matid(lg(cyc)-1);
3111 if (flag & nf_GEN) u1 = U;
3112 }
3113 y = bid_grp(nf, u1, cyc, gen, x, sarch);
3114 if (!(flag & nf_INIT)) return y;
3115 U = split_U(U, sprk);
3116 return mkvec5(mkvec2(x, arch), y, mkvec2(fa,fa2), mkvec2(sprk, sarch), U);
3117 }
3118 GEN
Idealstarmod(GEN nf,GEN ideal,long flag,GEN MOD)3119 Idealstarmod(GEN nf, GEN ideal, long flag, GEN MOD)
3120 {
3121 pari_sp av = avma;
3122 if (!nf) nf = nfinit(pol_x(0), DEFAULTPREC);
3123 return gerepilecopy(av, Idealstarmod_i(nf, ideal, flag, MOD));
3124 }
3125 GEN
Idealstar(GEN nf,GEN ideal,long flag)3126 Idealstar(GEN nf, GEN ideal, long flag) { return Idealstarmod(nf, ideal, flag, NULL); }
3127 GEN
Idealstarprk(GEN nf,GEN pr,long k,long flag)3128 Idealstarprk(GEN nf, GEN pr, long k, long flag)
3129 {
3130 pari_sp av = avma;
3131 GEN z = Idealstarmod_i(nf, mkmat2(mkcol(pr),mkcols(k)), flag, NULL);
3132 return gerepilecopy(av, z);
3133 }
3134
3135 /* FIXME: obsolete */
3136 GEN
zidealstarinitgen(GEN nf,GEN ideal)3137 zidealstarinitgen(GEN nf, GEN ideal)
3138 { return Idealstar(nf,ideal, nf_INIT|nf_GEN); }
3139 GEN
zidealstarinit(GEN nf,GEN ideal)3140 zidealstarinit(GEN nf, GEN ideal)
3141 { return Idealstar(nf,ideal, nf_INIT); }
3142 GEN
zidealstar(GEN nf,GEN ideal)3143 zidealstar(GEN nf, GEN ideal)
3144 { return Idealstar(nf,ideal, nf_GEN); }
3145
3146 GEN
idealstarmod(GEN nf,GEN ideal,long flag,GEN MOD)3147 idealstarmod(GEN nf, GEN ideal, long flag, GEN MOD)
3148 {
3149 switch(flag)
3150 {
3151 case 0: return Idealstarmod(nf,ideal, nf_GEN, MOD);
3152 case 1: return Idealstarmod(nf,ideal, nf_INIT, MOD);
3153 case 2: return Idealstarmod(nf,ideal, nf_INIT|nf_GEN, MOD);
3154 default: pari_err_FLAG("idealstar");
3155 }
3156 return NULL; /* LCOV_EXCL_LINE */
3157 }
3158 GEN
idealstar0(GEN nf,GEN ideal,long flag)3159 idealstar0(GEN nf, GEN ideal,long flag) { return idealstarmod(nf, ideal, flag, NULL); }
3160
3161 void
check_nfelt(GEN x,GEN * den)3162 check_nfelt(GEN x, GEN *den)
3163 {
3164 long l = lg(x), i;
3165 GEN t, d = NULL;
3166 if (typ(x) != t_COL) pari_err_TYPE("check_nfelt", x);
3167 for (i=1; i<l; i++)
3168 {
3169 t = gel(x,i);
3170 switch (typ(t))
3171 {
3172 case t_INT: break;
3173 case t_FRAC:
3174 if (!d) d = gel(t,2); else d = lcmii(d, gel(t,2));
3175 break;
3176 default: pari_err_TYPE("check_nfelt", x);
3177 }
3178 }
3179 *den = d;
3180 }
3181
3182 GEN
vecmodii(GEN x,GEN y)3183 vecmodii(GEN x, GEN y)
3184 { pari_APPLY_same(modii(gel(x,i), gel(y,i))) }
3185 GEN
ZV_snf_gcd(GEN x,GEN mod)3186 ZV_snf_gcd(GEN x, GEN mod)
3187 { pari_APPLY_same(gcdii(gel(x,i), mod)); }
3188
3189 GEN
vecmoduu(GEN x,GEN y)3190 vecmoduu(GEN x, GEN y)
3191 { pari_APPLY_ulong(uel(x,i) % uel(y,i)) }
3192
3193 /* assume a true bnf and bid */
3194 GEN
ideallog_units0(GEN bnf,GEN bid,GEN MOD)3195 ideallog_units0(GEN bnf, GEN bid, GEN MOD)
3196 {
3197 GEN nf = bnf_get_nf(bnf), D, y, C, cyc;
3198 long j, lU = lg(bnf_get_logfu(bnf)); /* r1+r2 */
3199 zlog_S S;
3200 init_zlog_mod(&S, bid, MOD);
3201 if (!S.hU) return zeromat(0,lU);
3202 cyc = bid_get_cyc(bid);
3203 if (MOD) cyc = ZV_snf_gcd(cyc, MOD);
3204 D = nfsign_fu(bnf, bid_get_archp(bid));
3205 y = cgetg(lU, t_MAT);
3206 if ((C = bnf_build_cheapfu(bnf)))
3207 { for (j = 1; j < lU; j++) gel(y,j) = zlog(nf, gel(C,j), gel(D,j), &S); }
3208 else
3209 {
3210 long i, l = lg(S.U), l0 = lg(S.sprk);
3211 GEN X, U;
3212 if (!(C = bnf_compactfu_mat(bnf))) bnf_build_units(bnf); /* error */
3213 X = gel(C,1); U = gel(C,2);
3214 for (j = 1; j < lU; j++) gel(y,j) = cgetg(l, t_COL);
3215 for (i = 1; i < l0; i++)
3216 {
3217 GEN sprk = gel(S.sprk, i);
3218 GEN Xi = sunits_makecoprime(X, sprk_get_pr(sprk), sprk_get_prk(sprk));
3219 for (j = 1; j < lU; j++)
3220 gcoeff(y,i,j) = famat_zlog_pr_coprime(nf, Xi, gel(U,j), sprk, MOD);
3221 }
3222 if (l0 != l)
3223 for (j = 1; j < lU; j++) gcoeff(y,l0,j) = Flc_to_ZC(gel(D,j));
3224 }
3225 y = vec_prepend(y, zlog(nf, bnf_get_tuU(bnf), nfsign_tu(bnf, S.archp), &S));
3226 for (j = 1; j <= lU; j++)
3227 gel(y,j) = vecmodii(ZMV_ZCV_mul(S.U, gel(y,j)), cyc);
3228 return y;
3229 }
3230 GEN
ideallog_units(GEN bnf,GEN bid)3231 ideallog_units(GEN bnf, GEN bid)
3232 { return ideallog_units0(bnf, bid, NULL); }
3233 GEN
log_prk_units(GEN nf,GEN D,GEN sprk)3234 log_prk_units(GEN nf, GEN D, GEN sprk)
3235 {
3236 GEN L, Ltu = log_prk(nf, gel(D,1), sprk, NULL);
3237 D = gel(D,2);
3238 if (lg(D) != 3 || typ(gel(D,2)) != t_MAT) L = veclog_prk(nf, D, sprk);
3239 else
3240 {
3241 GEN X = gel(D,1), U = gel(D,2);
3242 long j, lU = lg(U);
3243 X = sunits_makecoprime(X, sprk_get_pr(sprk), sprk_get_prk(sprk));
3244 L = cgetg(lU, t_MAT);
3245 for (j = 1; j < lU; j++)
3246 gel(L,j) = famat_zlog_pr_coprime(nf, X, gel(U,j), sprk, NULL);
3247 }
3248 return vec_prepend(L, Ltu);
3249 }
3250
3251 static GEN
ideallog_i(GEN nf,GEN x,zlog_S * S)3252 ideallog_i(GEN nf, GEN x, zlog_S *S)
3253 {
3254 pari_sp av = avma;
3255 GEN y;
3256 if (!S->hU) return cgetg(1, t_COL);
3257 if (typ(x) == t_MAT)
3258 y = famat_zlog(nf, x, NULL, S);
3259 else
3260 y = zlog(nf, x, NULL, S);
3261 y = ZMV_ZCV_mul(S->U, y);
3262 return gerepileupto(av, vecmodii(y, bid_get_cyc(S->bid)));
3263 }
3264 GEN
ideallogmod(GEN nf,GEN x,GEN bid,GEN mod)3265 ideallogmod(GEN nf, GEN x, GEN bid, GEN mod)
3266 {
3267 zlog_S S;
3268 if (!nf)
3269 {
3270 if (mod) pari_err_IMPL("Zideallogmod");
3271 return Zideallog(bid, x);
3272 }
3273 checkbid(bid); init_zlog_mod(&S, bid, mod);
3274 return ideallog_i(checknf(nf), x, &S);
3275 }
3276 GEN
ideallog(GEN nf,GEN x,GEN bid)3277 ideallog(GEN nf, GEN x, GEN bid) { return ideallogmod(nf, x, bid, NULL); }
3278
3279 /*************************************************************************/
3280 /** **/
3281 /** JOIN BID STRUCTURES, IDEAL LISTS **/
3282 /** **/
3283 /*************************************************************************/
3284 /* bid1, bid2: for coprime modules m1 and m2 (without arch. part).
3285 * Output: bid for m1 m2 */
3286 static GEN
join_bid(GEN nf,GEN bid1,GEN bid2)3287 join_bid(GEN nf, GEN bid1, GEN bid2)
3288 {
3289 pari_sp av = avma;
3290 long nbgen, l1,l2;
3291 GEN I1,I2, G1,G2, sprk1,sprk2, cyc1,cyc2, sarch;
3292 GEN sprk, fa,fa2, U, cyc, y, u1 = NULL, x, gen;
3293
3294 I1 = bid_get_ideal(bid1);
3295 I2 = bid_get_ideal(bid2);
3296 if (gequal1(gcoeff(I1,1,1))) return bid2; /* frequent trivial case */
3297 G1 = bid_get_grp(bid1);
3298 G2 = bid_get_grp(bid2);
3299 x = idealmul(nf, I1,I2);
3300 fa = famat_mul_shallow(bid_get_fact(bid1), bid_get_fact(bid2));
3301 fa2= famat_mul_shallow(bid_get_fact2(bid1), bid_get_fact2(bid2));
3302 sprk1 = bid_get_sprk(bid1);
3303 sprk2 = bid_get_sprk(bid2);
3304 sprk = shallowconcat(sprk1, sprk2);
3305
3306 cyc1 = abgrp_get_cyc(G1); l1 = lg(cyc1);
3307 cyc2 = abgrp_get_cyc(G2); l2 = lg(cyc2);
3308 gen = (lg(G1)>3 && lg(G2)>3)? gen_1: NULL;
3309 nbgen = l1+l2-2;
3310 cyc = ZV_snf_group(shallowconcat(cyc1,cyc2), &U, gen? &u1: NULL);
3311 if (nbgen)
3312 {
3313 GEN U1 = bid_get_U(bid1), U2 = bid_get_U(bid2);
3314 U1 = l1==1? const_vec(lg(sprk1), cgetg(1,t_MAT))
3315 : ZM_ZMV_mul(vecslice(U, 1, l1-1), U1);
3316 U2 = l2==1? const_vec(lg(sprk2), cgetg(1,t_MAT))
3317 : ZM_ZMV_mul(vecslice(U, l1, nbgen), U2);
3318 U = shallowconcat(U1, U2);
3319 }
3320 else
3321 U = const_vec(lg(sprk), cgetg(1,t_MAT));
3322
3323 if (gen)
3324 {
3325 GEN uv = zkchinese1init2(nf, I2, I1, x);
3326 gen = shallowconcat(zkVchinese1(gel(uv,1), abgrp_get_gen(G1)),
3327 zkVchinese1(gel(uv,2), abgrp_get_gen(G2)));
3328 }
3329 sarch = bid_get_sarch(bid1); /* trivial */
3330 y = bid_grp(nf, u1, cyc, gen, x, sarch);
3331 x = mkvec2(x, bid_get_arch(bid1));
3332 y = mkvec5(x, y, mkvec2(fa, fa2), mkvec2(sprk, sarch), U);
3333 return gerepilecopy(av,y);
3334 }
3335
3336 typedef struct _ideal_data {
3337 GEN nf, emb, L, pr, prL, sgnU, archp;
3338 } ideal_data;
3339
3340 /* z <- ( z | f(v[i])_{i=1..#v} ) */
3341 static void
concat_join(GEN * pz,GEN v,GEN (* f)(ideal_data *,GEN),ideal_data * data)3342 concat_join(GEN *pz, GEN v, GEN (*f)(ideal_data*,GEN), ideal_data *data)
3343 {
3344 long i, nz, lv = lg(v);
3345 GEN z, Z;
3346 if (lv == 1) return;
3347 z = *pz; nz = lg(z)-1;
3348 *pz = Z = cgetg(lv + nz, typ(z));
3349 for (i = 1; i <=nz; i++) gel(Z,i) = gel(z,i);
3350 Z += nz;
3351 for (i = 1; i < lv; i++) gel(Z,i) = f(data, gel(v,i));
3352 }
3353 static GEN
join_idealinit(ideal_data * D,GEN x)3354 join_idealinit(ideal_data *D, GEN x)
3355 { return join_bid(D->nf, x, D->prL); }
3356 static GEN
join_ideal(ideal_data * D,GEN x)3357 join_ideal(ideal_data *D, GEN x)
3358 { return idealmulpowprime(D->nf, x, D->pr, D->L); }
3359 static GEN
join_unit(ideal_data * D,GEN x)3360 join_unit(ideal_data *D, GEN x)
3361 {
3362 GEN bid = join_idealinit(D, gel(x,1)), u = gel(x,2), v = mkvec(D->emb);
3363 if (lg(u) != 1) v = shallowconcat(u, v);
3364 return mkvec2(bid, v);
3365 }
3366
3367 GEN
log_prk_units_init(GEN bnf)3368 log_prk_units_init(GEN bnf)
3369 {
3370 GEN U = bnf_has_fu(bnf);
3371 if (U) U = matalgtobasis(bnf_get_nf(bnf), U);
3372 else if (!(U = bnf_compactfu_mat(bnf))) (void)bnf_build_units(bnf);
3373 return mkvec2(bnf_get_tuU(bnf), U);
3374 }
3375 /* flag & nf_GEN : generators, otherwise no
3376 * flag &2 : units, otherwise no
3377 * flag &4 : ideals in HNF, otherwise bid
3378 * flag &8 : omit ideals which cannot be conductors (pr^1 with Npr=2) */
3379 static GEN
Ideallist(GEN bnf,ulong bound,long flag)3380 Ideallist(GEN bnf, ulong bound, long flag)
3381 {
3382 const long cond = flag & 8;
3383 const long do_units = flag & 2, big_id = !(flag & 4);
3384 const long istar_flag = (flag & nf_GEN) | nf_INIT;
3385 pari_sp av, av0 = avma;
3386 long i, j;
3387 GEN nf, z, p, fa, id, BOUND, U, empty = cgetg(1,t_VEC);
3388 forprime_t S;
3389 ideal_data ID;
3390 GEN (*join_z)(ideal_data*, GEN) =
3391 do_units? &join_unit
3392 : (big_id? &join_idealinit: &join_ideal);
3393
3394 if (do_units)
3395 {
3396 bnf = checkbnf(bnf);
3397 nf = bnf_get_nf(bnf);
3398 }
3399 else
3400 nf = checknf(bnf);
3401 if ((long)bound <= 0) return empty;
3402 id = matid(nf_get_degree(nf));
3403 if (big_id) id = Idealstar(nf,id, istar_flag);
3404
3405 /* z[i] will contain all "objects" of norm i. Depending on flag, this means
3406 * an ideal, a bid, or a couple [bid, log(units)]. Such objects are stored
3407 * in vectors, computed one primary component at a time; join_z
3408 * reconstructs the global object */
3409 BOUND = utoipos(bound);
3410 z = const_vec(bound, empty);
3411 U = do_units? log_prk_units_init(bnf): NULL;
3412 gel(z,1) = mkvec(U? mkvec2(id, cgetg(1,t_VEC)): id);
3413 ID.nf = nf;
3414
3415 p = cgetipos(3);
3416 u_forprime_init(&S, 2, bound);
3417 av = avma;
3418 while ((p[2] = u_forprime_next(&S)))
3419 {
3420 if (DEBUGLEVEL>1) err_printf("%ld ",p[2]);
3421 fa = idealprimedec_limit_norm(nf, p, BOUND);
3422 for (j = 1; j < lg(fa); j++)
3423 {
3424 GEN pr = gel(fa,j), z2 = leafcopy(z);
3425 ulong Q, q = upr_norm(pr);
3426 long l;
3427 ID.pr = ID.prL = pr;
3428 if (cond && q == 2) { l = 2; Q = 4; } else { l = 1; Q = q; }
3429 for (; Q <= bound; l++, Q *= q) /* add pr^l */
3430 {
3431 ulong iQ;
3432 ID.L = utoipos(l);
3433 if (big_id) {
3434 ID.prL = Idealstarprk(nf, pr, l, istar_flag);
3435 if (U)
3436 ID.emb = Q == 2? cgetg(1,t_VEC)
3437 : log_prk_units(nf, U, gel(bid_get_sprk(ID.prL),1));
3438 }
3439 for (iQ = Q,i = 1; iQ <= bound; iQ += Q,i++)
3440 concat_join(&gel(z,iQ), gel(z2,i), join_z, &ID);
3441 }
3442 }
3443 if (gc_needed(av,1))
3444 {
3445 if(DEBUGMEM>1) pari_warn(warnmem,"Ideallist");
3446 z = gerepilecopy(av, z);
3447 }
3448 }
3449 return gerepilecopy(av0, z);
3450 }
3451 GEN
ideallist0(GEN bnf,long bound,long flag)3452 ideallist0(GEN bnf,long bound, long flag) {
3453 if (flag<0 || flag>15) pari_err_FLAG("ideallist");
3454 return Ideallist(bnf,bound,flag);
3455 }
3456 GEN
ideallist(GEN bnf,long bound)3457 ideallist(GEN bnf,long bound) { return Ideallist(bnf,bound,4); }
3458
3459 /* bid = for module m (without arch. part), arch = archimedean part.
3460 * Output: bid for [m,arch] */
3461 static GEN
join_bid_arch(GEN nf,GEN bid,GEN archp)3462 join_bid_arch(GEN nf, GEN bid, GEN archp)
3463 {
3464 pari_sp av = avma;
3465 GEN G, U;
3466 GEN sprk, cyc, y, u1 = NULL, x, sarch, gen;
3467
3468 checkbid(bid);
3469 G = bid_get_grp(bid);
3470 x = bid_get_ideal(bid);
3471 sarch = nfarchstar(nf, bid_get_ideal(bid), archp);
3472 sprk = bid_get_sprk(bid);
3473
3474 gen = (lg(G)>3)? gel(G,3): NULL;
3475 cyc = diagonal_shallow(shallowconcat(gel(G,2), sarch_get_cyc(sarch)));
3476 cyc = ZM_snf_group(cyc, &U, gen? &u1: NULL);
3477 y = bid_grp(nf, u1, cyc, gen, x, sarch);
3478 U = split_U(U, sprk);
3479 y = mkvec5(mkvec2(x, archp), y, gel(bid,3), mkvec2(sprk, sarch), U);
3480 return gerepilecopy(av,y);
3481 }
3482 static GEN
join_arch(ideal_data * D,GEN x)3483 join_arch(ideal_data *D, GEN x) {
3484 return join_bid_arch(D->nf, x, D->archp);
3485 }
3486 static GEN
join_archunit(ideal_data * D,GEN x)3487 join_archunit(ideal_data *D, GEN x) {
3488 GEN bid = join_arch(D, gel(x,1)), u = gel(x,2), v = mkvec(D->emb);
3489 if (lg(u) != 1) v = shallowconcat(u, v);
3490 return mkvec2(bid, v);
3491 }
3492
3493 /* L from ideallist, add archimedean part */
3494 GEN
ideallistarch(GEN bnf,GEN L,GEN arch)3495 ideallistarch(GEN bnf, GEN L, GEN arch)
3496 {
3497 pari_sp av;
3498 long i, j, l = lg(L), lz;
3499 GEN v, z, V;
3500 ideal_data ID;
3501 GEN (*join_z)(ideal_data*, GEN);
3502
3503 if (typ(L) != t_VEC) pari_err_TYPE("ideallistarch",L);
3504 if (l == 1) return cgetg(1,t_VEC);
3505 z = gel(L,1);
3506 if (typ(z) != t_VEC) pari_err_TYPE("ideallistarch",z);
3507 z = gel(z,1); /* either a bid or [bid,U] */
3508 ID.nf = checknf(bnf);
3509 ID.archp = vec01_to_indices(arch);
3510 if (lg(z) == 3) { /* the latter: do units */
3511 if (typ(z) != t_VEC) pari_err_TYPE("ideallistarch",z);
3512 ID.emb = zm_to_ZM( rowpermute(nfsign_units(bnf,NULL,1), ID.archp) );
3513 join_z = &join_archunit;
3514 } else
3515 join_z = &join_arch;
3516 av = avma; V = cgetg(l, t_VEC);
3517 for (i = 1; i < l; i++)
3518 {
3519 z = gel(L,i); lz = lg(z);
3520 gel(V,i) = v = cgetg(lz,t_VEC);
3521 for (j=1; j<lz; j++) gel(v,j) = join_z(&ID, gel(z,j));
3522 }
3523 return gerepilecopy(av,V);
3524 }
3525