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