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