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