Annotation of OpenXM_contrib/pari/src/basemath/buch4.c, Revision 1.1
1.1 ! maekawa 1: /*******************************************************************/
! 2: /* */
! 3: /* S-CLASS GROUP AND NORM SYMBOLS */
! 4: /* (Denis Simon, desimon@math.u-bordeaux.fr) */
! 5: /* */
! 6: /*******************************************************************/
! 7: /* $Id: buch4.c,v 1.1.1.1 1999/09/16 13:47:29 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: static long
! 11: psquare(GEN a,GEN p)
! 12: {
! 13: long v;
! 14: GEN ap;
! 15:
! 16: if (gcmp0(a) || gcmp1(a)) return 1;
! 17:
! 18: if (!cmpis(p,2))
! 19: {
! 20: v=vali(a); if (v&1) return 0;
! 21: return (smodis(shifti(a,-v),8)==1);
! 22: }
! 23:
! 24: ap=stoi(1); v=pvaluation(a,p,&ap);
! 25: if (v&1) return 0;
! 26: return (kronecker(ap,p)==1);
! 27: }
! 28:
! 29: static long
! 30: lemma6(GEN pol,GEN p,long nu,GEN x)
! 31: {
! 32: long i,lambda,mu,ltop=avma;
! 33: GEN gx,gpx;
! 34:
! 35: for (i=lgef(pol)-2,gx=(GEN) pol[i+1]; i>1; i--)
! 36: gx=addii(mulii(gx,x),(GEN) pol[i]);
! 37: if (psquare(gx,p)) return 1;
! 38:
! 39: for (i=lgef(pol)-2,gpx=mulis((GEN) pol[i+1],i-1); i>2; i--)
! 40: gpx=addii(mulii(gpx,x),mulis((GEN) pol[i],i-2));
! 41:
! 42: lambda=pvaluation(gx,p,&gx);
! 43: if (gcmp0(gpx)) mu=BIGINT; else mu=pvaluation(gpx,p,&gpx);
! 44: avma=ltop;
! 45:
! 46: if (lambda>(mu<<1)) return 1;
! 47: if (lambda>=(nu<<1) && mu>=nu) return 0;
! 48: return -1;
! 49: }
! 50:
! 51: static long
! 52: lemma7(GEN pol,long nu,GEN x)
! 53: { long i,odd4,lambda,mu,mnl,ltop=avma;
! 54: GEN gx,gpx,oddgx;
! 55:
! 56: for (i=lgef(pol)-2,gx=(GEN) pol[i+1]; i>1; i--)
! 57: gx=addii(mulii(gx,x),(GEN) pol[i]);
! 58: if (psquare(gx,gdeux)) return 1;
! 59:
! 60: for (i=lgef(pol)-2,gpx=gmulgs((GEN) pol[i+1],i-1); i>2; i--)
! 61: gpx=gadd(gmul(gpx,x),gmulgs((GEN) pol[i],i-2));
! 62:
! 63: lambda=vali(gx);
! 64: if (gcmp0(gpx)) mu=BIGINT; else mu=vali(gpx);
! 65: oddgx=shifti(gx,-lambda);
! 66: mnl=mu+nu-lambda;
! 67: odd4=smodis(oddgx,4);
! 68: avma=ltop;
! 69: if (lambda>(mu<<1)) return 1;
! 70: if (nu > mu)
! 71: { if (mnl==1 && (lambda&1) == 0) return 1;
! 72: if (mnl==2 && (lambda&1) == 0 && odd4==1) return 1;
! 73: }
! 74: else
! 75: { if (lambda>=(nu<<1)) return 0;
! 76: if (lambda==((nu-1)<<1) && odd4==1) return 0;
! 77: }
! 78: return -1;
! 79: }
! 80:
! 81: static long
! 82: zpsol(GEN pol,GEN p,long nu,GEN pnu,GEN x0)
! 83: {
! 84: long i,result,ltop=avma;
! 85: GEN x,pnup;
! 86:
! 87: result = (cmpis(p,2)) ? lemma6(pol,p,nu,x0) : lemma7(pol,nu,x0);
! 88: if (result==+1) return 1; if (result==-1) return 0;
! 89: x=gcopy(x0); pnup=mulii(pnu,p);
! 90: for (i=0; i<itos(p); i++)
! 91: {
! 92: x=addii(x,pnu);
! 93: if (zpsol(pol,p,nu+1,pnup,x)) { avma=ltop; return 1; }
! 94: }
! 95: avma=ltop; return 0;
! 96: }
! 97:
! 98: /* vaut 1 si l'equation y^2=Pol(x) a une solution p-adique entiere
! 99: * 0 sinon. Les coefficients sont entiers.
! 100: */
! 101: long
! 102: zpsoluble(GEN pol,GEN p)
! 103: {
! 104: if ((typ(pol)!=t_POL && typ(pol)!=t_INT) || typ(p)!=t_INT )
! 105: err(typeer,"zpsoluble");
! 106: return zpsol(pol,p,0,gun,gzero);
! 107: }
! 108:
! 109: /* vaut 1 si l'equation y^2=Pol(x) a une solution p-adique rationnelle
! 110: * (eventuellement infinie), 0 sinon. Les coefficients sont entiers.
! 111: */
! 112: long
! 113: qpsoluble(GEN pol,GEN p)
! 114: {
! 115: if ((typ(pol)!=t_POL && typ(pol)!=t_INT) || typ(p)!=t_INT )
! 116: err(typeer,"qpsoluble");
! 117: if (zpsol(pol,p,0,gun,gzero)) return 1;
! 118: return (zpsol(polrecip(pol),p,1,p,gzero));
! 119: }
! 120:
! 121: /* p premier a 2. renvoie 1 si a est un carre dans ZK_p, 0 sinon */
! 122: static long
! 123: psquarenf(GEN nf,GEN a,GEN p)
! 124: {
! 125: long v,ltop=avma;
! 126: GEN norm,ap;
! 127:
! 128: if (gcmp0(a)) return 1;
! 129: v=idealval(nf,a,p); if (v&1) return 0;
! 130: ap=gdiv(a,gpuigs(basistoalg(nf,(GEN)p[2]),v));
! 131:
! 132: norm=gshift(idealnorm(nf,p),-1);
! 133: ap=gmul(ap,gmodulsg(1,(GEN)p[1]));
! 134: ap=gaddgs(gpui(ap,norm,0),-1);
! 135: if (gcmp0(ap)) { avma=ltop; return 1; }
! 136: ap=lift(lift(ap));
! 137: v = idealval(nf,ap,p); avma=ltop;
! 138: return (v>0);
! 139: }
! 140:
! 141: static long
! 142: check2(GEN nf, GEN a, GEN zinit)
! 143: {
! 144: GEN zlog=zideallog(nf,a,zinit), p1 = gmael(zinit,2,2);
! 145: long i;
! 146:
! 147: for (i=1; i<lg(p1); i++)
! 148: if (!mpodd((GEN)p1[i]) && mpodd((GEN)zlog[i])) return 0;
! 149: return 1;
! 150: }
! 151:
! 152: /* p divise 2. renvoie 1 si a est un carre dans ZK_p, 0 sinon */
! 153: static long
! 154: psquare2nf(GEN nf,GEN a,GEN p,GEN zinit)
! 155: {
! 156: long v, ltop=avma;
! 157:
! 158: if (gcmp0(a)) return 1;
! 159: v=idealval(nf,a,p); if (v&1) return 0;
! 160: a = gdiv(a,gmodulcp(gpuigs(gmul((GEN)nf[7],(GEN)p[2]),v),(GEN)nf[1]));
! 161: v = check2(nf,a,zinit); avma = ltop; return v;
! 162: }
! 163:
! 164: /* p divise 2, (a,p)=1. renvoie 1 si a est un carre de (ZK / p^q)*, 0 sinon. */
! 165: static long
! 166: psquare2qnf(GEN nf,GEN a,GEN p,long q)
! 167: {
! 168: long v, ltop=avma;
! 169: GEN zinit = zidealstarinit(nf,idealpows(nf,p,q));
! 170:
! 171: v = check2(nf,a,zinit); avma = ltop; return v;
! 172: }
! 173:
! 174: static long
! 175: lemma6nf(GEN nf,GEN pol,GEN p,long nu,GEN x)
! 176: {
! 177: long i,lambda,mu,ltop=avma;
! 178: GEN gx,gpx;
! 179:
! 180: for (i=lgef(pol)-2,gx=(GEN) pol[i+1]; i>1; i--)
! 181: gx = gadd(gmul(gx,x),(GEN) pol[i]);
! 182: if (psquarenf(nf,gx,p)) { avma=ltop; return 1; }
! 183: lambda = idealval(nf,gx,p);
! 184:
! 185: for (i=lgef(pol)-2,gpx=gmulgs((GEN) pol[i+1],i-1); i>2; i--)
! 186: gpx=gadd(gmul(gpx,x),gmulgs((GEN) pol[i],i-2));
! 187: mu = gcmp0(gpx)? BIGINT: idealval(nf,gpx,p);
! 188:
! 189: avma=ltop;
! 190: if (lambda > mu<<1) return 1;
! 191: if (lambda >= nu<<1 && mu >= nu) return 0;
! 192: return -1;
! 193: }
! 194:
! 195: static long
! 196: lemma7nf(GEN nf,GEN pol,GEN p,long nu,GEN x,GEN zinit)
! 197: {
! 198: long res,i,lambda,mu,q,ltop=avma;
! 199: GEN gx,gpx,p1;
! 200:
! 201: for (i=lgef(pol)-2, gx=(GEN) pol[i+1]; i>1; i--)
! 202: gx=gadd(gmul(gx,x),(GEN) pol[i]);
! 203: if (psquare2nf(nf,gx,p,zinit)) { avma=ltop; return 1; }
! 204: lambda=idealval(nf,gx,p);
! 205:
! 206: for (i=lgef(pol)-2,gpx=gmulgs((GEN) pol[i+1],i-1); i>2; i--)
! 207: gpx=gadd(gmul(gpx,x),gmulgs((GEN) pol[i],i-2));
! 208: if (!gcmp0(gpx)) mu=idealval(nf,gpx,p); else mu=BIGINT;
! 209:
! 210: if (lambda>(mu<<1)) { avma=ltop; return 1; }
! 211: if (nu > mu)
! 212: {
! 213: if (lambda&1) { avma=ltop; return -1; }
! 214: q=mu+nu-lambda; res=1;
! 215: }
! 216: else
! 217: {
! 218: if (lambda>=(nu<<1)) { avma=ltop; return 0; }
! 219: if (lambda&1) { avma=ltop; return -1; }
! 220: q=(nu<<1)-lambda; res=0;
! 221: }
! 222: if (q > itos((GEN) p[3])<<1) { avma=ltop; return -1; }
! 223: p1 = gmodulcp(gpuigs(gmul((GEN)nf[7],(GEN)p[2]), lambda), (GEN)nf[1]);
! 224: if (!psquare2qnf(nf,gdiv(gx,p1), p,q)) res = -1;
! 225: avma=ltop; return res;
! 226: }
! 227:
! 228: static long
! 229: zpsolnf(GEN nf,GEN pol,GEN p,long nu,GEN pnu,GEN x0,GEN repr,GEN zinit)
! 230: {
! 231: long i,result,ltop=avma;
! 232: GEN pnup;
! 233:
! 234: nf=checknf(nf);
! 235: if (cmpis((GEN) p[1],2))
! 236: result=lemma6nf(nf,pol,p,nu,x0);
! 237: else
! 238: result=lemma7nf(nf,pol,p,nu,x0,zinit);
! 239: if (result== 1) return 1;
! 240: if (result==-1) return 0;
! 241: pnup=gmul(pnu,gmodulcp(gmul((GEN) nf[7],(GEN) p[2]),(GEN) nf[1]));
! 242: nu++;
! 243: for (i=1; i<lg(repr); i++)
! 244: if (zpsolnf(nf,pol,p,nu,pnup,gadd(x0,gmul(pnu,(GEN)repr[i])),repr,zinit))
! 245: { avma=ltop; return 1; }
! 246: avma=ltop; return 0;
! 247: }
! 248:
! 249: /* calcule un systeme de representants Zk/p */
! 250: static GEN
! 251: repres(GEN nf,GEN p)
! 252: {
! 253: long i,j,k,f,pp,ppf,ppi;
! 254: GEN mat,fond,rep;
! 255:
! 256: fond=cgetg(1,t_VEC);
! 257: mat=idealhermite(nf,p);
! 258: for (i=1; i<lg(mat); i++)
! 259: if (!gcmp1(gmael(mat,i,i)))
! 260: fond = concatsp(fond,gmael(nf,7,i));
! 261: f=lg(fond)-1;
! 262: pp=itos((GEN) p[1]);
! 263: for (i=1,ppf=1; i<=f; i++) ppf*=pp;
! 264: rep=cgetg(ppf+1,t_VEC);
! 265: rep[1]=zero; ppi=1;
! 266: for (i=0; i<f; i++,ppi*=pp)
! 267: for (j=1; j<pp; j++)
! 268: for (k=1; k<=ppi; k++)
! 269: rep[j*ppi+k]=ladd((GEN) rep[k],gmulsg(j,(GEN) fond[i+1]));
! 270: return gmodulcp(rep,(GEN) nf[1]);
! 271: }
! 272:
! 273: /* =1 si l'equation y^2 = z^deg(pol) * pol(x/z) a une solution rationnelle
! 274: * p-adique (eventuellement (1,y,0) = oo)
! 275: * =0 sinon.
! 276: * Les coefficients de pol doivent etre des entiers de nf.
! 277: * p est un ideal premier sous forme primedec.
! 278: */
! 279: long
! 280: qpsolublenf(GEN nf,GEN pol,GEN p)
! 281: {
! 282: GEN repr,zinit,p1;
! 283: long ltop=avma;
! 284:
! 285: if (gcmp0(pol)) return 1;
! 286: if (typ(pol)!=t_POL) err(notpoler,"qpsolublenf");
! 287: if (typ(p)!=t_VEC || lg(p)!=6)
! 288: err(talker,"not a prime ideal in qpsolublenf");
! 289: nf=checknf(nf);
! 290:
! 291: if (cmpis((GEN) p[1],2))
! 292: {
! 293: if (psquarenf(nf,(GEN) pol[2],p)) return 1;
! 294: if (psquarenf(nf, leading_term(pol),p)) return 1;
! 295: zinit=gzero;
! 296: }
! 297: else
! 298: {
! 299: zinit=zidealstarinit(nf,idealpows(nf,p,1+2*idealval(nf,gdeux,p)));
! 300: if (psquare2nf(nf,(GEN) pol[2],p,zinit)) return 1;
! 301: if (psquare2nf(nf, leading_term(pol),p,zinit)) return 1;
! 302: }
! 303: repr=repres(nf,p);
! 304: if (zpsolnf(nf,pol,p,0,gun,gzero,repr,zinit)) { avma=ltop; return 1; }
! 305: p1 = gmodulcp(gmul((GEN) nf[7],(GEN) p[2]),(GEN) nf[1]);
! 306: if (zpsolnf(nf,polrecip(pol),p,1,p1,gzero,repr,zinit))
! 307: { avma=ltop; return 1; }
! 308:
! 309: avma=ltop; return 0;
! 310: }
! 311:
! 312: /* =1 si l'equation y^2 = pol(x) a une solution entiere p-adique
! 313: * =0 sinon.
! 314: * Les coefficients de pol doivent etre des entiers de nf.
! 315: * p est un ideal premier sous forme primedec.
! 316: */
! 317: long
! 318: zpsolublenf(GEN nf,GEN pol,GEN p)
! 319: {
! 320: GEN repr,zinit;
! 321: long ltop=avma;
! 322:
! 323: if (gcmp0(pol)) return 1;
! 324: if (typ(pol)!=t_POL) err(notpoler,"zpsolublenf");
! 325: if (typ(p)!=t_VEC || lg(p)!=6)
! 326: err(talker,"not a prime ideal in zpsolublenf");
! 327: nf=checknf(nf);
! 328:
! 329: if (cmpis((GEN)p[1],2))
! 330: {
! 331: if (psquarenf(nf,(GEN) pol[2],p)) return 1;
! 332: zinit=gzero;
! 333: }
! 334: else
! 335: {
! 336: zinit=zidealstarinit(nf,idealpows(nf,p,1+2*idealval(nf,gdeux,p)));
! 337: if (psquare2nf(nf,(GEN) pol[2],p,zinit)) return 1;
! 338: }
! 339: repr=repres(nf,p);
! 340: if (zpsolnf(nf,pol,p,0,gun,gzero,repr,zinit)) { avma=ltop; return 1; }
! 341: avma=ltop; return 0;
! 342: }
! 343:
! 344: static long
! 345: hilb2nf(GEN nf,GEN a,GEN b,GEN p)
! 346: {
! 347: GEN pol;
! 348: long ltop=avma;
! 349:
! 350: a=lift(a); b=lift(b);
! 351: pol=polx[0]; pol=gadd(gmul(a,gsqr(pol)),b);
! 352: if (qpsolublenf(nf,pol,p)) { avma=ltop; return 1; }
! 353: avma=ltop; return -1;
! 354: }
! 355:
! 356: /* pr doit etre sous la forme primedec */
! 357: static GEN
! 358: nfmodid2(GEN nf,GEN x,GEN pr)
! 359: {
! 360: x=lift(x);
! 361: x=gmod(x,lift(basistoalg(nf,(GEN)pr[2])));
! 362: return gmul(x,gmodulsg(1,(GEN)pr[1]));
! 363: }
! 364:
! 365: long
! 366: nfhilbertp(GEN nf,GEN a,GEN b,GEN p)
! 367: /* calcule le symbole de Hilbert quadratique local (a,b)_p
! 368: * en l'ideal premier p du corps nf,
! 369: * a et b sont des elements non nuls de nf, sous la forme
! 370: * de polymods ou de polynomes, et p renvoye par primedec.
! 371: */
! 372: {
! 373: GEN aux,aux2;
! 374: long ta=typ(a),tb=typ(b),alpha,beta,sign,rep,ltop=avma;
! 375:
! 376: if ((ta!=t_INT && ta!=t_POLMOD && ta!=t_POL)
! 377: || (tb!=t_INT && tb!=t_POLMOD && tb!=t_POL))
! 378: err (typeer,"nfhilbertp");
! 379: if (gcmp0(a) || gcmp0(b))
! 380: err (talker,"0 argument in nfhilbertp");
! 381: checkprimeid(p); nf=checknf(nf);
! 382:
! 383: if (!cmpis((GEN) p[1],2)) return hilb2nf(nf,a,b,p);
! 384:
! 385: if (ta != t_POLMOD) a=gmodulcp(a,(GEN)nf[1]);
! 386: if (tb != t_POLMOD) b=gmodulcp(b,(GEN)nf[1]);
! 387:
! 388: alpha=idealval(nf,a,p); beta=idealval(nf,b,p);
! 389: if ((alpha&1) == 0 && (beta&1) == 0) { avma=ltop; return 1; }
! 390: aux2=shifti(idealnorm(nf,p),-1);
! 391: if (alpha&1 && beta&1 && mpodd(aux2)) sign=1; else sign=-1;
! 392: aux=nfmodid2(nf,gdiv(gpuigs(a,beta),gpuigs(b,alpha)),p); /* ?????? */
! 393: aux=gaddgs(powgi(aux,aux2),sign);
! 394: aux=lift(lift(aux));
! 395: if (gcmp0(aux)) rep=1; else rep=(idealval(nf,aux,p)>=1);
! 396: avma=ltop; return rep? 1: -1;
! 397: }
! 398:
! 399: /* calcule le symbole de Hilbert quadratique global (a,b):
! 400: * = 1 si l'equation X^2-aY^2-bZ^2=0 a une solution non triviale,
! 401: * =-1 sinon,
! 402: * a et b doivent etre non nuls.
! 403: */
! 404: long
! 405: nfhilbert(GEN nf,GEN a,GEN b)
! 406: {
! 407: long ta=typ(a),tb=typ(b),r1,i,ltop=avma;
! 408: GEN S,al,bl;
! 409:
! 410: nf=checknf(nf);
! 411: if ((ta!=t_INT && ta!=t_POLMOD && ta!=t_POL)
! 412: || (tb!=t_INT && tb!=t_POLMOD && tb!=t_POL))
! 413: err (typeer,"nfhilbert");
! 414: if (gcmp0(a) || gcmp0(b))
! 415: err (talker,"0 argument in nfhilbert");
! 416:
! 417: al=lift(a); bl=lift(b);
! 418: /* solutions locales aux places infinies reelles */
! 419: r1=itos(gmael(nf,2,1));
! 420: for (i=1; i<=r1; i++)
! 421: if (signe(poleval(al,gmael(nf,6,i))) < 0 &&
! 422: signe(poleval(bl,gmael(nf,6,i))) < 0)
! 423: {
! 424: if (DEBUGLEVEL>=4)
! 425: {
! 426: fprintferr("nfhilbert not soluble at a real place\n");
! 427: flusherr();
! 428: }
! 429: avma=ltop; return -1;
! 430: }
! 431: if (ta!=t_POLMOD) a=gmodulcp(a,(GEN)nf[1]);
! 432: if (tb!=t_POLMOD) b=gmodulcp(b,(GEN)nf[1]);
! 433:
! 434: /* solutions locales aux places finies (celles qui divisent 2ab)*/
! 435:
! 436: S=(GEN) idealfactor(nf,gmul(gmulsg(2,a),b))[1];
! 437: /* formule du produit ==> on peut eviter un premier */
! 438: for (i=lg(S)-1; i>1; i--)
! 439: if (nfhilbertp(nf,a,b,(GEN) S[i])==-1)
! 440: {
! 441: if (DEBUGLEVEL >=4)
! 442: {
! 443: fprintferr("nfhilbert not soluble at finite place: ");
! 444: outerr((GEN)S[i]); flusherr();
! 445: }
! 446: avma=ltop; return -1;
! 447: }
! 448: avma=ltop; return 1;
! 449: }
! 450:
! 451: long
! 452: nfhilbert0(GEN nf,GEN a,GEN b,GEN p)
! 453: {
! 454: if (p) return nfhilbertp(nf,a,b,p);
! 455: return nfhilbert(nf,a,b);
! 456: }
! 457:
! 458: GEN vconcat(GEN Q1, GEN Q2);
! 459: GEN mathnfspec(GEN x, GEN *ptperm, GEN *ptdep, GEN *ptB, GEN *ptC);
! 460: /* S a list of prime ideal in primedec format. Return res:
! 461: * res[1] = generators of (S-units / units), as polynomials
! 462: * res[2] = [perm, HB, den], for bnfissunit
! 463: * res[3] = [] (was: log. embeddings of res[1])
! 464: * res[4] = S-regulator ( = R * det(res[2]) * \prod log(Norm(S[i])))
! 465: * res[5] = S class group
! 466: * res[6] = S
! 467: */
! 468: GEN
! 469: bnfsunit(GEN bnf,GEN S,long prec)
! 470: {
! 471: long i,j,ls,ltop=avma,lbot;
! 472: GEN p1,nf,classgp,gen,M,U,H;
! 473: GEN sunit,card,sreg,res,pow,fa = cgetg(3, t_MAT);
! 474:
! 475: if (typ(S) != t_VEC) err(typeer,"bnfsunit");
! 476: bnf = checkbnf(bnf); nf=(GEN)bnf[7];
! 477: classgp=gmael(bnf,8,1);
! 478: gen = (GEN)classgp[3];
! 479:
! 480: sreg = gmael(bnf,8,2);
! 481: res=cgetg(7,t_VEC);
! 482: res[1]=res[2]=res[3]=lgetg(1,t_VEC);
! 483: res[4]=(long)sreg;
! 484: res[5]=(long)classgp;
! 485: res[6]=(long)S; ls=lg(S);
! 486:
! 487: /* M = relation matrix for the S class group (in terms of the class group
! 488: * generators given by gen)
! 489: * 1) ideals in S
! 490: */
! 491: M = cgetg(ls,t_MAT);
! 492: for (i=1; i<ls; i++)
! 493: {
! 494: p1 = (GEN)S[i]; checkprimeid(p1);
! 495: M[i] = (long)isprincipal(bnf,p1);
! 496: }
! 497: /* 2) relations from bnf class group */
! 498: M = concatsp(M, diagonal((GEN) classgp[2]));
! 499:
! 500: /* S class group */
! 501: H = hnfall(M); U = (GEN)H[2]; H= (GEN)H[1];
! 502: card = gun;
! 503: if (lg(H) > 1)
! 504: { /* non trivial (rare!) */
! 505: GEN SNF, ClS = cgetg(4,t_VEC);
! 506:
! 507: SNF = smith2(H); p1 = (GEN)SNF[3];
! 508: card = dethnf_i(p1);
! 509: ClS[1] = (long)card; /* h */
! 510: for(i=1; i<lg(p1); i++)
! 511: if (gcmp1((GEN)p1[i])) break;
! 512: setlg(p1,i);
! 513: ClS[2]=(long)p1; /* cyc */
! 514:
! 515: p1=cgetg(i,t_VEC); pow=invmat((GEN)SNF[1]);
! 516: fa[1] = (long)gen;
! 517: for(i--; i; i--)
! 518: {
! 519: fa[2] = pow[i];
! 520: p1[i] = (long)factorback(fa, nf);
! 521: }
! 522: ClS[3]=(long)p1; /* gen */
! 523: res[5]=(long) ClS;
! 524: }
! 525:
! 526: /* S-units */
! 527: if (ls>1)
! 528: {
! 529: GEN den, Sperm, perm, dep, B, U1 = U;
! 530: long lH, lB;
! 531:
! 532: /* U1 = upper left corner of U, invertible. S * U1 = principal ideals
! 533: * whose generators generate the S-units */
! 534: setlg(U1,ls); p1 = cgetg(ls, t_MAT); /* p1 is junk for mathnfspec */
! 535: for (i=1; i<ls; i++) { setlg(U1[i],ls); p1[i] = lgetg(1,t_COL); }
! 536: H = mathnfspec(U1,&perm,&dep,&B,&p1);
! 537: lH = lg(H);
! 538: lB = lg(B);
! 539: if (lg(dep) > 1 && lg(dep[1]) > 1) err(bugparier,"bnfsunit");
! 540: /* [ H B ] [ H^-1 - H^-1 B ]
! 541: * perm o HNF(U1) = [ 0 Id ], inverse = [ 0 Id ]
! 542: * (permute the rows)
! 543: * S * HNF(U1) = _integral_ generators for S-units = sunit */
! 544: Sperm = cgetg(ls, t_VEC); sunit = cgetg(ls, t_VEC);
! 545: for (i=1; i<ls; i++) Sperm[i] = S[perm[i]]; /* S o perm */
! 546:
! 547: setlg(Sperm, lH); fa[1] = (long)Sperm;
! 548: for (i=1; i<lH; i++)
! 549: {
! 550: fa[2] = H[i];
! 551: sunit[i] = (long)factorback(fa, nf);
! 552: }
! 553: for (i=1; i<lB; i++)
! 554: {
! 555: fa[2] = B[i]; j = lH-1 + i;
! 556: sunit[j] = (long)idealmul(nf, (GEN)Sperm[j], factorback(fa, nf));
! 557: }
! 558: for (i=1; i<ls; i++)
! 559: sunit[i] = isprincipalgenforce(bnf, (GEN)sunit[i])[2];
! 560:
! 561: p1 = cgetg(4,t_VEC);
! 562: den = dethnf_i(H); H = gmul(den, invmat(H));
! 563: p1[1] = (long)perm;
! 564: p1[2] = (long)concatsp(H, gneg(gmul(H,B))); /* top part of inverse * den */
! 565: p1[3] = (long)den; /* keep denominator separately */
! 566: sunit = basistoalg(nf,sunit);
! 567: res[2] = (long)p1; /* HNF in split form perm + (H B) [0 Id missing] */
! 568: res[1] = (long)lift_intern(sunit);
! 569: }
! 570:
! 571: /* S-regulator */
! 572: sreg = gmul(sreg,card);
! 573: for (i=1; i<ls; i++)
! 574: {
! 575: GEN p = (GEN)S[i];
! 576: if (typ(p) == t_VEC) p = (GEN) p[1];
! 577: sreg = gmul(sreg,glog(p,prec));
! 578: }
! 579: res[4]=(long) sreg; lbot=avma;
! 580: return gerepile(ltop,lbot,gcopy(res));
! 581: }
! 582:
! 583: /* cette fonction est l'equivalent de isunit, sauf qu'elle donne le resultat
! 584: * avec des s-unites: si x n'est pas une s-unite alors issunit=[]~;
! 585: * si x est une s-unite alors
! 586: * x=\prod_{i=0}^r {e_i^issunit[i]}*prod{i=r+1}^{r+s} {s_i^issunit[i]}
! 587: * ou les e_i sont les unites du corps (comme dans isunit)
! 588: * et les s_i sont les s-unites calculees par sunit (dans le meme ordre).
! 589: */
! 590: GEN
! 591: bnfissunit(GEN bnf,GEN suni,GEN x)
! 592: {
! 593: long lB,cH,i,k,ls,tetpil, av = avma;
! 594: GEN den,gen,S,v,p1,xp,xm,xb,N,HB,perm;
! 595:
! 596: bnf = checkbnf(bnf);
! 597: if (typ(suni)!=t_VEC || lg(suni)!=7) err(typeer,"bnfissunit");
! 598: switch (typ(x))
! 599: {
! 600: case t_INT: case t_FRAC: case t_FRACN:
! 601: case t_POL: case t_COL:
! 602: x = basistoalg(bnf,x); break;
! 603: case t_POLMOD: break;
! 604: default: err(typeer,"bnfissunit");
! 605: }
! 606: if (gcmp0(x)) return cgetg(1,t_COL);
! 607:
! 608: S = (GEN) suni[6]; ls=lg(S);
! 609: if (ls==1) return isunit(bnf,x);
! 610:
! 611: p1 = (GEN)suni[2];
! 612: perm = (GEN)p1[1];
! 613: HB = (GEN)p1[2]; den = (GEN)p1[3];
! 614: cH = lg(HB[1]) - 1;
! 615: lB = lg(HB) - cH;
! 616: xb = algtobasis(bnf,x); p1 = denom(content(xb));
! 617: N = mulii(gnorm(gmul(x,p1)), p1); /* relevant primes divide N */
! 618: v = cgetg(ls, t_VECSMALL);
! 619: for (i=1; i<ls; i++)
! 620: {
! 621: GEN P = (GEN)S[i];
! 622: v[i] = (resii(N, (GEN)P[1]) == gzero)? element_val(bnf,xb,P): 0;
! 623: }
! 624: /* here, x = S v */
! 625: p1 = cgetg(ls, t_COL);
! 626: for (i=1; i<ls; i++) p1[i] = lstoi(v[perm[i]]); /* p1 = v o perm */
! 627: v = gmul(HB, p1);
! 628: for (i=1; i<=cH; i++)
! 629: {
! 630: GEN w = gdiv((GEN)v[i], den);
! 631: if (typ(w) != t_INT) { avma = av; return cgetg(1,t_COL); }
! 632: v[i] = (long)w;
! 633: }
! 634: p1 += cH;
! 635: p1[0] = evaltyp(t_COL) | evallg(lB);
! 636: v = concatsp(v, p1); /* append bottom of p1 (= [0 Id] part) */
! 637:
! 638: xp = gun; xm = gun; gen = (GEN)suni[1];
! 639: for (i=1; i<ls; i++)
! 640: {
! 641: k = -itos((GEN)v[i]); if (!k) continue;
! 642: p1 = basistoalg(bnf, (GEN)gen[i]);
! 643: if (k > 0) xp = gmul(xp, gpuigs(p1, k));
! 644: else xm = gmul(xm, gpuigs(p1,-k));
! 645: }
! 646: if (xp != gun) x = gmul(x,xp);
! 647: if (xm != gun) x = gdiv(x,xm);
! 648: p1 = isunit(bnf,x);
! 649: if (lg(p1)==1) { avma = av; return cgetg(1,t_COL); }
! 650: tetpil=avma; return gerepile(av,tetpil,concat(p1,v));
! 651: }
! 652:
! 653: static void
! 654: vecconcat(GEN bnf,GEN relnf,GEN vec,GEN *prod,GEN *S1,GEN *S2)
! 655: {
! 656: long i;
! 657:
! 658: for (i=1; i<lg(vec); i++)
! 659: if (signe(resii(*prod,(GEN)vec[i])))
! 660: {
! 661: *prod=mulii(*prod,(GEN)vec[i]);
! 662: *S1=concatsp(*S1,primedec(bnf,(GEN)vec[i]));
! 663: *S2=concatsp(*S2,primedec(relnf,(GEN)vec[i]));
! 664: }
! 665: }
! 666:
! 667: /* bnf est le corps de base (buchinitfu).
! 668: * ext definit l'extension relative:
! 669: * ext[1] est une equation relative du corps,
! 670: * telle qu'une de ses racines engendre le corps sur Q.
! 671: * ext[2] exprime le generateur (y) du corps de base,
! 672: * en fonction de la racine (x) de ext[1],
! 673: * ext[3] est le buchinitfu (sur Q) de l'extension.
! 674:
! 675: * si flag=0 c'est qu'on sait a l'avance que l'extension est galoisienne,
! 676: * et dans ce cas la reponse est exacte.
! 677: * si flag>0 alors on ajoue dans S tous les ideaux qui divisent p<=flag.
! 678: * si flag<0 alors on ajoute dans S tous les ideaux qui divisent -flag.
! 679:
! 680: * la reponse est un vecteur v a 2 composantes telles que
! 681: * x=N(v[1])*v[2].
! 682: * x est une norme ssi v[2]=1.
! 683: */
! 684: GEN
! 685: rnfisnorm(GEN bnf,GEN ext,GEN x,long flag,long PREC)
! 686: {
! 687: long lgsunitrelnf,i,lbot, ltop = avma;
! 688: GEN relnf,aux,vec,tors,xnf,H,Y,M,A,suni,sunitrelnf,sunitnormnf,prod;
! 689: GEN res = cgetg(3,t_VEC), S1,S2;
! 690:
! 691: if (typ(ext)!=t_VEC || lg(ext)!=4) err (typeer,"bnfisnorm");
! 692: bnf = checkbnf(bnf); relnf = (GEN)ext[3];
! 693: if (gcmp0(x) || gcmp1(x) || (gcmp_1(x) && (degree((GEN)ext[1])&1)))
! 694: {
! 695: res[1]=lcopy(x); res[2]=un; return res;
! 696: }
! 697:
! 698: /* construction de l'ensemble S des ideaux
! 699: qui interviennent dans les solutions */
! 700:
! 701: prod=gun; S1=S2=cgetg(1,t_VEC);
! 702: if (!gcmp1(gmael3(relnf,8,1,1)))
! 703: {
! 704: GEN genclass=gmael3(relnf,8,1,3);
! 705: vec=cgetg(1,t_VEC);
! 706: for(i=1;i<lg(genclass);i++)
! 707: if (!gcmp1(ggcd(gmael4(relnf,8,1,2,i), stoi(degree((GEN)ext[1])))))
! 708: vec=concatsp(vec,(GEN)factor(gmael3(genclass,i,1,1))[1]);
! 709: vecconcat(bnf,relnf,vec,&prod,&S1,&S2);
! 710: }
! 711:
! 712: if (flag>1)
! 713: {
! 714: for (i=2; i<=flag; i++)
! 715: if (isprime(stoi(i)) && signe(resis(prod,i)))
! 716: {
! 717: prod=mulis(prod,i);
! 718: S1=concatsp(S1,primedec(bnf,stoi(i)));
! 719: S2=concatsp(S2,primedec(relnf,stoi(i)));
! 720: }
! 721: }
! 722: else if (flag<0)
! 723: vecconcat(bnf,relnf,(GEN)factor(stoi(-flag))[1],&prod,&S1,&S2);
! 724:
! 725: if (flag)
! 726: {
! 727: GEN normdiscrel=divii(gmael(relnf,7,3),
! 728: gpuigs(gmael(bnf,7,3),lg(ext[1])-3));
! 729: vecconcat(bnf,relnf,(GEN) factor(absi(normdiscrel))[1],
! 730: &prod,&S1,&S2);
! 731: }
! 732: vec=(GEN) idealfactor(bnf,x)[1]; aux=cgetg(2,t_VEC);
! 733: for (i=1; i<lg(vec); i++)
! 734: if (signe(resii(prod,gmael(vec,i,1))))
! 735: {
! 736: aux[1]=vec[i]; S1=concatsp(S1,aux);
! 737: }
! 738: xnf=lift(x);
! 739: xnf=gsubst(xnf,varn(xnf),(GEN)ext[2]);
! 740: vec=(GEN) idealfactor(relnf,xnf)[1];
! 741: for (i=1; i<lg(vec); i++)
! 742: if (signe(resii(prod,gmael(vec,i,1))))
! 743: {
! 744: aux[1]=vec[i]; S2=concatsp(S2,aux);
! 745: }
! 746:
! 747: res[1]=un; res[2]=(long)x;
! 748: tors=cgetg(2,t_VEC); tors[1]=mael3(relnf,8,4,2);
! 749:
! 750: /* calcul sur les S-unites */
! 751:
! 752: suni=bnfsunit(bnf,S1,PREC);
! 753: A=lift(bnfissunit(bnf,suni,x));
! 754: sunitrelnf=(GEN) bnfsunit(relnf,S2,PREC)[1];
! 755: if (lg(sunitrelnf)>1)
! 756: {
! 757: sunitrelnf=lift(basistoalg(relnf,sunitrelnf));
! 758: sunitrelnf=concatsp(tors,sunitrelnf);
! 759: }
! 760: else sunitrelnf=tors;
! 761: aux=(GEN)relnf[8];
! 762: if (lg(aux)>=6) aux=(GEN)aux[5];
! 763: else
! 764: {
! 765: aux=buchfu(relnf);
! 766: if(gcmp0((GEN)aux[2]))
! 767: err(talker,"please increase precision to have units in bnfisnorm");
! 768: aux=(GEN)aux[1];
! 769: }
! 770: if (lg(aux)>1)
! 771: sunitrelnf=concatsp(gtrans(aux),sunitrelnf);
! 772: lgsunitrelnf=lg(sunitrelnf);
! 773: M=cgetg(lgsunitrelnf+1,t_MAT);
! 774: sunitnormnf=cgetg(lgsunitrelnf,t_VEC);
! 775: for (i=1; i<lgsunitrelnf; i++)
! 776: {
! 777: sunitnormnf[i]=lnorm(gmodulcp((GEN) sunitrelnf[i],(GEN)ext[1]));
! 778: M[i]=llift(bnfissunit(bnf,suni,(GEN) sunitnormnf[i]));
! 779: }
! 780: M[lgsunitrelnf]=lgetg(lg(A),t_COL);
! 781: for (i=1; i<lg(A); i++) mael(M,lgsunitrelnf,i)=zero;
! 782: mael(M,lgsunitrelnf,lg(mael(bnf,7,6))-1)=mael3(bnf,8,4,1);
! 783: H=hnfall(M); Y=inverseimage(gmul(M,(GEN) H[2]),A);
! 784: Y=gmul((GEN) H[2],Y);
! 785: for (aux=(GEN)res[1],i=1; i<lgsunitrelnf; i++)
! 786: aux=gmul(aux,gpuigs(gmodulcp((GEN) sunitrelnf[i],(GEN)ext[1]),
! 787: itos(gfloor((GEN)Y[i]))));
! 788: res[1]=(long)aux;
! 789: res[2]=ldiv(x,gnorm(gmodulcp(lift(aux),(GEN)ext[1])));
! 790:
! 791: lbot=avma; return gerepile(ltop,lbot,gcopy(res));
! 792: }
! 793:
! 794: GEN
! 795: bnfisnorm(GEN bnf,GEN x,long flag,long PREC)
! 796: {
! 797: long ltop = avma, lbot;
! 798: GEN ext = cgetg(4,t_VEC);
! 799:
! 800: bnf = checkbnf(bnf);
! 801: ext[1] = mael(bnf,7,1);
! 802: ext[2] = zero;
! 803: ext[3] = (long) bnf;
! 804: bnf = buchinitfu(polx[MAXVARN],NULL,NULL,0); lbot = avma;
! 805: return gerepile(ltop,lbot,rnfisnorm(bnf,ext,x,flag,PREC));
! 806: }
! 807:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to buch4.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