Annotation of OpenXM_contrib/pari-2.2/src/basemath/base3.c, Revision 1.2
1.2 ! noro 1: /* $Id: base3.c,v 1.87 2002/09/08 17:09:39 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /*******************************************************************/
17: /* */
18: /* BASIC NF OPERATIONS */
19: /* */
20: /*******************************************************************/
21: #include "pari.h"
1.2 ! noro 22: extern GEN u_FpM_deplin(GEN x, ulong p);
! 23: extern GEN u_FpM_inv(GEN a, ulong p);
1.1 noro 24: extern GEN famat_ideallog(GEN nf, GEN g, GEN e, GEN bid);
25: extern GEN famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid);
26: extern GEN vconcat(GEN A, GEN B);
27:
28: /*******************************************************************/
29: /* */
30: /* OPERATIONS OVER NUMBER FIELD ELEMENTS. */
31: /* These are always represented as column vectors over the */
32: /* integral basis nf[7] */
33: /* */
34: /*******************************************************************/
35:
36: int
37: isnfscalar(GEN x)
38: {
39: long lx=lg(x),i;
40:
1.2 ! noro 41: if (typ(x) != t_COL) return 0;
1.1 noro 42: for (i=2; i<lx; i++)
43: if (!gcmp0((GEN)x[i])) return 0;
44: return 1;
45: }
46:
47: static GEN
48: scal_mul(GEN nf, GEN x, GEN y, long ty)
49: {
1.2 ! noro 50: gpmem_t av=avma, tetpil;
1.1 noro 51: GEN p1;
52:
53: if (!is_extscalar_t(ty))
54: {
55: if (ty!=t_COL) err(typeer,"nfmul");
56: y = gmul((GEN)nf[7],y);
57: }
58: p1 = gmul(x,y); tetpil=avma;
59: return gerepile(av,tetpil,algtobasis(nf,p1));
60: }
61:
62: static GEN
63: get_tab(GEN nf, long *N)
64: {
65: GEN tab = (typ(nf) == t_MAT)? nf: (GEN)nf[9];
66: *N = lg(tab[1])-1; return tab;
67: }
68:
69: GEN
1.2 ! noro 70: mul_by_tab(GEN tab, GEN x, GEN y)
1.1 noro 71: {
1.2 ! noro 72: gpmem_t av;
! 73: long i,j,k,N;
! 74: GEN s,v,c,p1;
1.1 noro 75:
1.2 ! noro 76: N = lg(x)-1;
1.1 noro 77: v=cgetg(N+1,t_COL); av=avma;
78: for (k=1; k<=N; k++)
79: {
80: if (k == 1)
81: s = gmul((GEN)x[1],(GEN)y[1]);
82: else
83: s = gadd(gmul((GEN)x[1],(GEN)y[k]),
84: gmul((GEN)x[k],(GEN)y[1]));
85: for (i=2; i<=N; i++)
86: {
87: c=gcoeff(tab,k,(i-1)*N+i);
88: if (signe(c))
89: {
90: p1 = gmul((GEN)x[i],(GEN)y[i]);
91: if (!gcmp1(c)) p1 = gmul(p1,c);
92: s = gadd(s, p1);
93: }
94: for (j=i+1; j<=N; j++)
95: {
96: c=gcoeff(tab,k,(i-1)*N+j);
97: if (signe(c))
98: {
99: p1 = gadd(gmul((GEN)x[i],(GEN)y[j]),
100: gmul((GEN)x[j],(GEN)y[i]));
101: if (!gcmp1(c)) p1 = gmul(p1,c);
102: s = gadd(s, p1);
103: }
104: }
105: }
106: v[k]=(long)gerepileupto(av,s); av=avma;
107: }
108: return v;
109: }
110:
1.2 ! noro 111: /* product of x and y in nf */
! 112: GEN
! 113: element_mul(GEN nf, GEN x, GEN y)
! 114: {
! 115: long N,tx,ty;
! 116: GEN tab;
! 117:
! 118: if (x == y) return element_sqr(nf,x);
! 119:
! 120: tx=typ(x); ty=typ(y);
! 121: nf=checknf(nf);
! 122: if (tx==t_POLMOD) x=checknfelt_mod(nf,x,"element_mul");
! 123: if (ty==t_POLMOD) y=checknfelt_mod(nf,y,"element_mul");
! 124: if (is_extscalar_t(tx)) return scal_mul(nf,x,y,ty);
! 125: if (is_extscalar_t(ty)) return scal_mul(nf,y,x,tx);
! 126: if (tx != t_COL || ty != t_COL) err(typeer,"element_mul");
! 127: if (isnfscalar(x)) return gmul((GEN)x[1],y);
! 128: if (isnfscalar(y)) return gmul((GEN)y[1],x);
! 129:
! 130: tab = get_tab(nf, &N);
! 131: return mul_by_tab(tab,x,y);
! 132: }
! 133:
1.1 noro 134: /* inverse of x in nf */
135: GEN
136: element_inv(GEN nf, GEN x)
137: {
1.2 ! noro 138: gpmem_t av=avma;
! 139: long i,N,tx=typ(x);
1.1 noro 140: GEN p1,p;
141:
142: nf=checknf(nf); N=degpol(nf[1]);
143: if (is_extscalar_t(tx))
144: {
1.2 ! noro 145: if (tx==t_POLMOD) (void)checknfelt_mod(nf,x,"element_inv");
1.1 noro 146: else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
147: return gerepileupto(av, algtobasis(nf, ginv(x)));
148: }
149: if (isnfscalar(x))
150: {
151: p1=cgetg(N+1,t_COL); p1[1]=linv((GEN)x[1]);
152: for (i=2; i<=N; i++) p1[i]=lcopy((GEN)x[i]);
153: return p1;
154: }
1.2 ! noro 155: if (tx != t_COL) err(typeer,"element_inv");
1.1 noro 156: p = NULL;
157: for (i=1; i<=N; i++)
158: if (typ(x[i])==t_INTMOD)
159: {
160: p = gmael(x,i,1);
161: x = lift(x); break;
162: }
1.2 ! noro 163: p1 = QX_invmod(gmul((GEN)nf[7],x), (GEN)nf[1]);
! 164: p1 = algtobasis_i(nf,p1);
1.1 noro 165:
166: if (p) p1 = FpV(p1, p);
167: return gerepileupto(av,p1);
168: }
169:
170: /* quotient of x and y in nf */
171: GEN
172: element_div(GEN nf, GEN x, GEN y)
173: {
1.2 ! noro 174: gpmem_t av=avma;
! 175: long i,N,tx=typ(x),ty=typ(y);
1.1 noro 176: GEN p1,p;
177:
178: nf=checknf(nf); N=degpol(nf[1]);
1.2 ! noro 179: if (tx==t_POLMOD) (void)checknfelt_mod(nf,x,"element_div");
1.1 noro 180: else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
181:
1.2 ! noro 182: if (ty==t_POLMOD) (void)checknfelt_mod(nf,y,"element_div");
1.1 noro 183: else if (ty==t_POL) y=gmodulcp(y,(GEN)nf[1]);
184:
185: if (is_extscalar_t(tx))
186: {
187: if (is_extscalar_t(ty)) p1=gdiv(x,y);
188: else
189: {
190: if (ty!=t_COL) err(typeer,"nfdiv");
191: p1=gdiv(x,gmodulcp(gmul((GEN)nf[7],y),(GEN)nf[1]));
192: }
193: return gerepileupto(av, algtobasis(nf,p1));
194: }
195: if (is_extscalar_t(ty))
196: {
197: if (tx!=t_COL) err(typeer,"nfdiv");
198: p1=gdiv(gmodulcp(gmul((GEN)nf[7],x),(GEN)nf[1]),y);
199: return gerepileupto(av, algtobasis(nf,p1));
200: }
1.2 ! noro 201: if (tx != t_COL || ty != t_COL) err(typeer,"element_div");
1.1 noro 202:
203: if (isnfscalar(y)) return gdiv(x,(GEN)y[1]);
204: if (isnfscalar(x))
205: {
206: p1=element_inv(nf,y);
207: return gerepileupto(av, gmul((GEN)x[1],p1));
208: }
209:
210: p = NULL;
211: for (i=1; i<=N; i++)
212: if (typ(x[i])==t_INTMOD)
213: {
214: p = gmael(x,i,1);
215: x = lift(x); break;
216: }
217: for (i=1; i<=N; i++)
218: if (typ(y[i])==t_INTMOD)
219: {
220: p1 = gmael(y,i,1);
221: if (p && !egalii(p,p1))
222: err(talker,"inconsistant prime moduli in element_inv");
223: y = lift(y); break;
224: }
225:
1.2 ! noro 226: p1 = gmul(gmul((GEN)nf[7],x), QX_invmod(gmul((GEN)nf[7],y), (GEN)nf[1]));
! 227: p1 = algtobasis_i(nf, gres(p1, (GEN)nf[1]));
1.1 noro 228: if (p) p1 = FpV(p1,p);
229: return gerepileupto(av,p1);
230: }
231:
232: /* product of INTEGERS (i.e vectors with integral coeffs) x and y in nf */
233: GEN
234: element_muli(GEN nf, GEN x, GEN y)
235: {
1.2 ! noro 236: gpmem_t av;
! 237: long i,j,k,N;
1.1 noro 238: GEN p1,s,v,c,tab;
239:
240: tab = get_tab(nf, &N);
1.2 ! noro 241: if (typ(x) != t_COL || lg(x) != N+1
! 242: || typ(y) != t_COL || lg(y) != N+1) err(typeer,"element_muli");
1.1 noro 243: v=cgetg(N+1,t_COL); av=avma;
244: for (k=1; k<=N; k++)
245: {
246: if (k == 1)
247: s = mulii((GEN)x[1],(GEN)y[1]);
248: else
249: s = addii(mulii((GEN)x[1],(GEN)y[k]),
250: mulii((GEN)x[k],(GEN)y[1]));
251: for (i=2; i<=N; i++)
252: {
253: c=gcoeff(tab,k,(i-1)*N+i);
254: if (signe(c))
255: {
256: p1 = mulii((GEN)x[i],(GEN)y[i]);
257: if (!gcmp1(c)) p1 = mulii(p1,c);
258: s = addii(s,p1);
259: }
260: for (j=i+1; j<=N; j++)
261: {
262: c=gcoeff(tab,k,(i-1)*N+j);
263: if (signe(c))
264: {
265: p1 = addii(mulii((GEN)x[i],(GEN)y[j]),
266: mulii((GEN)x[j],(GEN)y[i]));
267: if (!gcmp1(c)) p1 = mulii(p1,c);
268: s = addii(s,p1);
269: }
270: }
271: }
1.2 ! noro 272: v[k] = lpileuptoint(av,s); av = avma;
1.1 noro 273: }
274: return v;
275: }
276:
277: /* product of INTEGERS (i.e vectors with integral coeffs) x and y in nf */
278: GEN
279: element_sqri(GEN nf, GEN x)
280: {
281: GEN p1,s,v,c,tab;
1.2 ! noro 282: gpmem_t av;
! 283: long i,j,k,N;
1.1 noro 284:
285: tab = get_tab(nf, &N);
286:
287: v=cgetg(N+1,t_COL); av=avma;
288: for (k=1; k<=N; k++)
289: {
290: if (k == 1)
291: s = sqri((GEN)x[1]);
292: else
293: s = shifti(mulii((GEN)x[1],(GEN)x[k]), 1);
294: for (i=2; i<=N; i++)
295: {
296: c=gcoeff(tab,k,(i-1)*N+i);
297: if (signe(c))
298: {
299: p1 = sqri((GEN)x[i]);
300: if (!gcmp1(c)) p1 = mulii(p1,c);
301: s = addii(s,p1);
302: }
303: for (j=i+1; j<=N; j++)
304: {
305: c=gcoeff(tab,k,(i-1)*N+j);
306: if (signe(c))
307: {
308: p1 = shifti(mulii((GEN)x[i],(GEN)x[j]),1);
309: if (!gcmp1(c)) p1 = mulii(p1,c);
310: s = addii(s,p1);
311: }
312: }
313: }
1.2 ! noro 314: v[k] = lpileuptoint(av,s); av = avma;
1.1 noro 315: }
316: return v;
317: }
318:
319: GEN
1.2 ! noro 320: sqr_by_tab(GEN tab, GEN x)
1.1 noro 321: {
1.2 ! noro 322: gpmem_t av = avma;
! 323: long i,j,k,N;
! 324: GEN p1,s,v,c;
1.1 noro 325:
1.2 ! noro 326: N = lg(x)-1;
1.1 noro 327: if (isnfscalar(x))
328: {
329: s=cgetg(N+1,t_COL); s[1]=lsqr((GEN)x[1]);
330: for (i=2; i<=N; i++) s[i]=lcopy((GEN)x[i]);
331: return s;
332: }
333: v=cgetg(N+1,t_COL); av = avma;
334: for (k=1; k<=N; k++)
335: {
336: if (k == 1)
337: s = gsqr((GEN)x[1]);
338: else
339: s = gmul2n(gmul((GEN)x[1],(GEN)x[k]), 1);
340: for (i=2; i<=N; i++)
341: {
342: c=gcoeff(tab,k,(i-1)*N+i);
343: if (signe(c))
344: {
345: p1 = gsqr((GEN)x[i]);
346: if (!gcmp1(c)) p1 = gmul(p1,c);
347: s = gadd(s,p1);
348: }
349: for (j=i+1; j<=N; j++)
350: {
351: c=gcoeff(tab,k,(i-1)*N+j);
352: if (signe(c))
353: {
354: p1 = gmul((GEN)x[i],(GEN)x[j]);
1.2 ! noro 355: p1 = gmul2n(p1,1);
! 356: if (!gcmp1(c)) p1 = gmul(p1, c);
1.1 noro 357: s = gadd(s,p1);
358: }
359: }
360: }
361: v[k]=(long)gerepileupto(av,s); av=avma;
362: }
363: return v;
364: }
365:
1.2 ! noro 366: /* square of x in nf */
! 367: GEN
! 368: element_sqr(GEN nf, GEN x)
! 369: {
! 370: long N, tx = typ(x);
! 371: GEN tab;
! 372:
! 373: nf = checknf(nf);
! 374: if (tx==t_POLMOD) x=checknfelt_mod(nf,x,"element_sqr");
! 375: if (is_extscalar_t(tx))
! 376: {
! 377: gpmem_t av = avma;
! 378: return gerepileupto(av, algtobasis(nf, gsqr(x)));
! 379: }
! 380: if (tx != t_COL) err(typeer,"element_sqr");
! 381:
! 382: tab = get_tab(nf, &N);
! 383: return sqr_by_tab(tab,x);
! 384: }
! 385:
! 386: static GEN
! 387: _mul(void *data, GEN x, GEN y) { return element_mul((GEN)data,x,y); }
! 388: static GEN
! 389: _sqr(void *data, GEN x) { return element_sqr((GEN)data,x); }
! 390:
1.1 noro 391: /* Compute x^n in nf, left-shift binary powering */
392: GEN
393: element_pow(GEN nf, GEN x, GEN n)
394: {
1.2 ! noro 395: gpmem_t av = avma;
! 396: long s,N;
! 397: GEN y, cx;
1.1 noro 398:
399: if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
400: nf=checknf(nf); N=degpol(nf[1]);
401: s=signe(n); if (!s) return gscalcol_i(gun,N);
1.2 ! noro 402: if (typ(x) != t_COL)
! 403: {
! 404: x = algtobasis(nf,x);
! 405: if (typ(x) != t_COL) err(typeer,"element_pow");
! 406: }
1.1 noro 407:
408: if (isnfscalar(x))
409: {
410: y = gscalcol_i(gun,N);
411: y[1] = (long)powgi((GEN)x[1],n); return y;
412: }
1.2 ! noro 413: x = primitive_part(x, &cx);
! 414: y = leftright_pow(x, n, (void*)nf, _sqr, _mul);
1.1 noro 415: if (s<0) y = element_inv(nf, y);
416: if (cx) y = gmul(y, powgi(cx, n));
417: return av==avma? gcopy(y): gerepileupto(av,y);
418: }
419:
1.2 ! noro 420: typedef struct {
! 421: GEN nf, p;
! 422: long I;
! 423: } eltmod_muldata;
! 424:
! 425: static GEN
! 426: _mulidmod(void *data, GEN x, GEN y)
! 427: {
! 428: eltmod_muldata *D = (eltmod_muldata*)data;
! 429: (void)y; /* ignored */
! 430: return FpV_red(element_mulid(D->nf, x, D->I), D->p);
! 431: }
! 432: static GEN
! 433: _sqrmod(void *data, GEN x)
1.1 noro 434: {
1.2 ! noro 435: eltmod_muldata *D = (eltmod_muldata*)data;
! 436: return FpV_red(element_sqri(D->nf, x), D->p);
1.1 noro 437: }
438:
439: /* x = I-th vector of the Z-basis of Z_K, in Z^n, compute lift(x^n mod p) */
440: GEN
441: element_powid_mod_p(GEN nf, long I, GEN n, GEN p)
442: {
1.2 ! noro 443: gpmem_t av = avma;
! 444: eltmod_muldata D;
! 445: long s,N;
! 446: GEN y;
1.1 noro 447:
448: if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
449: nf=checknf(nf); N=degpol(nf[1]);
450: s=signe(n);
451: if (!s || I == 1) return gscalcol_i(gun,N);
452: y = zerocol(N); y[I] = un;
1.2 ! noro 453: D.nf = nf;
! 454: D.p = p;
! 455: D.I = I;
! 456: y = leftright_pow(y,n, (void*)&D, &_sqrmod, &_mulidmod);
! 457: if (s < 0) y = FpV_red(element_inv(nf,y), p);
1.1 noro 458: return av==avma? gcopy(y): gerepileupto(av,y);
459: }
460:
1.2 ! noro 461: GEN
! 462: element_mulidid(GEN nf, long i, long j)
! 463: {
! 464: long N;
! 465: GEN tab = get_tab(nf, &N);
! 466: tab += (i-1)*N; return (GEN)tab[j];
! 467: }
! 468:
1.1 noro 469: /* Outputs x.w_i, where w_i is the i-th elt of the integral basis */
470: GEN
471: element_mulid(GEN nf, GEN x, long i)
472: {
1.2 ! noro 473: gpmem_t av;
! 474: long j,k,N;
1.1 noro 475: GEN s,v,c,p1, tab;
476:
477: if (i==1) return gcopy(x);
478: tab = get_tab(nf, &N);
479: if (typ(x) != t_COL || lg(x) != N+1) err(typeer,"element_mulid");
480: tab += (i-1)*N;
481: v=cgetg(N+1,t_COL); av=avma;
482: for (k=1; k<=N; k++)
483: {
484: s = gzero;
485: for (j=1; j<=N; j++)
486: if (signe(c = gcoeff(tab,k,j)) && !gcmp0(p1 = (GEN)x[j]))
487: {
488: if (!gcmp1(c)) p1 = gmul(p1,c);
489: s = gadd(s,p1);
490: }
491: v[k]=lpileupto(av,s); av=avma;
492: }
493: return v;
494: }
495:
1.2 ! noro 496: /* table of multiplication by wi in ZK = Z[w1,..., wN] */
! 497: GEN
! 498: eltmulid_get_table(GEN nf, long i)
! 499: {
! 500: long k,N;
! 501: GEN m, tab;
! 502:
! 503: tab = get_tab(nf, &N);
! 504: tab += (i-1)*N;
! 505: m = cgetg(N+1,t_COL);
! 506: for (k=1; k<=N; k++) m[k] = tab[k];
! 507: return m;
! 508: }
! 509:
! 510: /* table of multiplication by x in ZK = Z[w1,..., wN] */
! 511: GEN
! 512: eltmul_get_table(GEN nf, GEN x)
! 513: {
! 514: long i, N = degpol(nf[1]);
! 515: GEN mul = cgetg(N+1,t_MAT);
! 516: if (typ(x) != t_COL) x = algtobasis(nf,x);
! 517: mul[1] = (long)x; /* assume w_1 = 1 */
! 518: for (i=2; i<=N; i++) mul[i] = (long)element_mulid(nf,x,i);
! 519: return mul;
! 520: }
! 521:
1.1 noro 522: /* valuation of integer x, with resp. to prime ideal P above p.
523: * p.P^(-1) = b Z_K, v >= val_p(norm(x)), and N = deg(nf)
524: * [b may be given as the 'multiplication by b' matrix]
525: */
526: long
1.2 ! noro 527: int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx)
1.1 noro 528: {
529: long i,k,w, N = degpol(nf[1]);
530: GEN r,a,y,mul;
531:
1.2 ! noro 532: mul = (typ(b) == t_MAT)? b: eltmul_get_table(nf, b);
1.1 noro 533: y = cgetg(N+1, t_COL); /* will hold the new x */
534: x = dummycopy(x);
1.2 ! noro 535: for(w=0;; w++)
1.1 noro 536: {
537: for (i=1; i<=N; i++)
538: { /* compute (x.b)_i */
539: a = mulii((GEN)x[1], gcoeff(mul,i,1));
540: for (k=2; k<=N; k++) a = addii(a, mulii((GEN)x[k], gcoeff(mul,i,k)));
541: /* is it divisible by p ? */
542: y[i] = ldvmdii(a,p,&r);
1.2 ! noro 543: if (signe(r))
! 544: {
! 545: if (newx) *newx = x;
! 546: return w;
! 547: }
1.1 noro 548: }
549: r=x; x=y; y=r;
550: }
551: }
552:
553: long
554: element_val(GEN nf, GEN x, GEN vp)
555: {
1.2 ! noro 556: gpmem_t av = avma;
! 557: long N,w,vcx,e;
1.1 noro 558: GEN cx,p;
559:
560: if (gcmp0(x)) return VERYBIGINT;
561: nf=checknf(nf); N=degpol(nf[1]);
562: checkprimeid(vp); p=(GEN)vp[1]; e=itos((GEN)vp[3]);
563: switch(typ(x))
564: {
565: case t_INT: case t_FRAC: case t_FRACN:
566: return ggval(x,p)*e;
567: case t_POLMOD: x = (GEN)x[2]; /* fall through */
568: case t_POL:
1.2 ! noro 569: x = algtobasis_i(nf,x); break;
1.1 noro 570: case t_COL:
571: if (lg(x)==N+1) break;
572: default: err(typeer,"element_val");
573: }
574: if (isnfscalar(x)) return ggval((GEN)x[1],p)*e;
575:
576: cx = content(x);
577: if (gcmp1(cx)) vcx=0; else { x = gdiv(x,cx); vcx = ggval(cx,p); }
1.2 ! noro 578: w = int_elt_val(nf,x,p,(GEN)vp[5],NULL);
1.1 noro 579: avma=av; return w + vcx*e;
580: }
581:
582: /* polegal without comparing variables */
583: long
584: polegal_spec(GEN x, GEN y)
585: {
586: long i = lgef(x);
587:
588: if (i != lgef(y)) return 0;
589: for (i--; i > 1; i--)
590: if (!gegal((GEN)x[i],(GEN)y[i])) return 0;
591: return 1;
592: }
593:
594: GEN
595: basistoalg(GEN nf, GEN x)
596: {
1.2 ! noro 597: long tx=typ(x),lx=lg(x),i,j,l;
1.1 noro 598: GEN z;
599:
600: nf=checknf(nf);
601: switch(tx)
602: {
603: case t_COL:
604: for (i=1; i<lx; i++)
605: {
606: long t = typ(x[i]);
607: if (is_matvec_t(t)) break;
608: }
609: if (i==lx)
610: {
611: z = cgetg(3,t_POLMOD);
612: z[1] = lcopy((GEN)nf[1]);
613: z[2] = lmul((GEN)nf[7],x); return z;
614: }
615: /* fall through */
616:
1.2 ! noro 617: case t_VEC: z=cgetg(lx,tx);
1.1 noro 618: for (i=1; i<lx; i++) z[i]=(long)basistoalg(nf,(GEN)x[i]);
619: return z;
1.2 ! noro 620: case t_MAT: z=cgetg(lx,t_MAT);
! 621: if (lx == 1) return z;
! 622: l = lg(x[1]);
! 623: for (j=1; j<lx; j++)
! 624: {
! 625: z[j] = lgetg(l,t_COL);
! 626: for (i=1; i<l; i++) coeff(z,i,j) = (long)basistoalg(nf,gcoeff(x,i,j));
! 627: }
! 628: return z;
1.1 noro 629:
630: case t_POLMOD:
631: if (!polegal_spec((GEN)nf[1],(GEN)x[1]))
632: err(talker,"not the same number field in basistoalg");
633: return gcopy(x);
634: default: z=cgetg(3,t_POLMOD);
635: z[1]=lcopy((GEN)nf[1]);
636: z[2]=lmul(x,polun[varn(nf[1])]); return z;
637: }
638: }
639:
1.2 ! noro 640: /* gmul(A, pol_to_vec(x)), A t_MAT (or t_VEC) of compatible dimensions */
1.1 noro 641: GEN
642: mulmat_pol(GEN A, GEN x)
643: {
644: long i,l;
645: GEN z;
1.2 ! noro 646: if (typ(x) != t_POL) return gmul(x,(GEN)A[1]); /* scalar */
! 647: l=lgef(x)-1; if (l == 1) return typ(A)==t_VEC? gzero: zerocol(lg(A[1])-1);
1.1 noro 648: x++; z = gmul((GEN)x[1], (GEN)A[1]);
649: for (i=2; i<l ; i++)
650: if (!gcmp0((GEN)x[i])) z = gadd(z, gmul((GEN)x[i], (GEN)A[i]));
651: return z;
652: }
653:
654: /* valid for scalars and polynomial, degree less than N.
1.2 ! noro 655: * No garbage collecting. No check. */
1.1 noro 656: GEN
1.2 ! noro 657: algtobasis_i(GEN nf, GEN x)
1.1 noro 658: {
659: GEN P = (GEN)nf[1];
660: long tx = typ(x), N = degpol(P);
661:
662: if (tx==t_POLMOD) { x=(GEN)x[2]; tx=typ(x); }
663: if (tx==t_POL)
664: {
665: if (varn(x) != varn(P))
666: err(talker,"incompatible variables in algtobasis");
667: if (degpol(x) >= N) x = gres(x,P);
668: return mulmat_pol((GEN)nf[8], x);
669: }
670: return gscalcol(x,N);
671: }
672:
673: GEN
674: algtobasis(GEN nf, GEN x)
675: {
1.2 ! noro 676: long tx=typ(x),lx=lg(x),i,N;
! 677: gpmem_t av=avma;
1.1 noro 678: GEN z;
679:
680: nf=checknf(nf);
681: switch(tx)
682: {
683: case t_VEC: case t_COL: case t_MAT:
684: z=cgetg(lx,tx);
685: for (i=1; i<lx; i++) z[i]=(long)algtobasis(nf,(GEN)x[i]);
686: return z;
687: case t_POLMOD:
688: if (!polegal_spec((GEN)nf[1],(GEN)x[1]))
689: err(talker,"not the same number field in algtobasis");
690: x = (GEN)x[2]; /* fall through */
691: case t_POL:
1.2 ! noro 692: return gerepileupto(av,algtobasis_i(nf,x));
1.1 noro 693:
694: default: N=degpol(nf[1]); return gscalcol(x,N);
695: }
696: }
697:
698: /* Given a and b in nf, gives an algebraic integer y in nf such that a-b.y
699: * is "small"
700: */
701: GEN
702: nfdiveuc(GEN nf, GEN a, GEN b)
703: {
1.2 ! noro 704: gpmem_t av=avma, tetpil;
1.1 noro 705: a = element_div(nf,a,b); tetpil=avma;
706: return gerepile(av,tetpil,ground(a));
707: }
708:
709: /* Given a and b in nf, gives a "small" algebraic integer r in nf
710: * of the form a-b.y
711: */
712: GEN
713: nfmod(GEN nf, GEN a, GEN b)
714: {
1.2 ! noro 715: gpmem_t av=avma,tetpil;
1.1 noro 716: GEN p1=gneg_i(element_mul(nf,b,ground(element_div(nf,a,b))));
717: tetpil=avma; return gerepile(av,tetpil,gadd(a,p1));
718: }
719:
720: /* Given a and b in nf, gives a two-component vector [y,r] in nf such
721: * that r=a-b.y is "small".
722: */
723: GEN
724: nfdivres(GEN nf, GEN a, GEN b)
725: {
1.2 ! noro 726: gpmem_t av=avma,tetpil;
1.1 noro 727: GEN p1,z, y = ground(element_div(nf,a,b));
728:
729: p1=gneg_i(element_mul(nf,b,y)); tetpil=avma;
730: z=cgetg(3,t_VEC); z[1]=lcopy(y); z[2]=ladd(a,p1);
731: return gerepile(av,tetpil,z);
732: }
733:
734: /*************************************************************************/
735: /** **/
736: /** (Z_K/I)^* **/
737: /** **/
738: /*************************************************************************/
1.2 ! noro 739: /* return sign(sigma_k(x)), x t_COL (integral, primitive) */
! 740: static long
! 741: eval_sign(GEN M, GEN x, long k)
! 742: {
! 743: long i, l = lg(x);
! 744: GEN z = mpmul(gcoeff(M,k,1), (GEN)x[1]);
! 745: for (i = 2; i < l; i++)
! 746: z = mpadd(z, mpmul(gcoeff(M,k,i), (GEN)x[i]));
! 747: if (lg(z) < DEFAULTPREC) err(precer,"zsigne");
! 748: return signe(z);
! 749: }
! 750:
! 751: static void
! 752: const_sign(GEN v, GEN arch, GEN s)
! 753: {
! 754: long i,j, l = lg(v);
! 755: for (j=i=1; i<l; i++)
! 756: if (signe(arch[i])) v[j++] = (long)s;
! 757: setlg(v, j);
! 758: }
1.1 noro 759:
760: /* return (column) vector of R1 signatures of x (coeff modulo 2)
1.2 ! noro 761: * if arch = NULL, assume arch = [0,..0] */
1.1 noro 762: GEN
763: zsigne(GEN nf,GEN x,GEN arch)
764: {
1.2 ! noro 765: GEN _0, _1, M, vecsign;
! 766: long i,j,l,s;
! 767: gpmem_t av0, av;
1.1 noro 768:
769: if (!arch) return cgetg(1,t_COL);
770: av0 = avma;
771: l = lg(arch); vecsign = cgetg(l,t_COL);
772: _0 = gmodulss(0,2);
773: _1 = gmodulss(1,2); av = avma;
1.2 ! noro 774: nf = checknf(nf);
! 775: M = gmael(nf,5,1);
1.1 noro 776: switch(typ(x))
777: {
778: case t_MAT: /* factorisation */
779: {
1.2 ! noro 780: GEN g = (GEN)x[1], e = (GEN)x[2];
! 781: const_sign(vecsign, arch, _0);
! 782: if (lg(vecsign)==1) { avma = av0; return cgetg(1, t_COL); }
! 783: for (i=1; i<lg(g); i++)
! 784: if (mpodd((GEN)e[i]))
! 785: vecsign = gadd(vecsign, zsigne(nf,(GEN)g[i],arch));
1.1 noro 786: return gerepileupto(av0, vecsign);
787: }
1.2 ! noro 788: case t_POLMOD: x = (GEN)x[2]; /* fall through */
! 789: case t_POL: x = algtobasis(nf, x); /* fall through */
! 790: case t_COL: if (!isnfscalar(x)) break;
! 791: x = (GEN)x[1]; /* fall through */
! 792: case t_INT: case t_FRAC:
! 793: {
! 794: s = gsigne(x); if (!s) err(talker,"zero element in zsigne");
! 795: const_sign(vecsign, arch, (s > 0)? _0: _1);
! 796: return vecsign;
! 797: }
1.1 noro 798: }
1.2 ! noro 799: x = Q_primpart(x);
1.1 noro 800: for (j=i=1; i<l; i++)
801: if (signe(arch[i]))
1.2 ! noro 802: vecsign[j++] = (eval_sign(M, x, i) > 0)? (long)_0: (long)_1;
1.1 noro 803: avma = av; setlg(vecsign,j); return vecsign;
804: }
805:
1.2 ! noro 806: /* return the t_COL vector of signs of x; the matrix of such if x is a vector
! 807: * of nf elements */
! 808: GEN
! 809: zsigns(GEN nf, GEN x)
! 810: {
! 811: long r1, i, l;
! 812: GEN arch, S;
! 813:
! 814: nf = checknf(nf); r1 = nf_get_r1(nf);
! 815: arch = cgetg(r1+1,t_VEC); for (i=1; i<=r1; i++) arch[i] = un;
! 816: if (typ(x) != t_VEC) return zsigne(nf, x, arch);
! 817: l = lg(x); S = cgetg(l, t_MAT);
! 818: for (i=1; i<l; i++) S[i] = (long)zsigne(nf, (GEN)x[i], arch);
! 819: return S;
! 820: }
! 821:
1.1 noro 822: GEN
823: lllreducemodmatrix(GEN x,GEN y)
824: {
1.2 ! noro 825: gpmem_t av = avma;
! 826: GEN z = lllint_ip(y,4);
1.1 noro 827: return gerepileupto(av, reducemodinvertible(x, z));
828: }
829:
830: /* for internal use...reduce x modulo (invertible) y */
831: GEN
832: reducemodinvertible(GEN x, GEN y)
833: {
834: return gadd(x, gneg(gmul(y, ground(gauss(y,x)))));
835: }
836:
837: /* Reduce column x modulo y in HNF */
838: GEN
1.2 ! noro 839: colreducemodHNF(GEN x, GEN y, GEN *Q)
1.1 noro 840: {
841: long i, l = lg(x);
842: GEN q;
843:
844: if (Q) *Q = cgetg(l,t_COL);
845: if (l == 1) return cgetg(1,t_COL);
846: for (i = l-1; i>0; i--)
847: {
1.2 ! noro 848: q = negi(diviiround((GEN)x[i], gcoeff(y,i,i)));
1.1 noro 849: if (Q) (*Q)[i] = (long)q;
850: if (signe(q)) x = gadd(x, gmul(q, (GEN)y[i]));
851: }
852: return x;
853: }
854:
855: /* for internal use...Reduce matrix x modulo matrix y */
856: GEN
857: reducemodmatrix(GEN x, GEN y)
858: {
859: return reducemodHNF(x, hnfmod(y,detint(y)), NULL);
860: }
861:
862: /* x = y Q + R */
863: GEN
864: reducemodHNF(GEN x, GEN y, GEN *Q)
865: {
866: long lx = lg(x), i;
867: GEN R = cgetg(lx, t_MAT);
868: if (Q)
869: {
870: GEN q = cgetg(lx, t_MAT); *Q = q;
1.2 ! noro 871: for (i=1; i<lx; i++) R[i] = (long)colreducemodHNF((GEN)x[i],y,(GEN*)(q+i));
1.1 noro 872: }
873: else
874: for (i=1; i<lx; i++)
875: {
1.2 ! noro 876: gpmem_t av = avma;
! 877: R[i] = lpileupto(av, colreducemodHNF((GEN)x[i],y,NULL));
1.1 noro 878: }
879: return R;
880: }
881:
1.2 ! noro 882: /* For internal use. Reduce x modulo ideal (assumed non-zero, in HNF). We
! 883: * want a non-zero result */
! 884: GEN
! 885: nfreducemodideal_i(GEN x0,GEN ideal)
! 886: {
! 887: GEN x = colreducemodHNF(x0, ideal, NULL);
! 888: if (gcmp0(x)) return (GEN)ideal[1];
! 889: return x == x0? gcopy(x) : x;
! 890: }
1.1 noro 891: GEN
892: nfreducemodideal(GEN nf,GEN x0,GEN ideal)
893: {
1.2 ! noro 894: return nfreducemodideal_i(x0, idealhermite(nf,ideal));
1.1 noro 895: }
896:
897: /* multiply y by t = 1 mod^* f such that sign(x) = sign(y) at arch = idele[2].
898: * If x == NULL, make y >> 0 at sarch */
899: GEN
900: set_sign_mod_idele(GEN nf, GEN x, GEN y, GEN idele, GEN sarch)
901: {
902: GEN s,arch,gen;
903: long nba,i;
904: if (!sarch) return y;
905: gen = (GEN)sarch[2]; nba = lg(gen);
906: if (nba == 1) return y;
907:
908: arch = (GEN)idele[2];
909: s = zsigne(nf,y,arch);
910: if (x) s = gadd(s, zsigne(nf,x,arch));
911: s = lift_intern(gmul((GEN)sarch[3],s));
912: for (i=1; i<nba; i++)
913: if (signe(s[i])) y = element_mul(nf,y,(GEN)gen[i]);
914: return y;
915: }
916:
917: /* compute elt = x mod idele, with sign(elt) = sign(x) at arch */
918: GEN
919: nfreducemodidele(GEN nf,GEN x,GEN idele,GEN sarch)
920: {
921: GEN y;
922:
923: if (gcmp0(x)) return gcopy(x);
924: if (!sarch || typ(idele)!=t_VEC || lg(idele)!=3)
925: return nfreducemodideal(nf,x,idele);
926:
927: y = nfreducemodideal(nf,x,(GEN)idele[1]);
928: y = set_sign_mod_idele(nf, x, y, idele, sarch);
929: return (gexpo(y) > gexpo(x))? x: y;
930: }
931:
932: GEN
933: element_powmodideal(GEN nf,GEN x,GEN k,GEN ideal)
934: {
935: GEN y = gscalcol_i(gun,degpol(nf[1]));
936: for(;;)
937: {
938: if (mpodd(k)) y=element_mulmodideal(nf,x,y,ideal);
939: k=shifti(k,-1); if (!signe(k)) return y;
940: x = element_sqrmodideal(nf,x,ideal);
941: }
942: }
943:
944: GEN
945: element_powmodidele(GEN nf,GEN x,GEN k,GEN idele,GEN sarch)
946: {
947: GEN y = element_powmodideal(nf,x,k,idele);
948: if (mpodd(k))
949: { if (!gcmp1(k)) y = set_sign_mod_idele(nf, x, y, idele, sarch); }
950: else
951: { if (!gcmp0(k)) y = set_sign_mod_idele(nf, NULL, y, idele, sarch); }
952: return y;
953: }
954:
955: /* H relation matrix among row of generators g. Let URV = D its SNF, newU R
956: * newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of newD,
957: * newU and newUi such that 1/U = (newUi, ?).
958: * Rationale: let (G,0) = g Ui be the new generators then
959: * 0 = G U R --> G D = 0, g = G newU and G = g newUi */
960: GEN
961: smithrel(GEN H, GEN *newU, GEN *newUi)
962: {
1.2 ! noro 963: GEN U,Ui,D,p1;
1.1 noro 964: long l,c;
965:
1.2 ! noro 966: D = smithall(H, &U, NULL);
! 967: l = lg(D);
1.1 noro 968: for (c=1; c<l; c++)
969: {
970: p1 = gcoeff(D,c,c);
971: if (is_pm1(p1)) break;
972: }
973: if (newU) *newU = rowextract_i(U, 1, c-1);
974: if (newUi) {
975: Ui = ginv(U); setlg(Ui, c);
976: *newUi = reducemodHNF(Ui, H, NULL);
977: }
978: setlg(D, c); return mattodiagonal_i(D);
979: }
980:
981: /* given 2 integral ideals x, y in HNF s.t x | y | x^2, compute the quotient
982: (1+x)/(1+y) in the form [[cyc],[gen],ux^-1]. */
983: static GEN
984: zidealij(GEN x, GEN y)
985: {
986: GEN U,R,G,p1,cyc,z,xi;
987: long j,N;
988:
989: xi = ginv(x); R = gmul(xi, y); /* relations between the 1 + x_i (HNF) */
990: cyc = smithrel(R, &U, &G);
991: N = lg(cyc); G = gmul(x,G); settyp(G, t_VEC); /* new generators */
992: for (j=1; j<N; j++)
993: {
994: p1 = (GEN)G[j];
995: p1[1] = laddsi(1, (GEN)p1[1]); /* 1 + g_j */
996: }
997: z = cgetg(4,t_VEC);
998: z[1] = (long)cyc;
999: z[2] = (long)G;
1000: z[3] = lmul(U,xi); return z;
1001: }
1002:
1003: /* smallest integer n such that g0^n=x modulo p prime. Assume g0 has order q
1004: * (p-1 if q = NULL) */
1005: GEN
1006: Fp_shanks(GEN x,GEN g0,GEN p, GEN q)
1007: {
1.2 ! noro 1008: gpmem_t av=avma,av1,lim;
! 1009: long lbaby,i,k,c;
1.1 noro 1010: GEN p1,smalltable,giant,perm,v,g0inv;
1011:
1012: x = modii(x,p);
1013: if (is_pm1(x) || egalii(p,gdeux)) { avma = av; return gzero; }
1014: p1 = addsi(-1, p); if (!q) q = p1;
1015: if (egalii(p1,x)) { avma = av; return shifti(q,-1); }
1016: p1 = racine(q);
1017: if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in Fp_shanks");
1018: lbaby=itos(p1)+1; smalltable=cgetg(lbaby+1,t_VEC);
1019: g0inv = mpinvmod(g0, p); p1 = x;
1020:
1021: c = 3 * lgefint(p);
1022: for (i=1;;i++)
1023: {
1024: av1 = avma;
1025: if (is_pm1(p1)) { avma=av; return stoi(i-1); }
1026: smalltable[i]=(long)p1; if (i==lbaby) break;
1.2 ! noro 1027: (void)new_chunk(c); p1 = mulii(p1,g0inv); /* HACK */
1.1 noro 1028: avma = av1; p1 = resii(p1, p);
1029: }
1030: giant = resii(mulii(x, mpinvmod(p1,p)), p);
1031: p1=cgetg(lbaby+1,t_VEC);
1032: perm = gen_sort(smalltable, cmp_IND | cmp_C, cmpii);
1033: for (i=1; i<=lbaby; i++) p1[i]=smalltable[perm[i]];
1034: smalltable=p1; p1=giant;
1035:
1036: av1 = avma; lim=stack_lim(av1,2);
1037: for (k=1;;k++)
1038: {
1039: i=tablesearch(smalltable,p1,cmpii);
1040: if (i)
1041: {
1042: v=addis(mulss(lbaby-1,k),perm[i]);
1043: return gerepileuptoint(av,addsi(-1,v));
1044: }
1045: p1 = resii(mulii(p1,giant), p);
1046:
1047: if (low_stack(lim, stack_lim(av1,2)))
1048: {
1049: if(DEBUGMEM>1) err(warnmem,"Fp_shanks, k = %ld", k);
1050: p1 = gerepileuptoint(av1, p1);
1051: }
1052: }
1053: }
1054:
1055: /* Pohlig-Hellman */
1056: GEN
1057: Fp_PHlog(GEN a, GEN g, GEN p, GEN ord)
1058: {
1.2 ! noro 1059: gpmem_t av = avma;
1.1 noro 1060: GEN v,t0,a0,b,q,g_q,n_q,ginv0,qj,ginv;
1061: GEN fa, ex;
1062: long e,i,j,l;
1063:
1064: if (!ord) ord = subis(p,1);
1065: if (typ(ord) == t_MAT)
1066: {
1067: fa = ord;
1068: ord= factorback(fa,NULL);
1069: }
1070: else
1071: fa = decomp(ord);
1072: if (typ(g) == t_INTMOD) g = lift_intern(g);
1073: ex = (GEN)fa[2];
1074: fa = (GEN)fa[1];
1075: l = lg(fa);
1076: ginv = mpinvmod(g,p);
1077: v = cgetg(l, t_VEC);
1078: for (i=1; i<l; i++)
1079: {
1080: q = (GEN)fa[i];
1081: e = itos((GEN)ex[i]);
1082: if (DEBUGLEVEL>5)
1083: fprintferr("Pohlig-Hellman: DL mod %Z^%ld\n",q,e);
1084: qj = new_chunk(e+1); qj[0] = un;
1085: for (j=1; j<=e; j++) qj[j] = lmulii((GEN)qj[j-1], q);
1086: t0 = divii(ord, (GEN)qj[e]);
1087: a0 = powmodulo(a, t0, p);
1088: ginv0 = powmodulo(ginv, t0, p); /* order q^e */
1089: g_q = powmodulo(g, divii(ord,q), p); /* order q */
1090: n_q = gzero;
1091: for (j=0; j<e; j++)
1092: {
1093: b = modii(mulii(a0, powmodulo(ginv0, n_q, p)), p);
1094: b = powmodulo(b, (GEN)qj[e-1-j], p);
1095: b = Fp_shanks(b, g_q, p, q);
1096: n_q = addii(n_q, mulii(b, (GEN)qj[j]));
1097: }
1098: v[i] = lmodulcp(n_q, (GEN)qj[e]);
1099: }
1100: return gerepileuptoint(av, lift(chinese(v,NULL)));
1101: }
1102:
1.2 ! noro 1103: /* discrete log in Fq for a in Fp^*, g primitive root in Fq^* */
! 1104: static GEN
! 1105: ff_PHlog_Fp(GEN a, GEN g, GEN T, GEN p)
! 1106: {
! 1107: gpmem_t av = avma;
! 1108: GEN q,n_q,ord,ordp;
! 1109:
! 1110: if (gcmp1(a) || egalii(p, gdeux)) { avma = av; return gzero; }
! 1111:
! 1112: ordp = subis(p, 1);
! 1113: ord = T? subis(gpowgs(p,degpol(T)), 1): p;
! 1114: if (egalii(a, ordp)) /* -1 */
! 1115: return gerepileuptoint(av, shifti(ord,-1));
! 1116:
! 1117: if (!T) q = NULL;
! 1118: else
! 1119: { /* we want < g > = Fp^* */
! 1120: q = divii(ord,ordp);
! 1121: g = FpXQ_pow(g,q,T,p);
! 1122: if (typ(g) == t_POL) g = constant_term(g);
! 1123: }
! 1124: n_q = Fp_PHlog(a,g,p,NULL);
! 1125: if (q) n_q = mulii(q, n_q);
! 1126: return gerepileuptoint(av, n_q);
! 1127: }
! 1128:
1.1 noro 1129: /* smallest n >= 0 such that g0^n=x modulo pr, assume g0 reduced
1.2 ! noro 1130: * q = order of g0 is prime (and != p) */
1.1 noro 1131: static GEN
1.2 ! noro 1132: ffshanks(GEN x, GEN g0, GEN q, GEN T, GEN p)
1.1 noro 1133: {
1.2 ! noro 1134: gpmem_t av = avma, av1, lim;
! 1135: long lbaby,i,k;
! 1136: GEN p1,smalltable,giant,perm,v,g0inv;
! 1137:
! 1138: if (!degpol(x)) x = constant_term(x);
! 1139: if (typ(x) == t_INT)
! 1140: {
! 1141: if (!gcmp1(modii(p,q))) return gzero;
! 1142: /* g0 in Fp^*, order q | (p-1) */
! 1143: if (typ(g0) == t_POL) g0 = constant_term(g0);
! 1144: return Fp_PHlog(x,g0,p,q);
! 1145: }
! 1146:
! 1147: p1 = racine(q);
! 1148: if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in ffshanks");
! 1149: lbaby = itos(p1)+1; smalltable = cgetg(lbaby+1,t_VEC);
! 1150: g0inv = FpXQ_inv(g0,T,p);
1.1 noro 1151: p1 = x;
1152:
1153: for (i=1;;i++)
1154: {
1.2 ! noro 1155: if (gcmp1(p1)) { avma = av; return stoi(i-1); }
1.1 noro 1156:
1157: smalltable[i]=(long)p1; if (i==lbaby) break;
1.2 ! noro 1158: p1 = FpXQ_mul(p1,g0inv, T,p);
1.1 noro 1159: }
1.2 ! noro 1160: giant = FpXQ_div(x,p1,T,p);
! 1161: perm = gen_sort(smalltable, cmp_IND | cmp_C, cmp_pol);
! 1162: smalltable = vecextract_p(smalltable, perm);
! 1163: p1 = giant;
1.1 noro 1164:
1165: av1 = avma; lim=stack_lim(av1,2);
1166: for (k=1;;k++)
1167: {
1.2 ! noro 1168: i = tablesearch(smalltable,p1,cmp_pol);
1.1 noro 1169: if (i)
1170: {
1.2 ! noro 1171: v = addis(mulss(lbaby-1,k), perm[i]);
! 1172: return gerepileuptoint(av, addsi(-1,v));
1.1 noro 1173: }
1.2 ! noro 1174: p1 = FpXQ_mul(p1, giant, T,p);
1.1 noro 1175:
1176: if (low_stack(lim, stack_lim(av1,2)))
1177: {
1.2 ! noro 1178: if(DEBUGMEM>1) err(warnmem,"ffshanks");
1.1 noro 1179: p1 = gerepileupto(av1, p1);
1180: }
1181: }
1182: }
1183:
1.2 ! noro 1184: /* same in Fp[X]/T */
1.1 noro 1185: GEN
1.2 ! noro 1186: ff_PHlog(GEN a, GEN g, GEN T, GEN p)
1.1 noro 1187: {
1.2 ! noro 1188: gpmem_t av = avma;
1.1 noro 1189: GEN v,t0,a0,b,q,g_q,n_q,ginv0,qj,ginv,ord,fa,ex;
1190: long e,i,j,l;
1191:
1.2 ! noro 1192: if (typ(a) == t_INT)
! 1193: return gerepileuptoint(av, ff_PHlog_Fp(a,g,T,p));
! 1194: /* f > 1 ==> T != NULL */
! 1195: ord = subis(gpowgs(p, degpol(T)), 1);
! 1196: fa = factor(ord);
! 1197: ex = (GEN)fa[2];
1.1 noro 1198: fa = (GEN)fa[1];
1199: l = lg(fa);
1.2 ! noro 1200: ginv = FpXQ_inv(g,T,p);
1.1 noro 1201: v = cgetg(l, t_VEC);
1202: for (i=1; i<l; i++)
1203: {
1204: q = (GEN)fa[i];
1205: e = itos((GEN)ex[i]);
1.2 ! noro 1206: if (DEBUGLEVEL>5) fprintferr("nf_Pohlig-Hellman: DL mod %Z^%ld\n",q,e);
1.1 noro 1207: qj = new_chunk(e+1); qj[0] = un;
1208: for (j=1; j<=e; j++) qj[j] = lmulii((GEN)qj[j-1], q);
1209: t0 = divii(ord, (GEN)qj[e]);
1.2 ! noro 1210: a0 = FpXQ_pow(a, t0, T,p);
! 1211: ginv0 = FpXQ_pow(ginv, t0, T,p); /* order q^e */
! 1212: g_q = FpXQ_pow(g, divii(ord,q), T,p); /* order q */
! 1213:
1.1 noro 1214: n_q = gzero;
1215: for (j=0; j<e; j++)
1216: {
1.2 ! noro 1217: b = FpXQ_mul(a0, FpXQ_pow(ginv0, n_q, T,p), T,p);
! 1218: b = FpXQ_pow(b, (GEN)qj[e-1-j], T,p);
! 1219: b = ffshanks(b, g_q, q, T,p);
1.1 noro 1220: n_q = addii(n_q, mulii(b, (GEN)qj[j]));
1221: }
1222: v[i] = lmodulcp(n_q, (GEN)qj[e]);
1223: }
1224: return gerepileuptoint(av, lift(chinese(v,NULL)));
1225: }
1226:
1.2 ! noro 1227: /* same in nf.zk / pr */
! 1228: GEN
! 1229: nf_PHlog(GEN nf, GEN a, GEN g, GEN pr)
! 1230: {
! 1231: gpmem_t av = avma;
! 1232: GEN T,p, modpr = nf_to_ff_init(nf, &pr, &T, &p);
! 1233: GEN A = nf_to_ff(nf,a,modpr);
! 1234: GEN G = nf_to_ff(nf,g,modpr);
! 1235: return gerepileuptoint(av, ff_PHlog(A,G,T,p));
! 1236: }
! 1237:
1.1 noro 1238: GEN
1239: znlog(GEN x, GEN g0)
1240: {
1.2 ! noro 1241: gpmem_t av=avma;
1.1 noro 1242: GEN p = (GEN)g0[1];
1243: if (typ(g0) != t_INTMOD)
1244: err(talker,"not an element of (Z/pZ)* in znlog");
1245: switch(typ(x))
1246: {
1247: case t_INT: break;
1248: default: x = gmul(x, gmodulsg(1,p));
1249: if (typ(x) != t_INTMOD)
1250: err(talker,"not an element of (Z/pZ)* in znlog");
1251: case t_INTMOD: x = (GEN)x[2]; break;
1252: }
1253: return gerepileuptoint(av, Fp_PHlog(x,(GEN)g0[2],p,NULL));
1254: }
1255:
1256: GEN
1257: dethnf(GEN mat)
1258: {
1259: long i,l = lg(mat);
1.2 ! noro 1260: gpmem_t av;
1.1 noro 1261: GEN s;
1262:
1263: if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
1264: av = avma; s = gcoeff(mat,1,1);
1265: for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
1266: return av==avma? gcopy(s): gerepileupto(av,s);
1267: }
1268:
1269: GEN
1270: dethnf_i(GEN mat)
1271: {
1.2 ! noro 1272: gpmem_t av;
! 1273: long i,l = lg(mat);
1.1 noro 1274: GEN s;
1275:
1276: if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
1277: av = avma; s = gcoeff(mat,1,1);
1278: for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
1279: return gerepileuptoint(av,s);
1280: }
1281:
1282: /* as above with cyc = diagonal(mat) */
1283: GEN
1284: detcyc(GEN cyc)
1285: {
1.2 ! noro 1286: gpmem_t av;
! 1287: long i,l = lg(cyc);
1.1 noro 1288: GEN s;
1289:
1290: if (l<3) return l<2? gun: icopy((GEN)cyc[1]);
1291: av = avma; s = (GEN)cyc[1];
1292: for (i=2; i<l; i++) s = mulii(s,(GEN)cyc[i]);
1293: return gerepileuptoint(av,s);
1294: }
1295:
1296: /* (U,V) = 1. Return y = x mod U, = 1 mod V (uv: u + v = 1, u in U, v in V) */
1.2 ! noro 1297: GEN
1.1 noro 1298: makeprimetoideal(GEN nf,GEN UV,GEN uv,GEN x)
1299: {
1300: GEN y = gadd((GEN)uv[1], element_mul(nf,x,(GEN)uv[2]));
1.2 ! noro 1301: return nfreducemodideal_i(y,UV);
1.1 noro 1302: }
1303:
1304: static GEN
1305: makeprimetoidealvec(GEN nf,GEN UV,GEN uv,GEN gen)
1306: {
1307: long i, lx = lg(gen);
1308: GEN y = cgetg(lx,t_VEC);
1309: for (i=1; i<lx; i++)
1310: y[i] = (long)makeprimetoideal(nf,UV,uv,(GEN)gen[i]);
1311: return y;
1312: }
1313:
1.2 ! noro 1314: GEN
! 1315: FpXQ_gener(GEN T, GEN p)
! 1316: {
! 1317: long i,j, k, vT = varn(T), f = degpol(T);
! 1318: GEN g, list, pf_1 = subis(gpowgs(p, f), 1);
! 1319: gpmem_t av0 = avma, av;
! 1320:
! 1321: list = (GEN)factor(pf_1)[1];
! 1322: k = lg(list)-1;
! 1323:
! 1324: for (i=1; i<=k; i++) list[i] = (long)diviiexact(pf_1, (GEN)list[i]);
! 1325: for (av = avma;; avma = av)
! 1326: {
! 1327: g = FpX_rand(f, vT, p);
! 1328: if (degpol(g) < 1) continue;
! 1329: for (j=1; j<=k; j++)
! 1330: if (gcmp1(FpXQ_pow(g, (GEN)list[j], T, p))) break;
! 1331: if (j > k) return gerepilecopy(av0, g);
! 1332: }
! 1333: }
! 1334:
1.1 noro 1335: /* Given an ideal pr^ep, and an integral ideal x (in HNF form) compute a list
1336: * of vectors, each with 5 components as follows :
1337: * [[clh],[gen1],[gen2],[signat2],U.X^-1]. Each component corresponds to
1338: * d_i,g_i,g'_i,s_i. Generators g_i are not necessarily prime to x, the
1339: * generators g'_i are. signat2 is the (horizontal) vector made of the
1340: * signatures (column vectors) of the g'_i. If x = NULL, the original ideal
1341: * was a prime power
1342: */
1343: static GEN
1344: zprimestar(GEN nf,GEN pr,GEN ep,GEN x,GEN arch)
1345: {
1.2 ! noro 1346: gpmem_t av = avma, av1, tetpil;
! 1347: long N, f, i, e, a, b;
1.1 noro 1348: GEN prh,p,pf_1,list,v,p1,p3,p4,prk,uv,g0,newgen,pra,prb;
1349: GEN *gptr[2];
1350:
1351: if(DEBUGLEVEL>3) { fprintferr("treating pr = %Z ^ %Z\n",pr,ep); flusherr(); }
1.2 ! noro 1352: prh = prime_to_ideal(nf,pr); N = degpol(nf[1]);
! 1353: f = itos((GEN)pr[4]);
! 1354: p = (GEN)pr[1];
1.1 noro 1355:
1356: pf_1 = addis(gpowgs(p,f), -1);
1.2 ! noro 1357: if (f==1)
! 1358: {
! 1359: v = zerocol(N);
! 1360: v[1] = gener(p)[2];
! 1361: }
1.1 noro 1362: else
1363: {
1.2 ! noro 1364: GEN T, modpr = zk_to_ff_init(nf, &pr, &T, &p);
! 1365: v = ff_to_nf(FpXQ_gener(T,p), modpr);
! 1366: v = algtobasis(nf, v);
1.1 noro 1367: }
1368: /* v generates (Z_K / pr)^* */
1.2 ! noro 1369: if(DEBUGLEVEL>3) fprintferr("v computed\n");
1.1 noro 1370: e = itos(ep); prk=(e==1)? pr: idealpow(nf,pr,ep);
1.2 ! noro 1371: if(DEBUGLEVEL>3) fprintferr("prk computed\n");
1.1 noro 1372: g0 = v;
1373: uv = NULL; /* gcc -Wall */
1374: if (x)
1375: {
1.2 ! noro 1376: uv = idealaddtoone(nf,prk, idealdivpowprime(nf,x,pr,ep));
1.1 noro 1377: g0 = makeprimetoideal(nf,x,uv,v);
1.2 ! noro 1378: if(DEBUGLEVEL>3) fprintferr("g0 computed\n");
1.1 noro 1379: }
1380:
1381: p1 = cgetg(6,t_VEC); list = _vec(p1);
1382: p1[1] = (long)_vec(pf_1);
1383: p1[2] = (long)_vec(v);
1384: p1[3] = (long)_vec(g0);
1385: p1[4] = (long)_vec(zsigne(nf,g0,arch));
1386: p1[5] = un;
1387: if (e==1) return gerepilecopy(av, list);
1388:
1389: a=1; b=2; av1=avma;
1390: pra = prh; prb = (e==2)? prk: idealpow(nf,pr,gdeux);
1391: for(;;)
1392: {
1.2 ! noro 1393: if(DEBUGLEVEL>3) fprintferr(" treating a = %ld, b = %ld\n",a,b);
1.1 noro 1394: p1 = zidealij(pra,prb);
1395: newgen = dummycopy((GEN)p1[2]);
1396: p3 = cgetg(lg(newgen),t_VEC);
1.2 ! noro 1397: if(DEBUGLEVEL>3) fprintferr("zidealij done\n");
1.1 noro 1398: for (i=1; i<lg(newgen); i++)
1399: {
1400: if (x) newgen[i]=(long)makeprimetoideal(nf,x,uv,(GEN)newgen[i]);
1401: p3[i]=(long)zsigne(nf,(GEN)newgen[i],arch);
1402: }
1403: p4 = cgetg(6,t_VEC);
1404: p4[1] = p1[1];
1405: p4[2] = p1[2];
1406: p4[3] = (long)newgen;
1407: p4[4] = (long)p3;
1408: p4[5] = p1[3]; p4 = _vec(p4);
1409:
1410: a=b; b=min(e,b<<1); tetpil = avma;
1411: list = concat(list, p4);
1412: if (a==b) return gerepile(av,tetpil,list);
1413:
1414: pra = gcopy(prb);
1415: gptr[0]=&pra; gptr[1]=&list;
1416: gerepilemanysp(av1,tetpil,gptr,2);
1.2 ! noro 1417: prb = (b==(a<<1))? idealpows(nf,pr,b): prk;
1.1 noro 1418: }
1419: }
1420:
1.2 ! noro 1421: /* x integral ideal, compute elements in 1+x whose sign matrix is invertible */
1.1 noro 1422: GEN
1423: zarchstar(GEN nf,GEN x,GEN arch,long nba)
1424: {
1.2 ! noro 1425: long limr, N, i, j, k, r, rr, zk, lgmat;
! 1426: gpmem_t av;
1.1 noro 1427: GEN p1,y,bas,genarch,mat,lambda,nfun,vun;
1428:
1.2 ! noro 1429: if (nba < 0)
! 1430: {
! 1431: nba = 0;
! 1432: for (i=1; i<lg(arch); i++)
! 1433: if (signe(arch[i])) nba++;
! 1434: }
1.1 noro 1435: y = cgetg(4,t_VEC);
1436: if (!nba)
1437: {
1438: y[1]=lgetg(1,t_VEC);
1439: y[2]=lgetg(1,t_VEC);
1440: y[3]=lgetg(1,t_MAT); return y;
1441: }
1442: N = degpol(nf[1]);
1443: p1 = cgetg(nba+1,t_VEC); for (i=1; i<=nba; i++) p1[i] = deux;
1444: y[1] = (long)p1; av = avma;
1445: if (N == 1)
1446: {
1447: p1 = gerepileuptoint(av, subsi(1, shifti(gcoeff(x,1,1),1)));
1448: y[2] = (long)_vec(_col(p1));
1449: y[3] = (long)gscalmat(gun,1); return y;
1450: }
1.2 ! noro 1451: zk = gcmp1(gcoeff(x,1,1));
1.1 noro 1452: genarch = cgetg(nba+1,t_VEC);
1453: mat = cgetg(nba+1,t_MAT); setlg(mat,2);
1.2 ! noro 1454: for (i=1; i<=nba; i++) mat[i] = lgetg(nba+1, t_VECSMALL);
1.1 noro 1455: lambda = cgetg(N+1, t_VECSMALL);
1456:
1457: bas = dummycopy(gmael(nf,5,1)); r = lg(arch);
1458: for (k=1,i=1; i<r; i++)
1459: if (signe(arch[i]))
1460: {
1461: if (k < i)
1462: for (j=1; j<=N; j++) coeff(bas,k,j) = coeff(bas,i,j);
1463: k++;
1464: }
1465: r = nba+1; for (i=1; i<=N; i++) setlg(bas[i], r);
1466: if (!zk)
1467: {
1.2 ! noro 1468: x = lllint_ip(x,4);
1.1 noro 1469: bas = gmul(bas, x);
1470: }
1471:
1472: vun = cgetg(nba+1, t_COL);
1473: for (i=1; i<=nba; i++) vun[i] = un;
1474: nfun = gscalcol(gun, N); lgmat = 1;
1475: for (r=1, rr=3; ; r<<=1, rr=(r<<1)+1)
1476: {
1477: if (DEBUGLEVEL>5) fprintferr("zarchstar: r = %ld\n",r);
1478: p1 = gpowgs(stoi(rr),N);
1.2 ! noro 1479: limr = (cmpis(p1,BIGINT) >= 0)? BIGINT: p1[2];
1.1 noro 1480: limr = (limr-1)>>1;
1481: k = zk? rr: 1; /* to avoid 0 */
1482: for ( ; k<=limr; k++)
1483: {
1.2 ! noro 1484: gpmem_t av1=avma;
! 1485: long kk = k;
1.1 noro 1486: GEN alpha = vun;
1487: for (i=1; i<=N; i++)
1488: {
1489: lambda[i] = (kk+r)%rr - r; kk/=rr;
1490: if (lambda[i])
1491: alpha = gadd(alpha, gmulsg(lambda[i],(GEN)bas[i]));
1492: }
1493: p1 = (GEN)mat[lgmat];
1494: for (i=1; i<=nba; i++)
1.2 ! noro 1495: p1[i] = (gsigne((GEN)alpha[i]) < 0)? 1: 0;
1.1 noro 1496:
1.2 ! noro 1497: if (u_FpM_deplin(mat, 2)) avma = av1;
1.1 noro 1498: else
1499: { /* new vector indep. of previous ones */
1500: avma = av1; alpha = nfun;
1501: for (i=1; i<=N; i++)
1502: if (lambda[i])
1503: alpha = gadd(alpha, gmulsg(lambda[i],(GEN)x[i]));
1504: genarch[lgmat++] = (long)alpha;
1505: if (lgmat > nba)
1506: {
1507: GEN *gptr[2];
1.2 ! noro 1508: mat = small_to_mat( u_FpM_inv(mat, 2) );
! 1509: gptr[0]=&genarch; gptr[1]=&mat;
1.1 noro 1510: gerepilemany(av,gptr,2);
1511: y[2]=(long)genarch;
1512: y[3]=(long)mat; return y;
1513: }
1514: setlg(mat,lgmat+1);
1515: }
1516: }
1517: }
1518: }
1519:
1520: GEN
1521: zinternallog_pk(GEN nf, GEN a0, GEN y, GEN pr, GEN prk, GEN list, GEN *psigne)
1522: {
1523: GEN a = a0, L,e,p1,cyc,gen;
1524: long i,j, llist = lg(list)-1;
1525:
1526: e = gzero; /* gcc -Wall */
1527: for (j=1; j<=llist; j++)
1528: {
1529: L = (GEN)list[j]; cyc=(GEN)L[1]; gen=(GEN)L[2];
1530: if (j==1)
1.2 ! noro 1531: p1 = nf_PHlog(nf,a,(GEN)gen[1],pr);
1.1 noro 1532: else
1533: {
1534: p1 = (GEN)a[1]; a[1] = laddsi(-1,(GEN)a[1]);
1535: e = gmul((GEN)L[5],a); a[1] = (long)p1;
1536: /* here lg(e) == lg(cyc) */
1537: p1 = (GEN)e[1];
1538: }
1539: for (i=1; i<lg(cyc); i++, p1 = (GEN)e[i])
1540: {
1541: p1 = modii(negi(p1), (GEN)cyc[i]);
1542: *++y = lnegi(p1);
1543: if (!signe(p1)) continue;
1544:
1545: if (psigne && mpodd(p1)) *psigne = gadd(*psigne,gmael(L,4,i));
1546: if (j == llist) continue;
1547:
1548: if (DEBUGLEVEL>5) fprintferr("do element_powmodideal\n");
1549: p1 = element_powmodideal(nf,(GEN)gen[i],p1,prk);
1550: a = element_mulmodideal(nf,a,p1,prk);
1551: }
1552: }
1553: return y;
1554: }
1555:
1556: /* Retourne la decomposition de a sur les nbgen generateurs successifs
1557: * contenus dans list_set et si index !=0 on ne fait ce calcul que pour
1558: * l'ideal premier correspondant a cet index en mettant a 0 les autres
1559: * composantes
1560: */
1561: static GEN
1562: zinternallog(GEN nf,GEN a,GEN list_set,long nbgen,GEN arch,GEN fa,long index)
1563: {
1564: GEN prlist,ep,y0,y,ainit,list,pr,prk,cyc,psigne,p1,p2;
1.2 ! noro 1565: gpmem_t av;
! 1566: long nbp,i,j,k;
1.1 noro 1567:
1568: y0 = y = cgetg(nbgen+1,t_COL); av=avma;
1569: prlist=(GEN)fa[1]; ep=(GEN)fa[2]; nbp=lg(ep)-1;
1570: i=typ(a); if (is_extscalar_t(i)) a = algtobasis(nf,a);
1571: if (DEBUGLEVEL>3)
1572: {
1573: fprintferr("entering zinternallog, "); flusherr();
1574: if (DEBUGLEVEL>5) fprintferr("with a = %Z\n",a);
1575: }
1576: ainit = a; psigne = zsigne(nf,ainit,arch);
1577: p2 = NULL; /* gcc -Wall */
1578: for (k=1; k<=nbp; k++)
1579: {
1580: list = (GEN)list_set[k];
1581: if (index && index!=k)
1582: {
1583: for (j=1; j<lg(list); j++)
1584: {
1585: cyc = gmael(list,j,1);
1586: for (i=1; i<lg(cyc); i++) *++y = zero;
1587: }
1588: continue;
1589: }
1590: pr = (GEN)prlist[k]; prk = idealpow(nf,pr,(GEN)ep[k]);
1591: y = zinternallog_pk(nf, ainit, y, pr, prk, list, &psigne);
1592: }
1593: p1 = lift_intern(gmul(gmael(list_set,k,3), psigne));
1594: avma=av; for (i=1; i<lg(p1); i++) *++y = p1[i];
1595: if (DEBUGLEVEL>3) { fprintferr("leaving\n"); flusherr(); }
1596: for (i=1; i<=nbgen; i++) y0[i] = licopy((GEN)y0[i]);
1597: return y0;
1598: }
1599:
1600: static GEN
1601: compute_gen(GEN nf, GEN u1, GEN gen, GEN bid)
1602: {
1603: long i, c = lg(u1);
1604: GEN basecl = cgetg(c,t_VEC);
1605:
1606: for (i=1; i<c; i++)
1607: basecl[i] = (long)famat_to_nf_modidele(nf, gen, (GEN)u1[i], bid);
1608: return basecl;
1609: }
1610:
1611: /* Compute [[ideal,arch], [h,[cyc],[gen]], idealfact, [liste], U]
1612: gen not computed unless add_gen != 0 */
1613: GEN
1614: zidealstarinitall(GEN nf, GEN ideal,long add_gen)
1615: {
1.2 ! noro 1616: gpmem_t av = avma;
! 1617: long i,j,k,nba,nbp,R1,nbgen,cp;
1.1 noro 1618: GEN p1,y,h,cyc,U,u1,fa,fa2,ep,x,arch,list,sarch,gen;
1619:
1620: nf = checknf(nf); R1=itos(gmael(nf,2,1));
1621: if (typ(ideal)==t_VEC && lg(ideal)==3)
1622: {
1623: arch=(GEN)ideal[2]; ideal = (GEN)ideal[1];
1624: i = typ(arch);
1625: if (!is_vec_t(i) || lg(arch) != R1+1)
1626: err(talker,"incorrect archimedean component in zidealstarinit");
1627: for (nba=0,i=1; i<=R1; i++)
1628: if (signe(arch[i])) nba++;
1629: }
1630: else
1631: {
1632: arch=cgetg(R1+1,t_VEC);
1633: for (nba=0,i=1; i<=R1; i++) arch[i]=zero;
1634: }
1635: x = idealhermite(nf,ideal);
1636: if (!gcmp1(denom(x)))
1637: err(talker,"zidealstarinit needs an integral ideal: %Z",x);
1638: p1=cgetg(3,t_VEC); ideal=p1;
1639: p1[1]=(long)x;
1640: p1[2]=(long)arch;
1641:
1642: fa=idealfactor(nf,x); list=(GEN)fa[1]; ep=(GEN)fa[2];
1643: nbp=lg(list)-1; fa2=cgetg(nbp+2,t_VEC);
1644:
1645: gen = cgetg(1,t_VEC);
1646: p1 = (nbp==1)? (GEN)NULL: x;
1647: for (i=1; i<=nbp; i++)
1648: {
1649: GEN L = zprimestar(nf,(GEN)list[i],(GEN)ep[i],p1,arch);
1650: fa2[i] = (long)L;
1651: for (j=1; j<lg(L); j++) gen = concatsp(gen,gmael(L,j,3));
1652: }
1653: sarch = zarchstar(nf,x,arch,nba);
1654: fa2[i]=(long)sarch;
1655: gen = concatsp(gen,(GEN)sarch[2]);
1656:
1657: nbgen = lg(gen)-1;
1658: h=cgetg(nbgen+1,t_MAT); cp=0;
1659: for (i=1; i<=nbp; i++)
1660: {
1661: list = (GEN)fa2[i];
1662: for (j=1; j<lg(list); j++)
1663: {
1664: GEN a, L = (GEN)list[j], e = (GEN)L[1], g = (GEN)L[3];
1665: for (k=1; k<lg(g); k++)
1666: {
1667: a = element_powmodidele(nf,(GEN)g[k],(GEN)e[k],ideal,sarch);
1668: h[++cp] = lneg(zinternallog(nf,a,fa2,nbgen,arch,fa,i));
1669: coeff(h,cp,cp) = e[k];
1670: }
1671: }
1672: }
1673: for (j=1; j<=nba; j++) { h[++cp]=(long)zerocol(nbgen); coeff(h,cp,cp)=deux; }
1674: if (cp!=nbgen) err(talker,"bug in zidealstarinit");
1675: h = hnfall_i(h,NULL,0);
1676: cyc = smithrel(h, &U, add_gen? &u1: NULL);
1677: p1 = cgetg(add_gen? 4: 3, t_VEC);
1678: p1[1] = (long)detcyc(cyc);
1679: p1[2] = (long)cyc;
1680:
1681: y = cgetg(6,t_VEC);
1682: y[1] = (long)ideal;
1683: y[2] = (long)p1;
1684: y[3] = (long)fa;
1685: y[4] = (long)fa2;
1686: y[5] = (long)U;
1687: if (add_gen) p1[3] = (long)compute_gen(nf,u1,gen, y);
1688: return gerepilecopy(av, y);
1689: }
1690:
1691: GEN
1692: zidealstarinitgen(GEN nf, GEN ideal)
1693: {
1694: return zidealstarinitall(nf,ideal,1);
1695: }
1696:
1697: GEN
1698: zidealstarinit(GEN nf, GEN ideal)
1699: {
1700: return zidealstarinitall(nf,ideal,0);
1701: }
1702:
1703: GEN
1704: zidealstar(GEN nf, GEN ideal)
1705: {
1.2 ! noro 1706: gpmem_t av = avma;
1.1 noro 1707: GEN y = zidealstarinitall(nf,ideal,1);
1708: return gerepilecopy(av,(GEN)y[2]);
1709: }
1710:
1711: GEN
1712: idealstar0(GEN nf, GEN ideal,long flag)
1713: {
1714: switch(flag)
1715: {
1716: case 0: return zidealstar(nf,ideal);
1717: case 1: return zidealstarinit(nf,ideal);
1718: case 2: return zidealstarinitgen(nf,ideal);
1719: default: err(flagerr,"idealstar");
1720: }
1721: return NULL; /* not reached */
1722: }
1723:
1724: void
1725: check_nfelt(GEN x, GEN *den)
1726: {
1727: long l = lg(x), i;
1728: GEN t, d = NULL;
1729: if (typ(x) != t_COL) err(talker,"%Z not a nfelt", x);
1730: for (i=1; i<l; i++)
1731: {
1732: t = (GEN)x[i];
1733: switch (typ(t))
1734: {
1735: case t_INT: case t_INTMOD: break;
1736: case t_FRACN: case t_FRAC:
1737: if (!d) d = (GEN)t[2]; else d = mpppcm(d, (GEN)t[2]);
1738: break;
1739: default: err(talker,"%Z not a nfelt", x);
1740: }
1741: }
1742: *den = d;
1743: }
1744:
1745: /* Given x (not necessarily integral), and bid as output by zidealstarinit,
1746: * compute the vector of components on the generators bid[2].
1747: * Assume (x,bid) = 1 */
1748: GEN
1749: zideallog(GEN nf, GEN x, GEN bid)
1750: {
1.2 ! noro 1751: gpmem_t av;
! 1752: long l,i,N,c;
1.1 noro 1753: GEN fa,fa2,ideal,arch,den,p1,cyc,y;
1754:
1755: nf=checknf(nf); checkbid(bid);
1756: cyc=gmael(bid,2,2); c=lg(cyc);
1757: y=cgetg(c,t_COL); av=avma;
1758: N = degpol(nf[1]); ideal = (GEN) bid[1];
1759: if (typ(ideal)==t_VEC && lg(ideal)==3)
1760: arch = (GEN)ideal[2];
1761: else
1762: arch = NULL;
1763: switch(typ(x))
1764: {
1765: case t_INT: case t_FRAC: case t_FRACN:
1766: x = gscalcol_i(x,N); break;
1767: case t_POLMOD: case t_POL:
1768: x = algtobasis(nf,x); break;
1769: case t_MAT:
1.2 ! noro 1770: if (lg(x) == 1) return zerocol(c-1);
1.1 noro 1771: return famat_ideallog(nf,(GEN)x[1],(GEN)x[2],bid);
1772: case t_COL: break;
1773: default: err(talker,"not an element in zideallog");
1774: }
1775: if (lg(x) != N+1) err(talker,"not an element in zideallog");
1776: check_nfelt(x, &den);
1777: if (den)
1778: {
1779: GEN g = cgetg(3, t_COL);
1780: GEN e = cgetg(3, t_COL);
1.2 ! noro 1781: g[1] = (long)Q_muli_to_int(x,den); e[1] = un;
! 1782: g[2] = (long)den; e[2] = lstoi(-1);
1.1 noro 1783: p1 = famat_ideallog(nf, g, e, bid);
1784: }
1785: else
1786: {
1787: l=lg(bid[5])-1; fa=(GEN)bid[3]; fa2=(GEN)bid[4];
1788: p1 = zinternallog(nf,x,fa2,l,arch,fa,0);
1789: p1 = gmul((GEN)bid[5],p1); /* apply smith */
1790: }
1791: if (lg(p1)!=c) err(bugparier,"zideallog");
1792: for (i=1; i<c; i++)
1793: y[i] = lmodii((GEN)p1[i],(GEN)cyc[i]);
1794: avma=av; /* following line does a gerepile ! */
1795: for (i=1; i<c; i++)
1796: y[i] = (long)icopy((GEN)y[i]);
1797: return y;
1798: }
1799:
1800: /* bid1, bid2: output from 'zidealstarinit' for coprime modules m1 and m2
1801: * (without arch. part).
1802: * Output: zidealstarinit [[ideal,arch],[h,[cyc],[gen]],idealfact,[liste],U]
1803: * associated to m1 m2 */
1804: GEN
1805: zidealstarinitjoin(GEN nf, GEN bid1, GEN bid2, long add_gen)
1806: {
1.2 ! noro 1807: gpmem_t av=avma;
! 1808: long i,nbp,nbgen,lx1,lx2,l1,l2,lx;
1.1 noro 1809: GEN module1,module2,struct1,struct2,fact1,fact2,liste1,liste2;
1810: GEN module,liste,fact,U,U1,U2,cyc,cyc1,cyc2,uv;
1811: GEN p1,p2,y,u1,x,fa1,fa2, gen = add_gen? gun: NULL;
1812:
1813: nf = checknf(nf); checkbid(bid1); checkbid(bid2);
1814: module1 = (GEN)bid1[1]; struct1 = (GEN)bid1[2]; fact1 = (GEN)bid1[3];
1815: module2 = (GEN)bid2[1]; struct2 = (GEN)bid2[2]; fact2 = (GEN)bid2[3];
1816: x = idealmul(nf,(GEN)module1[1],(GEN)module2[1]);
1817: module = cgetg(3,t_VEC);
1818: module[1] = (long)x;
1819: module[2] = module1[2];
1820: if (!gcmp0((GEN)module1[2]) || !gcmp0((GEN)module2[2]))
1821: err(talker,"non-0 Archimedian components in zidealstarinitjoin");
1822:
1823: fa1 = (GEN)fact1[1]; lx1 = lg(fa1);
1824: fa2 = (GEN)fact2[1]; lx2 = lg(fa2);
1825: fact = cgetg(3,t_MAT);
1826: fact[1] = (long)concatsp(fa1,fa2);
1827: fact[2] = (long)concatsp((GEN)fact1[2],(GEN)fact2[2]);
1828: for (i=1; i<lx1; i++)
1829: if (isinvector(fa2,(GEN)fa1[i],lx2-1))
1830: err(talker,"noncoprime ideals in zidealstarinitjoin");
1831:
1832: liste1 = (GEN)bid1[4]; lx1 = lg(liste1);
1833: liste2 = (GEN)bid2[4]; lx2 = lg(liste2);
1834: /* concat (liste1 - last elt [zarchstar]) + liste2 */
1835: lx = lx1+lx2-2; liste = cgetg(lx,t_VEC);
1836: for (i=1; i<lx1-1; i++) liste[i] = liste1[i];
1837: for ( ; i<lx; i++) liste[i] = liste2[i-lx1+2];
1838:
1839: U1 = (GEN)bid1[5]; cyc1 = (GEN)struct1[2]; l1 = lg(cyc1);
1840: U2 = (GEN)bid2[5]; cyc2 = (GEN)struct2[2]; l2 = lg(cyc2);
1841: nbgen = l1+l2-2;
1842: cyc = diagonal(concatsp(cyc1,cyc2));
1843: cyc = smithrel(cyc, &U, gen? &u1: NULL);
1844:
1845: if (nbgen)
1846: U = concatsp(
1847: gmul(vecextract_i(U, 1, l1-1), U1) ,
1848: gmul(vecextract_i(U, l1, nbgen), U2)
1849: );
1850:
1851: if (gen)
1852: {
1853: if (lg(struct1)<=3 || lg(struct2)<=3)
1854: err(talker,"please apply idealstar(,,2) and not idealstar(,,1)");
1855:
1856: uv = idealaddtoone(nf,(GEN)module1[1],(GEN)module2[1]);
1857: p1 = makeprimetoidealvec(nf,x,uv,(GEN)struct1[3]);
1858: p2 = (GEN)uv[1]; uv[1] = uv[2]; uv[2] = (long)p2;
1859: p2 = makeprimetoidealvec(nf,x,uv,(GEN)struct2[3]);
1860: gen = concatsp(p1,p2);
1861: nbp = lg((GEN)fact[1])-1;
1862: }
1863: p1 = cgetg(gen? 4: 3,t_VEC);
1864: p1[1] = (long)detcyc(cyc);
1865: p1[2] = (long)cyc;
1866:
1867: y = cgetg(6,t_VEC);
1868: y[1] = (long)module;
1869: y[2] = (long)p1;
1870: y[3] = (long)fact;
1871: y[4] = (long)liste;
1872: y[5] = (long)U;
1873: if (gen) p1[3] = (long)compute_gen(nf,u1,gen, y);
1874: return gerepilecopy(av,y);
1875: }
1876:
1877: /* bid1: output from 'zidealstarinit' for module m1 (without arch. part)
1878: * arch: archimedean part.
1879: * Output: zidealstarinit [[ideal,arch],[h,[cyc],[gen]],idealfact,[liste],U]
1880: * associated to m1.arch */
1881: GEN
1882: zidealstarinitjoinarch(GEN nf, GEN bid1, GEN arch, long nba, long add_gen)
1883: {
1.2 ! noro 1884: gpmem_t av=avma;
! 1885: long i,nbp,lx1;
1.1 noro 1886: GEN module1,struct1,fact1,liste1,U1,U;
1887: GEN module,liste,cyc,p1,y,u1,x,sarch, gen = add_gen? gun: NULL;
1888:
1889: nf = checknf(nf); checkbid(bid1);
1890: module1 = (GEN)bid1[1]; struct1 = (GEN)bid1[2]; fact1 = (GEN)bid1[3];
1891: x = (GEN)module1[1];
1892: module = cgetg(3,t_VEC);
1893: module[1] = (long)x;
1894: module[2] = (long)arch;
1895: if (!gcmp0((GEN)module1[2]))
1896: err(talker,"non-0 Archimedian components in zidealstarinitjoinarch");
1897:
1898: sarch = zarchstar(nf,x,arch,nba);
1899: liste1 = (GEN)bid1[4]; lx1 = lg(liste1);
1900: liste = cgetg(lx1,t_VEC);
1901: for (i=1; i<lx1-1; i++) liste[i]=liste1[i];
1902: liste[i] = (long)sarch;
1903:
1904: U1 = (GEN)bid1[5]; lx1 = lg(U1);
1905: cyc = diagonal(concatsp((GEN)struct1[2],(GEN)sarch[1]));
1906: cyc = smithrel(cyc, &U, gen? &u1: NULL);
1907:
1908: if (gen)
1909: {
1910: if (lg(struct1)<=3)
1911: err(talker,"please apply idealstar(,,2) and not idealstar(,,1)");
1912: gen = concatsp((GEN)struct1[3],(GEN)sarch[2]);
1913: nbp = lg((GEN)fact1[1])-1;
1914: }
1915: p1 = cgetg(gen? 4: 3, t_VEC);
1916: p1[1] = (long)detcyc(cyc);
1917: p1[2] = (long)cyc;
1918:
1919: y = cgetg(6,t_VEC);
1920: y[1] = (long)module;
1921: y[2] = (long)p1;
1922: y[3] = (long)fact1;
1923: y[4] = (long)liste;
1924: y[5] = (long)U;
1925: if (gen) p1[3] = (long)compute_gen(nf,u1,gen, y);
1926: return gerepilecopy(av,y);
1927: }
1928:
1929: /* calcule la matrice des zinternallog des unites */
1930: GEN
1931: logunitmatrix(GEN nf,GEN funits,GEN racunit,GEN bid)
1932: {
1933: long R,j,sizeh;
1934: GEN m,fa2,fa,arch;
1935:
1936: R=lg(funits)-1; m=cgetg(R+2,t_MAT);
1937: fa2=(GEN)bid[4]; sizeh=lg(bid[5])-1; arch=gmael(bid,1,2);
1938: fa=(GEN)bid[3];
1939: m[1]=(long)zinternallog(nf,racunit,fa2,sizeh,arch,fa,0);
1940: for (j=2; j<=R+1; j++)
1941: m[j]=(long)zinternallog(nf,(GEN)funits[j-1],fa2,sizeh,arch,fa,0);
1942: return m;
1943: }
1944:
1945: /* calcule la matrice des zinternallog des unites */
1946: static GEN
1947: logunitmatrixarch(GEN nf,GEN funits,GEN racunit,GEN bid)
1948: {
1949: long R,j;
1950: GEN m,liste,structarch,arch;
1951:
1952: R=lg(funits)-1; m=cgetg(R+2,t_MAT); arch=gmael(bid,1,2);
1953: liste=(GEN)bid[4]; structarch=(GEN)liste[lg(liste)-1];
1954: m[1]=(long)zsigne(nf,racunit,arch);
1955: for (j=2; j<=R+1; j++)
1956: m[j]=(long)zsigne(nf,(GEN)funits[j-1],arch);
1957: return lift_intern(gmul((GEN)structarch[3],m));
1958: }
1959:
1960: static void
1961: init_units(GEN bnf, GEN *funits, GEN *racunit)
1962: {
1963: GEN p1;
1964: bnf = checkbnf(bnf); p1=(GEN)bnf[8];
1965: if (lg(p1)==5) *funits=(GEN)buchfu(bnf)[1];
1966: else
1967: {
1968: if (lg(p1)!=7) err(talker,"incorrect big number field");
1969: *funits=(GEN)p1[5];
1970: }
1971: *racunit=gmael(p1,4,2);
1972: }
1973:
1974: /* flag &1 : generateurs, sinon non
1975: * flag &2 : unites, sinon pas.
1976: * flag &4 : ideaux, sinon zidealstar.
1977: */
1978: static GEN
1979: ideallistzstarall(GEN bnf,long bound,long flag)
1980: {
1981: byteptr ptdif=diffptr;
1.2 ! noro 1982: gpmem_t lim,av0=avma,av;
! 1983: long i,j,k,l,q2,lp1,q;
1.1 noro 1984: long do_gen = flag & 1, do_units = flag & 2, big_id = !(flag & 4);
1985: GEN y,nf,p,z,z2,p1,p2,p3,fa,pr,ideal,lu,lu2,funits,racunit,embunit;
1986:
1987: nf = checknf(bnf);
1988: if (bound <= 0) return cgetg(1,t_VEC);
1989: z = cgetg(bound+1,t_VEC);
1990: for (i=2; i<=bound; i++) z[i] = lgetg(1,t_VEC);
1991:
1992: ideal = idmat(degpol(nf[1]));
1993: if (big_id) ideal = zidealstarinitall(nf,ideal,do_gen);
1994: z[1] = (long)_vec(ideal);
1995: if (do_units)
1996: {
1997: init_units(bnf,&funits,&racunit);
1998: lu = cgetg(bound+1,t_VEC);
1999: for (i=2; i<=bound; i++) lu[i]=lgetg(1,t_VEC);
2000: lu[1] = (long)_vec(logunitmatrix(nf,funits,racunit,ideal));
2001: }
2002:
2003: p=cgeti(3); p[1]=evalsigne(1) | evallgefint(3);
2004: av=avma; lim=stack_lim(av,1);
2005: lu2 = embunit = NULL; /* gcc -Wall */
2006: if (bound > (long)maxprime()) err(primer1);
2007: for (p[2]=0; p[2]<=bound; )
2008: {
1.2 ! noro 2009: NEXT_PRIME_VIADIFF(p[2], ptdif);
1.1 noro 2010: if (DEBUGLEVEL>1) { fprintferr("%ld ",p[2]); flusherr(); }
2011: fa = primedec(nf,p);
2012: for (j=1; j<lg(fa); j++)
2013: {
2014: pr = (GEN)fa[j]; p1 = powgi(p,(GEN)pr[4]);
2015: if (is_bigint(p1) || (q = itos(p1)) > bound) continue;
2016:
2017: q2=q; ideal=pr; z2=dummycopy(z);
2018: if (do_units) lu2=dummycopy(lu);
2019: for (l=2; ;l++)
2020: {
2021: if (big_id) ideal = zidealstarinitall(nf,ideal,do_gen);
2022: if (do_units) embunit = logunitmatrix(nf,funits,racunit,ideal);
2023: for (i=q; i<=bound; i+=q)
2024: {
2025: p1 = (GEN)z[i/q]; lp1 = lg(p1);
2026: if (lp1 == 1) continue;
2027:
2028: p2 = cgetg(lp1,t_VEC);
2029: for (k=1; k<lp1; k++)
2030: if (big_id)
2031: p2[k] = (long)zidealstarinitjoin(nf,(GEN)p1[k],ideal,do_gen);
2032: else
2033: p2[k] = (long)idealmul(nf,(GEN)p1[k],ideal);
2034: z2[i] = (long)concatsp((GEN)z2[i],p2);
2035: if (do_units)
2036: {
2037: p1 = (GEN)lu[i/q];
2038: p2 = cgetg(lp1,t_VEC);
2039: for (k=1; k<lp1; k++)
2040: p2[k] = (long)vconcat((GEN)p1[k],embunit);
2041: lu2[i] = (long)concatsp((GEN)lu2[i],p2);
2042: }
2043: }
2044: q *= q2; if ((ulong)q > (ulong)bound) break;
2045: ideal = idealpows(nf,pr,l);
2046: }
2047: z = z2; if (do_units) lu = lu2;
2048: }
2049: if (low_stack(lim, stack_lim(av,1)))
2050: {
2051: GEN *gptr[2]; gptr[0]=&z; gptr[1]=&lu;
2052: if(DEBUGMEM>1) err(warnmem,"ideallistzstarall");
2053: gerepilemany(av,gptr,do_units?2:1);
2054: }
2055: }
2056: if (!do_units) return gerepilecopy(av0, z);
2057: y = cgetg(3,t_VEC);
2058: y[1] = lcopy(z);
2059: lu2 = cgetg(lg(z),t_VEC);
2060: for (i=1; i<lg(z); i++)
2061: {
2062: p1=(GEN)z[i]; p2=(GEN)lu[i]; lp1=lg(p1);
2063: p3=cgetg(lp1,t_VEC); lu2[i]=(long)p3;
2064: for (j=1; j<lp1; j++) p3[j] = lmul(gmael(p1,j,5),(GEN)p2[j]);
2065: }
2066: y[2]=(long)lu2; return gerepileupto(av0, y);
2067: }
2068:
2069: GEN
2070: ideallist0(GEN bnf,long bound, long flag)
2071: {
2072: if (flag<0 || flag>4) err(flagerr,"ideallist");
2073: return ideallistzstarall(bnf,bound,flag);
2074: }
2075:
2076: GEN
2077: ideallistzstar(GEN nf,long bound)
2078: {
2079: return ideallistzstarall(nf,bound,0);
2080: }
2081:
2082: GEN
2083: ideallistzstargen(GEN nf,long bound)
2084: {
2085: return ideallistzstarall(nf,bound,1);
2086: }
2087:
2088: GEN
2089: ideallistunit(GEN nf,long bound)
2090: {
2091: return ideallistzstarall(nf,bound,2);
2092: }
2093:
2094: GEN
2095: ideallistunitgen(GEN nf,long bound)
2096: {
2097: return ideallistzstarall(nf,bound,3);
2098: }
2099:
2100: GEN
2101: ideallist(GEN bnf,long bound)
2102: {
2103: return ideallistzstarall(bnf,bound,4);
2104: }
2105:
2106: static GEN
2107: ideallist_arch(GEN nf,GEN list,GEN arch,long flun)
2108: {
2109: long nba,i,j,lx,ly;
2110: GEN p1,z,p2;
2111:
1.2 ! noro 2112: if (typ(arch) != t_VEC) err(typeer,"ideallistarch");
1.1 noro 2113: nba=0; for (i=1; i<lg(arch); i++) if (signe(arch[i])) nba++;
2114: lx=lg(list); z=cgetg(lx,t_VEC);
2115: for (i=1; i<lx; i++)
2116: {
2117: p2=(GEN)list[i]; checkbid(p2); ly=lg(p2);
2118: p1=cgetg(ly,t_VEC); z[i]=(long)p1;
2119: for (j=1; j<ly; j++)
2120: p1[j]=(long)zidealstarinitjoinarch(nf,(GEN)p2[j],arch,nba,flun);
2121: }
2122: return z;
2123: }
2124:
2125: static GEN
2126: ideallistarchall(GEN bnf,GEN list,GEN arch,long flag)
2127: {
1.2 ! noro 2128: gpmem_t av;
1.1 noro 2129: long i,j,lp1, do_units = flag & 2;
2130: GEN nf = checknf(bnf), p1,p2,p3,racunit,funits,lu2,lu,embunit,z,y;
2131:
2132: if (typ(list) != t_VEC || (do_units && lg(list) != 3))
2133: err(typeer, "ideallistarch");
2134: if (lg(list) == 1) return cgetg(1,t_VEC);
2135: if (do_units)
2136: {
2137: y = cgetg(3,t_VEC);
2138: z = (GEN)list[1];
2139: lu= (GEN)list[2]; if (typ(lu) != t_VEC) err(typeer, "ideallistarch");
2140: }
2141: else
2142: {
2143: z = list;
2144: y = lu = NULL; /* gcc -Wall */
2145: }
2146: if (typ(z) != t_VEC) err(typeer, "ideallistarch");
2147: z = ideallist_arch(nf,z,arch, flag & 1);
2148: if (!do_units) return z;
2149:
2150: y[1]=(long)z; av=avma;
2151: init_units(bnf,&funits,&racunit);
2152: lu2=cgetg(lg(z),t_VEC);
2153: for (i=1; i<lg(z); i++)
2154: {
2155: p1=(GEN)z[i]; p2=(GEN)lu[i]; lp1=lg(p1);
2156: p3=cgetg(lp1,t_VEC); lu2[i]=(long)p3;
2157: for (j=1; j<lp1; j++)
2158: {
2159: embunit = logunitmatrixarch(nf,funits,racunit,(GEN)p1[j]);
2160: p3[j] = lmul(gmael(p1,j,5), vconcat((GEN)p2[j],embunit));
2161: }
2162: }
2163: y[2]=lpilecopy(av,lu2); return y;
2164: }
2165:
2166: GEN
2167: ideallistarch(GEN nf, GEN list, GEN arch)
2168: {
2169: return ideallistarchall(nf,list,arch,0);
2170: }
2171:
2172: GEN
2173: ideallistarchgen(GEN nf, GEN list, GEN arch)
2174: {
2175: return ideallistarchall(nf,list,arch,1);
2176: }
2177:
2178: GEN
2179: ideallistunitarch(GEN bnf,GEN list,GEN arch)
2180: {
2181: return ideallistarchall(bnf,list,arch,2);
2182: }
2183:
2184: GEN
2185: ideallistunitarchgen(GEN bnf,GEN list,GEN arch)
2186: {
2187: return ideallistarchall(bnf,list,arch,3);
2188: }
2189:
2190: GEN
2191: ideallistarch0(GEN nf, GEN list, GEN arch,long flag)
2192: {
2193: if (!arch) arch=cgetg(1,t_VEC);
2194: if (flag<0 || flag>3) err(flagerr,"ideallistarch");
2195: return ideallistarchall(nf,list,arch,flag);
2196: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to base3.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath