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