Annotation of OpenXM_contrib/pari-2.2/src/basemath/base4.c, Revision 1.2
1.2 ! noro 1: /* $Id: base4.c,v 1.143 2002/09/08 10:41:21 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: /* (continued) */
20: /* */
21: /*******************************************************************/
22: #include "pari.h"
23: #include "parinf.h"
24:
1.2 ! noro 25: extern GEN makeprimetoideal(GEN nf,GEN UV,GEN uv,GEN x);
! 26: extern GEN gauss_triangle_i(GEN A, GEN B,GEN t);
! 27: extern GEN hnf_invimage(GEN A, GEN b);
! 28: extern int hnfdivide(GEN A, GEN B);
! 29: extern GEN colreducemodHNF(GEN x, GEN y, GEN *Q);
! 30: extern GEN zinternallog_pk(GEN nf,GEN a0,GEN y,GEN pr,GEN prk,GEN list,GEN *psigne);
! 31: extern GEN special_anti_uniformizer(GEN nf, GEN pr);
! 32: extern GEN set_sign_mod_idele(GEN nf, GEN x, GEN y, GEN idele, GEN sarch);
! 33: extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx);
1.1 noro 34:
1.2 ! noro 35: static GEN mat_ideal_two_elt2(GEN nf, GEN x, GEN a);
1.1 noro 36:
37: /*******************************************************************/
38: /* */
39: /* IDEAL OPERATIONS */
40: /* */
41: /*******************************************************************/
42:
43: /* A valid ideal is either principal (valid nf_element), or prime, or a matrix
44: * on the integer basis (preferably HNF).
45: * A prime ideal is of the form [p,a,e,f,b], where the ideal is p.Z_K+a.Z_K,
46: * p is a rational prime, a belongs to Z_K, e=e(P/p), f=f(P/p), and b
47: * Lenstra constant (p.P^(-1)= p Z_K + b Z_K).
48: *
49: * An idele is a couple[I,V] where I is a valid ideal and V a row vector
50: * with r1+r2 components (real or complex). For instance, if M=(a), V
51: * contains the complex logarithms of the first r1+r2 conjugates of a
52: * (depends on the chosen generator a). All subroutines work with either
53: * ideles or ideals (an omitted V is assumed to be 0).
54: *
1.2 ! noro 55: * All ideals are output in HNF form. */
1.1 noro 56:
57: /* types and conversions */
58:
59: static long
60: idealtyp(GEN *ideal, GEN *arch)
61: {
62: GEN x = *ideal;
63: long t,lx,tx = typ(x);
64:
65: if (tx==t_VEC && lg(x)==3)
66: { *arch = (GEN)x[2]; x = (GEN)x[1]; tx = typ(x); }
67: else
68: *arch = NULL;
69: switch(tx)
70: {
71: case t_MAT: lx = lg(x);
72: if (lx>2) t = id_MAT;
73: else
74: {
75: t = id_PRINCIPAL;
76: x = (lx==2)? (GEN)x[1]: gzero;
77: }
78: break;
79:
80: case t_VEC: if (lg(x)!=6) err(idealer2);
81: t = id_PRIME; break;
82:
83: case t_POL: case t_POLMOD: case t_COL:
84: t = id_PRINCIPAL; break;
85: default:
86: if (tx!=t_INT && !is_frac_t(tx)) err(idealer2);
87: t = id_PRINCIPAL;
88: }
89: *ideal = x; return t;
90: }
91:
92: static GEN
93: prime_to_ideal_aux(GEN nf, GEN vp)
94: {
1.2 ! noro 95: GEN m = eltmul_get_table(nf, (GEN)vp[2]);
! 96: return hnfmodid(m, (GEN)vp[1]);
1.1 noro 97: }
98:
1.2 ! noro 99: /* vp = [a,x,...], a in Z. Return (a,x) [HACK: vp need not be prime] */
1.1 noro 100: GEN
101: prime_to_ideal(GEN nf, GEN vp)
102: {
1.2 ! noro 103: gpmem_t av=avma;
1.1 noro 104: if (typ(vp) == t_INT) return gscalmat(vp, degpol(nf[1]));
105: return gerepileupto(av, prime_to_ideal_aux(nf,vp));
106: }
107:
108: /* x = ideal in matrix form. Put it in hnf. */
109: static GEN
110: idealmat_to_hnf(GEN nf, GEN x)
111: {
112: long rx,i,j,N;
1.2 ! noro 113: GEN m, cx;
1.1 noro 114:
115: N=degpol(nf[1]); rx=lg(x)-1;
116: if (!rx) return gscalmat(gzero,N);
117:
1.2 ! noro 118: x = Q_primitive_part(x, &cx);
1.1 noro 119: if (rx >= N) m = x;
120: else
121: {
122: m=cgetg(rx*N + 1,t_MAT);
123: for (i=1; i<=rx; i++)
124: for (j=1; j<=N; j++)
125: m[(i-1)*N + j] = (long) element_mulid(nf,(GEN)x[i],j);
126: }
1.2 ! noro 127: x = hnfmod(m, detint(m));
! 128: return cx? gmul(x,cx): x;
1.1 noro 129: }
130:
131: int
132: ishnfall(GEN x)
133: {
134: long i,j, lx = lg(x);
135: for (i=2; i<lx; i++)
136: {
137: if (gsigne(gcoeff(x,i,i)) <= 0) return 0;
138: for (j=1; j<i; j++)
139: if (!gcmp0(gcoeff(x,i,j))) return 0;
140: }
141: return (gsigne(gcoeff(x,1,1)) > 0);
142: }
1.2 ! noro 143: /* same x is known to be integral */
! 144: int
! 145: Z_ishnfall(GEN x)
! 146: {
! 147: long i,j, lx = lg(x);
! 148: for (i=2; i<lx; i++)
! 149: {
! 150: if (signe(gcoeff(x,i,i)) <= 0) return 0;
! 151: for (j=1; j<i; j++)
! 152: if (signe(gcoeff(x,i,j))) return 0;
! 153: }
! 154: return (signe(gcoeff(x,1,1)) > 0);
! 155: }
1.1 noro 156:
157: GEN
158: idealhermite_aux(GEN nf, GEN x)
159: {
160: long N,tx,lx;
1.2 ! noro 161: GEN z, cx;
1.1 noro 162:
163: tx = idealtyp(&x,&z);
164: if (tx == id_PRIME) return prime_to_ideal_aux(nf,x);
165: if (tx == id_PRINCIPAL)
166: {
1.2 ! noro 167: x = eltmul_get_table(nf, x);
1.1 noro 168: return idealmat_to_hnf(nf,x);
169: }
170: N=degpol(nf[1]); lx = lg(x);
171: if (lg(x[1]) != N+1) err(idealer2);
172:
173: if (lx == N+1 && ishnfall(x)) return x;
174: if (lx <= N) return idealmat_to_hnf(nf,x);
1.2 ! noro 175: x = Q_primitive_part(x, &cx);
1.1 noro 176: x = hnfmod(x,detint(x));
1.2 ! noro 177: return cx? gmul(x, cx): x;
1.1 noro 178: }
179:
180: GEN
181: idealhermite(GEN nf, GEN x)
182: {
1.2 ! noro 183: gpmem_t av=avma;
1.1 noro 184: GEN p1;
185: nf = checknf(nf); p1 = idealhermite_aux(nf,x);
186: if (p1==x || p1==(GEN)x[1]) return gcopy(p1);
187: return gerepileupto(av,p1);
188: }
189:
1.2 ! noro 190: GEN
! 191: principalideal(GEN nf, GEN x)
1.1 noro 192: {
193: GEN z = cgetg(2,t_MAT);
1.2 ! noro 194:
! 195: nf = checknf(nf);
1.1 noro 196: switch(typ(x))
197: {
198: case t_INT: case t_FRAC: case t_FRACN:
1.2 ! noro 199: x = gscalcol(x, degpol(nf[1])); break;
1.1 noro 200:
201: case t_POLMOD:
202: x = checknfelt_mod(nf,x,"principalideal");
203: /* fall through */
204: case t_POL:
1.2 ! noro 205: x = algtobasis(nf,x);
1.1 noro 206: break;
207:
208: case t_MAT:
209: if (lg(x)!=2) err(typeer,"principalideal");
210: x = (GEN)x[1];
211: case t_COL:
1.2 ! noro 212: if (lg(x)-1==degpol(nf[1])) { x = gcopy(x); break; }
1.1 noro 213: default: err(typeer,"principalideal");
214: }
215: z[1]=(long)x; return z;
216: }
217:
1.2 ! noro 218: static GEN
! 219: mylog(GEN x, long prec)
! 220: {
! 221: if (gcmp0(x)) err(precer,"get_arch");
! 222: return glog(x,prec);
! 223: }
! 224:
! 225: GEN get_arch(GEN nf,GEN x,long prec);
! 226:
! 227: static GEN
! 228: famat_get_arch(GEN nf, GEN x, long prec)
1.1 noro 229: {
1.2 ! noro 230: GEN A, a, g = (GEN)x[1], e = (GEN)x[2];
! 231: long i, l = lg(e);
! 232:
! 233: if (l <= 1)
! 234: return get_arch(nf, gun, prec);
! 235: A = NULL; /* -Wall */
! 236: for (i=1; i<l; i++)
! 237: {
! 238: a = get_arch(nf, (GEN)g[i], prec);
! 239: a = gmul((GEN)e[i], a);
! 240: A = (i == 1)? a: gadd(A, a);
! 241: }
! 242: return A;
1.1 noro 243: }
244:
245: static GEN
1.2 ! noro 246: scalar_get_arch(long R1, long RU, GEN x, long prec)
1.1 noro 247: {
1.2 ! noro 248: GEN v = cgetg(RU+1,t_VEC);
! 249: GEN p1 = glog(x, prec);
! 250: long i;
! 251:
! 252: for (i=1; i<=R1; i++) v[i] = (long)p1;
! 253: if (i <= RU)
! 254: {
! 255: p1 = gmul2n(p1,1);
! 256: for ( ; i<=RU; i++) v[i] = (long)p1;
! 257: }
! 258: return v;
1.1 noro 259: }
260:
1.2 ! noro 261: /* For internal use. Get archimedean components: [e_i log( sigma_i(x) )],
! 262: * with e_i = 1 (resp 2.) for i <= R1 (resp. > R1) */
1.1 noro 263: GEN
264: get_arch(GEN nf,GEN x,long prec)
265: {
1.2 ! noro 266: long i, RU, R1 = nf_get_r1(nf);
! 267: GEN v;
! 268:
! 269: RU = lg(nf[6]) - 1;
! 270: switch(typ(x))
! 271: {
! 272: case t_MAT: return famat_get_arch(nf,x,prec);
! 273: case t_POLMOD:
! 274: case t_POL: x = algtobasis_i(nf,x); /* fall through */
! 275: case t_COL: if (!isnfscalar(x)) break;
! 276: x = (GEN)x[1]; /* fall through */
! 277: default: /* rational number */
! 278: return scalar_get_arch(R1, RU, x, prec);
! 279: }
1.1 noro 280:
281: v = cgetg(RU+1,t_VEC);
282: if (isnfscalar(x)) /* rational number */
283: {
1.2 ! noro 284: GEN p1 = glog((GEN)x[1], prec);
! 285:
! 286: for (i=1; i<=R1; i++) v[i] = (long)p1;
! 287: if (i <= RU)
! 288: {
! 289: p1 = gmul2n(p1,1);
! 290: for ( ; i<=RU; i++) v[i] = (long)p1;
! 291: }
! 292: }
! 293: x = gmul(gmael(nf,5,1),x);
! 294: for (i=1; i<=R1; i++) v[i] = (long)mylog((GEN)x[i],prec);
! 295: for ( ; i<=RU; i++) v[i] = lmul2n(mylog((GEN)x[i],prec),1);
! 296: return v;
! 297: }
! 298:
! 299: GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);
! 300:
! 301: static GEN
! 302: famat_get_arch_real(GEN nf,GEN x,GEN *emb,long prec)
! 303: {
! 304: GEN A, T, a, t, g = (GEN)x[1], e = (GEN)x[2];
! 305: long i, l = lg(e);
! 306:
! 307: if (l <= 1)
! 308: return get_arch_real(nf, gun, emb, prec);
! 309: A = T = NULL; /* -Wall */
! 310: for (i=1; i<l; i++)
! 311: {
! 312: a = get_arch_real(nf, (GEN)g[i], &t, prec);
! 313: if (!a) return NULL;
! 314: a = gmul((GEN)e[i], a);
! 315: t = vecpow(t, (GEN)e[i]);
! 316: if (i == 1) { A = a; T = t; }
! 317: else { A = gadd(A, a); T = vecmul(T, t); }
1.1 noro 318: }
1.2 ! noro 319: *emb = T; return A;
! 320: }
! 321:
! 322: static GEN
! 323: scalar_get_arch_real(long R1, long RU, GEN u, GEN *emb, long prec)
! 324: {
! 325: GEN v, x, p1;
! 326: long i, s;
! 327:
! 328: s = gsigne(u);
! 329: if (!s) err(talker,"0 in get_arch_real");
! 330: x = cgetg(RU+1, t_COL);
! 331: for (i=1; i<=RU; i++) x[i] = (long)u;
! 332:
! 333: v = cgetg(RU+1, t_COL);
! 334: p1 = (s > 0)? glog(u,prec): gzero;
! 335: for (i=1; i<=R1; i++) v[i] = (long)p1;
! 336: if (i <= RU)
1.1 noro 337: {
1.2 ! noro 338: p1 = gmul2n(p1,1);
! 339: for ( ; i<=RU; i++) v[i] = (long)p1;
1.1 noro 340: }
1.2 ! noro 341: *emb = x; return v;
! 342: }
! 343:
! 344: static int
! 345: low_prec(GEN x)
! 346: {
! 347: return gcmp0(x) || (typ(x) == t_REAL && lg(x) == 3);
1.1 noro 348: }
349:
350: /* as above but return NULL if precision problem, and set *emb to the
351: * embeddings of x */
352: GEN
1.2 ! noro 353: get_arch_real(GEN nf, GEN x, GEN *emb, long prec)
1.1 noro 354: {
1.2 ! noro 355: long i, RU, R1 = nf_get_r1(nf);
! 356: GEN v, t;
! 357:
! 358: RU = lg(nf[6])-1;
! 359: switch(typ(x))
! 360: {
! 361: case t_MAT: return famat_get_arch_real(nf,x,emb,prec);
1.1 noro 362:
1.2 ! noro 363: case t_POLMOD:
! 364: case t_POL: x = algtobasis_i(nf,x); /* fall through */
! 365: case t_COL: if (!isnfscalar(x)) break;
! 366: x = (GEN)x[1]; /* fall through */
! 367: default: /* rational number */
! 368: return scalar_get_arch_real(R1, RU, x, emb, prec);
! 369: }
1.1 noro 370: v = cgetg(RU+1,t_COL);
1.2 ! noro 371: x = gmul(gmael(nf,5,1), x);
! 372: for (i=1; i<=R1; i++)
1.1 noro 373: {
1.2 ! noro 374: t = gabs((GEN)x[i],prec); if (low_prec(t)) return NULL;
! 375: v[i] = llog(t,prec);
1.1 noro 376: }
1.2 ! noro 377: for ( ; i<=RU; i++)
1.1 noro 378: {
1.2 ! noro 379: t = gnorm((GEN)x[i]); if (low_prec(t)) return NULL;
! 380: v[i] = llog(t,prec);
1.1 noro 381: }
382: *emb = x; return v;
383: }
384:
385: GEN
386: principalidele(GEN nf, GEN x, long prec)
387: {
1.2 ! noro 388: GEN p1, y = cgetg(3,t_VEC);
! 389: gpmem_t av;
1.1 noro 390:
1.2 ! noro 391: p1 = principalideal(nf,x);
1.1 noro 392: y[1] = (long)p1;
1.2 ! noro 393: av = avma; p1 = get_arch(nf,(GEN)p1[1],prec);
1.1 noro 394: y[2] = lpileupto(av,p1); return y;
395: }
396:
397: /* GP functions */
398:
399: GEN
400: ideal_two_elt0(GEN nf, GEN x, GEN a)
401: {
402: if (!a) return ideal_two_elt(nf,x);
403: return ideal_two_elt2(nf,x,a);
404: }
405:
406: GEN
407: idealpow0(GEN nf, GEN x, GEN n, long flag, long prec)
408: {
409: if (flag) return idealpowred(nf,x,n,prec);
410: return idealpow(nf,x,n);
411: }
412:
413: GEN
414: idealmul0(GEN nf, GEN x, GEN y, long flag, long prec)
415: {
416: if (flag) return idealmulred(nf,x,y,prec);
417: return idealmul(nf,x,y);
418: }
419:
420: GEN
421: idealdiv0(GEN nf, GEN x, GEN y, long flag)
422: {
423: switch(flag)
424: {
425: case 0: return idealdiv(nf,x,y);
426: case 1: return idealdivexact(nf,x,y);
427: default: err(flagerr,"idealdiv");
428: }
429: return NULL; /* not reached */
430: }
431:
432: GEN
433: idealaddtoone0(GEN nf, GEN arg1, GEN arg2)
434: {
435: if (!arg2) return idealaddmultoone(nf,arg1);
436: return idealaddtoone(nf,arg1,arg2);
437: }
438:
439: GEN
440: idealhnf0(GEN nf, GEN a, GEN b)
441: {
1.2 ! noro 442: gpmem_t av;
! 443: GEN x;
1.1 noro 444: if (!b) return idealhermite(nf,a);
445:
446: /* HNF of aZ_K+bZ_K */
1.2 ! noro 447: av = avma; nf = checknf(nf);
! 448: x = concatsp(eltmul_get_table(nf,a), eltmul_get_table(nf,b));
! 449: return gerepileupto(av, idealmat_to_hnf(nf, x));
1.1 noro 450: }
451:
452: GEN
453: idealhermite2(GEN nf, GEN a, GEN b)
454: {
455: return idealhnf0(nf,a,b);
456: }
457:
458: static int
459: ok_elt(GEN x, GEN xZ, GEN y)
460: {
1.2 ! noro 461: gpmem_t av = avma;
1.1 noro 462: int r = gegal(x, hnfmodid(y, xZ));
463: avma = av; return r;
464: }
465:
466: static GEN
467: addmul_col(GEN a, long s, GEN b)
468: {
469: long i,l;
470: if (!s) return a? dummycopy(a): a;
471: if (!a) return gmulsg(s,b);
472: l = lg(a);
473: for (i=1; i<l; i++)
474: if (signe(b[i])) a[i] = laddii((GEN)a[i], mulsi(s, (GEN)b[i]));
475: return a;
476: }
477:
478: /* a <-- a + s * b, all coeffs integers */
479: static GEN
480: addmul_mat(GEN a, long s, GEN b)
481: {
482: long j,l;
483: if (!s) return a? dummycopy(a): a; /* copy otherwise next call corrupts a */
484: if (!a) return gmulsg(s,b);
485: l = lg(a);
486: for (j=1; j<l; j++)
487: (void)addmul_col((GEN)a[j], s, (GEN)b[j]);
488: return a;
489: }
490:
491: /* if x square matrix, assume it is HNF */
492: static GEN
493: mat_ideal_two_elt(GEN nf, GEN x)
494: {
1.2 ! noro 495: GEN y,a,beta,cx,xZ,mul;
! 496: long i,lm, N = degpol(nf[1]);
! 497: gpmem_t av0,av,tetpil;
1.1 noro 498:
1.2 ! noro 499: y = cgetg(3,t_VEC); av = avma;
! 500: if (lg(x[1]) != N+1) err(typeer,"ideal_two_elt");
1.1 noro 501: if (N == 2)
502: {
503: y[1] = lcopy(gcoeff(x,1,1));
504: y[2] = lcopy((GEN)x[2]); return y;
505: }
506:
1.2 ! noro 507: x = Q_primitive_part(x, &cx); if (!cx) cx = gun;
1.1 noro 508: if (lg(x) != N+1) x = idealhermite_aux(nf,x);
1.2 ! noro 509:
1.1 noro 510: xZ = gcoeff(x,1,1);
511: if (gcmp1(xZ))
512: {
513: y[1] = lpilecopy(av,cx);
1.2 ! noro 514: y[2] = (long)gscalcol(cx, N); return y;
1.1 noro 515: }
1.2 ! noro 516: av0 = avma;
1.1 noro 517: a = NULL; /* gcc -Wall */
518: beta= cgetg(N+1, t_VEC);
519: mul = cgetg(N+1, t_VEC); lm = 1; /* = lg(mul) */
520: /* look for a in x such that a O/xZ = x O/xZ */
521: for (i=2; i<=N; i++)
522: {
1.2 ! noro 523: gpmem_t av1 = avma;
! 524: GEN t, y = eltmul_get_table(nf, (GEN)x[i]);
! 525: t = FpM_red(y, xZ);
! 526: if (gcmp0(t)) { avma = av1; continue; }
! 527: if (ok_elt(x,xZ, t)) { a = (GEN)x[i]; break; }
1.1 noro 528: beta[lm]= x[i];
1.2 ! noro 529: /* mul[i] = { canonical generators for x[i] O/xZ as Z-module } */
1.1 noro 530: mul[lm] = (long)t; lm++;
531: }
1.2 ! noro 532: if (i > N)
1.1 noro 533: {
534: GEN z = cgetg(lm, t_VECSMALL);
1.2 ! noro 535: gpmem_t av1;
! 536: ulong c = 0;
1.1 noro 537:
538: setlg(mul, lm);
539: setlg(beta,lm);
1.2 ! noro 540: if (DEBUGLEVEL>3) fprintferr("ideal_two_elt, hard case:\n");
1.1 noro 541: for(av1=avma;;avma=av1)
542: {
1.2 ! noro 543: if (++c == 100) {
! 544: if (DEBUGLEVEL>3) fprintferr("using approximation theorem\n");
! 545: a = mat_ideal_two_elt2(nf, x, xZ); goto END;
! 546: }
1.1 noro 547: for (a=NULL,i=1; i<lm; i++)
548: {
549: long t = (mymyrand() >> (BITS_IN_RANDOM-5)) - 7; /* in [-7,8] */
550: z[i] = t;
551: a = addmul_mat(a, t, (GEN)mul[i]);
552: }
553: /* a = matrix (NOT HNF) of ideal generated by beta.z in O/xZ */
554: if (a && ok_elt(x,xZ, a)) break;
555: }
556: for (a=NULL,i=1; i<lm; i++)
557: a = addmul_col(a, z[i], (GEN)beta[i]);
558: }
1.2 ! noro 559: END:
! 560: a = centermod(a, xZ);
! 561: tetpil = avma;
1.1 noro 562: y[1] = lmul(xZ,cx);
563: y[2] = lmul(a, cx);
564: gerepilemanyvec(av,tetpil,y+1,2); return y;
565: }
566:
1.2 ! noro 567: /* Given an ideal x, returns [a,alpha] such that a is in Q,
! 568: * x = a Z_K + alpha Z_K, alpha in K^*
! 569: * a = 0 or alpha = 0 are possible, but do not try to determine whether
! 570: * x is principal. */
1.1 noro 571: GEN
572: ideal_two_elt(GEN nf, GEN x)
573: {
574: GEN z;
575: long N, tx = idealtyp(&x,&z);
576:
577: nf=checknf(nf);
578: if (tx==id_MAT) return mat_ideal_two_elt(nf,x);
579:
580: N=degpol(nf[1]); z=cgetg(3,t_VEC);
581: if (tx == id_PRINCIPAL)
582: {
583: switch(typ(x))
584: {
585: case t_INT: case t_FRAC: case t_FRACN:
586: z[1]=lcopy(x);
587: z[2]=(long)zerocol(N); return z;
588:
589: case t_POLMOD:
590: x = checknfelt_mod(nf, x, "ideal_two_elt");
591: /* fall through */
592: case t_POL:
593: z[1]=zero; z[2]=(long)algtobasis(nf,x); return z;
594: case t_COL:
595: if (lg(x)==N+1) { z[1]=zero; z[2]=lcopy(x); return z; }
596: }
597: }
598: else if (tx == id_PRIME)
599: {
600: z[1]=lcopy((GEN)x[1]);
601: z[2]=lcopy((GEN)x[2]); return z;
602: }
603: err(typeer,"ideal_two_elt");
604: return NULL; /* not reached */
605: }
606:
607: /* factorization */
608:
609: /* x integral ideal in HNF, return v_p(Nx), *vz = v_p(x \cap Z)
610: * Use x[i,i] | x[1,1], i > 0 */
611: long
612: val_norm(GEN x, GEN p, long *vz)
613: {
614: long i,l = lg(x), v;
615: *vz = v = pvaluation(gcoeff(x,1,1), p, NULL);
616: if (!v) return 0;
617: for (i=2; i<l; i++)
618: v += pvaluation(gcoeff(x,i,i), p, NULL);
619: return v;
620: }
621:
622: /* return factorization of Nx */
623: GEN
624: factor_norm(GEN x)
625: {
626: GEN f = factor(gcoeff(x,1,1)), p = (GEN)f[1], e = (GEN)f[2];
627: long i,k, l = lg(p);
628: for (i=1; i<l; i++)
629: e[i] = (long)val_norm(x,(GEN)p[i], &k);
630: settyp(e, t_VECSMALL); return f;
631: }
632:
633: GEN
634: idealfactor(GEN nf, GEN x)
635: {
1.2 ! noro 636: gpmem_t av,tetpil;
! 637: long tx,i,j,k,lf,lc,N,l,v,vc,e;
1.1 noro 638: GEN f,f1,f2,c1,c2,y1,y2,y,p1,p2,cx,P;
639:
640: tx = idealtyp(&x,&y);
641: if (tx == id_PRIME)
642: {
643: y=cgetg(3,t_MAT);
1.2 ! noro 644: y[1]=(long)_col(gcopy(x));
! 645: y[2]=(long)_col(gun); return y;
1.1 noro 646: }
1.2 ! noro 647: av = avma;
! 648: nf = checknf(nf);
! 649: N = degpol(nf[1]);
! 650: if (tx == id_PRINCIPAL) x = idealhermite_aux(nf, x);
! 651: else
! 652: if (lg(x) != N+1) x = idealmat_to_hnf(nf,x);
1.1 noro 653: if (lg(x)==1) err(talker,"zero ideal in idealfactor");
1.2 ! noro 654: x = Q_primitive_part(x, &cx);
! 655: if (!cx)
1.1 noro 656: {
657: c1 = c2 = NULL; /* gcc -Wall */
1.2 ! noro 658: lc = 1;
1.1 noro 659: }
660: else
661: {
1.2 ! noro 662: f = factor(cx);
1.1 noro 663: c1 = (GEN)f[1];
664: c2 = (GEN)f[2]; lc = lg(c1);
665: }
666: f = factor_norm(x);
667: f1 = (GEN)f[1];
668: f2 = (GEN)f[2]; lf = lg(f1);
669: y1 = cgetg((lf+lc-2)*N+1,t_COL);
670: y2 = cgetg((lf+lc-2)*N+1,t_VECSMALL);
671: k = 1;
672: for (i=1; i<lf; i++)
673: {
674: p1 = primedec(nf,(GEN)f1[i]);
675: l = f2[i]; /* = v_p(Nx) */
1.2 ! noro 676: vc = cx? ggval(cx,(GEN)f1[i]): 0;
1.1 noro 677: for (j=1; j<lg(p1); j++)
678: {
679: P = (GEN)p1[j]; e = itos((GEN)P[3]);
680: v = idealval(nf,x,P);
681: l -= v*itos((GEN)P[4]);
682: v += vc*e; if (!v) continue;
683: y1[k] = (long)P;
684: y2[k] = v; k++;
685: if (l == 0) break; /* now only the content contributes */
686: }
687: if (vc == 0) continue;
688: for (j++; j<lg(p1); j++)
689: {
690: P = (GEN)p1[j]; e = itos((GEN)P[3]);
691: y1[k] = (long)P;
692: y2[k] = vc*e; k++;
693: }
694: }
695: for (i=1; i<lc; i++)
696: {
697: /* p | Nx already treated */
698: if (divise(gcoeff(x,1,1),(GEN)c1[i])) continue;
699: p1 = primedec(nf,(GEN)c1[i]);
700: vc = itos((GEN)c2[i]);
701: for (j=1; j<lg(p1); j++)
702: {
703: P = (GEN)p1[j]; e = itos((GEN)P[3]);
704: y1[k] = (long)P;
705: y2[k] = vc*e; k++;
706: }
707: }
708: tetpil=avma; y=cgetg(3,t_MAT);
709: p1=cgetg(k,t_COL); y[1]=(long)p1;
710: p2=cgetg(k,t_COL); y[2]=(long)p2;
711: for (i=1; i<k; i++) { p1[i]=lcopy((GEN)y1[i]); p2[i]=lstoi(y2[i]); }
712: y = gerepile(av,tetpil,y);
713: return sort_factor_gen(y, cmp_prime_ideal);
714: }
715:
716: /* P prime ideal in primedec format. Return valuation(ix) at P */
717: long
718: idealval(GEN nf, GEN ix, GEN P)
719: {
1.2 ! noro 720: gpmem_t av=avma,av1,lim;
! 721: long N,v,vd,w,e,f,i,j,k, tx = typ(ix);
1.1 noro 722: GEN mul,mat,a,x,y,r,bp,p,pk,cx;
723:
724: nf=checknf(nf); checkprimeid(P);
725: if (is_extscalar_t(tx) || tx==t_COL) return element_val(nf,ix,P);
726: p=(GEN)P[1]; N=degpol(nf[1]);
727: tx = idealtyp(&ix,&a);
1.2 ! noro 728: ix = Q_primitive_part(ix, &cx);
1.1 noro 729: if (tx != id_MAT)
730: ix = idealhermite_aux(nf,ix);
731: else
732: {
733: checkid(ix,N);
734: if (lg(ix) != N+1) ix=idealmat_to_hnf(nf,ix);
735: }
736: e = itos((GEN)P[3]);
737: f = itos((GEN)P[4]);
738: /* 0 <= v_P(ix) <= floor[v_p(Nix) / f] */
739: i = val_norm(ix,p, &k) / f;
740: /* 0 <= ceil[v_P(ix) / e] <= v_p(ix \cap Z) --> v_P <= e * v_p */
741: j = k * e;
742: v = min(i,j); /* v_P(ix) <= v */
1.2 ! noro 743: vd = cx? ggval(cx,p) * e: 0;
1.1 noro 744: if (!v) { avma = av; return vd; }
745:
746: mul = cgetg(N+1,t_MAT); bp=(GEN)P[5];
747: mat = cgetg(N+1,t_MAT);
748: for (j=1; j<=N; j++)
749: {
750: mul[j] = (long)element_mulid(nf,bp,j);
751: x = (GEN)ix[j];
752: y = cgetg(N+1, t_COL); mat[j] = (long)y;
753: for (i=1; i<=N; i++)
754: { /* compute (x.b)_i, ix in HNF ==> x[j+1..N] = 0 */
755: a = mulii((GEN)x[1], gcoeff(mul,i,1));
756: for (k=2; k<=j; k++) a = addii(a, mulii((GEN)x[k], gcoeff(mul,i,k)));
757: /* is it divisible by p ? */
758: y[i] = ldvmdii(a,p,&r);
759: if (signe(r)) { avma = av; return vd; }
760: }
761: }
762: pk = gpowgs(p, v-1);
763: av1 = avma; lim=stack_lim(av1,3);
764: y = cgetg(N+1,t_COL);
765: /* can compute mod p^(v-w) */
766: for (w=1; w<v; w++)
767: {
768: for (j=1; j<=N; j++)
769: {
770: x = (GEN)mat[j];
771: for (i=1; i<=N; i++)
772: { /* compute (x.b)_i */
773: a = mulii((GEN)x[1], gcoeff(mul,i,1));
774: for (k=2; k<=N; k++) a = addii(a, mulii((GEN)x[k], gcoeff(mul,i,k)));
775: /* is it divisible by p ? */
776: y[i] = ldvmdii(a,p,&r);
777: if (signe(r)) { avma = av; return w + vd; }
778: if (lgefint(y[i]) > lgefint(pk)) y[i] = lresii((GEN)y[i], pk);
779: }
780: r=x; mat[j]=(long)y; y=r;
781: if (low_stack(lim,stack_lim(av1,3)))
782: {
783: GEN *gptr[3]; gptr[0]=&y; gptr[1]=&mat; gptr[2]=&pk;
784: if(DEBUGMEM>1) err(warnmem,"idealval");
785: gerepilemany(av1,gptr,3);
786: }
787: }
788: pk = gdivexact(pk,p);
789: }
790: avma = av; return w + vd;
791: }
792:
793: /* gcd and generalized Bezout */
794:
795: GEN
796: idealadd(GEN nf, GEN x, GEN y)
797: {
1.2 ! noro 798: gpmem_t av=avma;
! 799: long N,tx,ty;
1.1 noro 800: GEN z,p1,dx,dy,dz;
801: int modid;
802:
803: tx = idealtyp(&x,&z);
804: ty = idealtyp(&y,&z);
805: nf=checknf(nf); N=degpol(nf[1]);
806: z = cgetg(N+1, t_MAT);
807: if (tx != id_MAT || lg(x)!=N+1) x = idealhermite_aux(nf,x);
808: if (ty != id_MAT || lg(y)!=N+1) y = idealhermite_aux(nf,y);
809: if (lg(x) == 1) return gerepileupto(av,y);
810: if (lg(y) == 1) return gerepileupto(av,x); /* check for 0 ideal */
1.2 ! noro 811: dx = denom(x);
! 812: dy = denom(y); dz = mulii(dx,dy);
! 813: if (gcmp1(dz)) dz = NULL; else {
! 814: x = Q_muli_to_int(x, dz);
! 815: y = Q_muli_to_int(y, dz);
! 816: }
1.1 noro 817: if (isnfscalar((GEN)x[1]) && isnfscalar((GEN)y[1]))
818: {
1.2 ! noro 819: p1 = mppgcd(gcoeff(x,1,1), gcoeff(y,1,1));
1.1 noro 820: modid = 1;
821: }
822: else
823: {
1.2 ! noro 824: p1 = mppgcd(detint(x), detint(y));
1.1 noro 825: modid = 0;
826: }
827: if (gcmp1(p1))
828: {
829: long i,j;
830: if (!dz) { avma=av; return idmat(N); }
1.2 ! noro 831: avma = (gpmem_t)dz; dz = gerepileupto((gpmem_t)z, ginv(dz));
1.1 noro 832: for (i=1; i<=N; i++)
833: {
834: z[i]=lgetg(N+1,t_COL);
835: for (j=1; j<=N; j++)
836: coeff(z,j,i) = (i==j)? (long)dz: zero;
837: }
838: return z;
839: }
840: z = concatsp(x,y);
841: z = modid? hnfmodid(z,p1): hnfmod(z, p1);
842: if (dz) z=gdiv(z,dz);
843: return gerepileupto(av,z);
844: }
845:
1.2 ! noro 846: /* Assume x,y integral non zero (presumably coprime, which is checked).
! 847: * Return a in x such that 1-a in y. If (x + y) \cap Z is 1, then addone_aux
! 848: * is more efficient. */
! 849: GEN
! 850: addone_aux2(GEN x, GEN y)
! 851: {
! 852: GEN U, H;
! 853: long i, l;
! 854:
! 855: H = hnfall_i(concatsp(x,y), &U, 1);
! 856: l = lg(H);
! 857: for (i=1; i<l; i++)
! 858: if (!gcmp1(gcoeff(H,i,i)))
! 859: err(talker,"ideals don't sum to Z_K in idealaddtoone");
! 860: U = (GEN)U[l]; setlg(U, lg(x));
! 861: return gmul(x, U);
! 862: }
! 863:
! 864: /* assume x,y integral, non zero in HNF */
1.1 noro 865: static GEN
1.2 ! noro 866: addone_aux(GEN x, GEN y)
1.1 noro 867: {
1.2 ! noro 868: GEN a, b, u, v;
1.1 noro 869:
1.2 ! noro 870: a = gcoeff(x,1,1);
! 871: b = gcoeff(y,1,1);
! 872: if (typ(a) != t_INT || typ(b) != t_INT)
! 873: err(talker,"ideals don't sum to Z_K in idealaddtoone");
! 874: if (gcmp1(bezout(a,b,&u,&v))) return gmul(u,(GEN)x[1]);
! 875: return addone_aux2(x, y);
1.1 noro 876: }
877:
878: GEN
879: idealaddtoone_i(GEN nf, GEN x, GEN y)
880: {
881: long t, fl = 1;
882: GEN p1,xh,yh;
883:
884: if (DEBUGLEVEL>4)
885: {
886: fprintferr(" entering idealaddtoone:\n");
887: fprintferr(" x = %Z\n",x);
888: fprintferr(" y = %Z\n",y);
889: }
890: t = idealtyp(&x,&p1);
1.2 ! noro 891: if (t != id_MAT || lg(x) == 1 || lg(x) != lg(x[1]))
1.1 noro 892: xh = idealhermite_aux(nf,x);
893: else
1.2 ! noro 894: { xh = x; fl = isnfscalar((GEN)x[1]); }
1.1 noro 895: t = idealtyp(&y,&p1);
896: if (t != id_MAT || lg(y) == 1 || lg(y) != lg(y[1]))
897: yh = idealhermite_aux(nf,y);
898: else
1.2 ! noro 899: { yh = y; if (fl) fl = isnfscalar((GEN)y[1]); }
1.1 noro 900: if (lg(xh) == 1)
901: {
902: if (lg(yh) == 1 || !gcmp1(gabs(gcoeff(yh,1,1),0)))
903: err(talker,"ideals don't sum to Z_K in idealaddtoone");
904: return algtobasis(nf, gzero);
905: }
906: if (lg(yh) == 1)
907: {
908: p1 = gcoeff(xh,1,1);
909: if (!gcmp1(gabs(p1,0)))
910: err(talker,"ideals don't sum to Z_K in idealaddtoone");
911: return algtobasis(nf, gneg(p1));
912: }
913:
1.2 ! noro 914: if (fl) p1 = addone_aux (xh,yh); /* xh[1], yh[1] scalar */
! 915: else p1 = addone_aux2(xh,yh);
1.1 noro 916: p1 = element_reduce(nf,p1, idealmullll(nf,x,y));
917: if (DEBUGLEVEL>4 && !gcmp0(p1))
918: fprintferr(" leaving idealaddtoone: %Z\n",p1);
919: return p1;
920: }
921:
922: /* ideal should be an idele (not mandatory). For internal use. */
923: GEN
924: ideleaddone_aux(GEN nf,GEN x,GEN ideal)
925: {
926: long i,nba,R1;
927: GEN p1,p2,p3,arch;
928:
929: (void)idealtyp(&ideal,&arch);
930: if (!arch) return idealaddtoone_i(nf,x,ideal);
931:
932: R1=itos(gmael(nf,2,1));
933: if (typ(arch)!=t_VEC && lg(arch)!=R1+1)
934: err(talker,"incorrect idele in idealaddtoone");
935: for (nba=0,i=1; i<lg(arch); i++)
936: if (signe(arch[i])) nba++;
937: if (!nba) return idealaddtoone_i(nf,x,ideal);
938:
939: p3 = idealaddtoone_i(nf,x,ideal);
940: if (gcmp0(p3)) p3=(GEN)idealhermite_aux(nf,x)[1];
941: p1=idealmullll(nf,x,ideal);
942:
943: p2=zarchstar(nf,p1,arch,nba);
944: p1=lift_intern(gmul((GEN)p2[3],zsigne(nf,p3,arch)));
945: p2=(GEN)p2[2]; nba=0;
946: for (i=1; i<lg(p1); i++)
947: if (signe(p1[i])) { nba=1; p3=element_mul(nf,p3,(GEN)p2[i]); }
948: if (gcmp0(p3)) return gcopy((GEN)x[1]); /* can happen if ideal = Z_K */
949: return nba? p3: gcopy(p3);
950: }
951:
1.2 ! noro 952: GEN
1.1 noro 953: unnf_minus_x(GEN x)
954: {
955: long i, N = lg(x);
956: GEN y = cgetg(N,t_COL);
957:
958: y[1] = lsub(gun,(GEN)x[1]);
959: for (i=2;i<N; i++) y[i] = lneg((GEN)x[i]);
960: return y;
961: }
962:
963: static GEN
964: addone(GEN f(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
965: {
1.2 ! noro 966: GEN z = cgetg(3,t_VEC), a;
! 967: gpmem_t av = avma;
! 968: nf = checknf(nf);
! 969: a = gerepileupto(av, f(nf,x,y));
! 970: z[1]=(long)a;
! 971: z[2]=(long)unnf_minus_x(a); return z;
! 972: }
1.1 noro 973:
1.2 ! noro 974: /* assume x,y HNF, non-zero */
! 975: static GEN
! 976: addone_nored(GEN x, GEN y)
! 977: {
! 978: GEN z = cgetg(3,t_VEC), a;
! 979: gpmem_t av = avma;
! 980: a = gerepileupto(av, addone_aux(x,y));
! 981: z[1] = (long)a;
! 982: z[2] = (long)unnf_minus_x(a); return z;
1.1 noro 983: }
984:
985: GEN
986: idealaddtoone(GEN nf, GEN x, GEN y)
987: {
988: return addone(idealaddtoone_i,nf,x,y);
989: }
990:
991: GEN
992: ideleaddone(GEN nf,GEN x,GEN idele)
993: {
994: return addone(ideleaddone_aux,nf,x,idele);
995: }
996:
997: /* given an element x in Z_K and an integral ideal y with x, y coprime,
998: outputs an element inverse of x modulo y */
999: GEN
1000: element_invmodideal(GEN nf, GEN x, GEN y)
1001: {
1.2 ! noro 1002: gpmem_t av=avma;
! 1003: long N,i, fl = 1;
1.1 noro 1004: GEN v,p1,xh,yh;
1005:
1006: nf=checknf(nf); N=degpol(nf[1]);
1.2 ! noro 1007: if (gcmp1(gcoeff(y,1,1))) return zerocol(N);
1.1 noro 1008: i = lg(y);
1009: if (typ(y)!=t_MAT || i==1 || i != lg(y[1])) yh=idealhermite_aux(nf,y);
1010: else
1011: { yh=y; fl = isnfscalar((GEN)y[1]); }
1012: switch (typ(x))
1013: {
1014: case t_POL: case t_POLMOD: case t_COL:
1015: xh = idealhermite_aux(nf,x); break;
1016: default: err(typeer,"element_invmodideal");
1017: return NULL; /* not reached */
1018: }
1.2 ! noro 1019: if (fl) p1 = addone_aux (xh,yh);
! 1020: else p1 = addone_aux2(xh,yh);
1.1 noro 1021: p1 = element_div(nf,p1,x);
1022: v = gerepileupto(av, nfreducemodideal(nf,p1,y));
1023: return v;
1024: }
1025:
1026: GEN
1027: idealaddmultoone(GEN nf, GEN list)
1028: {
1.2 ! noro 1029: gpmem_t av=avma,tetpil;
! 1030: long N,i,i1,j,k;
1.1 noro 1031: GEN z,v,v1,v2,v3,p1;
1032:
1033: nf=checknf(nf); N=degpol(nf[1]);
1034: if (DEBUGLEVEL>4)
1035: {
1036: fprintferr(" entree dans idealaddmultoone() :\n");
1037: fprintferr(" list = "); outerr(list);
1038: }
1039: if (typ(list)!=t_VEC && typ(list)!=t_COL)
1040: err(talker,"not a list in idealaddmultoone");
1041: k=lg(list); z=cgetg(1,t_MAT); list = dummycopy(list);
1042: if (k==1) err(talker,"ideals don't sum to Z_K in idealaddmultoone");
1043: for (i=1; i<k; i++)
1044: {
1045: p1=(GEN)list[i];
1046: if (typ(p1)!=t_MAT || lg(p1)!=lg(p1[1]))
1047: list[i] = (long)idealhermite_aux(nf,p1);
1048: z = concatsp(z,(GEN)list[i]);
1049: }
1050: v=hnfperm(z); v1=(GEN)v[1]; v2=(GEN)v[2]; v3=(GEN)v[3]; j=0;
1051: for (i=1; i<=N; i++)
1052: {
1053: if (!gcmp1(gcoeff(v1,i,i)))
1054: err(talker,"ideals don't sum to Z_K in idealaddmultoone");
1055: if (gcmp1((GEN)v3[i])) j=i;
1056: }
1057:
1058: v=(GEN)v2[(k-2)*N+j]; z=cgetg(k,t_VEC);
1059: for (i=1; i<k; i++)
1060: {
1061: p1=cgetg(N+1,t_COL); z[i]=(long)p1;
1062: for (i1=1; i1<=N; i1++) p1[i1]=v[(i-1)*N+i1];
1063: }
1064: tetpil=avma; v=cgetg(k,typ(list));
1065: for (i=1; i<k; i++) v[i]=lmul((GEN)list[i],(GEN)z[i]);
1066: if (DEBUGLEVEL>2)
1067: { fprintferr(" sortie de idealaddmultoone v = "); outerr(v); }
1068: return gerepile(av,tetpil,v);
1069: }
1070:
1071: /* multiplication */
1072:
1073: /* x integral ideal (without archimedean component) in HNF form
1.2 ! noro 1074: * y = [a,alpha] corresponds to the ideal aZ_K+alpha Z_K (a is a
1.1 noro 1075: * rational integer). Multiply them
1076: */
1077: static GEN
1.2 ! noro 1078: idealmulspec(GEN nf, GEN x, GEN y)
1.1 noro 1079: {
1080: long i, N=lg(x)-1;
1.2 ! noro 1081: GEN m, mod, a = (GEN)y[1], alpha = (GEN)y[2];
1.1 noro 1082:
1083: if (isnfscalar(alpha))
1.2 ! noro 1084: return gmul(mppgcd(a, (GEN)alpha[1]),x);
1.1 noro 1085: mod = mulii(a, gcoeff(x,1,1));
1086: m = cgetg((N<<1)+1,t_MAT);
1087: for (i=1; i<=N; i++) m[i]=(long)element_muli(nf,alpha,(GEN)x[i]);
1088: for (i=1; i<=N; i++) m[i+N]=lmul(a,(GEN)x[i]);
1089: return hnfmodid(m,mod);
1090: }
1091:
1092: /* x ideal (matrix form,maximal rank), vp prime ideal (primedec). Output the
1093: * product. Can be used for arbitrary vp of the form [p,a,e,f,b], IF vp
1094: * =pZ_K+aZ_K, p is an integer, and norm(vp) = p^f; e and b are not used.
1.2 ! noro 1095: * For internal use. */
1.1 noro 1096: GEN
1097: idealmulprime(GEN nf, GEN x, GEN vp)
1098: {
1.2 ! noro 1099: GEN cx;
! 1100: x = Q_primitive_part(x, &cx);
! 1101: x = idealmulspec(nf,x,vp);
! 1102: return cx? gmul(x,cx): x;
! 1103: }
1.1 noro 1104:
1.2 ! noro 1105: GEN arch_mul(GEN x, GEN y);
! 1106: GEN vecpow(GEN x, GEN n);
! 1107: GEN vecmul(GEN x, GEN y);
! 1108:
! 1109: static GEN
! 1110: mul(GEN nf, GEN x, GEN y)
! 1111: {
! 1112: GEN yZ = gcoeff(y,1,1);
! 1113: if (is_pm1(yZ)) return gcopy(x);
! 1114: y = mat_ideal_two_elt(nf,y);
! 1115: return idealmulspec(nf, x, y);
1.1 noro 1116: }
1117:
1118: /* Assume ix and iy are integral in HNF form (or ideles of the same form).
1119: * HACK: ideal in iy can be of the form [a,b], a in Z, b in Z_K
1120: * For internal use. */
1121: GEN
1122: idealmulh(GEN nf, GEN ix, GEN iy)
1123: {
1124: long f = 0;
1125: GEN res,x,y;
1126:
1127: if (typ(ix)==t_VEC) {f=1; x=(GEN)ix[1];} else x=ix;
1128: if (typ(iy)==t_VEC && typ(iy[1]) != t_INT) {f+=2; y=(GEN)iy[1];} else y=iy;
1129: res = f? cgetg(3,t_VEC): NULL;
1130:
1.2 ! noro 1131: if (typ(y) == t_VEC)
! 1132: y = idealmulspec(nf,x,y);
! 1133: else
! 1134: {
! 1135: GEN xZ = gcoeff(x,1,1);
! 1136: GEN yZ = gcoeff(x,1,1);
! 1137: y = (cmpii(xZ, yZ) < 0)? mul(nf,y,x): mul(nf,x,y);
! 1138: }
1.1 noro 1139: if (!f) return y;
1140:
1.2 ! noro 1141: res[1] = (long)y;
! 1142: if (f==3) y = arch_mul((GEN)ix[2], (GEN)iy[2]);
1.1 noro 1143: else
1144: {
1145: y = (f==2)? (GEN)iy[2]: (GEN)ix[2];
1146: y = gcopy(y);
1147: }
1.2 ! noro 1148: res[2] = (long)y; return res;
! 1149: }
! 1150:
! 1151: GEN
! 1152: mul_content(GEN cx, GEN cy)
! 1153: {
! 1154: if (!cx) return cy;
! 1155: if (!cy) return cx;
! 1156: return gmul(cx,cy);
1.1 noro 1157: }
1158:
1159: /* x and y are ideals in matrix form */
1160: static GEN
1161: idealmat_mul(GEN nf, GEN x, GEN y)
1162: {
1.2 ! noro 1163: long i,j, rx = lg(x)-1, ry = lg(y)-1;
! 1164: GEN cx, cy, m;
1.1 noro 1165:
1.2 ! noro 1166: x = Q_primitive_part(x, &cx);
! 1167: y = Q_primitive_part(y, &cy); cx = mul_content(cx,cy);
! 1168: if (rx <= 2 || ry <= 2)
1.1 noro 1169: {
1.2 ! noro 1170: m = cgetg(rx*ry+1,t_MAT);
1.1 noro 1171: for (i=1; i<=rx; i++)
1172: for (j=1; j<=ry; j++)
1.2 ! noro 1173: m[(i-1)*ry+j] = (long)element_muli(nf,(GEN)x[i],(GEN)y[j]);
1.1 noro 1174: if (isnfscalar((GEN)x[1]) && isnfscalar((GEN)y[1]))
1175: {
1176: GEN p1 = mulii(gcoeff(x,1,1),gcoeff(y,1,1));
1177: y = hnfmodid(m, p1);
1178: }
1179: else
1.2 ! noro 1180: y = hnfmod(m, detint(m));
1.1 noro 1181: }
1182: else
1183: {
1.2 ! noro 1184: if (!Z_ishnfall(x)) x = idealmat_to_hnf(nf,x);
! 1185: if (!Z_ishnfall(y)) y = idealmat_to_hnf(nf,y);
! 1186: y = idealmulh(nf,x,y);
1.1 noro 1187: }
1.2 ! noro 1188: return cx? gmul(y,cx): y;
1.1 noro 1189: }
1190:
1.2 ! noro 1191: /* operations on elements in factored form */
! 1192:
! 1193: GEN
! 1194: to_famat(GEN g, GEN e)
! 1195: {
! 1196: GEN h = cgetg(3,t_MAT);
! 1197: if (typ(g) != t_COL) { g = dummycopy(g); settyp(g, t_COL); }
! 1198: if (typ(e) != t_COL) { e = dummycopy(e); settyp(e, t_COL); }
! 1199: h[1] = (long)g;
! 1200: h[2] = (long)e; return h;
! 1201: }
! 1202:
! 1203: GEN
! 1204: to_famat_all(GEN x, GEN y) { return to_famat(_col(x), _col(y)); }
! 1205:
1.1 noro 1206: static GEN
1.2 ! noro 1207: append(GEN v, GEN x)
1.1 noro 1208: {
1.2 ! noro 1209: long i, l = lg(v);
! 1210: GEN w = cgetg(l+1, typ(v));
! 1211: for (i=1; i<l; i++) w[i] = lcopy((GEN)v[i]);
! 1212: w[i] = lcopy(x); return w;
1.1 noro 1213: }
1214:
1215: /* add x^1 to factorisation f */
1216: static GEN
1217: famat_add(GEN f, GEN x)
1218: {
1219: GEN t,h = cgetg(3,t_MAT);
1220: if (lg(f) == 1)
1221: {
1222: t=cgetg(2,t_COL); h[1]=(long)t; t[1]=lcopy(x);
1223: t=cgetg(2,t_COL); h[2]=(long)t; t[1]=un;
1224: }
1225: else
1226: {
1.2 ! noro 1227: h[1] = (long)append((GEN)f[1], x); /* x may be a t_COL */
1.1 noro 1228: h[2] = (long)concat((GEN)f[2], gun);
1229: }
1230: return h;
1231: }
1232:
1233: /* cf merge_factor_i */
1.2 ! noro 1234: GEN
1.1 noro 1235: famat_mul(GEN f, GEN g)
1236: {
1237: GEN h;
1238: if (typ(g) != t_MAT) return famat_add(f, g);
1239: if (lg(f) == 1) return gcopy(g);
1240: if (lg(g) == 1) return gcopy(f);
1241: h = cgetg(3,t_MAT);
1242: h[1] = (long)concat((GEN)f[1], (GEN)g[1]);
1243: h[2] = (long)concat((GEN)f[2], (GEN)g[2]);
1244: return h;
1245: }
1246:
1247: static GEN
1248: famat_sqr(GEN f)
1249: {
1250: GEN h;
1251: if (lg(f) == 1) return cgetg(1,t_MAT);
1252: h = cgetg(3,t_MAT);
1253: h[1] = lcopy((GEN)f[1]);
1254: h[2] = lmul2n((GEN)f[2],1);
1255: return h;
1256: }
1.2 ! noro 1257: GEN
1.1 noro 1258: famat_inv(GEN f)
1259: {
1260: GEN h;
1261: if (lg(f) == 1) return cgetg(1,t_MAT);
1262: h = cgetg(3,t_MAT);
1263: h[1] = lcopy((GEN)f[1]);
1264: h[2] = lneg((GEN)f[2]);
1265: return h;
1266: }
1.2 ! noro 1267: GEN
1.1 noro 1268: famat_pow(GEN f, GEN n)
1269: {
1270: GEN h;
1271: if (lg(f) == 1) return cgetg(1,t_MAT);
1.2 ! noro 1272: if (typ(f) != t_MAT) return to_famat_all(f,n);
1.1 noro 1273: h = cgetg(3,t_MAT);
1274: h[1] = lcopy((GEN)f[1]);
1275: h[2] = lmul((GEN)f[2],n);
1276: return h;
1277: }
1278:
1279: GEN
1280: famat_to_nf(GEN nf, GEN f)
1281: {
1282: GEN t, *x, *e;
1283: long i;
1284: if (lg(f) == 1) return gun;
1285:
1286: x = (GEN*)f[1];
1287: e = (GEN*)f[2];
1288: t = element_pow(nf, x[1], e[1]);
1289: for (i=lg(x)-1; i>1; i--)
1290: t = element_mul(nf, t, element_pow(nf, x[i], e[i]));
1291: return t;
1292: }
1293:
1.2 ! noro 1294: /* "compare" two nf elt. Goal is to quickly sort for uniqueness of
! 1295: * representation, not uniqueness of represented element ! */
! 1296: static int
! 1297: elt_cmp(GEN x, GEN y)
! 1298: {
! 1299: long tx = typ(x), ty = typ(y);
! 1300: if (ty == tx)
! 1301: return (tx == t_POL || tx == t_POLMOD)? cmp_pol(x,y): lexcmp(x,y);
! 1302: return tx - ty;
! 1303: }
! 1304: static int
! 1305: elt_egal(GEN x, GEN y)
1.1 noro 1306: {
1.2 ! noro 1307: if (typ(x) == typ(y)) return gegal(x,y);
! 1308: return 0;
1.1 noro 1309: }
1310:
1311: GEN
1.2 ! noro 1312: famat_reduce(GEN fa)
! 1313: {
! 1314: GEN E, F, G, L, g, e;
! 1315: long i, k, l;
! 1316:
! 1317: if (lg(fa) == 1) return fa;
! 1318: g = (GEN)fa[1]; l = lg(g);
! 1319: e = (GEN)fa[2];
! 1320: L = gen_sort(g, cmp_IND|cmp_C, &elt_cmp);
! 1321: G = cgetg(l, t_COL);
! 1322: E = cgetg(l, t_COL);
! 1323: /* merge */
! 1324: for (k=i=1; i<l; i++,k++)
! 1325: {
! 1326: G[k] = g[L[i]];
! 1327: E[k] = e[L[i]];
! 1328: if (k > 1 && elt_egal((GEN)G[k], (GEN)G[k-1]))
! 1329: {
! 1330: E[k-1] = laddii((GEN)E[k], (GEN)E[k-1]);
! 1331: k--;
! 1332: }
! 1333: }
! 1334: /* kill 0 exponents */
! 1335: l = k;
! 1336: for (k=i=1; i<l; i++)
! 1337: if (!gcmp0((GEN)E[i]))
! 1338: {
! 1339: G[k] = G[i];
! 1340: E[k] = E[i]; k++;
! 1341: }
! 1342: F = cgetg(3, t_MAT);
! 1343: setlg(G, k); F[1] = (long)G;
! 1344: setlg(E, k); F[2] = (long)E; return F;
! 1345: }
1.1 noro 1346:
1.2 ! noro 1347: /* assume (num(g[i]), id) = 1 for all i. Return prod g[i]^e[i] mod id */
1.1 noro 1348: GEN
1349: famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN id)
1350: {
1351: GEN t = NULL, ch,h,n,z,idZ = gcoeff(id,1,1);
1352: long i, lx = lg(g);
1353: if (is_pm1(idZ)) lx = 1; /* id = Z_K */
1354: for (i=1; i<lx; i++)
1355: {
1356: n = (GEN)e[i]; if (!signe(n)) continue;
1357: h = (GEN)g[i]; ch = denom(h);
1358: if (!is_pm1(ch))
1359: {
1360: h = gmul(h,ch); ch = mpinvmod(ch,idZ);
1361: h = gmod(gmul(h,ch), idZ);
1362: }
1363: z = element_powmodideal(nf, h, n, id);
1364: t = (t == NULL)? z: element_mulmodideal(nf, t, z, id);
1365: }
1366: return t? t: gscalcol(gun, lg(id)-1);
1367: }
1368:
1.2 ! noro 1369: /* assume pr has degree 1 and coprime to numerator(x) */
! 1370: static GEN
! 1371: nf_to_Fp_simple(GEN x, GEN modpr, GEN p)
1.1 noro 1372: {
1.2 ! noro 1373: GEN c, r = zk_to_ff(Q_primitive_part(x, &c), modpr);
! 1374: if (c) r = gmod(gmul(r, c), p);
! 1375: return r;
1.1 noro 1376: }
1377:
1.2 ! noro 1378: static GEN
! 1379: famat_to_Fp_simple(GEN nf, GEN x, GEN modpr, GEN p)
1.1 noro 1380: {
1.2 ! noro 1381: GEN h, n, t = gun, g = (GEN)x[1], e = (GEN)x[2], q = subis(p,1);
! 1382: long i, l = lg(g);
! 1383:
! 1384: for (i=1; i<l; i++)
1.1 noro 1385: {
1386: n = (GEN)e[i]; n = modii(n,q);
1387: if (!signe(n)) continue;
1.2 ! noro 1388:
! 1389: h = (GEN)g[i];
! 1390: switch(typ(h))
! 1391: {
! 1392: case t_POL: case t_POLMOD: h = algtobasis(nf, h); /* fall through */
! 1393: case t_COL: h = nf_to_Fp_simple(h, modpr, p); break;
! 1394: default: h = gmod(h, p);
! 1395: }
1.1 noro 1396: t = mulii(t, powmodulo(h, n, p)); /* not worth reducing */
1397: }
1398: return modii(t, p);
1399: }
1400:
1.2 ! noro 1401: /* cf famat_to_nf_modideal_coprime, but id is a prime of degree 1 (=pr) */
1.1 noro 1402: GEN
1.2 ! noro 1403: to_Fp_simple(GEN nf, GEN x, GEN pr)
1.1 noro 1404: {
1.2 ! noro 1405: GEN T, p, modpr = zk_to_ff_init(nf, &pr, &T, &p);
1.1 noro 1406: switch(typ(x))
1407: {
1.2 ! noro 1408: case t_COL: return nf_to_Fp_simple(x,modpr,p);
! 1409: case t_MAT: return famat_to_Fp_simple(nf,x,modpr,p);
1.1 noro 1410: default: err(impl,"generic conversion to finite field");
1411: }
1412: return NULL;
1413: }
1414:
1.2 ! noro 1415: /* Compute t = prod g[i]^e[i] mod pr^k, assuming (t, pr) = 1.
1.1 noro 1416: * Method: modify each g[i] so that it becomes coprime to pr :
1417: * x / (p^k u) --> x * (b/p)^v_pr(x) / z^k u, where z = b^e/p^(e-1)
1418: * b/p = vp^(-1) times something prime to p; both numerator and denominator
1419: * are integral and coprime to pr. Globally, we multiply by (b/p)^v_pr(t) = 1.
1420: *
1421: * EX = exponent of (O_K / pr^k)^* used to reduce the product in case the
1422: * e[i] are large */
1423: static GEN
1.2 ! noro 1424: famat_makecoprime(GEN nf, GEN g, GEN e, GEN pr, GEN prk, GEN EX)
1.1 noro 1425: {
1.2 ! noro 1426: long i, l = lg(g);
! 1427: GEN prkZ,cx,x,u, zpow = gzero, p = (GEN)pr[1], b = (GEN)pr[5];
! 1428: GEN mul = eltmul_get_table(nf, b);
1.1 noro 1429: GEN newg = cgetg(l+1, t_VEC); /* room for z */
1430:
1.2 ! noro 1431: prkZ = gcoeff(prk, 1,1);
1.1 noro 1432: for (i=1; i < l; i++)
1433: {
1434: x = (GEN)g[i];
1435: if (typ(x) != t_COL) x = algtobasis(nf, x);
1.2 ! noro 1436: x = Q_remove_denom(x, &cx);
! 1437: if (cx)
! 1438: {
! 1439: long k = pvaluation(cx, p, &u);
! 1440: if (!gcmp1(u)) /* could avoid the inversion, but prkZ is small--> cheap */
! 1441: x = gmul(x, mpinvmod(u, prkZ));
! 1442: if (k)
! 1443: zpow = addii(zpow, mulsi(k, (GEN)e[i]));
! 1444: }
! 1445: (void)int_elt_val(nf, x, p, mul, &x);
! 1446: newg[i] = (long)colreducemodHNF(x, prk, NULL);
1.1 noro 1447: }
1448: if (zpow == gzero) setlg(newg, l);
1449: else
1450: {
1.2 ! noro 1451: newg[i] = (long)FpV_red(special_anti_uniformizer(nf, pr), prkZ);
1.1 noro 1452: e = concatsp(e, negi(zpow));
1453: }
1454: e = gmod(e, EX);
1.2 ! noro 1455: return famat_to_nf_modideal_coprime(nf, newg, e, prk);
1.1 noro 1456: }
1457:
1458: GEN
1459: famat_ideallog(GEN nf, GEN g, GEN e, GEN bid)
1460: {
1.2 ! noro 1461: gpmem_t av = avma;
1.1 noro 1462: GEN vp = gmael(bid, 3,1), ep = gmael(bid, 3,2), arch = gmael(bid,1,2);
1463: GEN cyc = gmael(bid,2,2), list_set = (GEN)bid[4], U = (GEN)bid[5];
1464: GEN p1,y0,x,y, psigne;
1.2 ! noro 1465: long i, l;
1.1 noro 1466: if (lg(cyc) == 1) return cgetg(1,t_COL);
1467: y0 = y = cgetg(lg(U), t_COL);
1468: psigne = zsigne(nf, to_famat(g,e), arch);
1.2 ! noro 1469: l = lg(vp);
! 1470: for (i=1; i < l; i++)
1.1 noro 1471: {
1472: GEN pr = (GEN)vp[i], prk;
1.2 ! noro 1473: prk = (l==2)? gmael(bid,1,1): idealpow(nf, pr, (GEN)ep[i]);
1.1 noro 1474: /* TODO: FIX group exponent (should be mod prk, not f !) */
1475: x = famat_makecoprime(nf, g, e, pr, prk, (GEN)cyc[1]);
1476: y = zinternallog_pk(nf, x, y, pr, prk, (GEN)list_set[i], &psigne);
1477: }
1478: p1 = lift_intern(gmul(gmael(list_set,i,3), psigne));
1479: for (i=1; i<lg(p1); i++) *++y = p1[i];
1480: y = gmul(U,y0);
1481: avma = av; x = cgetg(lg(y), t_COL);
1482: for (i=1; i<lg(y); i++)
1483: x[i] = lmodii((GEN)y[i], (GEN)cyc[i]);
1484: return x;
1485: }
1486:
1487: /* prod g[i]^e[i] mod bid, assume (g[i], id) = 1 */
1488: GEN
1489: famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid)
1490: {
1491: GEN t,sarch,module,cyc,fa2;
1492: long lc;
1493: if (lg(g) == 1) return gscalcol_i(gun, degpol(nf[1])); /* 1 */
1494: module = (GEN)bid[1];
1495: fa2 = (GEN)bid[4]; sarch = (GEN)fa2[lg(fa2)-1];
1496: cyc = gmael(bid,2,2); lc = lg(cyc);
1497: t = NULL;
1498: if (lc != 1)
1499: {
1500: GEN EX = (GEN)cyc[1]; /* group exponent */
1501: GEN id = (GEN)module[1];
1502: t = famat_to_nf_modideal_coprime(nf,g, gmod(e,EX), id);
1503: }
1504: if (!t) t = gun;
1505: return set_sign_mod_idele(nf, to_famat(g,e), t, module, sarch);
1506: }
1507:
1508: GEN
1.2 ! noro 1509: famat_to_arch(GEN nf, GEN fa, long prec)
! 1510: {
! 1511: GEN g,e, y = NULL;
! 1512: long i,l;
! 1513:
! 1514: if (typ(fa) != t_MAT) return get_arch(nf, fa, prec);
! 1515: if (lg(fa) == 1) return zerovec(lg(nf[6])-1);
! 1516: g = (GEN)fa[1];
! 1517: e = (GEN)fa[2]; l = lg(e);
! 1518: for (i=1; i<l; i++)
! 1519: {
! 1520: GEN t, x = (GEN)g[i];
! 1521: if (typ(x) != t_COL) x = algtobasis(nf,x);
! 1522: x = Q_primpart(x);
! 1523: /* multiplicative arch would be better (save logs), but exponents overflow
! 1524: * [ could keep track of expo separately, but not worth it ] */
! 1525: t = gmul(get_arch(nf,x,prec), (GEN)e[i]);
! 1526: y = y? gadd(y,t): t;
! 1527: }
! 1528: return y;
! 1529: }
! 1530:
! 1531: GEN
1.1 noro 1532: vecmul(GEN x, GEN y)
1533: {
1534: long i,lx, tx = typ(x);
1535: GEN z;
1536: if (is_scalar_t(tx)) return gmul(x,y);
1537: lx = lg(x); z = cgetg(lx,tx);
1538: for (i=1; i<lx; i++) z[i] = (long)vecmul((GEN)x[i], (GEN)y[i]);
1539: return z;
1540: }
1541:
1542: GEN
1543: vecinv(GEN x)
1544: {
1545: long i,lx, tx = typ(x);
1546: GEN z;
1547: if (is_scalar_t(tx)) return ginv(x);
1548: lx = lg(x); z = cgetg(lx, tx);
1549: for (i=1; i<lx; i++) z[i] = (long)vecinv((GEN)x[i]);
1550: return z;
1551: }
1552:
1553: GEN
1.2 ! noro 1554: vecpow(GEN x, GEN n)
! 1555: {
! 1556: long i,lx, tx = typ(x);
! 1557: GEN z;
! 1558: if (is_scalar_t(tx)) return powgi(x,n);
! 1559: lx = lg(x); z = cgetg(lx, tx);
! 1560: for (i=1; i<lx; i++) z[i] = (long)vecpow((GEN)x[i], n);
! 1561: return z;
! 1562: }
! 1563:
! 1564: GEN
1.1 noro 1565: vecdiv(GEN x, GEN y) { return vecmul(x, vecinv(y)); }
1566:
1567: /* x,y assumed to be of the same type, either
1568: * t_VEC: logarithmic distance components
1569: * t_COL: multiplicative distance components [FIXME: find decent type]
1570: * t_POLMOD: nf elt
1571: * t_MAT: factorisation of nf elt */
1572: GEN
1573: arch_mul(GEN x, GEN y) {
1574: switch (typ(x)) {
1575: case t_POLMOD: return gmul(x, y);
1576: case t_COL: return vecmul(x, y);
1577: case t_MAT: return (x == y)? famat_sqr(x): famat_mul(x,y);
1578: default: return (x == y)? gmul2n(x,1): gadd(x,y); /* t_VEC */
1579: }
1580: }
1581:
1582: GEN
1583: arch_inv(GEN x) {
1584: switch (typ(x)) {
1585: case t_POLMOD: return ginv(x);
1.2 ! noro 1586: case t_COL: return vecinv(x);
1.1 noro 1587: case t_MAT: return famat_inv(x);
1.2 ! noro 1588: default: return gneg(x); /* t_VEC */
1.1 noro 1589: }
1590: }
1591:
1592: GEN
1593: arch_pow(GEN x, GEN n) {
1594: switch (typ(x)) {
1595: case t_POLMOD: return powgi(x,n);
1.2 ! noro 1596: case t_COL: return vecpow(x,n);
1.1 noro 1597: case t_MAT: return famat_pow(x,n);
1598: default: return gmul(n,x);
1599: }
1600: }
1601:
1602: /* output the ideal product ix.iy (don't reduce) */
1603: GEN
1604: idealmul(GEN nf, GEN x, GEN y)
1605: {
1.2 ! noro 1606: gpmem_t av;
! 1607: long tx,ty,f;
1.1 noro 1608: GEN res,ax,ay,p1;
1609:
1610: tx = idealtyp(&x,&ax);
1611: ty = idealtyp(&y,&ay);
1612: if (tx>ty) {
1613: res=ax; ax=ay; ay=res;
1614: res=x; x=y; y=res;
1615: f=tx; tx=ty; ty=f;
1616: }
1617: f = (ax||ay); res = f? cgetg(3,t_VEC): NULL; /* product is an idele */
1618: nf=checknf(nf); av=avma;
1619: switch(tx)
1620: {
1621: case id_PRINCIPAL:
1622: switch(ty)
1623: {
1624: case id_PRINCIPAL:
1625: p1 = idealhermite_aux(nf, element_mul(nf,x,y));
1626: break;
1627: case id_PRIME:
1.2 ! noro 1628: {
! 1629: GEN mx = eltmul_get_table(nf, x);
! 1630: GEN mpi= eltmul_get_table(nf, (GEN)y[2]);
! 1631: p1 = concatsp(gmul(mx,(GEN)y[1]), gmul(mx,mpi));
! 1632: p1 = idealmat_to_hnf(nf, p1);
1.1 noro 1633: break;
1.2 ! noro 1634: }
1.1 noro 1635: default: /* id_MAT */
1.2 ! noro 1636: p1 = idealmat_to_hnf(nf, gmul(eltmul_get_table(nf,x), y));
1.1 noro 1637: }break;
1638:
1639: case id_PRIME:
1640: p1 = (ty==id_PRIME)? prime_to_ideal_aux(nf,y)
1641: : idealmat_to_hnf(nf,y);
1642: p1 = idealmulprime(nf,p1,x); break;
1643:
1644: default: /* id_MAT */
1645: p1 = idealmat_mul(nf,x,y);
1646: }
1647: p1 = gerepileupto(av,p1);
1648: if (!f) return p1;
1649:
1650: if (ax && ay)
1651: ax = arch_mul(ax, ay);
1652: else
1653: ax = gcopy(ax? ax: ay);
1654: res[1]=(long)p1; res[2]=(long)ax; return res;
1655: }
1656:
1657: /* norm of an ideal */
1658: GEN
1659: idealnorm(GEN nf, GEN x)
1660: {
1.2 ! noro 1661: gpmem_t av = avma,tetpil;
1.1 noro 1662: GEN y;
1663:
1664: nf = checknf(nf);
1665: switch(idealtyp(&x,&y))
1666: {
1667: case id_PRIME:
1668: return powgi((GEN)x[1],(GEN)x[4]);
1669: case id_PRINCIPAL:
1670: x = gnorm(basistoalg(nf,x)); break;
1671: default:
1672: if (lg(x) != lgef(nf[1])-2) x = idealhermite_aux(nf,x);
1673: x = dethnf(x);
1674: }
1675: tetpil=avma; return gerepile(av,tetpil,gabs(x,0));
1676: }
1677:
1678: /* inverse */
1679:
1680: /* rewritten from original code by P.M & M.H.
1681: *
1682: * I^(-1) = { x \in K, Tr(x D^(-1) I) \in Z }, D different of K/Q
1683: *
1684: * nf[5][6] = d_K * D^(-1) is integral = d_K T^(-1), T = (Tr(wi wj))
1685: * nf[5][7] = same in 2-elt form */
1686: static GEN
1687: hnfideal_inv(GEN nf, GEN I)
1688: {
1.2 ! noro 1689: GEN J, dI, IZ,dual;
1.1 noro 1690:
1.2 ! noro 1691: I = Q_remove_denom(I, &dI);
1.1 noro 1692: if (lg(I)==1) err(talker, "cannot invert zero ideal");
1693: IZ = gcoeff(I,1,1); /* I \cap Z */
1694: if (!signe(IZ)) err(talker, "cannot invert zero ideal");
1695: J = idealmulh(nf,I, gmael(nf,5,7));
1696: /* I in HNF, hence easily inverted; multiply by IZ to get integer coeffs
1697: * d_K cancels while solving the linear equations. */
1698: dual = gtrans( gauss_triangle_i(J, gmael(nf,5,6), IZ) );
1699: dual = hnfmodid(dual, IZ);
1700: if (dI) IZ = gdiv(IZ,dI);
1701: return gdiv(dual,IZ);
1702: }
1703:
1704: /* return p * P^(-1) [integral] */
1705: GEN
1706: pidealprimeinv(GEN nf, GEN x)
1707: {
1708: GEN y=cgetg(6,t_VEC); y[1]=x[1]; y[2]=x[5];
1709: y[3]=y[5]=zero; y[4]=lsubsi(degpol(nf[1]), (GEN)x[4]);
1710: return prime_to_ideal_aux(nf,y);
1711: }
1712:
1713: /* Calcule le dual de mat_id pour la forme trace */
1714: GEN
1715: idealinv(GEN nf, GEN x)
1716: {
1717: GEN res,ax;
1.2 ! noro 1718: gpmem_t av=avma;
! 1719: long tx = idealtyp(&x,&ax);
1.1 noro 1720:
1721: res = ax? cgetg(3,t_VEC): NULL;
1722: nf=checknf(nf); av=avma;
1723: switch (tx)
1724: {
1725: case id_MAT:
1726: if (lg(x) != lg(x[1])) x = idealmat_to_hnf(nf,x);
1727: if (lg(x)-1 != degpol(nf[1])) err(consister,"idealinv");
1728: x = hnfideal_inv(nf,x); break;
1729: case id_PRINCIPAL: tx = typ(x);
1730: if (is_const_t(tx)) x = ginv(x);
1731: else
1732: {
1733: switch(tx)
1734: {
1735: case t_COL: x = gmul((GEN)nf[7],x); break;
1736: case t_POLMOD: x = (GEN)x[2]; break;
1737: }
1.2 ! noro 1738: x = QX_invmod(x,(GEN)nf[1]);
1.1 noro 1739: }
1740: x = idealhermite_aux(nf,x); break;
1741: case id_PRIME:
1742: x = gdiv(pidealprimeinv(nf,x), (GEN)x[1]);
1743: }
1744: x = gerepileupto(av,x); if (!ax) return x;
1745: res[1]=(long)x;
1746: res[2]=(long)arch_inv(ax); return res;
1747: }
1748:
1749: /* return x such that vp^n = x/d */
1750: static GEN
1751: idealpowprime_spec(GEN nf, GEN vp, GEN n, GEN *d)
1752: {
1753: GEN n1, x, r;
1754: long s = signe(n);
1755:
1756: if (s == 0) err(talker, "0th power in idealpowprime_spec");
1757: if (s < 0) n = negi(n);
1758: /* now n > 0 */
1759: x = dummycopy(vp);
1.2 ! noro 1760: if (is_pm1(n)) /* n = 1 special cased for efficiency */
1.1 noro 1761: {
1.2 ! noro 1762: if (s < 0)
! 1763: {
! 1764: x[2] = x[5];
! 1765: *d = (GEN)x[1];
! 1766: }
! 1767: else
! 1768: *d = NULL;
1.1 noro 1769: }
1770: else
1771: {
1.2 ! noro 1772: n1 = dvmdii(n, (GEN)x[3], &r);
! 1773: if (signe(r)) n1 = addis(n1,1); /* n1 = ceil(n/e) */
! 1774: x[1] = (long)powgi((GEN)x[1],n1);
! 1775: if (s < 0)
! 1776: {
! 1777: x[2] = ldiv(element_pow(nf,(GEN)x[5],n), powgi((GEN)vp[1],subii(n,n1)));
! 1778: *d = (GEN)x[1];
! 1779: }
! 1780: else
! 1781: {
! 1782: x[2] = (long)element_pow(nf,(GEN)x[2],n);
! 1783: *d = NULL;
! 1784: }
1.1 noro 1785: }
1786: return x;
1787: }
1788:
1789: static GEN
1790: idealpowprime(GEN nf, GEN vp, GEN n)
1791: {
1792: GEN x, d;
1793: long s = signe(n);
1794:
1795: nf = checknf(nf);
1796: if (s == 0) return idmat(degpol(nf[1]));
1797: x = idealpowprime_spec(nf, vp, n, &d);
1798: x = prime_to_ideal_aux(nf,x);
1799: if (d) x = gdiv(x, d);
1800: return x;
1801: }
1802:
1.2 ! noro 1803: /* x * vp^n. Assume x in HNF (possibly non-integral) */
1.1 noro 1804: GEN
1805: idealmulpowprime(GEN nf, GEN x, GEN vp, GEN n)
1806: {
1.2 ! noro 1807: GEN cx,y,dx;
1.1 noro 1808:
1809: if (!signe(n)) return x;
1810: nf = checknf(nf);
1.2 ! noro 1811:
! 1812: /* inert, special cased for efficiency */
! 1813: if (itos((GEN)vp[4]) == degpol(nf[1]))
! 1814: return gmul(x, powgi((GEN)vp[1], n));
! 1815:
! 1816: y = idealpowprime_spec(nf, vp, n, &dx);
! 1817: x = Q_primitive_part(x, &cx);
! 1818: if (cx && dx)
! 1819: {
! 1820: cx = gdiv(cx, dx);
! 1821: if (typ(cx) != t_FRAC) dx = NULL;
! 1822: else { dx = (GEN)cx[2]; cx = (GEN)cx[1]; }
! 1823: if (is_pm1(cx)) cx = NULL;
! 1824: }
! 1825: x = idealmulspec(nf,x,y);
! 1826: if (cx) x = gmul(x,cx);
! 1827: if (dx) x = gdiv(x,dx);
! 1828: return x;
! 1829: }
! 1830: GEN
! 1831: idealdivpowprime(GEN nf, GEN x, GEN vp, GEN n)
! 1832: {
! 1833: return idealmulpowprime(nf,x,vp, negi(n));
1.1 noro 1834: }
1835:
1836: /* raise the ideal x to the power n (in Z) */
1837: GEN
1838: idealpow(GEN nf, GEN x, GEN n)
1839: {
1.2 ! noro 1840: gpmem_t av;
! 1841: long tx,N,s;
1.1 noro 1842: GEN res,ax,m,cx,n1,a,alpha;
1843:
1844: if (typ(n) != t_INT) err(talker,"non-integral exponent in idealpow");
1845: tx = idealtyp(&x,&ax);
1846: res = ax? cgetg(3,t_VEC): NULL;
1847: nf = checknf(nf);
1848: av=avma; N=degpol(nf[1]); s=signe(n);
1849: if (!s) x = idmat(N);
1850: else
1851: switch(tx)
1852: {
1853: case id_PRINCIPAL: tx = typ(x);
1854: if (!is_const_t(tx))
1855: switch(tx)
1856: {
1857: case t_COL: x = gmul((GEN)nf[7],x);
1858: case t_POL: x = gmodulcp(x,(GEN)nf[1]);
1859: }
1860: x = powgi(x,n);
1861: x = idealhermite_aux(nf,x); break;
1862: case id_PRIME:
1863: x = idealpowprime(nf,x,n); break;
1864: default:
1865: n1 = (s<0)? negi(n): n;
1866:
1.2 ! noro 1867: x = Q_primitive_part(x, &cx);
1.1 noro 1868: a=ideal_two_elt(nf,x); alpha=(GEN)a[2]; a=(GEN)a[1];
1869: alpha = element_pow(nf,alpha,n1);
1.2 ! noro 1870: m = eltmul_get_table(nf, alpha);
! 1871: x = hnfmodid(m, powgi(a,n1));
1.1 noro 1872: if (s<0) x = hnfideal_inv(nf,x);
1873: if (cx) x = gmul(x, powgi(cx,n));
1874: }
1875: x = gerepileupto(av, x);
1876: if (!ax) return x;
1877: ax = arch_pow(ax, n);
1878: res[1]=(long)x;
1879: res[2]=(long)ax;
1880: return res;
1881: }
1882:
1883: /* Return ideal^e in number field nf. e is a C integer. */
1884: GEN
1885: idealpows(GEN nf, GEN ideal, long e)
1886: {
1.2 ! noro 1887: long court[] = {evaltyp(t_INT) | _evallg(3),0,0};
1.1 noro 1888: affsi(e,court); return idealpow(nf,ideal,court);
1889: }
1890:
1.2 ! noro 1891: static GEN
! 1892: _idealmulred(GEN nf, GEN x, GEN y, long prec)
! 1893: {
! 1894: return ideallllred(nf,idealmul(nf,x,y), NULL, prec);
! 1895: }
! 1896:
! 1897: typedef struct {
! 1898: GEN nf;
! 1899: long prec;
! 1900: } idealred_muldata;
! 1901:
! 1902: static GEN
! 1903: _mul(void *data, GEN x, GEN y)
! 1904: {
! 1905: idealred_muldata *D = (idealred_muldata *)data;
! 1906: return _idealmulred(D->nf,x,y,D->prec);
! 1907: }
! 1908: static GEN
! 1909: _sqr(void *data, GEN x)
1.1 noro 1910: {
1.2 ! noro 1911: return _mul(data,x,x);
1.1 noro 1912: }
1913:
1914: /* compute x^n (x ideal, n integer), reducing along the way */
1915: GEN
1916: idealpowred(GEN nf, GEN x, GEN n, long prec)
1917: {
1.2 ! noro 1918: gpmem_t av=avma;
! 1919: long s = signe(n);
! 1920: idealred_muldata D;
! 1921: GEN y;
1.1 noro 1922:
1923: if (typ(n) != t_INT) err(talker,"non-integral exponent in idealpowred");
1.2 ! noro 1924: if (s == 0) return idealpow(nf,x,n);
! 1925: D.nf = nf;
! 1926: D.prec= prec;
! 1927: y = leftright_pow(x, n, (void*)&D, &_sqr, &_mul);
! 1928:
1.1 noro 1929: if (s < 0) y = idealinv(nf,y);
1.2 ! noro 1930: if (s < 0 || is_pm1(n)) y = ideallllred(nf,y,NULL,prec);
1.1 noro 1931: return gerepileupto(av,y);
1932: }
1933:
1934: GEN
1935: idealmulred(GEN nf, GEN x, GEN y, long prec)
1936: {
1.2 ! noro 1937: gpmem_t av = avma;
! 1938: return gerepileupto(av, _idealmulred(nf,x,y,prec));
1.1 noro 1939: }
1940:
1941: long
1942: isideal(GEN nf,GEN x)
1943: {
1.2 ! noro 1944: gpmem_t av;
! 1945: long N,i,j,tx=typ(x),lx;
1.1 noro 1946:
1947: nf=checknf(nf); lx=lg(x);
1948: if (tx==t_VEC && lx==3) { x=(GEN)x[1]; tx=typ(x); lx=lg(x); }
1949: if (is_scalar_t(tx))
1950: return (tx==t_INT || tx==t_FRAC || tx==t_FRACN || tx==t_POL ||
1951: (tx==t_POLMOD && gegal((GEN)nf[1],(GEN)x[1])));
1952: if (typ(x)==t_VEC) return (lx==6);
1953: if (typ(x)!=t_MAT) return 0;
1954: if (lx == 1) return 1;
1.2 ! noro 1955: N = degpol(nf[1]); if (lg(x[1])-1 != N) return 0;
1.1 noro 1956:
1.2 ! noro 1957: av = avma;
! 1958: if (lx-1 != N) x = idealmat_to_hnf(nf,x);
! 1959: x = primpart(x);
! 1960:
! 1961: for (i=1; i<=N; i++)
! 1962: for (j=2; j<=N; j++)
! 1963: if (! hnf_invimage(x, element_mulid(nf,(GEN)x[i],j)))
! 1964: {
! 1965: avma = av; return 0;
! 1966: }
1.1 noro 1967: avma=av; return 1;
1968: }
1969:
1970: GEN
1971: idealdiv(GEN nf, GEN x, GEN y)
1972: {
1.2 ! noro 1973: gpmem_t av=avma,tetpil;
1.1 noro 1974: GEN z=idealinv(nf,y);
1975:
1976: tetpil=avma; return gerepile(av,tetpil,idealmul(nf,x,z));
1977: }
1978:
1979: /* This routine computes the quotient x/y of two ideals in the number field nf.
1980: * It assumes that the quotient is an integral ideal. The idea is to find an
1981: * ideal z dividing y such that gcd(Nx/Nz, Nz) = 1. Then
1982: *
1983: * x + (Nx/Nz) x
1984: * ----------- = ---
1985: * y + (Ny/Nz) y
1986: *
1987: * Proof: we can assume x and y are integral. Let p be any prime ideal
1988: *
1989: * If p | Nz, then it divides neither Nx/Nz nor Ny/Nz (since Nx/Nz is the
1990: * product of the integers N(x/y) and N(y/z)). Both the numerator and the
1991: * denominator on the left will be coprime to p. So will x/y, since x/y is
1992: * assumed integral and its norm N(x/y) is coprime to p.
1993: *
1994: * If instead p does not divide Nz, then v_p (Nx/Nz) = v_p (Nx) >= v_p(x).
1995: * Hence v_p (x + Nx/Nz) = v_p(x). Likewise for the denominators. QED.
1996: *
1997: * Peter Montgomery. July, 1994. */
1998: GEN
1999: idealdivexact(GEN nf, GEN x0, GEN y0)
2000: {
1.2 ! noro 2001: gpmem_t av = avma;
1.1 noro 2002: GEN x,y,Nx,Ny,Nz, cy = content(y0);
2003:
2004: nf = checknf(nf);
2005: if (gcmp0(cy)) err(talker, "cannot invert zero ideal");
2006:
2007: x = gdiv(x0,cy); Nx = idealnorm(nf,x);
2008: if (gcmp0(Nx)) { avma = av; return gcopy(x0); } /* numerator is zero */
2009:
2010: y = gdiv(y0,cy); Ny = idealnorm(nf,y);
2011: if (!gcmp1(denom(x)) || !divise(Nx,Ny))
2012: err(talker, "quotient not integral in idealdivexact");
2013: /* Find a norm Nz | Ny such that gcd(Nx/Nz, Nz) = 1 */
2014: for (Nz = Ny;;)
2015: {
2016: GEN p1 = mppgcd(Nz, divii(Nx,Nz));
2017: if (is_pm1(p1)) break;
2018: Nz = divii(Nz,p1);
2019: }
2020: /* Replace x/y by x+(Nx/Nz) / y+(Ny/Nz) */
2021: x = idealhermite_aux(nf, x);
2022: x = hnfmodid(x, divii(Nx,Nz));
2023: /* y reduced to unit ideal ? */
2024: if (Nz == Ny) return gerepileupto(av, x);
2025:
2026: y = idealhermite_aux(nf, y);
2027: y = hnfmodid(y, divii(Ny,Nz));
2028: y = hnfideal_inv(nf,y);
2029: return gerepileupto(av, idealmat_mul(nf,x,y));
2030: }
2031:
2032: GEN
2033: idealintersect(GEN nf, GEN x, GEN y)
2034: {
1.2 ! noro 2035: gpmem_t av=avma;
! 2036: long lz,i,N;
1.1 noro 2037: GEN z,dx,dy;
2038:
2039: nf=checknf(nf); N=degpol(nf[1]);
2040: if (idealtyp(&x,&z)!=t_MAT || lg(x)!=N+1) x=idealhermite_aux(nf,x);
2041: if (idealtyp(&y,&z)!=t_MAT || lg(y)!=N+1) y=idealhermite_aux(nf,y);
2042: if (lg(x) == 1 || lg(y) == 1) { avma = av; return cgetg(1, t_MAT); }
2043: dx=denom(x); if (!gcmp1(dx)) y = gmul(y,dx);
2044: dy=denom(y); if (!gcmp1(dy)) { x = gmul(x,dy); dx = mulii(dx,dy); }
2045: z = kerint(concatsp(x,y)); lz=lg(z);
2046: for (i=1; i<lz; i++) setlg(z[i], N+1);
2047: z = gmul(x,z);
2048: z = hnfmodid(z, glcm(gcoeff(x,1,1), gcoeff(y,1,1)));
2049: if (!gcmp1(dx)) z = gdiv(z,dx);
2050: return gerepileupto(av,z);
2051: }
2052:
1.2 ! noro 2053: /*******************************************************************/
! 2054: /* */
! 2055: /* T2-IDEAL REDUCTION */
! 2056: /* */
! 2057: /*******************************************************************/
! 2058:
! 2059: GEN
! 2060: computeGtwist(GEN nf, GEN vdir)
1.1 noro 2061: {
1.2 ! noro 2062: long i, j, k, l, lG, v, r1, r2;
! 2063: GEN p1, G = gmael(nf,5,2);
1.1 noro 2064:
1.2 ! noro 2065: if (!vdir) return G;
! 2066: l = lg(vdir); lG = lg(G);
! 2067: p1 = dummycopy(G);
! 2068: nf_get_sign(nf, &r1, &r2);
! 2069: for (i=1; i<l; i++)
! 2070: {
! 2071: v = vdir[i]; if (!v) continue;
! 2072: if (i <= r1) {
! 2073: for (j=1; j<lG; j++) coeff(p1,i,j) = lmul2n(gcoeff(p1,i,j), v);
! 2074: } else {
! 2075: k = (i<<1) - r1;
! 2076: for (j=1; j<lG; j++)
! 2077: {
! 2078: coeff(p1,k-1,j) = lmul2n(gcoeff(p1,k-1,j), v);
! 2079: coeff(p1,k ,j) = lmul2n(gcoeff(p1,k ,j), v);
! 2080: }
! 2081: }
1.1 noro 2082: }
1.2 ! noro 2083: return p1;
1.1 noro 2084: }
2085:
2086: static GEN
2087: chk_vdir(GEN nf, GEN vdir)
2088: {
1.2 ! noro 2089: long l, i, t;
! 2090: GEN v;
1.1 noro 2091: if (!vdir || gcmp0(vdir)) return NULL;
1.2 ! noro 2092: l = lg(vdir);
! 2093: if (l != lg(nf[6])) err(talker, "incorrect vector length in idealred");
! 2094: t = typ(vdir);
! 2095: if (t == t_VECSMALL) return vdir;
! 2096: if (t != t_VEC) err(talker, "not a vector in idealred");
! 2097: v = cgetg(l, t_VECSMALL);
! 2098: for (i=1; i<l; i++) v[i] = itos(gceil((GEN)vdir[i]));
! 2099: return v;
1.1 noro 2100: }
2101:
2102: /* assume I in NxN matrix form (not necessarily HNF) */
2103: static GEN
1.2 ! noro 2104: idealred_elt_i(GEN *ptnf, GEN I, GEN vdir, long *ptprec)
1.1 noro 2105: {
1.2 ! noro 2106: GEN G, u, y, nf = *ptnf;
1.1 noro 2107: long i, e, prec = *ptprec;
2108:
1.2 ! noro 2109: e = gexpo(I)>>TWOPOTBITS_IN_LONG;
! 2110: if (e < 0) e = 0;
1.1 noro 2111: for (i=1; ; i++)
2112: {
1.2 ! noro 2113: G = computeGtwist(nf,vdir);
! 2114: y = gmul(G, I);
! 2115: u = lllintern(y, 100, 1, prec);
1.1 noro 2116: if (u) break;
2117:
2118: if (i == MAXITERPOL) err(accurer,"ideallllred");
2119: prec = (prec<<1)-2;
2120: if (DEBUGLEVEL) err(warnprec,"ideallllred",prec);
1.2 ! noro 2121: nf = nfnewprec(nf, e + prec);
1.1 noro 2122: }
2123: *ptprec = prec;
2124: *ptnf = nf;
2125: return gmul(I, (GEN)u[1]); /* small elt in I */
2126: }
2127:
2128: GEN
1.2 ! noro 2129: ideallllred_elt(GEN nf, GEN I, GEN vdir)
1.1 noro 2130: {
2131: long prec = DEFAULTPREC;
1.2 ! noro 2132: return idealred_elt_i(&nf, I, vdir, &prec);
! 2133: }
! 2134: GEN
! 2135: idealred_elt(GEN nf, GEN I)
! 2136: {
! 2137: return ideallllred_elt(nf, I, NULL);
1.1 noro 2138: }
2139:
2140: GEN
2141: ideallllred(GEN nf, GEN I, GEN vdir, long prec)
2142: {
1.2 ! noro 2143: gpmem_t av = avma;
1.1 noro 2144: long N,i,nfprec;
2145: GEN J,I0,Ired,res,aI,y,x,Nx,b,c1,c,pol;
2146:
2147: nf = checknf(nf); nfprec = nfgetprec(nf);
2148: if (prec <= 0) prec = nfprec;
2149: pol = (GEN)nf[1]; N = degpol(pol);
2150: Nx = x = c = c1 = NULL;
2151: if (idealtyp(&I,&aI) == id_PRINCIPAL)
2152: {
2153: if (gcmp0(I)) { y=gun; I=cgetg(1,t_MAT); } else { y=I; I=idmat(N); }
2154: goto END;
2155: }
2156:
2157: if (DEBUGLEVEL>5) msgtimer("entering idealllred");
2158: I0 = I;
2159: if (typ(I) != id_MAT || lg(I) != N+1) I = idealhermite_aux(nf,I);
1.2 ! noro 2160: I = primitive_part(I, &c1);
1.1 noro 2161: if (2 * expi(gcoeff(I,1,1)) >= bit_accuracy(nfprec))
1.2 ! noro 2162: Ired = lllint_ip(I,4);
1.1 noro 2163: else
2164: Ired = I;
1.2 ! noro 2165: y = idealred_elt_i(&nf, Ired, chk_vdir(nf,vdir), &prec);
1.1 noro 2166:
2167: if (isnfscalar(y))
2168: { /* already reduced */
2169: if (!aI) I = gcopy(I);
2170: y = NULL; goto END;
2171: }
2172: if (DEBUGLEVEL>5) msgtimer("LLL reduction");
2173:
2174: x = gmul((GEN)nf[7], y); Nx = subres(pol,x);
1.2 ! noro 2175: b = gmul(Nx, QX_invmod(x,pol));
! 2176: b = algtobasis_i(nf,b);
1.1 noro 2177: J = cgetg(N+1,t_MAT); /* = I Nx / x integral */
2178: for (i=1; i<=N; i++)
2179: J[i] = (long)element_muli(nf,b,(GEN)Ired[i]);
1.2 ! noro 2180: J = Q_primitive_part(J, &c);
1.1 noro 2181: /* c = content (I Nx / x) = Nx / den(I/x) --> d = den(I/x) = Nx / c
2182: * J = (d I / x); I[1,1] = I \cap Z --> d I[1,1] belongs to J and Z */
2183: if (isnfscalar((GEN)I[1]))
1.2 ! noro 2184: b = mulii(gcoeff(I,1,1), c? diviiexact(Nx, c): Nx);
1.1 noro 2185: else
2186: b = detint(J);
2187: I = hnfmodid(J,b);
2188: if (DEBUGLEVEL>5) msgtimer("new ideal");
2189:
2190: END:
2191: if (!aI) return gerepileupto(av, I);
2192:
2193: switch(typ(aI))
2194: {
1.2 ! noro 2195: case t_POLMOD: /* compute y, I0 = J y */
1.1 noro 2196: if (!Nx) y = c1;
2197: else
2198: {
1.2 ! noro 2199: c = mul_content(c,c1);
! 2200: y = c? gmul(x, gdiv(c,Nx)): gdiv(x, Nx);
! 2201: }
! 2202: break;
! 2203:
! 2204: case t_MAT: /* compute y, I0 = J y */
! 2205: if (!Nx) y = c1;
! 2206: else
! 2207: {
! 2208: c = mul_content(c,c1);
! 2209: y = c? gmul(y, gdiv(c,Nx)): gdiv(y, Nx);
1.1 noro 2210: }
2211: break;
2212:
2213: case t_COL:
2214: if (y) y = vecinv(gmul(gmael(nf,5,1), y));
2215: break;
2216:
2217: default:
2218: if (y) y = gneg_i(get_arch(nf,y,prec));
2219: break;
2220: }
2221: if (y) aI = arch_mul(aI,y);
2222: res = cgetg(3,t_VEC);
2223: res[1] = (long)I;
2224: res[2] = (long)aI; return gerepilecopy(av, res);
2225: }
2226:
2227: GEN
2228: minideal(GEN nf, GEN x, GEN vdir, long prec)
2229: {
1.2 ! noro 2230: gpmem_t av = avma;
! 2231: long N, tx;
! 2232: GEN y;
1.1 noro 2233:
2234: nf = checknf(nf);
2235: vdir = chk_vdir(nf,vdir);
2236: N = degpol(nf[1]);
2237: tx = idealtyp(&x,&y);
2238: if (tx == id_PRINCIPAL) return gcopy(x);
2239: if (tx != id_MAT || lg(x) != N+1) x = idealhermite_aux(nf,x);
2240:
1.2 ! noro 2241: y = gmul(computeGtwist(nf,vdir), x);
! 2242: y = gmul(x, (GEN)lll(y,prec)[1]);
1.1 noro 2243: return gerepileupto(av, principalidele(nf,y,prec));
2244: }
1.2 ! noro 2245:
! 2246: /*******************************************************************/
! 2247: /* */
! 2248: /* APPROXIMATION THEOREM */
! 2249: /* */
! 2250: /*******************************************************************/
! 2251:
! 2252: /* write x = x1 x2, x2 maximal s.t. (x2,f) = 1, return x2 */
1.1 noro 2253: static GEN
1.2 ! noro 2254: coprime_part(GEN x, GEN f)
1.1 noro 2255: {
1.2 ! noro 2256: for (;;)
! 2257: {
! 2258: f = mppgcd(x, f); if (is_pm1(f)) break;
! 2259: x = diviiexact(x, f);
! 2260: }
! 2261: return x;
1.1 noro 2262: }
2263:
1.2 ! noro 2264: /* x t_INT, f ideal. Write x = x1 x2, sqf(x1) | f, (x2,f) = 1. Return x2 */
! 2265: static GEN
! 2266: nf_coprime_part(GEN nf, GEN x, GEN *listpr)
1.1 noro 2267: {
1.2 ! noro 2268: long v, j, lp = lg(listpr), N = degpol(nf[1]);
! 2269: GEN x1, x2, ex, p, pr;
1.1 noro 2270:
1.2 ! noro 2271: #if 0 /*1) via many gcds. Expensive ! */
! 2272: GEN f = idealprodprime(nf, (GEN)listpr);
! 2273: f = hnfmodid(f, x); /* first gcd is less expensive since x in Z */
! 2274: x = gscalmat(x, N);
! 2275: for (;;)
1.1 noro 2276: {
1.2 ! noro 2277: if (gcmp1(gcoeff(f,1,1))) break;
! 2278: x = idealdivexact(nf, x, f);
! 2279: f = hnfmodid(concatsp(f,x), gcoeff(x,1,1)); /* gcd(f,x) */
1.1 noro 2280: }
1.2 ! noro 2281: x2 = x;
! 2282: #else /*2) from prime decomposition */
! 2283: x1 = NULL;
! 2284: for (j=1; j<lp; j++)
1.1 noro 2285: {
1.2 ! noro 2286: pr = listpr[j]; p = (GEN)pr[1];
! 2287: v = ggval(x, p); if (!v) continue;
! 2288:
! 2289: ex = mulsi(v, (GEN)pr[3]); /* = v_pr(x) > 0 */
! 2290: x1 = x1? idealmulpowprime(nf, x1, pr, ex)
! 2291: : idealpow(nf, pr, ex);
1.1 noro 2292: }
1.2 ! noro 2293: x = gscalmat(x, N);
! 2294: x2 = x1? idealdivexact(nf, x, x1): x;
! 2295: #endif
! 2296: return x2;
! 2297: }
! 2298:
! 2299: /* assume L is a list of prime ideals. Return the product */
! 2300: GEN
! 2301: idealprodprime(GEN nf, GEN L)
! 2302: {
! 2303: long l = lg(L), i;
! 2304: GEN z;
! 2305:
! 2306: if (l == 1) return idmat(degpol(nf[1]));
! 2307: z = prime_to_ideal(nf, (GEN)L[1]);
! 2308: for (i=2; i<l; i++) z = idealmulprime(nf,z, (GEN)L[i]);
! 2309: return z;
! 2310: }
! 2311:
! 2312: /* assume L is a list of prime ideals. Return prod L[i]^e[i] */
! 2313: GEN
! 2314: factorbackprime(GEN nf, GEN L, GEN e)
! 2315: {
! 2316: long l = lg(L), i;
! 2317: GEN z;
! 2318:
! 2319: if (l == 1) return idmat(degpol(nf[1]));
! 2320: z = idealpow(nf, (GEN)L[1], (GEN)e[1]);
! 2321: for (i=2; i<l; i++)
! 2322: if (signe(e[i])) z = idealmulpowprime(nf,z, (GEN)L[i],(GEN)e[i]);
! 2323: return z;
! 2324: }
! 2325:
! 2326: /* compute anti-uniformizer for pr, coprime to f outside of pr, integral
! 2327: * outside of p below pr.
! 2328: * sqf = product or primes dividing f, NULL if f a prime power*/
! 2329: GEN
! 2330: anti_unif_mod_f(GEN nf, GEN pr, GEN sqf)
! 2331: {
! 2332: GEN U, V, UV, uv, cx, t, p = (GEN)pr[1];
! 2333: if (!sqf) return gdiv((GEN)pr[5], p);
1.1 noro 2334: else
2335: {
1.2 ! noro 2336: GEN d,d1,d2, sqfZ = gcoeff(sqf,1,1); /* = (sqf \cap Z) = (V \cap Z) * p */
! 2337: U = idealpow(nf,pr,gdeux);
! 2338: V = idealdivpowprime(nf,sqf,pr,gun);
! 2339: uv = addone_nored(U, V);
! 2340: UV = idealmul(nf,sqf,pr);
! 2341: t = (GEN)pr[2];
! 2342: t = makeprimetoideal(nf, UV, uv, t); /* cf unif_mod_f */
! 2343:
! 2344: t = element_inv(nf, t);
! 2345: t = primitive_part(t, &cx); /* p | denom(cx) */
! 2346: d = denom(cx);
! 2347: d2 = coprime_part(d, sqfZ);
! 2348: d1 = diviiexact(d, d2); /* p | d1 */
! 2349: cx = gmod(gmul(cx,d1), mulii(d1, gcoeff(V,1,1)));
! 2350: t = gdiv(colreducemodHNF(gmul(cx,t), gmul(d1,V), NULL), d1);
! 2351: return t; /* v_pr[i](t) = -1, v_pr[j](t) = 0 if i != j */
1.1 noro 2352: }
1.2 ! noro 2353: }
1.1 noro 2354:
1.2 ! noro 2355: /* compute integral uniformizer for pr, coprime to f outside of pr.
! 2356: * sqf = product or primes dividing f, NULL if f a prime power*/
! 2357: GEN
! 2358: unif_mod_f(GEN nf, GEN pr, GEN sqf)
! 2359: {
! 2360: GEN U, V, UV, uv, t = (GEN)pr[2];
! 2361: if (!sqf) return t;
! 2362: else
1.1 noro 2363: {
1.2 ! noro 2364: U = idealpow(nf,pr,gdeux);
! 2365: V = idealdivpowprime(nf,sqf,pr,gun);
! 2366: uv = addone_nored(U, V);
! 2367: UV = idealmulprime(nf,sqf,pr);
! 2368: return makeprimetoideal(nf, UV, uv, t);
1.1 noro 2369: }
1.2 ! noro 2370: }
! 2371:
! 2372: static GEN
! 2373: appr_reduce(GEN s, GEN y)
! 2374: {
! 2375: long i, N = lg(y)-1;
! 2376: GEN d, u, z = cgetg(N+2,t_MAT);
! 2377:
! 2378: s = gmod(s, gcoeff(y,1,1));
! 2379: y = lllint_ip(y,100);
! 2380: for (i=1; i<=N; i++) z[i] = y[i];
! 2381: z[N+1] = (long)s;
! 2382: u = (GEN)ker(z)[1];
! 2383: u = Q_remove_denom(u, NULL);
! 2384: d = (GEN)u[N+1]; setlg(u,N+1);
! 2385: for (i=1; i<=N; i++) u[i] = (long)diviiround((GEN)u[i],d);
! 2386: return gadd(s, gmul(y,u));
! 2387: }
! 2388:
! 2389: /* L0 in K^*.
! 2390: * If ptd1 == NULL, assume (L0,f) = 1
! 2391: * return L integral, L0 = L mod f
! 2392: *
! 2393: * Otherwise, assume v_pr(L0) <= 0 for all pr | f and set *ptd1 = d1
! 2394: * return L integral, L0 = L/d1 mod f, and such that
! 2395: * if (L*I,f) = 1 for some integral I, then d1 | L*I */
! 2396: GEN
! 2397: make_integral(GEN nf, GEN L0, GEN f, GEN *listpr, GEN *ptd1)
! 2398: {
! 2399: GEN fZ, t, L, D2, d1, d2, d;
! 2400:
! 2401: if (ptd1) *ptd1 = NULL;
! 2402: L = Q_remove_denom(L0, &d);
! 2403: if (!d) return L0;
! 2404:
! 2405: /* L0 = L / d, L integral */
! 2406: fZ = gcoeff(f,1,1);
! 2407: /* Kill denom part coprime to fZ */
! 2408: d2 = coprime_part(d, fZ);
! 2409: t = mpinvmod(d2, fZ); if (!is_pm1(t)) L = gmul(L,t);
! 2410: if (egalii(d, d2)) return L;
! 2411:
! 2412: d1 = diviiexact(d, d2);
! 2413: /* L0 = (L / d1) mod f. d1 not coprime to f
! 2414: * write (d1) = D1 D2, D2 minimal, (D2,f) = 1. */
! 2415: D2 = nf_coprime_part(nf, d1, listpr);
! 2416: t = idealaddtoone_i(nf, D2, f); /* in D2, 1 mod f */
! 2417: L = element_muli(nf,t,L);
! 2418:
! 2419: /* if (L0, f) = 1, then L in D1 ==> in D1 D2 = (d1) */
! 2420: if (!ptd1) return Q_div_to_int(L, d1); /* exact division */
! 2421:
! 2422: *ptd1 = d1; return L;
! 2423: }
! 2424:
! 2425: void
! 2426: check_listpr(GEN x)
! 2427: {
! 2428: long l = lg(x), i;
! 2429: for (i=1; i<l; i++) checkprimeid((GEN)x[i]);
! 2430: }
! 2431:
! 2432: /* Given a prime ideal factorization with possibly zero or negative
! 2433: * exponents, gives b such that v_p(b) = v_p(x) for all prime ideals pr | x
! 2434: * and v_pr(b)> = 0 for all other pr.
! 2435: * For optimal performance, all [anti-]uniformizers should be precomputed,
! 2436: * but no support for this yet.
! 2437: * No garbage collecting */
! 2438: GEN
! 2439: idealapprfact_i(GEN nf, GEN x)
! 2440: {
! 2441: gpmem_t av = avma;
! 2442: GEN tau, pi, z, d, sqf, sqfsafe, list, e, e2;
! 2443: long s, flag, i, r, N;
! 2444:
! 2445: nf = checknf(nf);
! 2446: if (typ(x) != t_MAT || lg(x) != 3)
! 2447: err(talker,"not a factorization in idealapprfact");
! 2448: check_listpr((GEN)x[1]);
! 2449: if (DEBUGLEVEL>4) {
! 2450: fprintferr(" entering idealapprfact():\n");
! 2451: fprintferr(" x = %Z\n", x);
! 2452: }
! 2453: list= (GEN)x[1];
! 2454: e = (GEN)x[2]; r = lg(e); N = degpol(nf[1]);
! 2455: if (r==1) return gscalcol_i(gun, N);
! 2456:
! 2457: sqf = (r == 2)? NULL: idealprodprime(nf, list);
! 2458: sqfsafe = NULL;
! 2459: z = gun; flag = 0;
1.1 noro 2460: for (i=1; i<r; i++)
2461: {
1.2 ! noro 2462: s = signe(e[i]);
! 2463: if (s < 0)
! 2464: {
! 2465: flag = 1;
! 2466: if (!sqfsafe) sqfsafe = sqf? sqf: prime_to_ideal(nf, (GEN)list[1]);
! 2467: tau = anti_unif_mod_f(nf, (GEN)list[i], sqf);
! 2468: tau = make_integral(nf, tau, sqfsafe, (GEN*)list, &d);
! 2469: tau = gdiv(tau, d);
! 2470: z = element_mul(nf, z, element_pow(nf, tau, negi((GEN)e[i])));
! 2471: }
! 2472: else if (s > 0)
! 2473: {
! 2474: pi = unif_mod_f(nf, (GEN)list[i], sqf);
! 2475: z = element_mul(nf, z, element_pow(nf, pi, (GEN)e[i]));
! 2476: }
! 2477: }
! 2478: if (z == gun) { avma = av; return gscalcol_i(gun, N); }
! 2479: e2 = cgetg(r, t_VEC);
! 2480: for (i=1; i<r; i++) e2[i] = laddis((GEN)e[i], 1);
! 2481: x = factorbackprime(nf, list,e2);
! 2482: if (flag) /* denominator */
! 2483: {
! 2484: z = Q_remove_denom(z, &d);
! 2485: x = Q_muli_to_int(x, d);
! 2486: }
! 2487: else d = NULL;
! 2488: z = appr_reduce(z, x);
! 2489: if (d) z = gdiv(z, d);
! 2490: return z;
! 2491: }
! 2492:
! 2493: GEN
! 2494: idealapprfact(GEN nf, GEN x) {
! 2495: gpmem_t av = avma;
! 2496: return gerepileupto(av, idealapprfact_i(nf, x));
! 2497: }
! 2498:
! 2499: GEN
! 2500: idealappr(GEN nf, GEN x) {
! 2501: gpmem_t av = avma;
! 2502: return gerepileupto(av, idealapprfact_i(nf, idealfactor(nf, x)));
! 2503: }
! 2504:
! 2505: GEN
! 2506: idealappr0(GEN nf, GEN x, long fl) {
! 2507: return fl? idealapprfact(nf, x): idealappr(nf, x);
1.1 noro 2508: }
2509:
2510: /* Given a prime ideal factorization x with possibly zero or negative exponents,
1.2 ! noro 2511: * and a vector w of elements of nf, gives a b such that
! 2512: * v_p(b-w_p)>=v_p(x) for all prime ideals p in the ideal factorization
1.1 noro 2513: * and v_p(b)>=0 for all other p, using the (standard) proof given in GTM 138.
1.2 ! noro 2514: * Certainly not the most efficient, but sure. */
1.1 noro 2515: GEN
1.2 ! noro 2516: idealchinese(GEN nf, GEN x, GEN w)
1.1 noro 2517: {
1.2 ! noro 2518: gpmem_t av = avma;
! 2519: long ty = typ(w), i,N,r;
! 2520: GEN L,e,z,t,y,v,s,den;
1.1 noro 2521:
2522: nf=checknf(nf); N=degpol(nf[1]);
1.2 ! noro 2523: if (typ(x) != t_MAT || lg(x) != 3)
1.1 noro 2524: err(talker,"not a prime ideal factorization in idealchinese");
1.2 ! noro 2525: L = (GEN)x[1]; r = lg(L);
! 2526: e = (GEN)x[2];
! 2527: if (!is_vec_t(ty) || lg(w) != r)
1.1 noro 2528: err(talker,"not a suitable vector of elements in idealchinese");
2529: if (r==1) return gscalcol_i(gun,N);
2530:
1.2 ! noro 2531: den = denom(w);
1.1 noro 2532: if (!gcmp1(den))
2533: {
1.2 ! noro 2534: GEN fa = famat_reduce(famat_mul(x, idealfactor(nf,den)));
! 2535: L = (GEN)fa[1]; r = lg(L);
! 2536: e = (GEN)fa[2];
1.1 noro 2537: }
2538: for (i=1; i<r; i++)
1.2 ! noro 2539: if (signe(e[i]) < 0) e[i] = zero;
! 2540:
! 2541: y = factorbackprime(nf, L, e);
! 2542: z = cgetg(r,t_VEC);
! 2543: for (i=1; i<r; i++)
! 2544: z[i] = (long)idealdivpowprime(nf,y, (GEN)L[i], (GEN)e[i]);
! 2545: v = idealaddmultoone(nf,z); s = NULL;
1.1 noro 2546: for (i=1; i<r; i++)
2547: {
1.2 ! noro 2548: t = element_mul(nf, (GEN)v[i], (GEN)w[i]);
! 2549: s = s? gadd(s,t): t;
1.1 noro 2550: }
1.2 ! noro 2551: return gerepileupto(av, appr_reduce(s,y));
! 2552: }
! 2553:
! 2554: static GEN
! 2555: mat_ideal_two_elt2(GEN nf, GEN x, GEN a)
! 2556: {
! 2557: GEN L, e, fact = idealfactor(nf,a);
! 2558: long i, v, r;
! 2559: L = (GEN)fact[1]; r = lg(L);
! 2560: e = (GEN)fact[2];
1.1 noro 2561: for (i=1; i<r; i++)
2562: {
1.2 ! noro 2563: v = idealval(nf,x,(GEN)L[i]);
! 2564: e[i] = lstoi(v);
1.1 noro 2565: }
1.2 ! noro 2566: return centermod(idealapprfact_i(nf,fact), gcoeff(x,1,1));
1.1 noro 2567: }
2568:
2569: /* Given an integral ideal x and a in x, gives a b such that
1.2 ! noro 2570: * x = aZ_K + bZ_K using the approximation theorem */
1.1 noro 2571: GEN
2572: ideal_two_elt2(GEN nf, GEN x, GEN a)
2573: {
1.2 ! noro 2574: gpmem_t av = avma;
! 2575: GEN cx, b;
1.1 noro 2576:
2577: nf = checknf(nf);
1.2 ! noro 2578: if (typ(a) != t_COL) a = algtobasis(nf, a);
1.1 noro 2579: x = idealhermite_aux(nf,x);
2580: if (gcmp0(x))
2581: {
2582: if (!gcmp0(a)) err(talker,"element not in ideal in ideal_two_elt2");
2583: avma=av; return gcopy(a);
2584: }
1.2 ! noro 2585: x = Q_primitive_part(x, &cx);
! 2586: if (cx) a = gdiv(a, cx);
! 2587: if (!hnf_invimage(x, a))
1.1 noro 2588: err(talker,"element does not belong to ideal in ideal_two_elt2");
2589:
1.2 ! noro 2590: b = mat_ideal_two_elt2(nf, x, a);
! 2591: b = cx? gmul(b,cx): gcopy(b);
! 2592: return gerepileupto(av, b);
1.1 noro 2593: }
2594:
1.2 ! noro 2595: /* Given 2 integral ideals x and y in nf, returns a beta in nf such that
! 2596: * beta * x is an integral ideal coprime to y */
1.1 noro 2597: GEN
2598: idealcoprime(GEN nf, GEN x, GEN y)
2599: {
1.2 ! noro 2600: gpmem_t av = avma;
! 2601: long i, r;
! 2602: GEN list, ep, fact = idealfactor(nf,y);
! 2603:
! 2604: list = (GEN)fact[1];
! 2605: ep = (GEN)fact[2]; r = lg(ep);
! 2606: for (i=1; i<r; i++) ep[i] = lstoi(-idealval(nf,x,(GEN)list[i]));
! 2607: return gerepileupto(av, idealapprfact_i(nf,fact));
! 2608: }
1.1 noro 2609:
1.2 ! noro 2610: /*******************************************************************/
! 2611: /* */
! 2612: /* LINEAR ALGEBRA OVER Z_K (HNF,SNF) */
! 2613: /* */
! 2614: /*******************************************************************/
1.1 noro 2615:
2616: /* returns the first index i<=n such that x=v[i] if it exits, 0 otherwise */
2617: long
2618: isinvector(GEN v, GEN x, long n)
2619: {
2620: long i;
2621:
2622: for (i=1; i<=n; i++)
2623: if (gegal((GEN)v[i],x)) return i;
2624: return 0;
2625: }
2626:
2627: GEN
2628: element_mulvec(GEN nf, GEN x, GEN v)
2629: {
1.2 ! noro 2630: long lv = lg(v), i;
1.1 noro 2631: GEN y = cgetg(lv,t_COL);
2632:
2633: if (typ(x) == t_COL)
2634: {
1.2 ! noro 2635: GEN mul = eltmul_get_table(nf,x);
! 2636: for (i=1; i<lv; i++) y[i] = lmul(mul,(GEN)v[i]);
1.1 noro 2637: }
2638: else
2639: { /* scalar */
1.2 ! noro 2640: for (i=1; i<lv; i++) y[i] = lmul(x, (GEN)v[i]);
1.1 noro 2641: }
2642: return y;
2643: }
2644:
2645: static GEN
2646: element_mulvecrow(GEN nf, GEN x, GEN m, long i, long lim)
2647: {
2648: long lv,j;
1.2 ! noro 2649: GEN y, mul = eltmul_get_table(nf,x);
1.1 noro 2650:
2651: lv=min(lg(m),lim+1); y=cgetg(lv,t_VEC);
1.2 ! noro 2652: for (j=1; j<lv; j++) y[j] = lmul(mul,gcoeff(m,i,j));
1.1 noro 2653: return y;
2654: }
2655:
2656: /* Given an element x and an ideal in matrix form (not necessarily HNF),
1.2 ! noro 2657: * gives an r such that x-r is in ideal and r is small. No checks */
1.1 noro 2658: GEN
2659: element_reduce(GEN nf, GEN x, GEN ideal)
2660: {
1.2 ! noro 2661: long tx = typ(x), N, i;
! 2662: gpmem_t av=avma;
! 2663: GEN u,d;
1.1 noro 2664:
2665: if (is_extscalar_t(tx))
1.2 ! noro 2666: x = algtobasis_i(checknf(nf), x);
1.1 noro 2667: N = lg(x);
2668: if (typ(ideal) != t_MAT || lg(ideal) != N) err(typeer,"element_reduce");
1.2 ! noro 2669: u = ker( concatsp(ideal,x) ); /* do NOT use deplin. Much, much slower -- KB */
! 2670: u = (GEN)u[1];
! 2671: d = (GEN)u[N]; setlg(u,N);
! 2672: for (i=1; i<N; i++) u[i] = (long)gdivround((GEN)u[i],d);
! 2673: return gerepileupto(av, gadd(x, gmul(ideal,u)));
1.1 noro 2674: }
2675:
2676: /* A torsion-free module M over Z_K will be given by a row vector [A,I] with
2677: * two components. I=[\a_1,...,\a_k] is a row vector of k fractional ideals
2678: * given in HNF. A is an nxk matrix (same k and n the rank of the module)
2679: * such that if A_j is the j-th column of A then M=\a_1A_1+...\a_kA_k. We say
2680: * that [A,I] is a pseudo-basis if k=n
2681: */
2682:
1.2 ! noro 2683: #define swap(x,y) { long _t=x; x=y; y=_t; }
! 2684:
! 2685: /* Given a torsion-free module x as above outputs a pseudo-basis for x in HNF */
1.1 noro 2686: GEN
2687: nfhermite(GEN nf, GEN x)
2688: {
1.2 ! noro 2689: long i, j, def, k, m;
! 2690: gpmem_t av0 = avma, av, lim;
1.1 noro 2691: GEN p1,p2,y,A,I,J;
2692:
1.2 ! noro 2693: nf = checknf(nf);
1.1 noro 2694: if (typ(x)!=t_VEC || lg(x)!=3) err(talker,"not a module in nfhermite");
1.2 ! noro 2695: A = (GEN)x[1];
! 2696: I = (GEN)x[2]; k = lg(A)-1;
1.1 noro 2697: if (typ(A)!=t_MAT) err(talker,"not a matrix in nfhermite");
2698: if (typ(I)!=t_VEC || lg(I)!=k+1)
2699: err(talker,"not a correct ideal list in nfhermite");
2700: if (!k) err(talker,"not a matrix of maximal rank in nfhermite");
1.2 ! noro 2701: m = lg(A[1])-1;
! 2702: if (k < m) err(talker,"not a matrix of maximal rank in nfhermite");
1.1 noro 2703:
2704: av = avma; lim = stack_lim(av, 1);
1.2 ! noro 2705: A = matalgtobasis(nf,A);
! 2706: I = dummycopy(I);
! 2707: J = zerovec(k); def = k+1;
1.1 noro 2708: for (i=m; i>=1; i--)
2709: {
2710: GEN den,p4,p5,p6,u,v,newid, invnewid = NULL;
2711:
2712: def--; j=def; while (j>=1 && gcmp0(gcoeff(A,i,j))) j--;
2713: if (!j) err(talker,"not a matrix of maximal rank in nfhermite");
1.2 ! noro 2714: if (j==def) j--; else { swap(A[j], A[def]); swap(I[j], I[def]); }
! 2715: p1 = gcoeff(A,i,def);
! 2716: p2 = element_inv(nf, p1);
! 2717: A[def] = (long)element_mulvec(nf,p2,(GEN)A[def]);
! 2718: I[def] = (long)idealmul(nf,p1,(GEN)I[def]);
1.1 noro 2719: for ( ; j; j--)
2720: {
1.2 ! noro 2721: p1 = gcoeff(A,i,j);
1.1 noro 2722: if (gcmp0(p1)) continue;
2723:
1.2 ! noro 2724: p2 = idealmul(nf,p1,(GEN)I[j]);
1.1 noro 2725: newid = idealadd(nf,p2,(GEN)I[def]);
1.2 ! noro 2726:
1.1 noro 2727: invnewid = hnfideal_inv(nf,newid);
2728: p4 = idealmul(nf, p2, invnewid);
2729: p5 = idealmul(nf,(GEN)I[def],invnewid);
2730: y = idealaddtoone(nf,p4,p5);
1.2 ! noro 2731: u = (GEN)y[1]; u = element_div(nf, u, p1);
! 2732: v = (GEN)y[2];
! 2733: p6 = gsub((GEN)A[j], element_mulvec(nf, p1,(GEN)A[def]));
! 2734: A[def] = ladd(element_mulvec(nf, u, (GEN)A[j]),
! 2735: element_mulvec(nf, v, (GEN)A[def]));
! 2736: A[j] = (long)p6;
! 2737: I[j] = (long)idealmul(nf,idealmul(nf,(GEN)I[j],(GEN)I[def]),invnewid);
! 2738: I[def] = (long)newid; den = denom((GEN)I[j]);
1.1 noro 2739: if (!gcmp1(den))
2740: {
1.2 ! noro 2741: I[j] = lmul(den,(GEN)I[j]);
! 2742: A[j] = ldiv((GEN)A[j],den);
1.1 noro 2743: }
2744: }
2745: if (!invnewid) invnewid = hnfideal_inv(nf,(GEN)I[def]);
1.2 ! noro 2746: J[def] = (long)invnewid;
! 2747: p1 = (GEN)I[def];
1.1 noro 2748: for (j=def+1; j<=k; j++)
2749: {
1.2 ! noro 2750: GEN c = gcoeff(A,i,j);
! 2751: p2 = gsub(element_reduce(nf,c, idealmul(nf,p1,(GEN)J[j])), c);
! 2752: A[j] = ladd((GEN)A[j], element_mulvec(nf, p2,(GEN)A[def]));
1.1 noro 2753: }
2754: if (low_stack(lim, stack_lim(av1,1)))
2755: {
1.2 ! noro 2756: GEN *gptr[3]; gptr[0]=&A; gptr[1]=&I; gptr[2]=&J;
1.1 noro 2757: if(DEBUGMEM>1) err(warnmem,"nfhermite, i = %ld", i);
1.2 ! noro 2758: gerepilemany(av,gptr,3);
1.1 noro 2759: }
2760: }
2761: y = cgetg(3,t_VEC);
1.2 ! noro 2762: A += k-m; A[0] = evaltyp(t_MAT)|evallg(m+1); y[1] = (long)A;
! 2763: I += k-m; I[0] = evaltyp(t_VEC)|evallg(m+1); y[2] = (long)I;
! 2764: return gerepilecopy(av0, y);
! 2765: }
! 2766:
! 2767: static GEN
! 2768: idealmulelt(GEN nf, GEN elt, GEN x)
! 2769: {
! 2770: long t = typ(elt);
! 2771: GEN z;
! 2772: if (t == t_POL || t == t_POLMOD) elt = algtobasis(nf,elt);
! 2773: if (isnfscalar(elt)) elt = (GEN)elt[1];
! 2774: z = element_mulvec(nf, elt, x);
! 2775: settyp(z, t_MAT); return z;
! 2776: }
! 2777:
! 2778: static GEN
! 2779: zero_nfbezout(GEN nf,GEN b, GEN A,GEN B,GEN *u,GEN *v,GEN *w,GEN *di)
! 2780: {
! 2781: gpmem_t av, tetpil;
! 2782: GEN pab,d;
! 2783:
! 2784: d=idealmulelt(nf,b,B); *di=idealinv(nf,idealmat_to_hnf(nf,d));
! 2785: av=avma; pab=idealmul(nf,A,B); tetpil=avma;
! 2786: *w=gerepile(av,tetpil, idealmul(nf,pab,*di));
! 2787: *v=element_inv(nf,b);
! 2788: *u=gzero; return d;
! 2789: }
! 2790:
! 2791: /* Given elements a,b and ideals A, B, outputs d = a.A+b.B and gives
! 2792: * di=d^-1, w=A.B.di, u, v such that au+bv=1 and u in A.di, v in
! 2793: * B.di. Assume A, B non-zero, but a or b can be zero (not both)
! 2794: */
! 2795: static GEN
! 2796: nfbezout(GEN nf,GEN a,GEN b, GEN A,GEN B, GEN *u,GEN *v,GEN *w,GEN *di)
! 2797: {
! 2798: GEN pab,pu,pv,pw,uv,d,dinv,pa,pb,pa1,pb1, *gptr[5];
! 2799: gpmem_t av, tetpil;
! 2800:
! 2801: if (gcmp0(a))
! 2802: {
! 2803: if (gcmp0(b)) err(talker,"both elements zero in nfbezout");
! 2804: return zero_nfbezout(nf,b,A,B,u,v,w,di);
! 2805: }
! 2806: if (gcmp0(b))
! 2807: return zero_nfbezout(nf,a,B,A,v,u,w,di);
! 2808:
! 2809: av = avma;
! 2810: pa = idealmulelt(nf,a,A);
! 2811: pb = idealmulelt(nf,b,B);
! 2812:
! 2813: d=idealadd(nf,pa,pb); dinv=idealinv(nf,d);
! 2814: pa1=idealmullll(nf,pa,dinv);
! 2815: pb1=idealmullll(nf,pb,dinv);
! 2816: uv=idealaddtoone(nf,pa1,pb1);
! 2817: pab=idealmul(nf,A,B); tetpil=avma;
! 2818:
! 2819: pu=element_div(nf,(GEN)uv[1],a);
! 2820: pv=element_div(nf,(GEN)uv[2],b);
! 2821: d=gcopy(d); dinv=gcopy(dinv);
! 2822: pw=idealmul(nf,pab,dinv);
! 2823:
! 2824: *u=pu; *v=pv; *w=pw; *di=dinv;
! 2825: gptr[0]=u; gptr[1]=v; gptr[2]=w; gptr[3]=di;
! 2826: gptr[4]=&d; gerepilemanysp(av,tetpil,gptr,5);
! 2827: return d;
1.1 noro 2828: }
2829:
2830: /* A torsion module M over Z_K will be given by a row vector [A,I,J] with
2831: * three components. I=[b_1,...,b_n] is a row vector of k fractional ideals
2832: * given in HNF, J=[a_1,...,a_n] is a row vector of n fractional ideals in
2833: * HNF. A is an nxn matrix (same n) such that if A_j is the j-th column of A
2834: * and e_n is the canonical basis of K^n, then
2835: * M=(b_1e_1+...+b_ne_n)/(a_1A_1+...a_nA_n)
2836: */
2837:
2838: /* We input a torsion module x=[A,I,J] as above, and output the
2839: * smith normal form as K=[\c_1,...,\c_n] such that x=Z_K/\c_1+...+Z_K/\c_n.
2840: */
2841: GEN
2842: nfsmith(GEN nf, GEN x)
2843: {
1.2 ! noro 2844: long i, j, k, l, c, n, m, N;
! 2845: gpmem_t av, tetpil, lim;
1.1 noro 2846: GEN p1,p2,p3,p4,z,b,u,v,w,d,dinv,unnf,A,I,J;
2847:
2848: nf=checknf(nf); N=degpol(nf[1]);
2849: if (typ(x)!=t_VEC || lg(x)!=4) err(talker,"not a module in nfsmith");
2850: A=(GEN)x[1]; I=(GEN)x[2]; J=(GEN)x[3];
2851: if (typ(A)!=t_MAT) err(talker,"not a matrix in nfsmith");
2852: n=lg(A)-1;
2853: if (typ(I)!=t_VEC || lg(I)!=n+1 || typ(J)!=t_VEC || lg(J)!=n+1)
2854: err(talker,"not a correct ideal list in nfsmith");
2855: if (!n) err(talker,"not a matrix of maximal rank in nfsmith");
2856: m=lg(A[1])-1;
2857: if (n<m) err(talker,"not a matrix of maximal rank in nfsmith");
2858: if (n>m) err(impl,"nfsmith for non square matrices");
2859:
2860: av=avma; lim=stack_lim(av,1);
2861: p1 = cgetg(n+1,t_MAT); for (j=1; j<=n; j++) p1[j]=A[j];
2862: A = p1; I = dummycopy(I); J=dummycopy(J);
2863: for (j=1; j<=n; j++)
2864: if (typ(I[j])!=t_MAT) I[j]=(long)idealhermite_aux(nf,(GEN)I[j]);
2865: for (j=1; j<=n; j++)
2866: if (typ(J[j])!=t_MAT) J[j]=(long)idealhermite_aux(nf,(GEN)J[j]);
2867: for (i=n; i>=2; i--)
2868: {
2869: do
2870: {
1.2 ! noro 2871: c = 0;
1.1 noro 2872: for (j=i-1; j>=1; j--)
2873: {
1.2 ! noro 2874: p1=gcoeff(A,i,j); if (gcmp0(p1)) continue;
! 2875:
! 2876: p2 = gcoeff(A,i,i);
! 2877: d = nfbezout(nf,p2,p1,(GEN)J[i],(GEN)J[j],&u,&v,&w,&dinv);
! 2878: if (gcmp0(u)) b = element_mulvec(nf,v,(GEN)A[j]);
! 2879: else
! 2880: {
! 2881: if (gcmp0(v)) b = element_mulvec(nf,u,(GEN)A[i]);
! 2882: else
! 2883: b = gadd(element_mulvec(nf,u,(GEN)A[i]),
! 2884: element_mulvec(nf,v,(GEN)A[j]));
! 2885: }
! 2886: A[j] = lsub(element_mulvec(nf,p2,(GEN)A[j]),
! 2887: element_mulvec(nf,p1,(GEN)A[i]));
! 2888: A[i] = (long)b; J[j] = (long)w; J[i] = (long)d;
1.1 noro 2889: }
2890: for (j=i-1; j>=1; j--)
2891: {
1.2 ! noro 2892: p1 = gcoeff(A,j,i); if (gcmp0(p1)) continue;
! 2893:
! 2894: p2 = gcoeff(A,i,i);
! 2895: d = nfbezout(nf,p2,p1,(GEN)I[i],(GEN)I[j],&u,&v,&w,&dinv);
! 2896: if (gcmp0(u)) b = element_mulvecrow(nf,v,A,j,i);
! 2897: else
! 2898: {
! 2899: if (gcmp0(v))
! 2900: b = element_mulvecrow(nf,u,A,i,i);
! 2901: else
! 2902: b = gadd(element_mulvecrow(nf,u,A,i,i),
! 2903: element_mulvecrow(nf,v,A,j,i));
! 2904: }
! 2905: p3 = gsub(element_mulvecrow(nf,p2,A,j,i),
! 2906: element_mulvecrow(nf,p1,A,i,i));
! 2907: for (k=1; k<=i; k++) { coeff(A,j,k) = p3[k]; coeff(A,i,k) = b[k]; }
! 2908: I[j] = (long)w; I[i] = (long)d; c = 1;
! 2909: }
! 2910: if (c) continue;
! 2911:
! 2912: b=gcoeff(A,i,i); if (gcmp0(b)) break;
! 2913:
! 2914: b=idealmul(nf,b,idealmul(nf,(GEN)J[i],(GEN)I[i]));
! 2915: for (k=1; k<i; k++)
! 2916: for (l=1; l<i; l++)
! 2917: {
! 2918: p3 = gcoeff(A,k,l);
! 2919: if (gcmp0(p3)) continue;
! 2920:
! 2921: p4 = idealmul(nf,p3,idealmul(nf,(GEN)J[l],(GEN)I[k]));
! 2922: if (hnfdivide(b, p4)) continue;
! 2923:
! 2924: b = idealdiv(nf,(GEN)I[k],(GEN)I[i]);
! 2925: p1 = idealdiv(nf,(GEN)J[i], idealmul(nf,p3,(GEN)J[l]));
! 2926: p4 = gauss(p4, b);
! 2927: l=1; while (l<=N && gcmp1(denom((GEN)p4[l]))) l++;
! 2928: if (l>N) err(talker,"bug2 in nfsmith");
! 2929: p3 = element_mulvecrow(nf,(GEN)b[l],A,k,i);
! 2930: for (l=1; l<=i; l++)
! 2931: coeff(A,i,l) = ladd(gcoeff(A,i,l),(GEN)p3[l]);
1.1 noro 2932:
1.2 ! noro 2933: k = l = i; c = 1;
! 2934: }
1.1 noro 2935: if (low_stack(lim, stack_lim(av,1)))
2936: {
1.2 ! noro 2937: GEN *gptr[3]; gptr[0]=&A; gptr[1]=&I; gptr[2]=&J;
1.1 noro 2938: if(DEBUGMEM>1) err(warnmem,"nfsmith");
1.2 ! noro 2939: gerepilemany(av,gptr,3);
1.1 noro 2940: }
2941: }
2942: while (c);
2943: }
2944: unnf=gscalcol_i(gun,N);
2945: p1=gcoeff(A,1,1); coeff(A,1,1)=(long)unnf;
2946: J[1]=(long)idealmul(nf,p1,(GEN)J[1]);
2947: for (i=2; i<=n; i++)
2948: if (!gegal(gcoeff(A,i,i),unnf)) err(talker,"bug in nfsmith");
2949: tetpil=avma; z=cgetg(n+1,t_VEC);
2950: for (i=1; i<=n; i++) z[i]=(long)idealmul(nf,(GEN)I[i],(GEN)J[i]);
2951: return gerepile(av,tetpil,z);
2952: }
2953:
2954: #define trivlift(x) ((typ(x)==t_POLMOD)? (GEN)x[2]: lift_intern(x))
2955:
2956: GEN
1.2 ! noro 2957: element_mulmodpr(GEN nf, GEN x, GEN y, GEN modpr)
1.1 noro 2958: {
1.2 ! noro 2959: gpmem_t av=avma;
1.1 noro 2960: GEN p1;
2961:
1.2 ! noro 2962: nf=checknf(nf); checkmodpr(modpr);
1.1 noro 2963: p1 = element_mul(nf,x,y);
1.2 ! noro 2964: return gerepileupto(av,nfreducemodpr(nf,p1,modpr));
1.1 noro 2965: }
2966:
2967: GEN
1.2 ! noro 2968: element_divmodpr(GEN nf, GEN x, GEN y, GEN modpr)
1.1 noro 2969: {
1.2 ! noro 2970: gpmem_t av=avma;
1.1 noro 2971: GEN p1;
2972:
1.2 ! noro 2973: nf=checknf(nf); checkmodpr(modpr);
1.1 noro 2974: p1=lift_intern(gdiv(gmodulcp(gmul((GEN)nf[7],trivlift(x)), (GEN)nf[1]),
2975: gmodulcp(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1])));
1.2 ! noro 2976: p1=algtobasis_i(nf,p1);
! 2977: return gerepileupto(av,nfreducemodpr(nf,p1,modpr));
1.1 noro 2978: }
2979:
2980: GEN
1.2 ! noro 2981: element_invmodpr(GEN nf, GEN y, GEN modpr)
1.1 noro 2982: {
1.2 ! noro 2983: gpmem_t av=avma;
1.1 noro 2984: GEN p1;
2985:
1.2 ! noro 2986: p1=QX_invmod(gmul((GEN)nf[7],trivlift(y)), (GEN)nf[1]);
! 2987: p1=algtobasis_i(nf,p1);
! 2988: return gerepileupto(av,nfreducemodpr(nf,p1,modpr));
1.1 noro 2989: }
2990:
2991: GEN
1.2 ! noro 2992: element_powmodpr(GEN nf,GEN x,GEN k,GEN pr)
1.1 noro 2993: {
1.2 ! noro 2994: gpmem_t av=avma;
! 2995: GEN z,T,p,modpr;
1.1 noro 2996:
1.2 ! noro 2997: nf = checknf(nf);
! 2998: modpr = nf_to_ff_init(nf,&pr,&T,&p);
! 2999: z = nf_to_ff(nf,x,modpr);
! 3000: z = FpXQ_pow(z,k,T,p);
! 3001: z = ff_to_nf(z,modpr);
! 3002: if (typ(z) != t_COL) z = algtobasis(nf, z);
! 3003: return gerepilecopy(av, z);
1.1 noro 3004: }
3005:
3006: GEN
1.2 ! noro 3007: nfkermodpr(GEN nf, GEN x, GEN pr)
1.1 noro 3008: {
1.2 ! noro 3009: gpmem_t av = avma;
! 3010: GEN T,p,modpr;
1.1 noro 3011:
1.2 ! noro 3012: nf = checknf(nf);
1.1 noro 3013: if (typ(x)!=t_MAT) err(typeer,"nfkermodpr");
1.2 ! noro 3014: modpr = nf_to_ff_init(nf, &pr,&T,&p);
! 3015: x = modprM(lift(x), nf, modpr);
! 3016: return gerepilecopy(av, modprM_lift(FqM_ker(x,T,p), modpr));
1.1 noro 3017: }
3018:
3019: /* a.x=b ou b est un vecteur */
3020: GEN
1.2 ! noro 3021: nfsolvemodpr(GEN nf, GEN a, GEN b, GEN pr)
1.1 noro 3022: {
1.2 ! noro 3023: gpmem_t av = avma;
! 3024: GEN T,p,modpr;
1.1 noro 3025:
1.2 ! noro 3026: nf = checknf(nf);
1.1 noro 3027: if (typ(a)!=t_MAT) err(typeer,"nfsolvemodpr");
1.2 ! noro 3028: modpr = nf_to_ff_init(nf, &pr,&T,&p);
! 3029: a = modprM(lift(a), nf, modpr);
! 3030: b = modprM(lift(b), nf, modpr);
! 3031: return gerepilecopy(av, modprM_lift(FqM_gauss(a,b,T,p), modpr));
1.1 noro 3032: }
3033:
1.2 ! noro 3034: #if 0
1.1 noro 3035: GEN
1.2 ! noro 3036: nfsuppl(GEN nf, GEN x, GEN pr)
1.1 noro 3037: {
1.2 ! noro 3038: gpmem_t av = avma;
! 3039: GEN T,p,modpr;
1.1 noro 3040:
1.2 ! noro 3041: nf = checknf(nf);
! 3042: if (typ(x)!=t_MAT) err(typeer,"nfkermodpr");
! 3043: modpr = nf_to_ff_init(nf, &pr,&T,&p);
! 3044: x = modprM(lift(x), nf, modpr);
! 3045: return gerepilecopy(av, modprM_lift(FqM_suppl(x,T,p), modpr));
1.1 noro 3046: }
1.2 ! noro 3047: #endif
1.1 noro 3048:
3049: /* Given a pseudo-basis pseudo, outputs a multiple of its ideal determinant */
3050: GEN
3051: nfdetint(GEN nf,GEN pseudo)
3052: {
3053: GEN pass,c,v,det1,piv,pivprec,vi,p1,x,I,unnf,zeronf,id,idprod;
1.2 ! noro 3054: long i, j, k, rg, n, n1, m, m1, cm=0, N;
! 3055: gpmem_t av=avma, av1, tetpil, lim;
1.1 noro 3056:
3057: nf=checknf(nf); N=degpol(nf[1]);
3058: if (typ(pseudo)!=t_VEC || lg(pseudo)!=3)
3059: err(talker,"not a module in nfdetint");
3060: x=(GEN)pseudo[1]; I=(GEN)pseudo[2];
3061: if (typ(x)!=t_MAT) err(talker,"not a matrix in nfdetint");
3062: n1=lg(x); n=n1-1; if (!n) return gun;
3063:
3064: m1=lg(x[1]); m=m1-1;
3065: if (typ(I)!=t_VEC || lg(I)!=n1)
3066: err(talker,"not a correct ideal list in nfdetint");
3067:
3068: unnf=gscalcol_i(gun,N); zeronf=zerocol(N);
3069: id=idmat(N); c=new_chunk(m1); for (k=1; k<=m; k++) c[k]=0;
3070: piv = pivprec = unnf;
3071:
3072: av1=avma; lim=stack_lim(av1,1);
3073: det1=idprod=gzero; /* dummy for gerepilemany */
3074: pass=cgetg(m1,t_MAT); v=cgetg(m1,t_COL);
3075: for (j=1; j<=m; j++)
3076: {
3077: v[j] = zero; /* dummy */
3078: p1=cgetg(m1,t_COL); pass[j]=(long)p1;
3079: for (i=1; i<=m; i++) p1[i]=(long)zeronf;
3080: }
3081: for (rg=0,k=1; k<=n; k++)
3082: {
3083: long t = 0;
3084: for (i=1; i<=m; i++)
3085: if (!c[i])
3086: {
3087: vi=element_mul(nf,piv,gcoeff(x,i,k));
3088: for (j=1; j<=m; j++)
3089: if (c[j]) vi=gadd(vi,element_mul(nf,gcoeff(pass,i,j),gcoeff(x,j,k)));
3090: v[i]=(long)vi; if (!t && !gcmp0(vi)) t=i;
3091: }
3092: if (t)
3093: {
3094: pivprec = piv;
3095: if (rg == m-1)
3096: {
3097: if (!cm)
3098: {
3099: cm=1; idprod = id;
3100: for (i=1; i<=m; i++)
3101: if (i!=t)
3102: idprod = (idprod==id)? (GEN)I[c[i]]
3103: : idealmul(nf,idprod,(GEN)I[c[i]]);
3104: }
3105: p1 = idealmul(nf,(GEN)v[t],(GEN)I[k]); c[t]=0;
3106: det1 = (typ(det1)==t_INT)? p1: idealadd(nf,p1,det1);
3107: }
3108: else
3109: {
3110: rg++; piv=(GEN)v[t]; c[t]=k;
3111: for (i=1; i<=m; i++)
3112: if (!c[i])
3113: {
3114: for (j=1; j<=m; j++)
3115: if (c[j] && j!=t)
3116: {
3117: p1=gsub(element_mul(nf,piv,gcoeff(pass,i,j)),
3118: element_mul(nf,(GEN)v[i],gcoeff(pass,t,j)));
3119: coeff(pass,i,j) = rg>1? (long) element_div(nf,p1,pivprec)
3120: : (long) p1;
3121: }
3122: coeff(pass,i,t)=lneg((GEN)v[i]);
3123: }
3124: }
3125: }
3126: if (low_stack(lim, stack_lim(av1,1)))
3127: {
3128: GEN *gptr[6];
3129: if(DEBUGMEM>1) err(warnmem,"nfdetint");
3130: gptr[0]=&det1; gptr[1]=ϖ gptr[2]=&pivprec; gptr[3]=&pass;
3131: gptr[4]=&v; gptr[5]=&idprod; gerepilemany(av1,gptr,6);
3132: }
3133: }
3134: if (!cm) { avma=av; return gscalmat(gzero,N); }
3135: tetpil=avma; return gerepile(av,tetpil,idealmul(nf,idprod,det1));
3136: }
3137:
3138: /* clean in place (destroy x) */
3139: static void
3140: nfcleanmod(GEN nf, GEN x, long lim, GEN detmat)
3141: {
3142: long lx=lg(x),i;
3143:
3144: if (lim<=0 || lim>=lx) lim=lx-1;
3145: for (i=1; i<=lim; i++)
3146: x[i]=(long)element_reduce(nf,(GEN)x[i],detmat);
3147: }
3148:
3149: /* A usage interne. Pas de verifs ni gestion de pile */
3150: GEN
3151: idealoplll(GEN op(GEN,GEN,GEN), GEN nf, GEN x, GEN y)
3152: {
3153: GEN z = op(nf,x,y), den = denom(z);
3154:
3155: if (gcmp1(den)) den = NULL; else z=gmul(den,z);
3156: z=gmul(z,lllintpartial(z));
3157: return den? gdiv(z,den): z;
3158: }
3159:
3160: GEN
3161: nfhermitemod(GEN nf, GEN pseudo, GEN detmat)
3162: {
1.2 ! noro 3163: long li, co, i, j, jm1, def, ldef, N;
! 3164: gpmem_t av0=avma, av, lim;
! 3165: GEN b,q,w,p1,d,u,v,den,x,I,J,dinv,wh,unnf;
1.1 noro 3166:
3167: nf=checknf(nf); N=degpol(nf[1]);
3168: if (typ(pseudo)!=t_VEC || lg(pseudo)!=3)
3169: err(talker,"not a module in nfhermitemod");
1.2 ! noro 3170: x = (GEN)pseudo[1];
! 3171: I = (GEN)pseudo[2];
! 3172: if (typ(x) != t_MAT) err(talker,"not a matrix in nfhermitemod");
! 3173: co = lg(x);
! 3174: if (typ(I) != t_VEC || lg(I) != co)
1.1 noro 3175: err(talker,"not a correct ideal list in nfhermitemod");
3176: if (co==1) return cgetg(1,t_MAT);
3177:
1.2 ! noro 3178: li= lg(x[1]);
! 3179: x = matalgtobasis(nf, x);
! 3180: I = dummycopy(I);
! 3181: unnf = gscalcol(gun,N);
1.1 noro 3182:
1.2 ! noro 3183: den = denom(detmat); if (!gcmp1(den)) detmat = gmul(den,detmat);
! 3184: detmat = lllint_ip(detmat,100);
1.1 noro 3185:
1.2 ! noro 3186: av = avma; lim = stack_lim(av,1);
! 3187: def = co; ldef = (li>co)? li-co+1: 1;
1.1 noro 3188: for (i=li-1; i>=ldef; i--)
3189: {
3190: def--; j=def-1; while (j && gcmp0(gcoeff(x,i,j))) j--;
3191: while (j)
3192: {
3193: jm1=j-1; if (!jm1) jm1=def;
3194: d=nfbezout(nf,gcoeff(x,i,j),gcoeff(x,i,jm1),(GEN)I[j],(GEN)I[jm1],
3195: &u,&v,&w,&dinv);
3196: if (gcmp0(u))
3197: p1 = element_mulvec(nf,v,(GEN)x[jm1]);
3198: else
3199: {
3200: p1 = element_mulvec(nf,u,(GEN)x[j]);
3201: if (!gcmp0(v)) p1=gadd(p1, element_mulvec(nf,v,(GEN)x[jm1]));
3202: }
3203: x[j]=lsub(element_mulvec(nf,gcoeff(x,i,j),(GEN)x[jm1]),
3204: element_mulvec(nf,gcoeff(x,i,jm1),(GEN)x[j]));
3205: nfcleanmod(nf,(GEN)x[j],i,idealdivlll(nf,detmat,w));
3206: nfcleanmod(nf,p1,i,idealmullll(nf,detmat,dinv));
3207: x[jm1]=(long)p1; I[j]=(long)w; I[jm1]=(long)d;
3208: j--; while (j && gcmp0(gcoeff(x,i,j))) j--;
3209: }
3210: if (low_stack(lim, stack_lim(av,1)))
3211: {
1.2 ! noro 3212: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&I;
1.1 noro 3213: if(DEBUGMEM>1) err(warnmem,"[1]: nfhermitemod");
1.2 ! noro 3214: gerepilemany(av,gptr,2);
1.1 noro 3215: }
3216: }
1.2 ! noro 3217: b = detmat; wh = cgetg(li,t_MAT); def--;
1.1 noro 3218: for (i=li-1; i>=1; i--)
3219: {
3220: d = nfbezout(nf,gcoeff(x,i,i+def),unnf,(GEN)I[i+def],b,&u,&v,&w,&dinv);
3221: p1 = element_mulvec(nf,u,(GEN)x[i+def]);
3222: nfcleanmod(nf,p1,i,idealmullll(nf,b,dinv));
1.2 ! noro 3223: wh[i] = (long)p1; coeff(wh,i,i) = (long)unnf;
! 3224: I[i+def] = (long)d;
! 3225: if (i>1) b = idealmul(nf,b,dinv);
1.1 noro 3226: }
3227: J=cgetg(li,t_VEC); J[1]=zero;
1.2 ! noro 3228: for (j=2; j<li; j++) J[j] = (long)idealinv(nf,(GEN)I[j+def]);
1.1 noro 3229: for (i=li-2; i>=1; i--)
3230: {
3231: for (j=i+1; j<li; j++)
3232: {
1.2 ! noro 3233: q = idealmul(nf,(GEN)I[i+def],(GEN)J[j]);
! 3234: p1 = gsub(element_reduce(nf,gcoeff(wh,i,j),q), gcoeff(wh,i,j));
! 3235: wh[j] = ladd((GEN)wh[j],element_mulvec(nf,p1,(GEN)wh[i]));
1.1 noro 3236: }
3237: if (low_stack(lim, stack_lim(av,1)))
3238: {
1.2 ! noro 3239: GEN *gptr[3]; gptr[0]=&wh; gptr[1]=&I; gptr[2]=&J;
1.1 noro 3240: if(DEBUGMEM>1) err(warnmem,"[2]: nfhermitemod");
1.2 ! noro 3241: gerepilemany(av,gptr,3);
1.1 noro 3242: }
3243: }
1.2 ! noro 3244: I += def; I[0] = evaltyp(t_VEC)|evallg(li);
! 3245: p1=cgetg(3,t_VEC);
! 3246: p1[1] = (long)wh;
! 3247: p1[2] = (long)I; return gerepilecopy(av0, p1);
1.1 noro 3248: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to base4.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