Annotation of OpenXM_contrib/pari-2.2/src/basemath/base2.c, Revision 1.2
1.2 ! noro 1: /* $Id: base2.c,v 1.181 2002/09/11 02:31:22 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: /* MAXIMAL ORDERS */
19: /* */
20: /*******************************************************************/
21: #include "pari.h"
1.2 ! noro 22: #include "parinf.h"
1.1 noro 23:
1.2 ! noro 24: #define RXQX_rem(x,y,T) RXQX_divrem((x),(y),(T),ONLY_REM)
! 25: #define FpX_rem(x,y,p) FpX_divres((x),(y),(p),ONLY_REM)
! 26: extern GEN addone_aux2(GEN x, GEN y);
! 27: extern GEN addshiftw(GEN x, GEN y, long d);
! 28: extern GEN gmul_mat_smallvec(GEN x, GEN y);
! 29: extern GEN hnf_invimage(GEN A, GEN b);
! 30: extern GEN norm_by_embed(long r1, GEN x);
! 31: extern GEN ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound);
! 32: extern GEN polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N);
! 33: extern GEN FqX_gcd(GEN P, GEN Q, GEN T, GEN p);
! 34: extern GEN FqX_factor(GEN x, GEN T, GEN p);
! 35: extern GEN DDF2(GEN x, long hint);
! 36: extern GEN eltabstorel(GEN x, GEN T, GEN pol, GEN k);
! 37: extern GEN element_powid_mod_p(GEN nf, long I, GEN n, GEN p);
! 38: extern GEN FpVQX_red(GEN z, GEN T, GEN p);
! 39: extern GEN Fp_factor_irred(GEN P,GEN l, GEN Q);
! 40: extern GEN RXQX_divrem(GEN x, GEN y, GEN T, GEN *pr);
! 41: extern GEN RXQX_mul(GEN x, GEN y, GEN T);
! 42: extern GEN ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LPRS);
! 43: extern GEN _ei(long n, long i);
! 44: extern GEN col_to_ff(GEN x, long v);
! 45: extern GEN element_mulidid(GEN nf, long i, long j);
! 46: extern GEN eltmulid_get_table(GEN nf, long i);
! 47: extern GEN idealaddtoone_i(GEN nf, GEN x, GEN y);
1.1 noro 48: extern GEN mat_to_vecpol(GEN x, long v);
1.2 ! noro 49: extern GEN merge_factor_i(GEN f, GEN g);
! 50: extern GEN mulmat_pol(GEN A, GEN x);
! 51: extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);
! 52: extern GEN pidealprimeinv(GEN nf, GEN x);
1.1 noro 53: extern GEN pol_to_monic(GEN pol, GEN *lead);
54: extern GEN quicktrace(GEN x, GEN sym);
1.2 ! noro 55: extern GEN sqr_by_tab(GEN tab, GEN x);
! 56: extern GEN to_polmod(GEN x, GEN mod);
! 57: extern GEN unnf_minus_x(GEN x);
! 58: extern int isrational(GEN x);
! 59: extern long int_elt_val(GEN nf, GEN x, GEN p, GEN bp, GEN *t);
! 60:
! 61: /* FIXME: backward compatibility. Should use the proper nf_* equivalents */
! 62: #define compat_PARTIAL 1
! 63: #define compat_ROUND2 2
1.1 noro 64:
65: static void
1.2 ! noro 66: allbase_check_args(GEN f, long flag, GEN *dx, GEN *ptw)
1.1 noro 67: {
1.2 ! noro 68: GEN w = *ptw;
! 69:
! 70: if (DEBUGLEVEL) (void)timer2();
1.1 noro 71: if (typ(f)!=t_POL) err(notpoler,"allbase");
72: if (degpol(f) <= 0) err(constpoler,"allbase");
1.2 ! noro 73:
! 74: *dx = w? factorback(w, NULL): ZX_disc(f);
! 75: if (!signe(*dx)) err(talker,"reducible polynomial in allbase");
! 76: if (!w) *ptw = auxdecomp(absi(*dx), (flag & nf_PARTIALFACT)? 0: 1);
1.1 noro 77: if (DEBUGLEVEL) msgtimer("disc. factorisation");
78: }
79:
80: /*******************************************************************/
81: /* */
82: /* ROUND 2 */
83: /* */
84: /*******************************************************************/
85: /* companion matrix of unitary polynomial x */
86: static GEN
87: companion(GEN x) /* cf assmat */
88: {
89: long i,j,l;
90: GEN y;
91:
92: l=degpol(x)+1; y=cgetg(l,t_MAT);
93: for (j=1; j<l; j++)
94: {
95: y[j] = lgetg(l,t_COL);
96: for (i=1; i<l-1; i++)
97: coeff(y,i,j)=(i+1==j)? un: zero;
98: coeff(y,i,j) = lneg((GEN)x[j+1]);
99: }
100: return y;
101: }
102:
103: /* assume x, y are square integer matrices of same dim. Multiply them */
104: static GEN
105: mulmati(GEN x, GEN y)
106: {
1.2 ! noro 107: gpmem_t av;
! 108: long n = lg(x),i,j,k;
1.1 noro 109: GEN z = cgetg(n,t_MAT),p1,p2;
110:
111: for (j=1; j<n; j++)
112: {
113: z[j] = lgetg(n,t_COL);
114: for (i=1; i<n; i++)
115: {
116: p1=gzero; av=avma;
117: for (k=1; k<n; k++)
118: {
119: p2=mulii(gcoeff(x,i,k),gcoeff(y,k,j));
120: if (p2 != gzero) p1=addii(p1,p2);
121: }
122: coeff(z,i,j)=lpileupto(av,p1);
123: }
124: }
125: return z;
126: }
127:
128: static GEN
1.2 ! noro 129: _mulmati(void *data /*ignored*/, GEN x, GEN y) {
! 130: (void)data; return mulmati(x,y);
! 131: }
! 132: static GEN
! 133: _sqrmati(void *data /*ignored*/, GEN x) {
! 134: (void)data; return mulmati(x,x);
! 135: }
! 136:
! 137: static GEN
! 138: powmati(GEN x, GEN n)
1.1 noro 139: {
1.2 ! noro 140: gpmem_t av = avma;
! 141: GEN y = leftright_pow(x, n, NULL, &_sqrmati, &_mulmati);
1.1 noro 142: return gerepileupto(av,y);
143: }
144:
145: static GEN
146: rtran(GEN v, GEN w, GEN q)
147: {
1.2 ! noro 148: gpmem_t av,tetpil;
1.1 noro 149: GEN p1;
150:
151: if (signe(q))
152: {
153: av=avma; p1=gneg(gmul(q,w)); tetpil=avma;
154: return gerepile(av,tetpil,gadd(v,p1));
155: }
156: return v;
157: }
158:
159: /* return (v - qw) mod m (only compute entries k0,..,n)
160: * v and w are expected to have entries smaller than m */
161: static GEN
1.2 ! noro 162: mtran(GEN v, GEN w, GEN q, GEN m, GEN mo2, long k0)
1.1 noro 163: {
1.2 ! noro 164: long k;
1.1 noro 165: GEN p1;
166:
167: if (signe(q))
1.2 ! noro 168: for (k=lg(v)-1; k >= k0; k--)
1.1 noro 169: {
1.2 ! noro 170: gpmem_t av = avma;
1.1 noro 171: p1 = subii((GEN)v[k], mulii(q,(GEN)w[k]));
1.2 ! noro 172: p1 = centermodii(p1, m, mo2);
! 173: v[k] = lpileuptoint(av, p1);
1.1 noro 174: }
175: return v;
176: }
177:
178: /* entries of v and w are C small integers */
179: static GEN
180: mtran_long(GEN v, GEN w, long q, long m, long k0)
181: {
182: long k, p1;
183:
184: if (q)
185: {
186: for (k=lg(v)-1; k>= k0; k--)
187: {
188: p1 = v[k] - q * w[k];
189: v[k] = p1 % m;
190: }
191: }
192: return v;
193: }
194:
195: /* coeffs of a are C-long integers */
196: static void
197: rowred_long(GEN a, long rmod)
198: {
199: long q,j,k,pro, c = lg(a), r = lg(a[1]);
200:
201: for (j=1; j<r; j++)
202: {
203: for (k=j+1; k<c; k++)
204: while (coeff(a,j,k))
205: {
206: q = coeff(a,j,j) / coeff(a,j,k);
207: pro=(long)mtran_long((GEN)a[j],(GEN)a[k],q,rmod, j);
208: a[j]=a[k]; a[k]=pro;
209: }
210: if (coeff(a,j,j) < 0)
211: for (k=j; k<r; k++) coeff(a,k,j)=-coeff(a,k,j);
212: for (k=1; k<j; k++)
213: {
214: q = coeff(a,j,k) / coeff(a,j,j);
215: a[k]=(long)mtran_long((GEN)a[k],(GEN)a[j],q,rmod, k);
216: }
217: }
218: /* don't update the 0s in the last columns */
219: for (j=1; j<r; j++)
220: for (k=1; k<r; k++) coeff(a,j,k) = lstoi(coeff(a,j,k));
221: }
222:
223: static void
224: rowred(GEN a, GEN rmod)
225: {
226: long j,k,pro, c = lg(a), r = lg(a[1]);
1.2 ! noro 227: gpmem_t av=avma, lim=stack_lim(av,1);
! 228: GEN q, rmodo2 = shifti(rmod,-1);
1.1 noro 229:
230: for (j=1; j<r; j++)
231: {
232: for (k=j+1; k<c; k++)
233: while (signe(gcoeff(a,j,k)))
234: {
1.2 ! noro 235: q=diviiround(gcoeff(a,j,j),gcoeff(a,j,k));
! 236: pro=(long)mtran((GEN)a[j],(GEN)a[k],q,rmod,rmodo2, j);
1.1 noro 237: a[j]=a[k]; a[k]=pro;
238: }
239: if (signe(gcoeff(a,j,j)) < 0)
240: for (k=j; k<r; k++) coeff(a,k,j)=lnegi(gcoeff(a,k,j));
241: for (k=1; k<j; k++)
242: {
1.2 ! noro 243: q=diviiround(gcoeff(a,j,k),gcoeff(a,j,j));
! 244: a[k]=(long)mtran((GEN)a[k],(GEN)a[j],q,rmod,rmodo2, k);
1.1 noro 245: }
246: if (low_stack(lim, stack_lim(av,1)))
247: {
248: long j1,k1;
249: GEN p1 = a;
250: if(DEBUGMEM>1) err(warnmem,"rowred j=%ld", j);
251: p1 = gerepilecopy(av,a);
252: for (j1=1; j1<r; j1++)
253: for (k1=1; k1<c; k1++) coeff(a,j1,k1) = coeff(p1,j1,k1);
254: }
255: }
256: }
257:
258: /* Calcule d/x ou d est entier et x matrice triangulaire inferieure
1.2 ! noro 259: * entiere dont les coeff diagonaux divisent d (resultat entier). */
1.1 noro 260: static GEN
1.2 ! noro 261: matinv(GEN x, GEN d)
1.1 noro 262: {
1.2 ! noro 263: gpmem_t av,av1;
! 264: long i,j,k, n = lg(x[1]); /* Warning: lg(x) from ordmax is bogus */
1.1 noro 265: GEN y,h;
266:
1.2 ! noro 267: y = idmat(n-1);
! 268: for (i=1; i<n; i++)
! 269: coeff(y,i,i) = (long)diviiexact(d,gcoeff(x,i,i));
1.1 noro 270: av=avma;
1.2 ! noro 271: for (i=2; i<n; i++)
1.1 noro 272: for (j=i-1; j; j--)
273: {
274: for (h=gzero,k=j+1; k<=i; k++)
275: {
276: GEN p1 = mulii(gcoeff(y,i,k),gcoeff(x,k,j));
277: if (p1 != gzero) h=addii(h,p1);
278: }
279: setsigne(h,-signe(h)); av1=avma;
280: coeff(y,i,j) = lpile(av,av1,divii(h,gcoeff(x,j,j)));
281: av = avma;
282: }
283: return y;
284: }
285:
286: static GEN
287: ordmax(GEN *cf, GEN p, long epsilon, GEN *ptdelta)
288: {
1.2 ! noro 289: long sp,i,n=lg(cf)-1;
! 290: gpmem_t av=avma, av2,limit;
! 291: GEN T,T2,Tn,m,v,delta,hard_case_exponent, *w;
1.1 noro 292: const GEN pp = sqri(p);
1.2 ! noro 293: const GEN ppo2 = shifti(pp,-1);
1.1 noro 294: const long pps = (2*expi(pp)+2<BITS_IN_LONG)? pp[2]: 0;
295:
296: if (cmpis(p,n) > 0)
297: {
1.2 ! noro 298: hard_case_exponent = NULL;
1.1 noro 299: sp = 0; /* gcc -Wall */
300: }
301: else
302: {
303: long k;
304: k = sp = itos(p);
305: i=1; while (k < n) { k *= sp; i++; }
1.2 ! noro 306: hard_case_exponent = stoi(i);
1.1 noro 307: }
308: T=cgetg(n+1,t_MAT); for (i=1; i<=n; i++) T[i]=lgetg(n+1,t_COL);
309: T2=cgetg(2*n+1,t_MAT); for (i=1; i<=2*n; i++) T2[i]=lgetg(n+1,t_COL);
310: Tn=cgetg(n*n+1,t_MAT); for (i=1; i<=n*n; i++) Tn[i]=lgetg(n+1,t_COL);
311: v = new_chunk(n+1);
1.2 ! noro 312: w = (GEN*)new_chunk(n+1);
1.1 noro 313:
314: av2 = avma; limit = stack_lim(av2,1);
315: delta=gun; m=idmat(n);
316:
317: for(;;)
318: {
1.2 ! noro 319: long j, k, h;
! 320: gpmem_t av0 = avma;
1.1 noro 321: GEN t,b,jp,hh,index,p1, dd = sqri(delta), ppdd = mulii(dd,pp);
1.2 ! noro 322: GEN ppddo2 = shifti(ppdd,-1);
1.1 noro 323:
324: if (DEBUGLEVEL > 3)
325: fprintferr("ROUND2: epsilon = %ld\tavma = %ld\n",epsilon,avma);
326:
1.2 ! noro 327: b=matinv(m,delta);
1.1 noro 328: for (i=1; i<=n; i++)
329: {
330: for (j=1; j<=n; j++)
331: for (k=1; k<=n; k++)
332: {
333: p1 = j==k? gcoeff(m,i,1): gzero;
334: for (h=2; h<=n; h++)
335: {
336: GEN p2 = mulii(gcoeff(m,i,h),gcoeff(cf[h],j,k));
337: if (p2!=gzero) p1 = addii(p1,p2);
338: }
1.2 ! noro 339: coeff(T,j,k) = (long)centermodii(p1, ppdd, ppddo2);
1.1 noro 340: }
341: p1 = mulmati(m, mulmati(T,b));
342: for (j=1; j<=n; j++)
343: for (k=1; k<=n; k++)
1.2 ! noro 344: coeff(p1,j,k)=(long)centermodii(divii(gcoeff(p1,j,k),dd),pp,ppo2);
1.1 noro 345: w[i] = p1;
346: }
347:
348: if (hard_case_exponent)
349: {
350: for (j=1; j<=n; j++)
351: {
352: for (i=1; i<=n; i++) coeff(T,i,j) = coeff(w[j],1,i);
353: /* ici la boucle en k calcule la puissance p mod p de w[j] */
354: for (k=1; k<sp; k++)
355: {
356: for (i=1; i<=n; i++)
357: {
358: p1 = gzero;
359: for (h=1; h<=n; h++)
360: {
361: GEN p2=mulii(gcoeff(T,h,j),gcoeff(w[j],h,i));
362: if (p2!=gzero) p1 = addii(p1,p2);
363: }
364: v[i] = lmodii(p1, p);
365: }
366: for (i=1; i<=n; i++) coeff(T,i,j)=v[i];
367: }
368: }
369: t = powmati(T, hard_case_exponent);
370: }
371: else
372: {
373: for (i=1; i<=n; i++)
374: for (j=1; j<=n; j++)
375: {
1.2 ! noro 376: gpmem_t av1 = avma;
1.1 noro 377: p1 = gzero;
378: for (k=1; k<=n; k++)
379: for (h=1; h<=n; h++)
380: {
381: const GEN r=modii(gcoeff(w[i],k,h),p);
382: const GEN s=modii(gcoeff(w[j],h,k),p);
1.2 ! noro 383: const GEN p2 = mulii(r,s);
1.1 noro 384: if (p2!=gzero) p1 = addii(p1,p2);
385: }
386: coeff(T,i,j) = lpileupto(av1,p1);
387: }
388: t = T;
389: }
390:
391: if (pps)
392: {
393: long ps = p[2];
394: for (i=1; i<=n; i++)
395: for (j=1; j<=n; j++)
396: {
397: coeff(T2,j,i)=(i==j)? ps: 0;
398: coeff(T2,j,n+i)=smodis(gcoeff(t,i,j),ps);
399: }
1.2 ! noro 400: rowred_long(T2,pps);
1.1 noro 401: }
402: else
403: {
404: for (i=1; i<=n; i++)
405: for (j=1; j<=n; j++)
406: {
407: coeff(T2,j,i)=(i==j)? (long)p: zero;
408: coeff(T2,j,n+i)=lmodii(gcoeff(t,i,j),p);
409: }
410: rowred(T2,pp);
411: }
1.2 ! noro 412: jp=matinv(T2,p);
1.1 noro 413: if (pps)
414: {
415: for (k=1; k<=n; k++)
416: {
1.2 ! noro 417: gpmem_t av1=avma;
1.1 noro 418: t = mulmati(mulmati(jp,w[k]), T2);
419: for (h=i=1; i<=n; i++)
420: for (j=1; j<=n; j++)
421: { coeff(Tn,k,h) = itos(divii(gcoeff(t,i,j), p)) % pps; h++; }
422: avma=av1;
423: }
424: avma = av0;
425: rowred_long(Tn,pps);
426: }
427: else
428: {
429: for (k=1; k<=n; k++)
430: {
431: t = mulmati(mulmati(jp,w[k]), T2);
432: for (h=i=1; i<=n; i++)
433: for (j=1; j<=n; j++)
434: { coeff(Tn,k,h) = ldivii(gcoeff(t,i,j), p); h++; }
435: }
436: rowred(Tn,pp);
437: }
438: for (index=gun,i=1; i<=n; i++)
439: index = mulii(index,gcoeff(Tn,i,i));
440: if (gcmp1(index)) break;
441:
1.2 ! noro 442: m = mulmati(matinv(Tn,index), m);
1.1 noro 443: hh = delta = mulii(index,delta);
444: for (i=1; i<=n; i++)
445: for (j=1; j<=n; j++)
446: hh = mppgcd(gcoeff(m,i,j),hh);
447: if (!is_pm1(hh))
448: {
449: m = gdiv(m,hh);
450: delta = divii(delta,hh);
451: }
452: epsilon -= 2 * ggval(index,p);
453: if (epsilon < 2) break;
454: if (low_stack(limit,stack_lim(av2,1)))
455: {
456: GEN *gptr[3]; gptr[0]=&m; gptr[1]=δ
457: if(DEBUGMEM>1) err(warnmem,"ordmax");
458: gerepilemany(av2, gptr,2);
459: }
460: }
461: {
462: GEN *gptr[2]; gptr[0]=&m; gptr[1]=δ
463: gerepilemany(av,gptr,2);
464: }
465: *ptdelta=delta; return m;
466: }
467:
468: /* Input:
469: * x normalized integral polynomial of degree n, defining K=Q(theta).
470: *
471: * code 0, 1 or (long)p if we want base, smallbase ou factoredbase (resp.).
472: * y is GEN *, which will receive the discriminant of K.
473: *
474: * Output
475: * 1) A t_COL whose n components are rationnal polynomials (with degree
476: * 0,1...n-1) : integral basis for K (putting x=theta).
477: * Rem: common denominator is in da.
478: *
479: * 2) discriminant of K (in *y).
480: */
1.2 ! noro 481: static GEN
! 482: allbase2(GEN f, int flag, GEN *dx, GEN *dK, GEN *ptw)
1.1 noro 483: {
1.2 ! noro 484: GEN w,w1,w2,a,pro,at,bt,b,da,db,q, *cf,*gptr[2];
! 485: gpmem_t av=avma,tetpil;
! 486: long n,h,j,i,k,r,s,t,v,mf;
! 487:
! 488: w = ptw? *ptw: NULL;
! 489: allbase_check_args(f,flag,dx, &w);
! 490: w1 = (GEN)w[1];
! 491: w2 = (GEN)w[2];
1.1 noro 492: v = varn(f); n = degpol(f); h = lg(w1)-1;
493: cf = (GEN*)cgetg(n+1,t_VEC);
494: cf[2]=companion(f);
495: for (i=3; i<=n; i++) cf[i]=mulmati(cf[2],cf[i-1]);
496:
497: a=idmat(n); da=gun;
498: for (i=1; i<=h; i++)
499: {
1.2 ! noro 500: gpmem_t av1 = avma;
1.1 noro 501: mf=itos((GEN)w2[i]); if (mf==1) continue;
502: if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
503:
504: b=ordmax(cf,(GEN)w1[i],mf,&db);
505: a=gmul(db,a); b=gmul(da,b);
506: da=mulii(db,da);
507: at=gtrans(a); bt=gtrans(b);
508: for (r=n; r; r--)
509: for (s=r; s; s--)
510: while (signe(gcoeff(bt,s,r)))
511: {
1.2 ! noro 512: q=diviiround(gcoeff(at,s,s),gcoeff(bt,s,r));
1.1 noro 513: pro=rtran((GEN)at[s],(GEN)bt[r],q);
514: for (t=s-1; t; t--)
515: {
1.2 ! noro 516: q=diviiround(gcoeff(at,t,s),gcoeff(at,t,t));
1.1 noro 517: pro=rtran(pro,(GEN)at[t],q);
518: }
519: at[s]=bt[r]; bt[r]=(long)pro;
520: }
521: for (j=n; j; j--)
522: {
523: for (k=1; k<j; k++)
524: {
525: while (signe(gcoeff(at,j,k)))
526: {
1.2 ! noro 527: q=diviiround(gcoeff(at,j,j),gcoeff(at,j,k));
1.1 noro 528: pro=rtran((GEN)at[j],(GEN)at[k],q);
529: at[j]=at[k]; at[k]=(long)pro;
530: }
531: }
532: if (signe(gcoeff(at,j,j))<0)
533: for (k=1; k<=j; k++) coeff(at,k,j)=lnegi(gcoeff(at,k,j));
534: for (k=j+1; k<=n; k++)
535: {
1.2 ! noro 536: q=diviiround(gcoeff(at,j,k),gcoeff(at,j,j));
1.1 noro 537: at[k]=(long)rtran((GEN)at[k],(GEN)at[j],q);
538: }
539: }
540: for (j=2; j<=n; j++)
541: if (egalii(gcoeff(at,j,j), gcoeff(at,j-1,j-1)))
542: {
543: coeff(at,1,j)=zero;
544: for (k=2; k<=j; k++) coeff(at,k,j)=coeff(at,k-1,j-1);
545: }
546: tetpil=avma; a=gtrans(at);
547: {
1.2 ! noro 548: GEN *gptr[2];
1.1 noro 549: da = icopy(da); gptr[0]=&a; gptr[1]=&da;
550: gerepilemanysp(av1,tetpil,gptr,2);
551: }
552: }
1.2 ! noro 553: *dK = *dx;
1.1 noro 554: for (j=1; j<=n; j++)
1.2 ! noro 555: *dK = diviiexact(mulii(*dK,sqri(gcoeff(a,j,j))), sqri(da));
! 556: tetpil=avma; *dK = icopy(*dK);
1.1 noro 557: at=cgetg(n+1,t_VEC); v=varn(f);
558: for (k=1; k<=n; k++)
559: {
560: q=cgetg(k+2,t_POL); at[k]=(long)q;
561: q[1] = evalsigne(1) | evallgef(2+k) | evalvarn(v);
1.2 ! noro 562: for (j=1; j<=k; j++) q[j+1] = ldiv(gcoeff(a,k,j),da);
1.1 noro 563: }
1.2 ! noro 564: gptr[0] = &at; gptr[1] = dK;
1.1 noro 565: gerepilemanysp(av,tetpil,gptr,2);
566: return at;
567: }
568:
569: GEN
1.2 ! noro 570: base2(GEN x, GEN *pdK) { return nfbasis(x, pdK, compat_ROUND2, NULL); }
1.1 noro 571:
572: GEN
1.2 ! noro 573: discf2(GEN x) { return nfdiscf0(x, compat_ROUND2, NULL); }
1.1 noro 574:
575: /*******************************************************************/
576: /* */
577: /* ROUND 4 */
578: /* */
579: /*******************************************************************/
580:
581: GEN nilord(GEN p,GEN fx,long mf,GEN gx,long flag);
582: GEN Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long r);
583: static GEN dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U);
584: static GEN maxord(GEN p,GEN f,long mf);
585: static GEN testb2(GEN p,GEN fa,long Fa,GEN theta,GEN pmf,long Ft,GEN ns);
586: static GEN testc2(GEN p,GEN fa,GEN pmr,GEN pmf,GEN alph2,
587: long Ea,GEN thet2,long Et,GEN ns);
588:
589: static int
590: fnz(GEN x,long j)
591: {
592: long i;
593: for (i=1; i<j; i++)
594: if (signe(x[i])) return 0;
595: return 1;
596: }
597:
1.2 ! noro 598: /* return list u[i], 2 by 2 coprime with the same prime divisors as ab */
! 599: static GEN
! 600: get_coprimes(GEN a, GEN b)
1.1 noro 601: {
1.2 ! noro 602: long i, k = 1;
! 603: GEN u = cgetg(3, t_VEC);
! 604: u[1] = (long)a;
! 605: u[2] = (long)b;
! 606: /* u1,..., uk 2 by 2 coprime */
! 607: while (k+1 < lg(u))
! 608: {
! 609: GEN d, c = (GEN)u[k+1];
! 610: if (is_pm1(c)) { k++; continue; }
! 611: for (i=1; i<=k; i++)
! 612: {
! 613: if (is_pm1(u[i])) continue;
! 614: d = mppgcd(c, (GEN)u[i]);
! 615: if (d == gun) continue;
! 616: c = diviiexact(c, d);
! 617: u[i] = (long)diviiexact((GEN)u[i], d);
! 618: u = concatsp(u, d);
! 619: }
! 620: u[++k] = (long)c;
! 621: }
! 622: for (i = k = 1; i < lg(u); i++)
! 623: if (!is_pm1(u[i])) u[k++] = u[i];
! 624: setlg(u, k); return u;
! 625: }
! 626:
! 627: /* return integer basis. Set dK = disc(K), dx = disc(f), w (possibly partial)
! 628: * factorization of dK. *ptw can be set by the caller, in which case it is
! 629: * taken to be the factorization of disc(f), then overwritten
! 630: * [no consistency check] */
! 631: GEN
! 632: allbase(GEN f, int flag, GEN *dx, GEN *dK, GEN *index, GEN *ptw)
! 633: {
! 634: GEN w, w1, w2, a, da, p1, ordmax;
! 635: long n, lw, i, j, k, l;
! 636:
! 637: if (flag & nf_ROUND2) return allbase2(f,flag,dx,dK,ptw);
! 638: w = ptw? *ptw: NULL;
! 639: allbase_check_args(f, flag, dx, &w);
! 640: w1 = (GEN)w[1];
! 641: w2 = (GEN)w[2];
! 642: n = degpol(f); lw = lg(w1);
! 643: ordmax = cgetg(1, t_VEC);
! 644: /* "complete" factorization first */
! 645: for (i=1; i<lw; i++)
! 646: {
! 647: long mf = itos((GEN)w2[i]);
! 648: if (mf == 1) { ordmax = concatsp(ordmax, gun); continue; }
! 649:
! 650: CATCH(invmoder) { /* caught false prime, update factorization */
! 651: GEN x = (GEN)global_err_data;
! 652: GEN p = mppgcd((GEN)x[1], (GEN)x[2]);
! 653: GEN N, u;
! 654: if (DEBUGLEVEL) err(warner,"impossible inverse: %Z", x);
! 655:
! 656: u = get_coprimes(p, diviiexact((GEN)x[1],p));
! 657: l = lg(u);
! 658: /* no small factors, but often a prime power */
! 659: for (k = 1; k < l; k++) u[k] = coeff(auxdecomp((GEN)u[k], 2),1,1);
! 660:
! 661: w1[i] = u[1];
! 662: w1 = concatsp(w1, vecextract_i(u, 2, l-1));
! 663: N = *dx;
! 664: w2[i] = lstoi(pvaluation(N, (GEN)w1[i], &N));
! 665: k = lw;
! 666: lw = lg(w1);
! 667: for ( ; k < lw; k++) w2[k] = lstoi(pvaluation(N, (GEN)w1[k], &N));
! 668: } RETRY {
! 669: if (DEBUGLEVEL) fprintferr("Treating p^k = %Z^%ld\n",w1[i],mf);
! 670: ordmax = concatsp(ordmax, _vec( maxord((GEN)w1[i],f,mf) ));
! 671: } ENDCATCH;
! 672: }
1.1 noro 673:
674: a = NULL; /* gcc -Wall */
1.2 ! noro 675: da = NULL;
! 676: for (i=1; i<lw; i++)
1.1 noro 677: {
1.2 ! noro 678: GEN db, b = (GEN)ordmax[i];
! 679: if (b == gun) continue;
! 680: db = gun;
1.1 noro 681: for (j=1; j<=n; j++)
682: {
683: p1 = denom(gcoeff(b,j,j));
684: if (cmpii(p1,db) > 0) db = p1;
685: }
1.2 ! noro 686: if (db == gun) continue;
! 687:
! 688: /* db = denom(diag(b)), (da,db) = 1 */
! 689: b = Q_muli_to_int(b,db);
! 690: if (!da) { da = db; a = b; }
! 691: else
! 692: {
! 693: j=1; while (j<=n && fnz((GEN)a[j],j) && fnz((GEN)b[j],j)) j++;
! 694: b = gmul(da,b);
! 695: a = gmul(db,a); da = mulii(da,db);
! 696: k = j-1; p1 = cgetg(2*n-k+1,t_MAT);
! 697: for (j=1; j<=k; j++)
! 698: {
! 699: p1[j] = a[j];
! 700: coeff(p1,j,j) = lmppgcd(gcoeff(a,j,j),gcoeff(b,j,j));
! 701: }
! 702: for ( ; j<=n; j++) p1[j] = a[j];
! 703: for ( ; j<=2*n-k; j++) p1[j] = b[j+k-n];
! 704: a = hnfmodid(p1, da);
1.1 noro 705: }
1.2 ! noro 706: if (DEBUGLEVEL>5) fprintferr("Result for prime %Z is:\n%Z\n",w1[i],b);
1.1 noro 707: }
1.2 ! noro 708: *dK = *dx;
1.1 noro 709: if (da)
710: {
1.2 ! noro 711: *index = diviiexact(da, gcoeff(a,1,1));
! 712: for (j=2; j<=n; j++) *index = mulii(*index, diviiexact(da, gcoeff(a,j,j)));
! 713: *dK = diviiexact(*dK, sqri(*index));
1.1 noro 714: for (j=n-1; j; j--)
715: if (cmpis(gcoeff(a,j,j),2) > 0)
716: {
1.2 ! noro 717: p1 = shifti(gcoeff(a,j,j),-1);
1.1 noro 718: for (k=j+1; k<=n; k++)
719: if (cmpii(gcoeff(a,j,k),p1) > 0)
720: for (l=1; l<=j; l++)
1.2 ! noro 721: coeff(a,l,k) = lsubii(gcoeff(a,l,k),gcoeff(a,l,j));
1.1 noro 722: }
1.2 ! noro 723: a = gdiv(a, da);
1.1 noro 724: }
1.2 ! noro 725: else
1.1 noro 726: {
1.2 ! noro 727: *index = gun;
! 728: a = idmat(n);
1.1 noro 729: }
1.2 ! noro 730:
! 731: if (ptw)
1.1 noro 732: {
1.2 ! noro 733: long lfa = 1;
! 734: GEN W1, W2, D = *dK;
! 735:
! 736: w = cgetg(3,t_MAT);
! 737: W1 = cgetg(lw, t_COL); w[1] = (long)W1;
! 738: W2 = cgetg(lw, t_COL); w[2] = (long)W2;
! 739: for (j=1; j<lw; j++)
1.1 noro 740: {
1.2 ! noro 741: k = pvaluation(D, (GEN)w1[j], &D);
! 742: if (k) { W1[lfa] = w1[j]; W2[lfa] = lstoi(k); lfa++; }
1.1 noro 743: }
1.2 ! noro 744: setlg(W1, lfa);
! 745: setlg(W2, lfa); *ptw = w;
1.1 noro 746: }
1.2 ! noro 747: return mat_to_vecpol(a, varn(f));
1.1 noro 748: }
749:
750: static GEN
751: update_fact(GEN x, GEN f)
752: {
753: GEN e,q,d = ZX_disc(x), g = cgetg(3, t_MAT), p = (GEN)f[1];
754: long iq,i,k,l;
755: if (typ(f)!=t_MAT || lg(f)!=3)
756: err(talker,"not a factorisation in nfbasis");
757: l = lg(p);
758: q = cgetg(l,t_COL); g[1]=(long)q;
759: e = cgetg(l,t_COL); g[2]=(long)e; iq = 1;
760: for (i=1; i<l; i++)
761: {
762: k = pvaluation(d, (GEN)p[i], &d);
1.2 ! noro 763: if (k) { q[iq] = p[i]; e[iq] = lstoi(k); iq++; }
1.1 noro 764: }
765: setlg(q,iq); setlg(e,iq);
766: return merge_factor_i(decomp(d), g);
767: }
768:
769: static GEN
1.2 ! noro 770: unscale_vecpol(GEN v, GEN h)
1.1 noro 771: {
1.2 ! noro 772: long i, l;
! 773: GEN w;
! 774: if (!h) return v;
! 775: l = lg(v); w = cgetg(l, typ(v));
! 776: for (i=1; i<l; i++) w[i] = (long)unscale_pol((GEN)v[i], h);
! 777: return w;
! 778: }
1.1 noro 779:
1.2 ! noro 780: /* FIXME: have to deal with compatibility flags */
! 781: static void
! 782: _nfbasis(GEN x0, long flag, GEN fa, GEN *pbas, GEN *pdK)
! 783: {
! 784: GEN x, dx, dK, basis, lead, index;
! 785: long fl = 0;
1.1 noro 786:
1.2 ! noro 787: if (typ(x0)!=t_POL) err(typeer,"nfbasis");
! 788: if (!degpol(x0)) err(zeropoler,"nfbasis");
! 789: check_pol_int(x0, "nfbasis");
1.1 noro 790:
1.2 ! noro 791: x = pol_to_monic(x0, &lead);
! 792: if (fa && gcmp0(fa)) fa = NULL; /* compatibility. NULL is the proper arg */
! 793: if (fa && lead) fa = update_fact(x, fa);
! 794: if (flag & compat_PARTIAL) fl |= nf_PARTIALFACT;
! 795: if (flag & compat_ROUND2) fl |= nf_ROUND2;
! 796: basis = allbase(x, fl, &dx, &dK, &index, &fa);
! 797: if (pbas) *pbas = unscale_vecpol(basis, lead);
! 798: if (pdK) *pdK = dK;
1.1 noro 799: }
800:
801: GEN
1.2 ! noro 802: nfbasis(GEN x, GEN *pdK, long flag, GEN fa)
1.1 noro 803: {
1.2 ! noro 804: gpmem_t av = avma;
! 805: GEN bas; _nfbasis(x, flag, fa, &bas, pdK);
! 806: gerepileall(av, 2, &bas, pdK); return bas;
1.1 noro 807: }
808:
809: GEN
1.2 ! noro 810: nfbasis0(GEN x, long flag, GEN fa)
1.1 noro 811: {
1.2 ! noro 812: gpmem_t av = avma;
! 813: GEN bas; _nfbasis(x, flag, fa, &bas, NULL);
! 814: return gerepilecopy(av, bas);
1.1 noro 815: }
816:
817: GEN
1.2 ! noro 818: nfdiscf0(GEN x, long flag, GEN fa)
1.1 noro 819: {
1.2 ! noro 820: gpmem_t av = avma;
! 821: GEN dK; _nfbasis(x, flag, fa, NULL, &dK);
! 822: return gerepilecopy(av, dK);
1.1 noro 823: }
824:
825: GEN
1.2 ! noro 826: base(GEN x, GEN *pdK) { return nfbasis(x, pdK, 0, NULL); }
1.1 noro 827:
828: GEN
1.2 ! noro 829: smallbase(GEN x, GEN *pdK) { return nfbasis(x, pdK, compat_PARTIAL, NULL); }
1.1 noro 830:
831: GEN
1.2 ! noro 832: factoredbase(GEN x, GEN fa, GEN *pdK) { return nfbasis(x, pdK, 0, fa); }
1.1 noro 833:
834: GEN
1.2 ! noro 835: discf(GEN x) { return nfdiscf0(x, 0, NULL); }
1.1 noro 836:
837: GEN
1.2 ! noro 838: smalldiscf(GEN x) { return nfdiscf0(x, nf_PARTIALFACT, NULL); }
1.1 noro 839:
840: GEN
1.2 ! noro 841: factoreddiscf(GEN x, GEN fa) { return nfdiscf0(x, 0, fa); }
1.1 noro 842:
843:
844: /* return U if Z[alpha] is not maximal or 2*dU < m-1; else return NULL */
845: static GEN
846: dedek(GEN f, long mf, GEN p,GEN g)
847: {
848: GEN k,h;
849: long dk;
850:
851: h = FpX_div(f,g,p);
852: k = gdivexact(gadd(f, gneg_i(gmul(g,h))), p);
853: k = FpX_gcd(k, FpX_gcd(g,h, p), p);
854:
855: dk = degpol(k);
1.2 ! noro 856: if (DEBUGLEVEL>2)
! 857: {
! 858: fprintferr(" dedek: gcd has degree %ld\n", dk);
! 859: if (DEBUGLEVEL>5) fprintferr("initial parameters p=%Z,\n f=%Z\n",p,f);
! 860: }
1.1 noro 861: if (2*dk >= mf-1) return FpX_div(f,k,p);
862: return dk? (GEN)NULL: f;
863: }
864:
865: /* p-maximal order of Af; mf = v_p(Disc(f)) */
866: static GEN
867: maxord(GEN p,GEN f,long mf)
868: {
1.2 ! noro 869: const gpmem_t av = avma;
! 870: long j,r, flw = (cmpsi(degpol(f),p) < 0);
1.1 noro 871: GEN w,g,h,res;
872:
873: if (flw)
874: {
875: h = NULL; r = 0; /* gcc -Wall */
876: g = FpX_div(f, FpX_gcd(f,derivpol(f), p), p);
877: }
878: else
879: {
880: w=(GEN)factmod(f,p)[1]; r=lg(w)-1;
881: g = h = lift_intern((GEN)w[r]); /* largest factor */
882: for (j=1; j<r; j++) g = FpX_red(gmul(g, lift_intern((GEN)w[j])), p);
883: }
884: res = dedek(f,mf,p,g);
885: if (res)
886: res = dbasis(p,f,mf,polx[varn(f)],res);
887: else
888: {
889: if (flw) { w=(GEN)factmod(f,p)[1]; r=lg(w)-1; h=lift_intern((GEN)w[r]); }
890: res = (r==1)? nilord(p,f,mf,h,0): Decomp(p,f,mf,polx[varn(f)],f,h,0);
891: }
892: return gerepileupto(av,res);
893: }
894:
895: /* do a centermod on integer or rational number */
896: static GEN
897: polmodiaux(GEN x, GEN y, GEN ys2)
898: {
1.2 ! noro 899: if (typ(x)!=t_INT) x = mulii((GEN)x[1], mpinvmod((GEN)x[2],y));
! 900: return centermodii(x,y,ys2);
1.1 noro 901: }
902:
903: /* x polynomial with integer or rational coeff. Reduce them mod y IN PLACE */
904: GEN
905: polmodi(GEN x, GEN y)
906: {
907: long lx=lgef(x), i;
908: GEN ys2 = shifti(y,-1);
909: for (i=2; i<lx; i++) x[i]=(long)polmodiaux((GEN)x[i],y,ys2);
910: return normalizepol_i(x, lx);
911: }
912:
913: /* same but not in place */
914: GEN
915: polmodi_keep(GEN x, GEN y)
916: {
917: long lx=lgef(x), i;
918: GEN ys2 = shifti(y,-1);
919: GEN z = cgetg(lx,t_POL);
920: for (i=2; i<lx; i++) z[i]=(long)polmodiaux((GEN)x[i],y,ys2);
921: z[1]=x[1]; return normalizepol_i(z, lx);
922: }
923:
924: static GEN
1.2 ! noro 925: shiftpol(GEN x, long v)
! 926: {
! 927: x = addshiftw(x, zeropol(v), 1);
! 928: setvarn(x,v); return x;
! 929: }
! 930:
! 931: /* Sylvester's matrix, mod p^m (assumes f1 monic) */
! 932: static GEN
! 933: sylpm(GEN f1,GEN f2,GEN pm)
! 934: {
! 935: long n,j,v=varn(f1);
! 936: GEN a,h;
! 937:
! 938: n=degpol(f1); a=cgetg(n+1,t_MAT);
! 939: h = FpX_res(f2,f1,pm);
! 940: for (j=1;; j++)
! 941: {
! 942: a[j] = (long)pol_to_vec(h,n);
! 943: if (j == n) break;
! 944: h = FpX_res(shiftpol(h,v),f1,pm);
! 945: }
! 946: return hnfmodid(a,pm);
! 947: }
! 948:
! 949: /* polynomial gcd mod p^m (assumes f1 monic) */
! 950: GEN
! 951: gcdpm(GEN f1,GEN f2,GEN pm)
! 952: {
! 953: gpmem_t av=avma,tetpil;
! 954: long n,c,v=varn(f1);
! 955: GEN a,col;
! 956:
! 957: n=degpol(f1); a=sylpm(f1,f2,pm);
! 958: for (c=1; c<=n; c++)
! 959: if (signe(resii(gcoeff(a,c,c),pm))) break;
! 960: if (c > n) { avma=av; return zeropol(v); }
! 961:
! 962: col = gdiv((GEN)a[c], gcoeff(a,c,c)); tetpil=avma;
! 963: return gerepile(av,tetpil, gtopolyrev(col,v));
! 964: }
! 965:
! 966: /* reduced resultant mod p^m (assumes x monic) */
! 967: GEN
! 968: respm(GEN x,GEN y,GEN pm)
! 969: {
! 970: gpmem_t av = avma;
! 971: GEN p1 = sylpm(x,y,pm);
! 972:
! 973: p1 = gcoeff(p1,1,1);
! 974: if (egalii(p1,pm)) { avma = av; return gzero; }
! 975: return gerepileuptoint(av, icopy(p1));
! 976: }
! 977:
! 978: static GEN
1.1 noro 979: dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U)
980: {
1.2 ! noro 981: long n=degpol(f),dU,i;
1.1 noro 982: GEN b,ha,pd,pdp;
983:
984: if (n == 1) return gscalmat(gun, 1);
1.2 ! noro 985: if (DEBUGLEVEL>5)
1.1 noro 986: {
1.2 ! noro 987: fprintferr(" entering Dedekind Basis with parameters p=%Z\n",p);
! 988: fprintferr(" f = %Z,\n alpha = %Z\n",f,alpha);
1.1 noro 989: }
1.2 ! noro 990: ha = pd = gpowgs(p,mf/2); pdp = mulii(pd,p);
1.1 noro 991: dU = typ(U)==t_POL? degpol(U): 0;
992: b = cgetg(n,t_MAT); /* Z[a] + U/p Z[a] is maximal */
993: /* skip first column = gscalcol(pd,n) */
1.2 ! noro 994: for (i=1; i<n; i++)
1.1 noro 995: {
1.2 ! noro 996: if (i == dU)
1.1 noro 997: {
1.2 ! noro 998: ha = gdiv(gmul(pd,RX_RXQ_compo(U,alpha,f)),p);
1.1 noro 999: ha = polmodi(ha,pdp);
1000: }
1001: else
1002: {
1.2 ! noro 1003: GEN d, mod;
1.1 noro 1004: ha = gmul(ha,alpha);
1.2 ! noro 1005: ha = Q_remove_denom(ha, &d);
! 1006: mod = d? mulii(pdp,d): pdp;
1.1 noro 1007: ha = FpX_res(ha, f, mod);
1.2 ! noro 1008: if (d) ha = gdivexact(ha,d);
1.1 noro 1009: }
1.2 ! noro 1010: b[i] = (long)pol_to_vec(ha,n);
1.1 noro 1011: }
1012: b = hnfmodid(b,pd);
1013: if (DEBUGLEVEL>5) fprintferr(" new order: %Z\n",b);
1.2 ! noro 1014: return gdiv(b, pd);
1.1 noro 1015: }
1016:
1017: static GEN
1018: get_partial_order_as_pols(GEN p, GEN f)
1019: {
1.2 ! noro 1020: GEN b = maxord(p,f, ggval(ZX_disc(f),p));
! 1021: return mat_to_vecpol(b, varn(f));
! 1022: }
! 1023:
! 1024: long
! 1025: FpX_val(GEN x0, GEN t, GEN p, GEN *py)
! 1026: {
! 1027: long k;
! 1028: GEN r, y, x = x0;
1.1 noro 1029:
1.2 ! noro 1030: for (k=0; ; k++)
! 1031: {
! 1032: y = FpX_divres(x, t, p, &r);
! 1033: if (signe(r)) break;
! 1034: x = y;
1.1 noro 1035: }
1.2 ! noro 1036: *py = x; return k;
! 1037: }
! 1038:
! 1039: /* e in Qp, f i Zp. Return r = e mod (f, pk) */
! 1040: static GEN
! 1041: QpX_mod(GEN e, GEN f, GEN pk)
! 1042: {
! 1043: GEN mod, d;
! 1044: e = Q_remove_denom(e, &d);
! 1045: mod = d? mulii(pk,d): pk;
! 1046: e = FpX_res(e, centermod(f, mod), mod);
! 1047: e = FpX_center(e, mod);
! 1048: if (d) e = gdiv(e, d);
! 1049: return e;
1.1 noro 1050: }
1051:
1052: /* if flag != 0, factorization to precision r (maximal order otherwise) */
1053: GEN
1054: Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long flag)
1055: {
1.2 ! noro 1056: GEN fred,res,pr,pk,ph,pdr,b1,b2,a,e,f1,f2;
1.1 noro 1057:
1058: if (DEBUGLEVEL>2)
1059: {
1060: fprintferr(" entering Decomp ");
1061: if (DEBUGLEVEL>5)
1062: {
1063: fprintferr("with parameters: p=%Z, expo=%ld\n",p,mf);
1064: if (flag) fprintferr("precision = %ld\n",flag);
1065: fprintferr(" f=%Z",f);
1066: }
1067: fprintferr("\n");
1068: }
1.2 ! noro 1069:
! 1070: (void)FpX_val(chi, nu, p, &b1); /* nu irreducible mod p */
! 1071: b2 = FpX_div(chi, b1, p);
! 1072: a = FpX_mul(FpXQ_inv(b2, b1, p), b2, p);
! 1073: pdr = respm(f, derivpol(f), gpowgs(p,mf+1));
! 1074: e = RX_RXQ_compo(a, theta, f);
1.1 noro 1075: e = gdiv(polmodi(gmul(pdr,e), mulii(pdr,p)),pdr);
1076:
1077: pr = flag? gpowgs(p,flag): mulii(p,sqri(pdr));
1.2 ! noro 1078: pk = p; ph = mulii(pdr,pr);
1.1 noro 1079: /* E(t) - e(t) belongs to p^k Op, which is contained in p^(k-df)*Zp[xi] */
1080: while (cmpii(pk,ph) < 0)
1.2 ! noro 1081: { /* e <-- e^2(3-2e) mod p^2k */
! 1082: pk = sqri(pk);
1.1 noro 1083: e = gmul(gsqr(e), gsubsg(3,gmul2n(e,1)));
1.2 ! noro 1084: e = QpX_mod(e, f, pk);
1.1 noro 1085: }
1.2 ! noro 1086: fred = centermod(f, ph);
! 1087: f1 = gcdpm(fred, gmul(pdr,gsubsg(1,e)), ph);
! 1088: fred = centermod(fred, pr); /* pr | ph */
! 1089: f1 = centermod(f1, pr);
! 1090: f2 = FpX_div(fred,f1, pr);
! 1091: f2 = FpX_center(f2, pr);
1.1 noro 1092:
1.2 ! noro 1093: if (DEBUGLEVEL>5)
1.1 noro 1094: {
1095: fprintferr(" leaving Decomp");
1.2 ! noro 1096: fprintferr(" with parameters: f1 = %Z\nf2 = %Z\ne = %Z\n\n", f1,f2,e);
1.1 noro 1097: }
1098:
1099: if (flag)
1100: {
1.2 ! noro 1101: b1 = factorpadic4(f1,p,flag);
! 1102: b2 = factorpadic4(f2,p,flag); res = cgetg(3,t_MAT);
! 1103: res[1] = (long)concatsp((GEN)b1[1],(GEN)b2[1]);
! 1104: res[2] = (long)concatsp((GEN)b1[2],(GEN)b2[2]); return res;
1.1 noro 1105: }
1106: else
1107: {
1.2 ! noro 1108: GEN ib1, ib2;
! 1109: long n, n1, n2, i;
! 1110: ib1 = get_partial_order_as_pols(p,f1); n1 = lg(ib1)-1;
! 1111: ib2 = get_partial_order_as_pols(p,f2); n2 = lg(ib2)-1;
! 1112: n = n1+n2;
! 1113: res = cgetg(n+1, t_VEC);
1.1 noro 1114: for (i=1; i<=n1; i++)
1.2 ! noro 1115: res[i] = (long)QpX_mod(gmul(gmul(pdr,(GEN)ib1[i]),e), f, pdr);
! 1116: e = gsubsg(1,e); ib2 -= n1;
! 1117: for ( ; i<=n; i++)
! 1118: res[i] = (long)QpX_mod(gmul(gmul(pdr,(GEN)ib2[i]),e), f, pdr);
! 1119: res = vecpol_to_mat(res, n);
! 1120: return gdiv(hnfmodid(res,pdr), pdr); /* normalized integral basis */
1.1 noro 1121: }
1122: }
1123:
1124: /* minimum extension valuation: res[0]/res[1] (both are longs) */
1.2 ! noro 1125: static void
! 1126: vstar(GEN p,GEN h, long *L, long *E)
1.1 noro 1127: {
1128: long m,first,j,k,v,w;
1129:
1130: m=degpol(h); first=1; k=1; v=0;
1131: for (j=1; j<=m; j++)
1132: if (! gcmp0((GEN)h[m-j+2]))
1133: {
1134: w = ggval((GEN)h[m-j+2],p);
1135: if (first || w*k < v*j) { v=w; k=j; }
1136: first=0;
1137: }
1138: m = cgcd(v,k);
1.2 ! noro 1139: *L = v/m;
! 1140: *E = k/m;
1.1 noro 1141: }
1142:
1143: /* reduce the element elt modulo rd, taking first care of the denominators */
1144: static GEN
1145: redelt(GEN elt, GEN rd, GEN pd)
1146: {
1147: GEN den, nelt, nrd, relt;
1148:
1.2 ! noro 1149: den = ggcd(Q_denom(elt), pd);
1.1 noro 1150: nelt = gmul(den, elt);
1151: nrd = gmul(den, rd);
1152:
1153: if (typ(elt) == t_POL)
1154: relt = polmodi(nelt, nrd);
1155: else
1156: relt = centermod(nelt, nrd);
1157:
1158: return gdiv(relt, den);
1159: }
1160:
1161: /* compute the Newton sums of g(x) mod pp from its coefficients */
1162: GEN
1163: polsymmodpp(GEN g, GEN pp)
1164: {
1.2 ! noro 1165: gpmem_t av1, av2;
! 1166: long d = degpol(g), i, k;
1.1 noro 1167: GEN s , y;
1168:
1.2 ! noro 1169: y = cgetg(d + 1, t_COL);
1.1 noro 1170: y[1] = lstoi(d);
1171: for (k = 1; k < d; k++)
1172: {
1.2 ! noro 1173: av1 = avma;
1.1 noro 1174: s = centermod(gmulsg(k, polcoeff0(g,d-k,-1)), pp);
1175: for (i = 1; i < k; i++)
1176: s = gadd(s, gmul((GEN)y[k-i+1], polcoeff0(g,d-i,-1)));
1.2 ! noro 1177: av2 = avma;
1.1 noro 1178: y[k + 1] = lpile(av1, av2, centermod(gneg(s), pp));
1179: }
1180:
1181: return y;
1182: }
1183:
1184: /* no GC */
1185: static GEN
1186: manage_cache(GEN chi, GEN pp, GEN ns)
1187: {
1188: long j, n = degpol(chi);
1189: GEN ns2, npp = (GEN)ns[n+1];
1190:
1191: if (gcmp(pp, npp) > 0)
1192: {
1193: if (DEBUGLEVEL > 4)
1194: fprintferr("newtonsums: result too large to fit in cache\n");
1.2 ! noro 1195: return polsymmodpp(chi, pp);
1.1 noro 1196: }
1197:
1198: if (!signe((GEN)ns[1]))
1199: {
1.2 ! noro 1200: ns2 = polsymmodpp(chi, pp);
! 1201: for (j = 1; j <= n; j++)
1.1 noro 1202: affii((GEN)ns2[j], (GEN)ns[j]);
1203: }
1.2 ! noro 1204:
1.1 noro 1205: return ns;
1206: }
1207:
1.2 ! noro 1208: /* compute the Newton sums modulo pp of the characteristic
1.1 noro 1209: polynomial of a(x) mod g(x) */
1210: static GEN
1211: newtonsums(GEN a, GEN chi, GEN pp, GEN ns)
1212: {
1213: GEN va, pa, s, ns2;
1.2 ! noro 1214: long j, k, n = degpol(chi);
! 1215: gpmem_t av2, lim;
1.1 noro 1216:
1217: ns2 = manage_cache(chi, pp, ns);
1218:
1219: av2 = avma;
1220: lim = stack_lim(av2, 1);
1221:
1222: pa = gun;
1223: va = zerovec(n);
1224:
1225: for (j = 1; j <= n; j++)
1226: {
1227: pa = gmul(pa, a);
1.2 ! noro 1228: pa = polmodi(pa, pp);
1.1 noro 1229: pa = gmod(pa, chi);
1.2 ! noro 1230: pa = polmodi(pa, pp);
1.1 noro 1231:
1232: s = gzero;
1233:
1234: for (k = 0; k <= n-1; k++)
1235: s = addii(s, mulii(polcoeff0(pa, k, -1), (GEN)ns2[k+1]));
1.2 ! noro 1236:
! 1237: va[j] = (long)centermod(s, pp);
! 1238:
1.1 noro 1239: if (low_stack(lim, stack_lim(av2, 1)))
1240: {
1241: GEN *gptr[2];
1242: gptr[0]=&pa; gptr[1]=&va;
1243: if(DEBUGMEM>1) err(warnmem, "newtonsums");
1244: gerepilemany(av2, gptr, 2);
1245: }
1246: }
1247:
1248: return va;
1249: }
1250:
1251: /* compute the characteristic polynomial of a mod g
1252: to a precision of pp using Newton sums */
1253: static GEN
1254: newtoncharpoly(GEN a, GEN chi, GEN pp, GEN ns)
1255: {
1256: GEN v, c, s, t;
1.2 ! noro 1257: long n = degpol(chi), j, k, vn = varn(chi);
! 1258: gpmem_t av = avma, av2, lim;
1.1 noro 1259:
1260: v = newtonsums(a, chi, pp, ns);
1261: av2 = avma;
1262: lim = stack_lim(av2, 1);
1263: c = cgetg(n + 2, t_VEC);
1264: c[1] = un;
1265: if (n%2) c[1] = lneg((GEN)c[1]);
1266: for (k = 2; k <= n+1; k++) c[k] = zero;
1267:
1268: for (k = 2; k <= n+1; k++)
1269: {
1270: s = gzero;
1271: for (j = 1; j < k; j++)
1272: {
1273: t = gmul((GEN)v[j], (GEN)c[k-j]);
1274: if (!(j%2)) t = gneg(t);
1275: s = gadd(s, t);
1276: }
1277: c[k] = ldiv(s, stoi(k - 1));
1278:
1279: if (low_stack(lim, stack_lim(av2, 1)))
1280: {
1281: if(DEBUGMEM>1) err(warnmem, "newtoncharpoly");
1282: c = gerepilecopy(av2, c);
1283: }
1284: }
1.2 ! noro 1285:
1.1 noro 1286: k = (n%2)? 1: 2;
1.2 ! noro 1287: for ( ; k <= n+1; k += 2)
1.1 noro 1288: c[k] = lneg((GEN)c[k]);
1289:
1290: return gerepileupto(av, gtopoly(c, vn));
1291: }
1292:
1.2 ! noro 1293: /* return v_p(n!) */
! 1294: static long
! 1295: val_fact(long n, GEN pp)
! 1296: {
! 1297: long p, q, v;
! 1298: if (is_bigint(pp)) return 0;
! 1299: q = p = itos(pp); v = 0;
! 1300: do { v += n/q; q *= p; } while (n >= q);
! 1301: return v;
! 1302: }
! 1303:
! 1304: static GEN
1.1 noro 1305: mycaract(GEN f, GEN beta, GEN p, GEN pp, GEN ns)
1306: {
1.2 ! noro 1307: GEN p1, chi, npp;
! 1308: long v = varn(f), n = degpol(f);
1.1 noro 1309:
1310: if (gcmp0(beta)) return zeropol(v);
1311:
1.2 ! noro 1312: beta = Q_primitive_part(beta,&p1);
1.1 noro 1313: if (!pp)
1.2 ! noro 1314: chi = ZX_caract(f, beta, v);
1.1 noro 1315: else
1316: {
1.2 ! noro 1317: npp = mulii(pp, gpowgs(p, val_fact(n, p)));
! 1318: if (p1) npp = mulii(npp, gpowgs(denom(p1), n));
1.1 noro 1319:
1320: chi = newtoncharpoly(beta, f, npp, ns);
1321: }
1322:
1.2 ! noro 1323: if (p1) chi = rescale_pol(chi,p1);
! 1324:
! 1325: if (!pp) return chi;
1.1 noro 1326:
1327: /* this can happen only if gamma is incorrect (not an integer) */
1.2 ! noro 1328: if (divise(Q_denom(chi), p)) return NULL;
! 1329:
! 1330: return redelt(chi, pp, pp);
1.1 noro 1331: }
1332:
1333: /* Factor characteristic polynomial of beta mod p */
1334: static GEN
1335: factcp(GEN p, GEN f, GEN beta, GEN pp, GEN ns)
1336: {
1.2 ! noro 1337: gpmem_t av;
1.1 noro 1338: GEN chi,nu, b = cgetg(4,t_VEC);
1.2 ! noro 1339: long l;
1.1 noro 1340:
1341: chi=mycaract(f,beta,p,pp,ns);
1342: av=avma; nu=(GEN)factmod(chi,p)[1]; l=lg(nu)-1;
1343: nu=lift_intern((GEN)nu[1]);
1344: b[1]=(long)chi;
1345: b[2]=lpilecopy(av,nu);
1346: b[3]=lstoi(l); return b;
1347: }
1348:
1349: /* return the prime element in Zp[phi] */
1350: static GEN
1351: getprime(GEN p, GEN chi, GEN phi, GEN chip, GEN nup, long *Lp, long *Ep)
1352: {
1.2 ! noro 1353: long v = varn(chi), L, E, r, s;
! 1354: GEN chin, pip, pp;
1.1 noro 1355:
1356: if (gegal(nup, polx[v]))
1357: chin = chip;
1358: else
1359: chin = mycaract(chip, nup, p, NULL, NULL);
1360:
1.2 ! noro 1361: vstar(p, chin, &L, &E);
! 1362: (void)cbezout(L, -E, &r, &s);
! 1363: if (r <= 0)
! 1364: {
! 1365: long q = 1 + ((-r) / E);
! 1366: r += q*E;
! 1367: s += q*L;
1.1 noro 1368: }
1369:
1.2 ! noro 1370: pip = RX_RXQ_compo(nup, phi, chi);
! 1371: pip = lift_intern(gpowgs(gmodulcp(pip, chi), r));
! 1372: pp = gpowgs(p, s);
1.1 noro 1373:
1.2 ! noro 1374: *Lp = L;
! 1375: *Ep = E; return gdiv(pip, pp);
1.1 noro 1376: }
1377:
1378: static GEN
1.2 ! noro 1379: update_alpha(GEN p, GEN fx, GEN alph, GEN chi, GEN pmr, GEN pmf, long mf,
1.1 noro 1380: GEN ns)
1381: {
1382: long l, v = varn(fx);
1383: GEN nalph = NULL, nchi, w, nnu, pdr, npmr, rep;
1384:
1385: affii(gzero, (GEN)ns[1]); /* kill cache */
1386:
1387: if (!chi)
1388: nchi = mycaract(fx, alph, p, pmf, ns);
1389: else
1390: {
1391: nchi = chi;
1392: nalph = alph;
1393: }
1394:
1395: pdr = respm(nchi, derivpol(nchi), pmr);
1396: for (;;)
1397: {
1398: if (signe(pdr)) break;
1.2 ! noro 1399: if (!nalph) nalph = gadd(alph, gmul(p, polx[v]));
! 1400: /* nchi was too reduced at this point; try a larger precision */
! 1401: pmf = sqri(pmf);
! 1402: nchi = mycaract(fx, nalph, p, pmf, ns);
! 1403: pdr = respm(nchi, derivpol(nchi), pmf);
1.1 noro 1404: if (signe(pdr)) break;
1405: if (DEBUGLEVEL >= 6)
1406: fprintferr(" non separable polynomial in update_alpha!\n");
1.2 ! noro 1407: /* at this point, we assume that chi is not square-free */
1.1 noro 1408: nalph = gadd(nalph, gmul(p, polx[v]));
1409: w = factcp(p, fx, nalph, NULL, ns);
1.2 ! noro 1410: nchi = (GEN)w[1];
! 1411: nnu = (GEN)w[2];
1.1 noro 1412: l = itos((GEN)w[3]);
1413: if (l > 1) return Decomp(p, fx, mf, nalph, nchi, nnu, 0);
1414: pdr = respm(nchi, derivpol(nchi), pmr);
1415: }
1416:
1417: if (is_pm1(pdr))
1418: npmr = gun;
1419: else
1420: {
1421: npmr = mulii(sqri(pdr), p);
1422: nchi = polmodi(nchi, npmr);
1.2 ! noro 1423: nalph = redelt(nalph? nalph: alph, npmr, pmf);
1.1 noro 1424: }
1425:
1426: affii(gzero, (GEN)ns[1]); /* kill cache again (contains data for fx) */
1427:
1428: rep = cgetg(5, t_VEC);
1429: rep[1] = (long)nalph;
1430: rep[2] = (long)nchi;
1431: rep[3] = (long)npmr;
1432: rep[4] = lmulii(p, pdr);
1433:
1434: return rep;
1435: }
1436:
1.2 ! noro 1437: /* flag != 0 iff we're looking for the p-adic factorization,
1.1 noro 1438: in which case it is the p-adic precision we want */
1439: GEN
1440: nilord(GEN p, GEN fx, long mf, GEN gx, long flag)
1441: {
1.2 ! noro 1442: long L, E, Fa, La, Ea, oE, Fg, eq, er, v = varn(fx), i, nv, Le, Ee, N, l, vn;
1.1 noro 1443: GEN p1, alph, chi, nu, w, phi, pmf, pdr, pmr, kapp, pie, chib, ns;
1.2 ! noro 1444: GEN gamm, chig, nug, delt, beta, eta, chie, nue, pia, opa;
1.1 noro 1445:
1.2 ! noro 1446: if (DEBUGLEVEL>2)
1.1 noro 1447: {
1.2 ! noro 1448: fprintferr(" entering Nilord");
! 1449: if (DEBUGLEVEL>4)
1.1 noro 1450: {
1451: fprintferr(" with parameters: p = %Z, expo = %ld\n", p, mf);
1452: fprintferr(" fx = %Z, gx = %Z", fx, gx);
1453: }
1454: fprintferr("\n");
1455: }
1.2 ! noro 1456:
1.1 noro 1457: pmf = gpowgs(p, mf + 1);
1458: pdr = respm(fx, derivpol(fx), pmf);
1459: pmr = mulii(sqri(pdr), p);
1460: pdr = mulii(p, pdr);
1461: chi = polmodi_keep(fx, pmr);
1462:
1463: alph = polx[v];
1.2 ! noro 1464: nu = gx;
1.1 noro 1465: N = degpol(fx);
1466: oE = 0;
1467: opa = NULL;
1468:
1469: /* used to cache the newton sums of chi */
1470: ns = cgetg(N + 2, t_COL);
1471: p1 = powgi(p, gceil(gdivsg(N, mulii(p, subis(p, 1))))); /* p^(N/(p(p-1))) */
1472: p1 = mulii(p1, mulii(pmf, gpowgs(pmr, N)));
1473: l = lg(p1); /* enough in general... */
1474: for (i = 1; i <= N + 1; i++) ns[i] = lgeti(l);
1475: ns[N+1] = (long)p1;
1476: affii(gzero, (GEN)ns[1]); /* zero means: need to be computed */
1477:
1478: for(;;)
1479: {
1.2 ! noro 1480: /* kappa needs to be recomputed */
1.1 noro 1481: kapp = NULL;
1482: Fa = degpol(nu);
1483: /* the prime element in Zp[alpha] */
1484: pia = getprime(p, chi, polx[v], chi, nu, &La, &Ea);
1485: pia = redelt(pia, pmr, pmf);
1486:
1487: if (Ea < oE)
1488: {
1489: alph = gadd(alph, opa);
1490: w = update_alpha(p, fx, alph, NULL, pmr, pmf, mf, ns);
1491: alph = (GEN)w[1];
1492: chi = (GEN)w[2];
1493: pmr = (GEN)w[3];
1494: pdr = (GEN)w[4];
1495: kapp = NULL;
1496: pia = getprime(p, chi, polx[v], chi, nu, &La, &Ea);
1497: pia = redelt(pia, pmr, pmf);
1498: }
1.2 ! noro 1499:
! 1500: oE = Ea; opa = RX_RXQ_compo(pia, alph, fx);
1.1 noro 1501:
1502: if (DEBUGLEVEL >= 5)
1503: fprintferr(" Fa = %ld and Ea = %ld \n", Fa, Ea);
1.2 ! noro 1504:
1.1 noro 1505: /* we change alpha such that nu = pia */
1506: if (La > 1)
1507: {
1.2 ! noro 1508: alph = gadd(alph, RX_RXQ_compo(pia, alph, fx));
1.1 noro 1509:
1510: w = update_alpha(p, fx, alph, NULL, pmr, pmf, mf, ns);
1511: alph = (GEN)w[1];
1512: chi = (GEN)w[2];
1513: pmr = (GEN)w[3];
1514: pdr = (GEN)w[4];
1515: }
1516:
1517: /* if Ea*Fa == N then O = Zp[alpha] */
1.2 ! noro 1518: if (Ea*Fa == N)
1.1 noro 1519: {
1520: if (flag) return NULL;
1521:
1522: alph = redelt(alph, sqri(p), pmf);
1523: return dbasis(p, fx, mf, alph, p);
1524: }
1525:
1526: /* during the process beta tends to a factor of chi */
1527: beta = lift_intern(gpowgs(gmodulcp(nu, chi), Ea));
1528:
1529: for (;;)
1530: {
1531: if (DEBUGLEVEL >= 5)
1532: fprintferr(" beta = %Z\n", beta);
1533:
1534: /* if pmf divides norm(beta) then it's useless */
1535: p1 = gmod(gnorm(gmodulcp(beta, chi)), pmf);
1536: if (signe(p1))
1537: {
1538: chib = NULL;
1539: vn = ggval(p1, p);
1540: eq = (long)(vn / N);
1541: er = (long)(vn*Ea/N - eq*Ea);
1542: }
1.2 ! noro 1543: else
1.1 noro 1544: {
1545: chib = mycaract(chi, beta, NULL, NULL, ns);
1.2 ! noro 1546: vstar(p, chib, &L, &E);
! 1547: eq = (long)(L / E);
! 1548: er = (long)(L*Ea / E - eq*Ea);
1.1 noro 1549: }
1.2 ! noro 1550:
! 1551: /* eq and er are such that gamma = beta.p^-eq.nu^-er is a unit */
1.1 noro 1552: if (eq) gamm = gdiv(beta, gpowgs(p, eq));
1553: else gamm = beta;
1554:
1555: if (er)
1556: {
1557: /* kappa = nu^-1 in Zp[alpha] */
1558: if (!kapp)
1559: {
1.2 ! noro 1560: kapp = QX_invmod(nu, chi);
1.1 noro 1561: kapp = redelt(kapp, pmr, pmf);
1562: kapp = gmodulcp(kapp, chi);
1563: }
1564: gamm = lift(gmul(gamm, gpowgs(kapp, er)));
1565: gamm = redelt(gamm, p, pmf);
1566: }
1567:
1568: if (DEBUGLEVEL >= 6)
1569: fprintferr(" gamma = %Z\n", gamm);
1570:
1571: if (er || !chib)
1.2 ! noro 1572: chig = mycaract(chi, gamm, p, pmf, ns);
1.1 noro 1573: else
1574: {
1575: chig = poleval(chib, gmul(polx[v], gpowgs(p, eq)));
1576: chig = gdiv(chig, gpowgs(p, N*eq));
1577: chig = polmodi(chig, pmf);
1578: }
1579:
1.2 ! noro 1580: if (!chig || !gcmp1(Q_denom(chig)))
1.1 noro 1581: {
1.2 ! noro 1582: /* Valuation of beta was wrong. This means that
! 1583: gamma fails the v*-test */
! 1584: chib = mycaract(chi, beta, p, NULL, ns);
! 1585: vstar(p, chib, &L, &E);
! 1586: eq = (long)(-L / E);
! 1587: er = (long)(-L*Ea / E - eq*Ea);
! 1588:
! 1589: if (eq) gamm = gmul(beta, gpowgs(p, eq));
! 1590: else gamm = beta;
! 1591: if (er)
1.1 noro 1592: {
1593: gamm = gmul(gamm, gpowgs(nu, er));
1594: gamm = gmod(gamm, chi);
1595: gamm = redelt(gamm, p, pmr);
1596: }
1.2 ! noro 1597: chig = mycaract(chi, gamm, p, pmf, ns);
1.1 noro 1598: }
1599:
1600: nug = (GEN)factmod(chig, p)[1];
1601: l = lg(nug) - 1;
1602: nug = lift((GEN)nug[l]);
1603:
1604: if (l > 1)
1605: {
1606: /* there are at least 2 factors mod. p => chi can be split */
1.2 ! noro 1607: phi = RX_RXQ_compo(gamm, alph, fx);
1.1 noro 1608: phi = redelt(phi, p, pmf);
1609: if (flag) mf += 3;
1610: return Decomp(p, fx, mf, phi, chig, nug, flag);
1611: }
1612:
1613: Fg = degpol(nug);
1614: if (Fa%Fg)
1615: {
1616: if (DEBUGLEVEL >= 5)
1617: fprintferr(" Increasing Fa\n");
1618: /* we compute a new element such F = lcm(Fa, Fg) */
1619: w = testb2(p, chi, Fa, gamm, pmf, Fg, ns);
1620: if (gcmp1((GEN)w[1]))
1621: {
1622: /* there are at least 2 factors mod. p => chi can be split */
1.2 ! noro 1623: phi = RX_RXQ_compo((GEN)w[2], alph, fx);
1.1 noro 1624: phi = redelt(phi, p, pmf);
1625: if (flag) mf += 3;
1626: return Decomp(p, fx, mf, phi, (GEN)w[3], (GEN)w[4], flag);
1627: }
1628: break;
1629: }
1630:
1631: /* we look for a root delta of nug in Fp[alpha] such that
1632: vp(gamma - delta) > 0. This root can then be used to
1633: improved the approximation given by beta */
1634: nv = fetch_var();
1635: w = Fp_factor_irred(nug, p, gsubst(nu, varn(nu), polx[nv]));
1636: if (degpol(w[1]) != 1)
1637: err(talker,"bug in nilord (no root). Is p a prime ?");
1638:
1639: for (i = 1;; i++)
1640: {
1641: if (i >= lg(w))
1.2 ! noro 1642: err(talker, "bug in nilord (no suitable root), is p = %Z a prime?",
1.1 noro 1643: (long)p);
1644: delt = gneg_i(gsubst(gcoeff(w, 2, i), nv, polx[v]));
1.2 ! noro 1645: eta = gsub(gamm, delt);
1.1 noro 1646: if (typ(delt) == t_INT)
1647: {
1648: chie = poleval(chig, gadd(polx[v], delt));
1649: nue = (GEN)factmod(chie, p)[1];
1650: l = lg(nue) - 1;
1651: nue = lift((GEN)nue[l]);
1652: }
1653: else
1.2 ! noro 1654: {
1.1 noro 1655: p1 = factcp(p, chi, eta, pmf, ns);
1656: chie = (GEN)p1[1];
1.2 ! noro 1657: nue = (GEN)p1[2];
1.1 noro 1658: l = itos((GEN)p1[3]);
1659: }
1660: if (l > 1)
1661: {
1662: /* there are at least 2 factors mod. p => chi can be split */
1.2 ! noro 1663: (void)delete_var();
! 1664: phi = RX_RXQ_compo(eta, alph, fx);
1.1 noro 1665: phi = redelt(phi, p, pmf);
1666: if (flag) mf += 3;
1667: return Decomp(p, fx, mf, phi, chie, nue, flag);
1668: }
1.2 ! noro 1669:
1.1 noro 1670: /* if vp(eta) = vp(gamma - delta) > 0 */
1671: if (gegal(nue, polx[v])) break;
1672: }
1.2 ! noro 1673: (void)delete_var();
1.1 noro 1674:
1.2 ! noro 1675: if (!signe(modii(constant_term(chie), pmr)))
1.1 noro 1676: chie = mycaract(chi, eta, p, pmf, ns);
1677:
1678: pie = getprime(p, chi, eta, chie, nue, &Le, &Ee);
1679: if (Ea%Ee)
1680: {
1681: if (DEBUGLEVEL >= 5)
1682: fprintferr(" Increasing Ea\n");
1683: pie = redelt(pie, p, pmf);
1684: /* we compute a new element such E = lcm(Ea, Ee) */
1685: w = testc2(p, chi, pmr, pmf, nu, Ea, pie, Ee, ns);
1686: if (gcmp1((GEN)w[1]))
1687: {
1688: /* there are at least 2 factors mod. p => chi can be split */
1.2 ! noro 1689: phi = RX_RXQ_compo((GEN)w[2], alph, fx);
1.1 noro 1690: phi = redelt(phi, p, pmf);
1691: if (flag) mf += 3;
1692: return Decomp(p, fx, mf, phi, (GEN)w[3], (GEN)w[4], flag);
1693: }
1694: break;
1695: }
1.2 ! noro 1696:
1.1 noro 1697: if (eq) delt = gmul(delt, gpowgs(p, eq));
1698: if (er) delt = gmul(delt, gpowgs(nu, er));
1699: beta = gsub(beta, delt);
1700: }
1701:
1702: /* we replace alpha by a new alpha with a larger F or E */
1.2 ! noro 1703: alph = RX_RXQ_compo((GEN)w[2], alph, fx);
1.1 noro 1704: chi = (GEN)w[3];
1705: nu = (GEN)w[4];
1706:
1707: w = update_alpha(p, fx, alph, chi, pmr, pmf, mf, ns);
1708: alph = (GEN)w[1];
1709: chi = (GEN)w[2];
1710: pmr = (GEN)w[3];
1711: pdr = (GEN)w[4];
1712:
1713: /* that can happen if p does not divide the field discriminant! */
1714: if (is_pm1(pmr))
1715: {
1716: if (flag)
1717: {
1718: p1 = lift((GEN)factmod(chi, p)[1]);
1719: l = lg(p1) - 1;
1720: if (l == 1) return NULL;
1721: phi = redelt(alph, p, pmf);
1722: return Decomp(p, fx, ggval(pmf, p), phi, chi, (GEN)p1[l], flag);
1723: }
1724: else
1725: return dbasis(p, fx, mf, alph, chi);
1726: }
1727: }
1728: }
1729:
1730: /* Returns [1,phi,chi,nu] if phi non-primary
1.2 ! noro 1731: * [2,phi,chi,nu] with F_phi = lcm (F_alpha, F_theta)
! 1732: * and E_phi = E_alpha
1.1 noro 1733: */
1734: static GEN
1735: testb2(GEN p, GEN fa, long Fa, GEN theta, GEN pmf, long Ft, GEN ns)
1736: {
1737: long m, Dat, t, v = varn(fa);
1738: GEN b, w, phi, h;
1739:
1.2 ! noro 1740: Dat = clcm(Fa, Ft) + 3;
1.1 noro 1741: b = cgetg(5, t_VEC);
1.2 ! noro 1742: m = p[2];
1.1 noro 1743: if (degpol(p) > 0 || m < 0) m = 0;
1744:
1745: for (t = 1;; t++)
1746: {
1747: h = m? stopoly(t, m, v): scalarpol(stoi(t), v);
1748: phi = gadd(theta, gmod(h, fa));
1749: w = factcp(p, fa, phi, pmf, ns);
1750: h = (GEN)w[3];
1751: if (h[2] > 1) { b[1] = un; break; }
1752: if (lgef(w[2]) == Dat) { b[1] = deux; break; }
1753: }
1.2 ! noro 1754:
! 1755: b[2] = (long)phi;
! 1756: b[3] = w[1];
! 1757: b[4] = w[2];
1.1 noro 1758: return b;
1759: }
1760:
1761: /* Returns [1, phi, chi, nu] if phi non-primary
1762: * [2, phi, chi, nu] if E_phi = lcm (E_alpha, E_theta)
1763: */
1764: static GEN
1.2 ! noro 1765: testc2(GEN p, GEN fa, GEN pmr, GEN pmf, GEN alph2, long Ea, GEN thet2,
1.1 noro 1766: long Et, GEN ns)
1767: {
1768: GEN b, c1, c2, c3, psi, phi, w, h;
1769: long r, s, t, v = varn(fa);
1770:
1771: b=cgetg(5, t_VEC);
1772:
1.2 ! noro 1773: (void)cbezout(Ea, Et, &r, &s); t = 0;
1.1 noro 1774: while (r < 0) { r = r + Et; t++; }
1775: while (s < 0) { s = s + Ea; t++; }
1776:
1.2 ! noro 1777: c1 = lift_intern(gpowgs(gmodulcp(alph2, fa), s));
! 1778: c2 = lift_intern(gpowgs(gmodulcp(thet2, fa), r));
! 1779: c3 = gdiv(gmod(gmul(c1, c2), fa), gpowgs(p, t));
1.1 noro 1780:
1781: psi = redelt(c3, pmr, pmr);
1782: phi = gadd(polx[v], psi);
1783:
1784: w = factcp(p, fa, phi, pmf, ns);
1785: h = (GEN)w[3];
1786: b[1] = (h[2] > 1)? un: deux;
1.2 ! noro 1787: b[2] = (long)phi;
! 1788: b[3] = w[1];
! 1789: b[4] = w[2];
1.1 noro 1790: return b;
1791: }
1792:
1.2 ! noro 1793: /*******************************************************************/
! 1794: /* */
! 1795: /* 2-ELT REPRESENTATION FOR PRIME IDEALS (dividing index) */
! 1796: /* */
! 1797: /*******************************************************************/
! 1798: /* beta = generators of prime ideal (vectors). Return "small" random elt in P */
! 1799: static GEN
! 1800: random_elt_in_P(GEN beta, long small)
1.1 noro 1801: {
1.2 ! noro 1802: long z, i, j, la, lbeta = lg(beta);
! 1803: GEN a = NULL, B;
1.1 noro 1804:
1.2 ! noro 1805: /* compute sum random(-7..8) * beta[i] */
! 1806: if (small)
1.1 noro 1807: {
1.2 ! noro 1808: la = lg(beta[1]);
! 1809: if (small > 7) small = 0;
! 1810: for (i=1; i<lbeta; i++)
! 1811: { /* Warning: change definition of 'small' if you change this */
! 1812: z = mymyrand() >> (BITS_IN_RANDOM-5); /* in [0,15] */
! 1813: if (z >= 9) z -= 7;
! 1814: if (small) z %= small;
! 1815: if (!z) continue;
! 1816: B = (GEN)beta[i];
! 1817: if (a)
! 1818: for (j = 1; j < la; j++) a[j] += z * B[j];
! 1819: else {
! 1820: a = cgetg(la, t_VECSMALL);
! 1821: for (j = 1; j <la; j++) a[j] = z * B[j];
! 1822: }
! 1823: }
1.1 noro 1824: }
1.2 ! noro 1825: else
! 1826: for (i=1; i<lbeta; i++)
! 1827: {
! 1828: z = mymyrand() >> (BITS_IN_RANDOM-5);
! 1829: if (z >= 9) z -= 7;
! 1830: if (!z) continue;
! 1831: B = gmulsg(z, (GEN)beta[i]);
! 1832: a = a? gadd(a, B): B;
! 1833: }
! 1834: return a;
1.1 noro 1835: }
1836:
1.2 ! noro 1837: /* routines to check uniformizer algebraically (as t_POL) */
! 1838:
! 1839: /* a t_POL, D t_INT (or NULL if 1), q = p^(f+1). a/D in pr | p, norm(pr) = pf.
! 1840: * Return 1 if (a/D,p) = pr, and 0 otherwise */
! 1841: static int
! 1842: is_uniformizer(GEN D, GEN a, GEN T, GEN q)
! 1843: {
! 1844: GEN N = ZX_resultant_all(T, a, D, 0); /* norm(a) */
! 1845: return (resii(N, q) != gzero);
! 1846: }
1.1 noro 1847:
1.2 ! noro 1848: /* a/D or a/D + p uniformizer ? */
1.1 noro 1849: static GEN
1.2 ! noro 1850: prime_check_elt(GEN D, GEN Dp, GEN a, GEN T, GEN q, int ramif)
! 1851: {
! 1852: if (is_uniformizer(D,a,T,q)) return a;
! 1853: /* can't succeed if e > 1 */
! 1854: if (!ramif && is_uniformizer(D,gadd(a,Dp),T,q)) return a;
! 1855: return NULL;
! 1856: }
! 1857:
! 1858: /* P given by an Fp basis */
! 1859: static GEN
! 1860: compute_pr2(GEN nf, GEN P, GEN p, int *ramif)
! 1861: {
! 1862: long i, j, m = lg(P)-1;
! 1863: GEN mul, P2 = gmul(p, P);
! 1864: for (i=1; i<=m; i++)
! 1865: {
! 1866: mul = eltmul_get_table(nf, (GEN)P[i]);
! 1867: for (j=1; j<=i; j++)
! 1868: P2 = concatsp(P2, gmul(mul, (GEN)P[j]));
! 1869: }
! 1870: P2 = hnfmodid(P2, sqri(p)); /* HNF of pr^2 */
! 1871: *ramif = egalii(gcoeff(P2,1,1), p); /* pr/p ramified ? */
! 1872: return P2;
! 1873: }
! 1874:
! 1875: static GEN
! 1876: random_unif_loop_pol(GEN nf, GEN P, GEN D, GEN Dp, GEN beta, GEN pol,
! 1877: GEN p, GEN q)
! 1878: {
! 1879: long small, i, c = 0, m = lg(P)-1, N = degpol(nf[1]), keep = getrand();
! 1880: int ramif;
! 1881: GEN a, P2;
! 1882: gpmem_t av;
! 1883:
! 1884: for(i=1; i<=m; i++)
! 1885: if ((a = prime_check_elt(D,Dp,(GEN)beta[i],pol,q,0))) return a;
! 1886:
! 1887: (void)setrand(1);
! 1888: if (DEBUGLEVEL) fprintferr("uniformizer_loop, hard case: ");
! 1889: P2 = compute_pr2(nf, P, p, &ramif);
! 1890:
! 1891: a = mulis(Dp, 8*degpol(nf[1]));
! 1892: if (is_bigint(a)) small = 0;
! 1893: else
! 1894: {
! 1895: ulong mod = itou(Dp);
! 1896: small = itos(p);
! 1897: for (i=1; i<=m; i++)
! 1898: beta[i] = (long)u_Fp_FpV(pol_to_vec((GEN)beta[i], N), mod);
! 1899: }
! 1900:
! 1901: for(av = avma; ;avma = av)
! 1902: {
! 1903: if (DEBUGLEVEL && (++c & 0x3f) == 0) fprintferr("%d ", c);
! 1904: a = random_elt_in_P(beta, small);
! 1905: if (!a) continue;
! 1906: if (small) a = vec_to_pol(small_to_vec(a), varn(pol));
! 1907: a = centermod(a, Dp);
! 1908: if ((a = prime_check_elt(D,Dp,a,pol,q,ramif)))
! 1909: {
! 1910: if (DEBUGLEVEL) fprintferr("\n");
! 1911: (void)setrand(keep); return a;
! 1912: }
! 1913: }
! 1914: }
! 1915:
! 1916: /* routines to check uniformizer numerically (as t_COL) */
! 1917:
! 1918: static int
! 1919: vec_is_uniformizer(GEN a, GEN q, long r1)
! 1920: {
! 1921: long e;
! 1922: GEN N = grndtoi( norm_by_embed(r1, a), &e );
! 1923: if (e > -5) err(precer, "vec_is_uniformizer");
! 1924: return (resii(N, q) != gzero);
! 1925: }
! 1926:
! 1927: static GEN
! 1928: prime_check_eltvec(GEN a, GEN M, GEN p, GEN q, long r1, long small, int ramif)
1.1 noro 1929: {
1.2 ! noro 1930: GEN ar;
! 1931: long i,l;
! 1932: a = centermod(a, p);
! 1933: if (small) ar = gmul_mat_smallvec(M, a);
! 1934: else ar = gmul(M, a);
! 1935:
! 1936: if (vec_is_uniformizer(ar, q, r1)) return small? small_to_col(a): a;
! 1937: /* can't succeed if e > 1 */
! 1938: if (ramif) return NULL;
! 1939: l = lg(ar);
! 1940: for (i=1; i<l; i++) ar[i] = ladd((GEN)ar[i], p);
! 1941: if (!vec_is_uniformizer(ar, q, r1)) return NULL;
! 1942: if (small) a = small_to_col(a);
! 1943: a[1] = laddii((GEN)a[1],p); return a;
1.1 noro 1944: }
1945:
1946: static GEN
1.2 ! noro 1947: random_unif_loop_vec(GEN nf, GEN P, GEN p, GEN q)
1.1 noro 1948: {
1.2 ! noro 1949: long small, r1, i, c = 0, m = lg(P)-1, keep = getrand();
! 1950: int ramif;
! 1951: GEN a, P2, beta, M = gmael(nf,5,1);
! 1952: gpmem_t av;
! 1953:
! 1954: r1 = nf_get_r1(nf);
! 1955: a = mulis(p, 8*degpol(nf[1]));
! 1956: if (is_bigint(a)) { small = 0; beta = P; }
! 1957: else
! 1958: {
! 1959: small = itos(p); beta = cgetg(m+1, t_MAT);
! 1960: for (i=1; i<=m; i++)
! 1961: beta[i] = (long)u_Fp_FpV((GEN)P[i], itou(p));
! 1962: }
! 1963: for(i=1; i<=m; i++)
! 1964: if ((a = prime_check_eltvec((GEN)beta[i],M,p,q,r1,small,0))) return a;
! 1965:
! 1966: (void)setrand(1);
! 1967: if (DEBUGLEVEL) fprintferr("uniformizer_loop, hard case: ");
! 1968: P2 = compute_pr2(nf, P, p, &ramif);
1.1 noro 1969:
1.2 ! noro 1970: for(av = avma; ;avma = av)
1.1 noro 1971: {
1.2 ! noro 1972: if (DEBUGLEVEL && (++c & 0x3f) == 0) fprintferr("%d ", c);
! 1973: a = random_elt_in_P(beta, small);
! 1974: if (!a) continue;
! 1975: if ((a = prime_check_eltvec(a,M,p,q,r1,small,0)))
! 1976: {
! 1977: if (DEBUGLEVEL) fprintferr("\n");
! 1978: (void)setrand(keep); return a;
! 1979: }
1.1 noro 1980: }
1981: }
1982:
1.2 ! noro 1983: /* L = Fp-bases for prime ideals of O_K/p. Return a p-uniformizer for
! 1984: * L[i] using the approximation theorem */
! 1985: static GEN
! 1986: uniformizer_appr(GEN nf, GEN L, long i, GEN p)
1.1 noro 1987: {
1.2 ! noro 1988: long j, l;
! 1989: GEN inter = NULL, P = (GEN)L[i], P2, Q, u, v, pi;
! 1990: int ramif;
1.1 noro 1991:
1.2 ! noro 1992: l = lg(L)-1;
! 1993: for (j=l; j; j--)
! 1994: {
! 1995: if (j == i) continue;
! 1996: Q = (GEN)L[j]; if (typ(Q) != t_MAT) break;
! 1997: inter = inter? FpM_intersect(inter, Q, p): Q;
! 1998: }
! 1999: if (!inter)
! 2000: {
! 2001: setlg(L, l);
! 2002: inter = idealprodprime(nf, L);
! 2003: }
! 2004: else
! 2005: {
! 2006: inter = hnfmodid(inter, p);
! 2007: for ( ; j; j--)
! 2008: {
! 2009: Q = (GEN)L[j];
! 2010: inter = idealmul(nf, inter, Q);
! 2011: }
! 2012: }
1.1 noro 2013:
1.2 ! noro 2014: /* inter = prod L[j], j != i */
! 2015: P2 = compute_pr2(nf, P, p, &ramif);
! 2016: u = addone_aux2(P2, inter);
! 2017: v = unnf_minus_x(u);
! 2018: if (!ramif) pi = gmul(p, v);
! 2019: else
! 2020: {
! 2021: l = lg(P)-1;
! 2022: for (j=l; j; j--)
! 2023: if (!hnf_invimage(P2, (GEN)P[j])) break;
! 2024: pi = element_muli(nf,(GEN)P[j],v);
! 2025: }
! 2026: pi = gadd(pi, u);
! 2027: pi = centermod(pi, p);
! 2028: if (!ramif && hnf_invimage(P2, pi)) pi[1] = laddii((GEN)pi[1], p);
! 2029: return pi;
1.1 noro 2030: }
2031:
1.2 ! noro 2032: /* Input: an ideal mod p, P != Z_K (Fp-basis, in matrix form)
! 2033: * Output: x such that P = (p,x) */
! 2034: static GEN
! 2035: uniformizer(GEN nf, GEN P, GEN p)
1.1 noro 2036: {
1.2 ! noro 2037: GEN q, D, Dp, w, beta, a, T = (GEN)nf[1];
! 2038: long i, f, N = degpol(T), m = lg(P)-1;
! 2039:
! 2040: if (!m) return gscalcol_i(p,N);
1.1 noro 2041:
1.2 ! noro 2042: /* we want v_p(Norm(beta)) = p^f, f = N-m */
! 2043: f = N-m;
! 2044: P = centermod(P, p);
! 2045: q = mulii(gpowgs(p,f), p);
! 2046:
! 2047: if (typ(nf[5]) == t_VEC) /* dummy nf from padicff */
! 2048: {
! 2049: GEN M = gmael(nf,5,1);
! 2050: long ex, prec = gprecision(M);
! 2051: ex = gexpo(M) + gexpo(mulsi(8 * N, p));
! 2052: /* enough prec to use norm_by_embed */
! 2053: if (N * ex <= bit_accuracy(prec))
! 2054: return random_unif_loop_vec(nf, P, p, q);
! 2055: }
! 2056:
! 2057: w = Q_remove_denom((GEN)nf[7], &D);
! 2058: if (D)
! 2059: {
! 2060: long v = pvaluation(D, p, &a);
! 2061: D = gpowgs(p, v);
! 2062: Dp = mulii(D, p);
! 2063: } else {
! 2064: w = dummycopy(w);
! 2065: Dp = p;
! 2066: }
! 2067: for (i=1; i<=N; i++)
! 2068: w[i] = (long)FpX_red((GEN)w[i], Dp); /* w[i] / D defined mod p */
! 2069: beta = gmul(w, P);
! 2070: a = random_unif_loop_pol(nf,P,D,Dp,beta,T,p,q);
! 2071: a = algtobasis_i(nf,a);
! 2072: if (D) a = gdivexact(a,D);
! 2073: a = centermod(a, p);
! 2074: if (!is_uniformizer(D,gmul(w,a), T,q)) a[1] = laddii((GEN)a[1],p);
! 2075: return a;
1.1 noro 2076: }
2077:
1.2 ! noro 2078: /* Assuming P = (p,u) prime, return tau such that p Z + tau Z = p P^(-1)*/
1.1 noro 2079: static GEN
1.2 ! noro 2080: anti_uniformizer(GEN nf, GEN p, GEN u)
1.1 noro 2081: {
1.2 ! noro 2082: gpmem_t av = avma;
! 2083: GEN mat = eltmul_get_table(nf, u);
! 2084: return gerepileupto(av, FpM_deplin(mat,p));
1.1 noro 2085: }
2086:
2087: /*******************************************************************/
2088: /* */
2089: /* BUCHMANN-LENSTRA ALGORITHM */
2090: /* */
2091: /*******************************************************************/
2092:
1.2 ! noro 2093: /* pr = (p,u) of ramification index e */
! 2094: GEN
! 2095: apply_kummer(GEN nf,GEN u,GEN e,GEN p)
! 2096: {
! 2097: GEN t, T = (GEN)nf[1], pr = cgetg(6,t_VEC);
! 2098: long f = degpol(u), N = degpol(T);
! 2099:
! 2100: pr[1] = (long)p;
! 2101: pr[3] = (long)e;
! 2102: pr[4] = lstoi(f);
! 2103: if (f == N) /* inert */
! 2104: {
! 2105: pr[2] = (long)gscalcol_i(p,N);
! 2106: pr[5] = (long)gscalcol_i(gun,N);
! 2107: }
! 2108: else
! 2109: { /* make sure v_pr(u) = 1 (automatic if e>1) */
! 2110: if (is_pm1(e) && !is_uniformizer(NULL,u, T, gpowgs(p,f+1)))
! 2111: u[2] = laddii((GEN)u[2], p);
! 2112: t = algtobasis_i(nf, FpX_div(T,u,p));
! 2113: pr[2] = (long)algtobasis_i(nf,u);
! 2114: pr[5] = (long)centermod(t, p);
! 2115: }
! 2116: return pr;
! 2117: }
! 2118:
! 2119: /* return a Z basis of Z_K's p-radical, phi = x--> x^p-x */
1.1 noro 2120: static GEN
1.2 ! noro 2121: pradical(GEN nf, GEN p, GEN *phi)
1.1 noro 2122: {
2123: long i,N = degpol(nf[1]);
1.2 ! noro 2124: GEN q,m,frob,rad;
1.1 noro 2125:
1.2 ! noro 2126: /* matrix of Frob: x->x^p over Z_K/p */
1.1 noro 2127: frob = cgetg(N+1,t_MAT);
2128: for (i=1; i<=N; i++)
1.2 ! noro 2129: frob[i] = (long)element_powid_mod_p(nf,i,p, p);
1.1 noro 2130:
1.2 ! noro 2131: m = frob; q = p;
! 2132: while (cmpis(q,N) < 0) { q = mulii(q,p); m = FpM_mul(m, frob, p); }
! 2133: rad = FpM_ker(m, p); /* m = Frob^k, s.t p^k >= N */
! 2134: for (i=1; i<=N; i++)
! 2135: coeff(frob,i,i) = lsubis(gcoeff(frob,i,i), 1);
! 2136: *phi = frob; return rad;
! 2137: }
! 2138:
! 2139: /* return powers of a: a^0, ... , a^d, d = dim A */
! 2140: static GEN
! 2141: get_powers(GEN mul, GEN p)
! 2142: {
! 2143: long i, d = lg(mul[1]);
! 2144: GEN z, pow = cgetg(d+2,t_MAT), P = pow+1;
! 2145:
! 2146: P[0] = (long)gscalcol_i(gun, d-1);
! 2147: z = (GEN)mul[1];
! 2148: for (i=1; i<=d; i++)
1.1 noro 2149: {
1.2 ! noro 2150: P[i] = (long)z; /* a^i */
! 2151: if (i!=d) z = FpM_FpV_mul(mul, z, p);
1.1 noro 2152: }
1.2 ! noro 2153: return pow;
1.1 noro 2154: }
2155:
1.2 ! noro 2156: /* minimal polynomial of a in A (dim A = d).
! 2157: * mul = multiplication table by a in A */
1.1 noro 2158: static GEN
1.2 ! noro 2159: pol_min(GEN mul, GEN p)
1.1 noro 2160: {
1.2 ! noro 2161: gpmem_t av = avma;
! 2162: GEN z, pow = get_powers(mul, p);
! 2163: z = FpM_deplin(pow, p);
! 2164: return gerepileupto(av, gtopolyrev(z,0));
1.1 noro 2165: }
2166:
2167: static GEN
1.2 ! noro 2168: get_pr(GEN nf, GEN L, long i, GEN p, int appr)
1.1 noro 2169: {
1.2 ! noro 2170: GEN pr, u, t, P = (GEN)L[i];
! 2171: long e, f;
! 2172: gpmem_t av;
1.1 noro 2173:
1.2 ! noro 2174: if (typ(P) == t_VEC) return P; /* already done (Kummer) */
! 2175:
! 2176: av = avma;
! 2177: if (appr) u = uniformizer_appr(nf, L, i, p);
! 2178: else u = uniformizer(nf, P, p);
! 2179: u = gerepilecopy(av, u);
! 2180: t = anti_uniformizer(nf,p,u);
! 2181: av = avma;
! 2182: e = 1 + int_elt_val(nf,t,p,t,NULL);
! 2183: f = degpol(nf[1]) - (lg(P)-1);
! 2184: avma = av;
! 2185: pr = cgetg(6,t_VEC);
! 2186: pr[1] = (long)p;
! 2187: pr[2] = (long)u;
! 2188: pr[3] = lstoi(e);
! 2189: pr[4] = lstoi(f);
! 2190: pr[5] = (long)t; return pr;
1.1 noro 2191: }
2192:
1.2 ! noro 2193: static int
! 2194: use_appr(GEN L, GEN pp, long N)
1.1 noro 2195: {
1.2 ! noro 2196: long i, f, l = lg(L);
! 2197: double prod = 1., NP;
! 2198: GEN P;
! 2199: gpmem_t av = avma;
1.1 noro 2200:
1.2 ! noro 2201: for (i=1; i<l; i++)
1.1 noro 2202: {
1.2 ! noro 2203: P = (GEN)L[i];
! 2204: f = typ(P)==t_VEC? itos((GEN)P[4]): N - (lg(P)-1);
! 2205: NP = gtodouble(gpowgs(pp, f));
! 2206: prod *= (1 - 1./NP);
1.1 noro 2207: }
1.2 ! noro 2208: avma = av;
! 2209: return (prod * N * 10 < 1);
1.1 noro 2210: }
2211:
1.2 ! noro 2212: /* prime ideal decomposition of p */
1.1 noro 2213: static GEN
1.2 ! noro 2214: _primedec(GEN nf, GEN p)
1.1 noro 2215: {
1.2 ! noro 2216: long i,k,c,iL,N;
! 2217: GEN ex,F,L,Lpr,Ip,H,phi,mat1,T,f,g,h,p1,UN;
! 2218: int appr;
! 2219:
! 2220: if (DEBUGLEVEL>2) (void)timer2();
! 2221: nf = checknf(nf); T = (GEN)nf[1];
! 2222: F = factmod(T,p);
! 2223: ex = (GEN)F[2];
! 2224: F = (GEN)F[1]; F = centerlift(F);
! 2225: if (DEBUGLEVEL>5) msgtimer("factmod");
1.1 noro 2226:
1.2 ! noro 2227: k = lg(F);
! 2228: if (signe(modii((GEN)nf[4],p))) /* p doesn't divide index */
1.1 noro 2229: {
1.2 ! noro 2230: L = cgetg(k,t_VEC);
! 2231: for (i=1; i<k; i++)
! 2232: L[i] = (long)apply_kummer(nf,(GEN)F[i],(GEN)ex[i],p);
! 2233: if (DEBUGLEVEL>5) msgtimer("simple primedec");
! 2234: return L;
! 2235: }
! 2236:
! 2237: g = (GEN)F[1];
! 2238: for (i=2; i<k; i++) g = FpX_mul(g,(GEN)F[i], p);
! 2239: h = FpX_div(T,g,p);
! 2240: f = FpX_red(gdivexact(gsub(gmul(g,h), T), p), p);
! 2241:
! 2242: N = degpol(T);
! 2243: L = cgetg(N+1,t_VEC); iL = 1;
! 2244: for (i=1; i<k; i++)
! 2245: if (is_pm1(ex[i]) || signe(FpX_res(f,(GEN)F[i],p)))
! 2246: L[iL++] = (long)apply_kummer(nf,(GEN)F[i],(GEN)ex[i],p);
! 2247: else /* F[i] | (f,g,h), happens at least once by Dedekind criterion */
! 2248: ex[i] = 0;
! 2249: if (DEBUGLEVEL>2) msgtimer("%ld Kummer factors", iL-1);
! 2250:
! 2251: /* phi matrix of x -> x^p - x in algebra Z_K/p */
! 2252: Ip = pradical(nf,p,&phi);
! 2253: if (DEBUGLEVEL>2) msgtimer("pradical");
! 2254:
! 2255: /* split etale algebra Z_K / (p,Ip) */
! 2256: h = cgetg(N+1,t_VEC);
! 2257: if (iL > 1)
! 2258: { /* split off Kummer factors */
! 2259: GEN mulbeta, beta = NULL;
! 2260: for (i=1; i<k; i++)
! 2261: if (!ex[i]) beta = beta? FpX_mul(beta, (GEN)F[i], p): (GEN)F[i];
! 2262: beta = FpV_red(algtobasis_i(nf,beta), p);
! 2263:
! 2264: mulbeta = FpM_red(eltmul_get_table(nf, beta), p);
! 2265: p1 = concatsp(mulbeta, Ip);
! 2266: /* Fp-base of ideal (Ip, beta) in ZK/p */
! 2267: h[1] = (long)FpM_image(p1, p);
1.1 noro 2268: }
1.2 ! noro 2269: else
! 2270: h[1] = (long)Ip;
! 2271:
! 2272: UN = gscalcol(gun, N);
! 2273: for (c=1; c; c--)
! 2274: { /* Let A:= (Z_K/p) / Ip; try to split A2 := A / Im H ~ Im M2
! 2275: H * ? + M2 * Mi2 = Id_N ==> M2 * Mi2 projector A --> A2 */
! 2276: GEN M, Mi, M2, Mi2, phi2;
! 2277: long dim;
! 2278:
! 2279: H = (GEN)h[c]; k = lg(H)-1;
! 2280: M = FpM_suppl(concatsp(H,UN), p);
! 2281: Mi = FpM_inv(M, p);
! 2282: M2 = vecextract_i(M, k+1,N); /* M = (H|M2) invertible */
! 2283: Mi2 = rowextract_i(Mi,k+1,N);
! 2284: /* FIXME: FpM_mul(,M2) could be done with vecextract_p */
! 2285: phi2 = FpM_mul(Mi2, FpM_mul(phi,M2, p), p);
! 2286: mat1 = FpM_ker(phi2, p);
! 2287: dim = lg(mat1)-1; /* A2 product of 'dim' fields */
! 2288: if (dim > 1)
! 2289: { /* phi2 v = 0 <==> a = M2 v in Ker phi */
! 2290: GEN I,R,r,a,mula,mul2, v = (GEN)mat1[2];
! 2291: long n;
! 2292:
! 2293: a = FpM_FpV_mul(M2,v, p);
! 2294: mula = FpM_red(eltmul_get_table(nf, a), p);
! 2295: mul2 = FpM_mul(Mi2, FpM_mul(mula,M2, p), p);
! 2296: R = rootmod(pol_min(mul2,p), p); /* totally split mod p */
! 2297:
! 2298: n = lg(R)-1;
! 2299: for (i=1; i<=n; i++)
! 2300: {
! 2301: r = lift_intern((GEN)R[i]);
! 2302: I = gaddmat_i(negi(r), mula);
! 2303: h[c++] = (long)FpM_image(concatsp(H, I), p);
! 2304: }
! 2305: if (n == dim)
! 2306: for (i=1; i<=n; i++) { H = (GEN)h[--c]; L[iL++] = (long)H; }
! 2307: }
! 2308: else /* A2 field ==> H maximal, f = N-k = dim(A2) */
! 2309: L[iL++] = (long)H;
! 2310: }
! 2311: setlg(L, iL);
! 2312: Lpr = cgetg(iL, t_VEC);
! 2313: if (DEBUGLEVEL>2) msgtimer("splitting %ld factors",iL-1);
! 2314: appr = use_appr(L,p,N);
! 2315: if (DEBUGLEVEL>2 && appr) fprintferr("using the approximation theorem\n");
! 2316: for (i=1; i<iL; i++) Lpr[i] = (long)get_pr(nf, L, i, p, appr);
! 2317: if (DEBUGLEVEL>2) msgtimer("finding uniformizers");
! 2318: return Lpr;
! 2319: }
! 2320:
! 2321: GEN
! 2322: primedec(GEN nf, GEN p)
! 2323: {
! 2324: gpmem_t av = avma;
! 2325: return gerepileupto(av, gen_sort(_primedec(nf,p), 0, cmp_prime_over_p));
1.1 noro 2326: }
2327:
1.2 ! noro 2328: /* return [Fp[x]: Fp] */
! 2329: static long
! 2330: ffdegree(GEN x, GEN frob, GEN p)
1.1 noro 2331: {
1.2 ! noro 2332: gpmem_t av = avma;
! 2333: long d, f = lg(frob)-1;
! 2334: GEN y = x;
1.1 noro 2335:
1.2 ! noro 2336: for (d=1; d < f; d++)
1.1 noro 2337: {
1.2 ! noro 2338: y = FpM_FpV_mul(frob, y, p);
! 2339: if (gegal(y, x)) break;
1.1 noro 2340: }
1.2 ! noro 2341: avma = av; return d;
1.1 noro 2342: }
2343:
2344: static GEN
1.2 ! noro 2345: lift_to_zk(GEN v, GEN c, long N)
1.1 noro 2346: {
1.2 ! noro 2347: GEN w = zerocol(N);
! 2348: long i, l = lg(c);
! 2349: for (i=1; i<l; i++) w[c[i]] = v[i];
! 2350: return w;
1.1 noro 2351: }
2352:
1.2 ! noro 2353: /* return integral x = 0 mod p/pr^e, (x,pr) = 1.
! 2354: * Don't reduce mod p here: caller may need result mod pr^k */
! 2355: GEN
! 2356: special_anti_uniformizer(GEN nf, GEN pr)
! 2357: {
! 2358: GEN p = (GEN)pr[1], e = (GEN)pr[3];
! 2359: return gdivexact(element_pow(nf,(GEN)pr[5],e), gpowgs(p,itos(e)-1));
! 2360: }
1.1 noro 2361:
1.2 ! noro 2362: /* return t = 1 mod pr, t = 0 mod p / pr^e(pr/p) */
1.1 noro 2363: static GEN
1.2 ! noro 2364: anti_uniformizer2(GEN nf, GEN pr)
1.1 noro 2365: {
1.2 ! noro 2366: GEN p = (GEN)pr[1], z;
! 2367: z = gmod(special_anti_uniformizer(nf, pr), p);
! 2368: z = eltmul_get_table(nf, z);
! 2369: z = hnfmodid(z, p);
! 2370: z = idealaddtoone_i(nf,pr,z);
! 2371: return unnf_minus_x(z);
! 2372: }
1.1 noro 2373:
1.2 ! noro 2374: #define mpr_TAU 1
! 2375: #define mpr_FFP 2
! 2376: #define mpr_PR 3
! 2377: #define mpr_T 4
! 2378: #define mpr_NFP 5
! 2379: #define SMALLMODPR 4
! 2380: #define LARGEMODPR 6
! 2381: static GEN
! 2382: modpr_TAU(GEN modpr)
! 2383: {
! 2384: GEN tau = (GEN)modpr[mpr_TAU];
! 2385: if (typ(tau) == t_INT && signe(tau) == 0) tau = NULL;
! 2386: return tau;
! 2387: }
1.1 noro 2388:
1.2 ! noro 2389: /* prh = HNF matrix, which is identity but for the first line. Return a
! 2390: * projector to ZK / prh ~ Z/prh[1,1] */
! 2391: GEN
! 2392: dim1proj(GEN prh)
! 2393: {
! 2394: long i, N = lg(prh)-1;
! 2395: GEN ffproj = cgetg(N+1, t_VEC);
! 2396: GEN x, q = gcoeff(prh,1,1);
! 2397: ffproj[1] = un;
! 2398: for (i=2; i<=N; i++)
1.1 noro 2399: {
1.2 ! noro 2400: x = gcoeff(prh,1,i);
! 2401: if (signe(x)) x = subii(q,x);
! 2402: ffproj[i] = (long)x;
1.1 noro 2403: }
1.2 ! noro 2404: return ffproj;
1.1 noro 2405: }
1.2 ! noro 2406:
1.1 noro 2407: static GEN
1.2 ! noro 2408: modprinit(GEN nf, GEN pr, int zk)
1.1 noro 2409: {
1.2 ! noro 2410: gpmem_t av = avma;
! 2411: GEN res, tau, mul, x, p, T, pow, ffproj, nfproj, prh, c, gf;
! 2412: long N, i, k, f;
! 2413:
! 2414: nf = checknf(nf); checkprimeid(pr);
! 2415: gf = (GEN)pr[4];
! 2416: f = itos(gf);
! 2417: N = degpol(nf[1]);
! 2418: prh = prime_to_ideal(nf, pr);
! 2419: tau = zk? gzero: anti_uniformizer2(nf, pr);
! 2420: p = (GEN)pr[1];
1.1 noro 2421:
1.2 ! noro 2422: if (f == 1)
! 2423: {
! 2424: res = cgetg(SMALLMODPR, t_COL);
! 2425: res[mpr_TAU] = (long)tau;
! 2426: res[mpr_FFP] = (long)dim1proj(prh);
! 2427: res[mpr_PR ] = (long)pr; return gerepilecopy(av, res);
! 2428: }
1.1 noro 2429:
1.2 ! noro 2430: c = cgetg(f+1, t_VECSMALL);
! 2431: ffproj = cgetg(N+1, t_MAT);
! 2432: for (k=i=1; i<=N; i++)
! 2433: {
! 2434: x = gcoeff(prh, i,i);
! 2435: if (!is_pm1(x)) { c[k] = i; ffproj[i] = (long)_ei(N, i); k++; }
! 2436: else
! 2437: ffproj[i] = lneg((GEN)prh[i]);
! 2438: }
! 2439: ffproj = rowextract_p(ffproj, c);
! 2440: if (! divise((GEN)nf[4], p))
! 2441: {
! 2442: GEN basis, proj_modT;
! 2443: if (N == f) T = (GEN)nf[1]; /* pr inert */
! 2444: else
! 2445: T = primpart(gmul((GEN)nf[7], (GEN)pr[2]));
1.1 noro 2446:
1.2 ! noro 2447: basis = (GEN)nf[7];
! 2448: proj_modT = cgetg(f+1, t_MAT);
! 2449: for (i=1; i<=f; i++)
! 2450: {
! 2451: GEN cx, w = Q_primitive_part((GEN)basis[c[i]], &cx);
! 2452: w = FpX_rem(w, T, p);
! 2453: if (cx)
1.1 noro 2454: {
1.2 ! noro 2455: cx = gmod(cx, p);
! 2456: w = FpX_Fp_mul(w, cx, p);
! 2457: }
! 2458: proj_modT[i] = (long)pol_to_vec(w, f); /* w_i mod (T,p) */
1.1 noro 2459: }
1.2 ! noro 2460: ffproj = FpM_mul(proj_modT,ffproj,p);
! 2461:
! 2462: res = cgetg(SMALLMODPR+1, t_COL);
! 2463: res[mpr_TAU] = (long)tau;
! 2464: res[mpr_FFP] = (long)ffproj;
! 2465: res[mpr_PR ] = (long)pr;
! 2466: res[mpr_T ] = (long)T; return gerepilecopy(av, res);
1.1 noro 2467: }
2468:
1.2 ! noro 2469: if (isprime(gf))
! 2470: {
! 2471: mul = eltmulid_get_table(nf, c[2]);
! 2472: mul = vecextract_p(mul, c);
! 2473: }
! 2474: else
! 2475: {
! 2476: GEN v, u, u2, frob;
! 2477: long deg,deg1,deg2;
1.1 noro 2478:
1.2 ! noro 2479: /* matrix of Frob: x->x^p over Z_K/pr = < w[c1], ..., w[cf] > over Fp */
! 2480: frob = cgetg(f+1, t_MAT);
! 2481: for (i=1; i<=f; i++)
1.1 noro 2482: {
1.2 ! noro 2483: x = element_powid_mod_p(nf,c[i],p, p);
! 2484: frob[i] = (long)FpM_FpV_mul(ffproj, x, p);
1.1 noro 2485: }
1.2 ! noro 2486: u = _ei(f,2); k = 2;
! 2487: deg1 = ffdegree(u, frob, p);
! 2488: while (deg1 < f)
1.1 noro 2489: {
1.2 ! noro 2490: k++; u2 = _ei(f, k);
! 2491: deg2 = ffdegree(u2, frob, p);
! 2492: deg = clcm(deg1,deg2);
! 2493: if (deg == deg1) continue;
! 2494: if (deg == deg2) { deg1 = deg2; u = u2; continue; }
! 2495: u = gadd(u, u2);
! 2496: while (ffdegree(u, frob, p) < deg) u = gadd(u, u2);
! 2497: deg1 = deg;
1.1 noro 2498: }
1.2 ! noro 2499: v = lift_to_zk(u,c,N);
! 2500:
! 2501: mul = cgetg(f+1,t_MAT);
! 2502: mul[1] = (long)v; /* assume w_1 = 1 */
! 2503: for (i=2; i<=f; i++) mul[i] = (long)element_mulid(nf,v,c[i]);
1.1 noro 2504: }
2505:
1.2 ! noro 2506: /* Z_K/pr = Fp(v), mul = mul by v */
! 2507: mul = FpM_red(mul, p);
! 2508: mul = FpM_mul(ffproj, mul, p);
! 2509:
! 2510: pow = get_powers(mul, p);
! 2511: T = gtopolyrev(FpM_deplin(pow, p), varn(nf[1]));
! 2512: nfproj = cgetg(f+1, t_MAT);
! 2513: for (i=1; i<=f; i++) nfproj[i] = (long)lift_to_zk((GEN)pow[i], c, N);
! 2514: nfproj = gmul((GEN)nf[7], nfproj);
1.1 noro 2515:
1.2 ! noro 2516: setlg(pow, f+1);
! 2517: ffproj = FpM_mul(FpM_inv(pow, p), ffproj, p);
1.1 noro 2518:
1.2 ! noro 2519: res = cgetg(LARGEMODPR, t_COL);
! 2520: res[mpr_TAU] = (long)tau;
! 2521: res[mpr_FFP] = (long)ffproj;
! 2522: res[mpr_PR ] = (long)pr;
! 2523: res[mpr_T ] = (long)T;
! 2524: res[mpr_NFP] = (long)nfproj; return gerepilecopy(av, res);
! 2525: }
1.1 noro 2526:
1.2 ! noro 2527: GEN
! 2528: nfmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 0); }
! 2529: GEN
! 2530: zkmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 1); }
1.1 noro 2531:
1.2 ! noro 2532: void
! 2533: checkmodpr(GEN modpr)
! 2534: {
! 2535: if (typ(modpr) != t_COL || lg(modpr) < SMALLMODPR)
! 2536: err(talker,"incorrect modpr format");
! 2537: checkprimeid((GEN)modpr[mpr_PR]);
1.1 noro 2538: }
2539:
1.2 ! noro 2540:
1.1 noro 2541: static GEN
1.2 ! noro 2542: to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p, int zk)
1.1 noro 2543: {
1.2 ! noro 2544: GEN modpr = (typ(*pr) == t_COL)? *pr: modprinit(nf, *pr, zk);
! 2545: *T = lg(modpr)==SMALLMODPR? NULL: (GEN)modpr[mpr_T];
! 2546: *pr = (GEN)modpr[mpr_PR];
! 2547: *p = (GEN)(*pr)[1]; return modpr;
! 2548: }
! 2549: GEN
! 2550: nf_to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
! 2551: GEN modpr = to_ff_init(nf,pr,T,p,0);
! 2552: GEN tau = modpr_TAU(modpr);
! 2553: if (!tau) modpr[mpr_TAU] = (long)anti_uniformizer2(nf, *pr);
! 2554: return modpr;
! 2555: }
! 2556: GEN
! 2557: zk_to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
! 2558: return to_ff_init(nf,pr,T,p,1);
1.1 noro 2559: }
2560:
1.2 ! noro 2561: /* assume x in 'basis' form (t_COL) */
1.1 noro 2562: GEN
1.2 ! noro 2563: zk_to_ff(GEN x, GEN modpr)
1.1 noro 2564: {
1.2 ! noro 2565: GEN pr = (GEN)modpr[mpr_PR];
! 2566: GEN p = (GEN)pr[1];
! 2567: GEN y = gmul((GEN)modpr[mpr_FFP], x);
! 2568: if (lg(modpr) == SMALLMODPR) return modii(y,p);
! 2569: y = FpV_red(y, p);
! 2570: return col_to_ff(y, varn(modpr[mpr_T]));
! 2571: }
! 2572:
! 2573: /* REDUCTION Modulo a prime ideal */
1.1 noro 2574:
1.2 ! noro 2575: /* assume x in t_COL form, v_pr(x) >= 0 */
! 2576: static GEN
! 2577: kill_denom(GEN x, GEN nf, GEN p, GEN modpr)
! 2578: {
! 2579: GEN cx, den = denom(x);
! 2580: long v;
! 2581: if (gcmp1(den)) return x;
1.1 noro 2582:
1.2 ! noro 2583: v = ggval(den,p);
! 2584: if (v)
1.1 noro 2585: {
1.2 ! noro 2586: GEN tau = modpr_TAU(modpr);
! 2587: if (!tau) err(talker,"modpr initialized for integers only!");
! 2588: x = element_mul(nf,x, element_pow(nf,tau,stoi(v)));
1.1 noro 2589: }
1.2 ! noro 2590: x = Q_primitive_part(x, &cx);
! 2591: x = FpV_red(x, p);
! 2592: if (cx) x = FpV_red(gmul(gmod(cx, p), x), p);
! 2593: return x;
1.1 noro 2594: }
2595:
2596: /* x integral, reduce mod prh in HNF */
2597: GEN
2598: nfreducemodpr_i(GEN x, GEN prh)
2599: {
2600: GEN p = gcoeff(prh,1,1);
2601: long i,j;
2602:
2603: x = dummycopy(x);
2604: for (i=lg(x)-1; i>=2; i--)
2605: {
2606: GEN t = (GEN)prh[i], p1 = resii((GEN)x[i], p);
2607: x[i] = (long)p1;
2608: if (signe(p1) && is_pm1(t[i]))
2609: {
2610: for (j=1; j<i; j++)
2611: x[j] = lsubii((GEN)x[j], mulii(p1, (GEN)t[j]));
2612: x[i] = zero;
2613: }
2614: }
2615: x[1] = lresii((GEN)x[1], p); return x;
2616: }
2617:
2618: GEN
1.2 ! noro 2619: nfreducemodpr(GEN nf, GEN x, GEN modpr)
1.1 noro 2620: {
1.2 ! noro 2621: gpmem_t av = avma;
! 2622: long i;
! 2623: GEN pr, p;
1.1 noro 2624:
1.2 ! noro 2625: checkmodpr(modpr);
! 2626: pr = (GEN)modpr[mpr_PR];
! 2627: p = (GEN)pr[1];
! 2628: if (typ(x) != t_COL) x = algtobasis(nf,x);
1.1 noro 2629: for (i=lg(x)-1; i>0; i--)
1.2 ! noro 2630: if (typ(x[i]) == t_INTMOD) { x = lift(x); break; }
! 2631: x = kill_denom(x, nf, p, modpr);
! 2632: x = ff_to_nf(zk_to_ff(x,modpr), modpr);
! 2633: if (typ(x) != t_COL) x = algtobasis(nf, x);
! 2634: return gerepileupto(av, FpV(x, p));
! 2635: }
! 2636:
! 2637: GEN
! 2638: nf_to_ff(GEN nf, GEN x, GEN modpr)
! 2639: {
! 2640: gpmem_t av = avma;
! 2641: GEN pr = (GEN)modpr[mpr_PR];
! 2642: GEN p = (GEN)pr[1];
! 2643: long t = typ(x);
! 2644:
! 2645: if (t == t_POLMOD) { x = (GEN)x[2]; t = typ(x); }
! 2646: switch(t)
! 2647: {
! 2648: case t_INT: return modii(x, p);
! 2649: case t_FRAC: return gmod(x, p);
! 2650: case t_POL: x = algtobasis(nf, x); break;
! 2651: case t_COL: break;
! 2652: default: err(typeer,"nf_to_ff");
! 2653: }
! 2654: x = kill_denom(x, nf, p, modpr);
! 2655: return gerepilecopy(av, zk_to_ff(x, modpr));
! 2656: }
! 2657:
! 2658: GEN
! 2659: ff_to_nf(GEN x, GEN modpr)
! 2660: {
! 2661: if (lg(modpr) < LARGEMODPR) return x;
! 2662: return mulmat_pol((GEN)modpr[mpr_NFP], x);
! 2663: }
! 2664: GEN
! 2665: modprM_lift(GEN x, GEN modpr)
! 2666: {
! 2667: long i,j,h,l = lg(x);
! 2668: GEN y = cgetg(l, t_MAT);
! 2669: if (l == 1) return y;
! 2670:
! 2671: h = lg(x[1]);
! 2672: for (j=1; j<l; j++)
1.1 noro 2673: {
1.2 ! noro 2674: GEN p1 = cgetg(h,t_COL); y[j] = (long)p1;
! 2675: for (i=1; i<h; i++) p1[i]=(long)ff_to_nf(gcoeff(x,i,j), modpr);
1.1 noro 2676: }
1.2 ! noro 2677: return y;
1.1 noro 2678: }
1.2 ! noro 2679: GEN
! 2680: modprX_lift(GEN x, GEN modpr)
! 2681: {
! 2682: long i, l;
! 2683: GEN z;
! 2684:
! 2685: if (typ(x)!=t_POL) return gcopy(x); /* scalar */
! 2686: l = lgef(x);
! 2687: z = cgetg(l, t_POL); z[1] = x[1];
! 2688: for (i=2; i<l; i++) z[i] = (long)ff_to_nf((GEN)x[i], modpr);
! 2689: return z;
! 2690: }
! 2691:
! 2692: /* reduce the coefficients of pol modulo modpr */
! 2693: GEN
! 2694: modprX(GEN x, GEN nf,GEN modpr)
! 2695: {
! 2696: long i, l;
! 2697: GEN z;
1.1 noro 2698:
1.2 ! noro 2699: if (typ(x)!=t_POL) return nf_to_ff(nf,x,modpr);
! 2700: l = lgef(x);
! 2701: z = cgetg(l,t_POL); z[1] = x[1];
! 2702: for (i=2; i<l; i++) z[i] = (long)nf_to_ff(nf,(GEN)x[i],modpr);
! 2703: return normalizepol(z);
! 2704: }
! 2705: GEN
! 2706: modprV(GEN z, GEN nf,GEN modpr)
! 2707: {
! 2708: long i,l = lg(z);
! 2709: GEN x;
! 2710: x = cgetg(l,typ(z));
! 2711: for (i=1; i<l; i++) x[i] = (long)nf_to_ff(nf,(GEN)z[i], modpr);
! 2712: return x;
! 2713: }
! 2714: /* assume z a t_VEC/t_COL/t_MAT */
1.1 noro 2715: GEN
1.2 ! noro 2716: modprM(GEN z, GEN nf,GEN modpr)
1.1 noro 2717: {
1.2 ! noro 2718: long i,l = lg(z), m;
! 2719: GEN x;
! 2720:
! 2721: if (typ(z) != t_MAT) return modprV(z,nf,modpr);
! 2722: x = cgetg(l,t_MAT); if (l==1) return x;
! 2723: m = lg((GEN)z[1]);
! 2724: for (i=1; i<l; i++) x[i] = (long)modprV((GEN)z[i],nf,modpr);
! 2725: return x;
1.1 noro 2726: }
2727:
2728: /* relative ROUND 2
2729: *
2730: * input: nf = base field K
2731: * x monic polynomial, coefficients in Z_K, degree n defining a relative
2732: * extension L=K(theta).
2733: * One MUST have varn(x) < varn(nf[1]).
2734: * output: a pseudo-basis [A,I] of Z_L, where A is in M_n(K) in HNF form and
2735: * I a vector of n ideals.
2736: */
2737:
2738: /* given MODULES x and y by their pseudo-bases in HNF, gives a
1.2 ! noro 2739: * pseudo-basis of the module generated by x and y. */
1.1 noro 2740: static GEN
2741: rnfjoinmodules(GEN nf, GEN x, GEN y)
2742: {
2743: long i,lx,ly;
2744: GEN p1,p2,z,Hx,Hy,Ix,Iy;
2745:
2746: if (x == NULL) return y;
1.2 ! noro 2747: Hx = (GEN)x[1]; lx = lg(Hx); Ix = (GEN)x[2];
! 2748: Hy = (GEN)y[1]; ly = lg(Hy); Iy = (GEN)y[2];
1.1 noro 2749: i = lx+ly-1;
2750: z = (GEN)gpmalloc(sizeof(long*)*(3+2*i));
2751: *z = evaltyp(t_VEC)|evallg(3);
2752: p1 = z+3; z[1]=(long)p1; *p1 = evaltyp(t_MAT)|evallg(i);
2753: p2 = p1+i; z[2]=(long)p2; *p2 = evaltyp(t_VEC)|evallg(i);
2754:
1.2 ! noro 2755: for (i=1; i<lx; i++) { p1[i] = Hx[i]; p2[i] = Ix[i]; }
! 2756: p1 += lx-1;
! 2757: p2 += lx-1;
! 2758: for (i=1; i<ly; i++) { p1[i] = Hy[i]; p2[i] = Iy[i]; }
1.1 noro 2759: x = nfhermite(nf,z); free(z); return x;
2760: }
2761:
1.2 ! noro 2762: typedef struct {
! 2763: GEN nf, multab, modpr,T,p;
! 2764: long h;
! 2765: } rnfeltmod_muldata;
! 2766:
1.1 noro 2767: static GEN
1.2 ! noro 2768: _mul(void *data, GEN x, GEN y/* base; ignored */)
1.1 noro 2769: {
1.2 ! noro 2770: rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
! 2771: GEN z = x? element_mulid(D->multab,x,D->h)
! 2772: : element_mulidid(D->multab,D->h,D->h);
! 2773: (void)y;
! 2774: return FpVQX_red(z,D->T,D->p);
1.1 noro 2775: }
2776: static GEN
1.2 ! noro 2777: _sqr(void *data, GEN x)
1.1 noro 2778: {
1.2 ! noro 2779: rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
! 2780: GEN z = x? sqr_by_tab(D->multab,x)
! 2781: : element_mulidid(D->multab,D->h,D->h);
! 2782: return FpVQX_red(z,D->T,D->p);
1.1 noro 2783: }
2784:
1.2 ! noro 2785: /* Compute x^n mod pr in the extension, assume n >= 0 */
1.1 noro 2786: static GEN
1.2 ! noro 2787: rnfelementid_powmod(GEN multab, long h, GEN n, GEN T, GEN p)
1.1 noro 2788: {
1.2 ! noro 2789: gpmem_t av = avma;
! 2790: GEN y;
! 2791: rnfeltmod_muldata D;
! 2792:
! 2793: if (!signe(n)) return gun;
1.1 noro 2794:
1.2 ! noro 2795: D.multab = multab;
! 2796: D.h = h;
! 2797: D.T = T;
! 2798: D.p = p;
! 2799: y = leftright_pow(NULL, n, (void*)&D, &_sqr, &_mul);
1.1 noro 2800: return gerepilecopy(av, y);
2801: }
2802:
1.2 ! noro 2803: #if 0
1.1 noro 2804: GEN
2805: mymod(GEN x, GEN p)
2806: {
2807: long i,lx, tx = typ(x);
2808: GEN y,p1;
2809:
2810: if (tx == t_INT) return resii(x,p);
2811: if (tx == t_FRAC)
2812: {
2813: p1 = resii((GEN)x[2], p);
2814: if (p1 != gzero) { cgiv(p1); return gmod(x,p); }
2815: return x;
2816: }
2817: if (!is_matvec_t(tx))
2818: err(bugparier, "mymod (missing type)");
2819: lx = lg(x); y = cgetg(lx,tx);
2820: for (i=1; i<lx; i++) y[i] = (long)mymod((GEN)x[i],p);
2821: return y;
2822: }
1.2 ! noro 2823: #endif
! 2824:
! 2825: #define FqX_div(x,y,T,p) FpXQX_divres((x),(y),(T),(p),NULL)
! 2826: #define FqX_mul FpXQX_mul
! 2827:
! 2828: /* relative Dedekind criterion over nf, applied to the order defined by a
! 2829: * root of irreducible polynomial P, modulo the prime ideal pr. Returns
! 2830: * [flag,basis,val], where basis is a pseudo-basis of the enlarged order,
! 2831: * flag is 1 iff this order is pr-maximal, and val is the valuation at pr of
! 2832: * the order discriminant */
! 2833: static GEN
! 2834: rnfdedekind_i(GEN nf,GEN P,GEN pr, GEN disc)
! 2835: {
! 2836: long vt, r, d, n, m, i, j;
! 2837: gpmem_t av = avma;
! 2838: GEN Prd,p1,p2,p,tau,g,matid;
! 2839: GEN modpr,res,h,k,base,nfT,T,gzk,hzk;
! 2840:
! 2841: if (!disc) disc = discsr(P);
! 2842: P = lift(P);
! 2843: nf = checknf(nf); nfT = (GEN)nf[1];
! 2844: res = cgetg(4,t_VEC);
! 2845: modpr = nf_to_ff_init(nf,&pr, &T, &p);
! 2846: tau = gmul((GEN)nf[7], (GEN)pr[5]);
! 2847: n = degpol(nfT);
! 2848: m = degpol(P);
! 2849:
! 2850: Prd = modprX(P, nf, modpr);
! 2851: p1 = (GEN)FqX_factor(Prd,T,p)[1];
! 2852: r = lg(p1); if (r < 2) err(constpoler,"rnfdedekind");
! 2853: g = (GEN)p1[1];
! 2854: for (i=2; i<r; i++) g = FqX_mul(g, (GEN)p1[i], T, p);
! 2855: h = FqX_div(Prd,g, T, p);
! 2856: gzk = modprX_lift(g, modpr);
! 2857: hzk = modprX_lift(h, modpr);
! 2858:
! 2859: k = gsub(P, RXQX_mul(gzk,hzk, nfT));
! 2860: k = gdiv(RXQX_mul(tau,k,nfT), p);
! 2861: k = modprX(k, nf, modpr);
! 2862:
! 2863: p2 = FqX_gcd(g,h, T,p);
! 2864: k = FqX_gcd(p2,k, T,p); d = degpol(k); /* <= m */
! 2865:
! 2866: vt = idealval(nf,disc,pr) - 2*d;
! 2867: res[3] = lstoi(vt);
! 2868: res[1] = (!d || vt<=1)? un: zero;
! 2869:
! 2870: base = cgetg(3,t_VEC);
! 2871: p1 = cgetg(m+d+1,t_MAT); base[1] = (long)p1;
! 2872: p2 = cgetg(m+d+1,t_VEC); base[2] = (long)p2;
! 2873: /* if d > 0, base[2] temporarily multiplied by p, for the final nfhermitemod
! 2874: * [ which requires integral ideals ] */
! 2875: matid = gscalmat(d? p: gun, n);
! 2876: for (j=1; j<=m; j++)
! 2877: {
! 2878: p1[j] = (long)_ei(m, j);
! 2879: p2[j] = (long)matid;
! 2880: }
! 2881: if (d)
! 2882: {
! 2883: GEN prinvp = pidealprimeinv(nf,pr); /* again multiplied by p */
! 2884: GEN X = polx[varn(P)], pal;
! 2885:
! 2886: pal = FqX_div(Prd,k, T,p);
! 2887: pal = modprX_lift(pal, modpr);
! 2888: for ( ; j<=m+d; j++)
! 2889: {
! 2890: p1[j] = (long)pol_to_vec(pal,m);
! 2891: p2[j] = (long)prinvp;
! 2892: pal = RXQX_rem(RXQX_mul(pal,X,nfT),P,nfT);
! 2893: }
! 2894: /* the modulus is integral */
! 2895: base = nfhermitemod(nf,base, gmul(gpowgs(p, m-d),
! 2896: idealpows(nf, prinvp, d)));
! 2897: base[2] = ldiv((GEN)base[2], p); /* cancel the factor p */
! 2898: }
! 2899: res[2] = (long)base; return gerepilecopy(av, res);
! 2900: }
! 2901:
! 2902: GEN
! 2903: rnfdedekind(GEN nf,GEN P0,GEN pr)
! 2904: {
! 2905: return rnfdedekind_i(nf,P0,pr,NULL);
! 2906: }
1.1 noro 2907:
2908: static GEN
1.2 ! noro 2909: rnfordmax(GEN nf, GEN pol, GEN pr, GEN disc)
1.1 noro 2910: {
1.2 ! noro 2911: gpmem_t av=avma,av1,lim;
! 2912: long i,j,k,n,v1,v2,vpol,m,cmpt,sep;
! 2913: GEN p,T,q,q1,modpr,A,Aa,Aaa,A1,I,R,p1,p2,p3,multab,multabmod,Aainv;
1.1 noro 2914: GEN pip,baseIp,baseOp,alpha,matprod,alphainv,matC,matG,vecpro,matH;
1.2 ! noro 2915: GEN neworder,H,Hid,alphalistinv,alphalist,prhinv,nfT,id,rnfId;
1.1 noro 2916:
2917: if (DEBUGLEVEL>1) fprintferr(" treating %Z\n",pr);
1.2 ! noro 2918: modpr = nf_to_ff_init(nf,&pr,&T,&p);
! 2919: p1 = rnfdedekind_i(nf,pol,modpr,disc);
1.1 noro 2920: if (gcmp1((GEN)p1[1])) return gerepilecopy(av,(GEN)p1[2]);
2921:
1.2 ! noro 2922: sep = itos((GEN)p1[3]);
! 2923: A = gmael(p1,2,1);
! 2924: I = gmael(p1,2,2);
! 2925:
! 2926: pip = basistoalg(nf, (GEN)pr[2]);
! 2927: nfT = (GEN)nf[1];
! 2928: n = degpol(pol); vpol = varn(pol);
! 2929: q = T? gpowgs(p,degpol(T)): p;
! 2930: q1 = q; while (cmpis(q1,n) < 0) q1 = mulii(q1,q);
! 2931: rnfId = idmat(degpol(pol));
! 2932: id = idmat(degpol(nfT));
! 2933:
! 2934: multab = cgetg(n*n+1,t_MAT);
! 2935: for (j=1; j<=n*n; j++) multab[j] = lgetg(n+1,t_COL);
! 2936: prhinv = idealinv(nf, pr);
! 2937: alphalistinv = cgetg(n+1,t_VEC);
! 2938: alphalist = cgetg(n+1,t_VEC);
! 2939: A1 = cgetg(n+1,t_MAT);
! 2940: av1 = avma; lim = stack_lim(av1,1);
1.1 noro 2941: for(cmpt=1; ; cmpt++)
2942: {
1.2 ! noro 2943: if (DEBUGLEVEL>1) fprintferr(" %ld%s pass\n",cmpt,eng_ord(cmpt));
! 2944: for (j=1; j<=n; j++)
1.1 noro 2945: {
1.2 ! noro 2946: if (gegal((GEN)I[j],id))
! 2947: {
! 2948: alphalist[j] = alphalistinv[j] = un;
! 2949: A1[j] = A[j]; continue;
! 2950: }
! 2951:
! 2952: p1 = ideal_two_elt(nf,(GEN)I[j]);
! 2953: if (gcmp0((GEN)p1[1]))
! 2954: alpha = (GEN)p1[2];
1.1 noro 2955: else
2956: {
1.2 ! noro 2957: v1 = element_val(nf,(GEN)p1[1],pr);
! 2958: v2 = element_val(nf,(GEN)p1[2],pr);
! 2959: alpha = (v1>v2)? (GEN)p1[2]: (GEN)p1[1];
1.1 noro 2960: }
1.2 ! noro 2961: alphalist[j] = (long)alpha;
! 2962: alphalistinv[j] = (long)element_inv(nf,alpha);
! 2963:
! 2964: p1 = cgetg(n+1,t_COL); A1[j] = (long)p1;
1.1 noro 2965: for (i=1; i<=n; i++)
1.2 ! noro 2966: p1[i] = (long)element_mul(nf,gcoeff(A,i,j),alpha);
1.1 noro 2967: }
1.2 ! noro 2968: Aa = basistoalg(nf,A1);
! 2969: Aainv = lift_intern(ginv(Aa));
! 2970: Aaa = lift_intern(mat_to_vecpol(Aa,vpol));
! 2971: for (i=1; i<=n; i++)
! 2972: for (j=i; j<=n; j++)
1.1 noro 2973: {
1.2 ! noro 2974: long tp;
! 2975: p1 = RXQX_rem(gmul((GEN)Aaa[i],(GEN)Aaa[j]), pol, nfT);
! 2976: tp = typ(p1);
! 2977: if (is_scalar_t(tp) || (tp==t_POL && varn(p1)>vpol))
! 2978: p2 = gmul(p1, (GEN)Aainv[1]);
! 2979: else
! 2980: p2 = mulmat_pol(Aainv, p1);
! 2981: for (k=1; k<=n; k++)
! 2982: {
! 2983: GEN c = gres((GEN)p2[k], nfT);
! 2984: coeff(multab,k,(i-1)*n+j) = (long)c;
! 2985: coeff(multab,k,(j-1)*n+i) = (long)c;
! 2986: }
1.1 noro 2987: }
1.2 ! noro 2988: multabmod = modprM(multab,nf,modpr);
! 2989: R = cgetg(n+1,t_MAT);
! 2990: R[1] = rnfId[1];
! 2991: for (j=2; j<=n; j++)
! 2992: R[j] = (long) rnfelementid_powmod(multabmod,j,q1,T,p);
! 2993: baseIp = FqM_ker(R,T,p);
! 2994: baseOp = lg(baseIp)==1? rnfId: FqM_suppl(baseIp,T,p);
! 2995: alpha = modprM_lift(baseOp, modpr);
! 2996: for (j=1; j<lg(baseIp); j++)
! 2997: {
! 2998: p1 = (GEN)alpha[j];
! 2999: for (i=1; i<=n; i++) p1[i] = (long)to_polmod((GEN)p1[i], nfT);
1.1 noro 3000: }
3001: for ( ; j<=n; j++)
3002: {
1.2 ! noro 3003: p1 = (GEN)alpha[j];
! 3004: for (i=1; i<=n; i++) p1[i] = lmul(pip, (GEN)p1[i]);
1.1 noro 3005: }
1.2 ! noro 3006: alphainv = lift_intern(ginv(alpha));
! 3007: alpha = lift_intern(alpha);
! 3008:
! 3009: matprod = cgetg(n+1,t_MAT);
1.1 noro 3010: for (j=1; j<=n; j++)
3011: {
1.2 ! noro 3012: p1 = cgetg(n+1,t_COL); matprod[j] = (long)p1;
1.1 noro 3013: for (i=1; i<=n; i++)
1.2 ! noro 3014: p1[i] = (long)gmod(element_mulid(multab, (GEN)alpha[i],j), nfT);
1.1 noro 3015: }
3016: matC = cgetg(n+1,t_MAT);
3017: for (j=1; j<=n; j++)
3018: {
1.2 ! noro 3019: p1 = cgetg(n*n+1,t_COL); matC[j] = (long)p1;
1.1 noro 3020: for (i=1; i<=n; i++)
3021: {
3022: p2 = gmul(alphainv, gcoeff(matprod,i,j));
3023: for (k=1; k<=n; k++)
1.2 ! noro 3024: {
! 3025: GEN c = gres((GEN)p2[k], nfT);
! 3026: p1[(i-1)*n+k] = (long)nf_to_ff(nf,c,modpr);
! 3027: }
1.1 noro 3028: }
3029: }
1.2 ! noro 3030: matG = FqM_ker(matC,T,p); m = lg(matG)-1;
! 3031: matG = modprM_lift(matG, modpr);
! 3032: vecpro = cgetg(3,t_VEC);
! 3033: vecpro[1] = (long)concatsp(matG,rnfId);
! 3034:
! 3035: p2 = cgetg(n+m+1,t_VEC);
! 3036: vecpro[2] = (long)p2;
! 3037: for (j=1; j<=m; j++) p2[j] = (long)prhinv;
1.1 noro 3038: p2 += m;
3039: for (j=1; j<=n; j++)
3040: p2[j] = (long)idealmul(nf,(GEN)I[j],(GEN)alphalistinv[j]);
3041:
1.2 ! noro 3042: matH = nfhermite(nf,vecpro);
! 3043: p1 = algtobasis(nf, gmul(Aa, basistoalg(nf,(GEN)matH[1])));
! 3044: p2 = (GEN)matH[2];
! 3045:
! 3046: neworder = cgetg(3,t_VEC);
! 3047: H = cgetg(n+1,t_MAT);
! 3048: Hid = cgetg(n+1,t_VEC);
1.1 noro 3049: for (j=1; j<=n; j++)
3050: {
1.2 ! noro 3051: alphainv = (GEN)alphalistinv[j];
! 3052: if (alphainv == gun) { H[j] = p1[j]; Hid[j] = p2[j]; continue; }
! 3053:
1.1 noro 3054: p3=cgetg(n+1,t_COL); H[j]=(long)p3;
3055: for (i=1; i<=n; i++)
1.2 ! noro 3056: p3[i] = (long)element_mul(nf,gcoeff(p1,i,j),alphainv);
! 3057: Hid[j] = (long)idealmul(nf,(GEN)p2[j],(GEN)alphalist[j]);
1.1 noro 3058: }
1.2 ! noro 3059: if (DEBUGLEVEL>3) { fprintferr(" new order:\n"); outerr(H); outerr(Hid); }
1.1 noro 3060: if (sep == 2 || gegal(I,Hid))
3061: {
1.2 ! noro 3062: neworder[1] = (long)H;
! 3063: neworder[2] = (long)Hid;
! 3064: return gerepilecopy(av, neworder);
1.1 noro 3065: }
3066:
1.2 ! noro 3067: A = H; I = Hid;
1.1 noro 3068: if (low_stack(lim, stack_lim(av1,1)) || (cmpt & 3) == 0)
3069: {
3070: GEN *gptr[2]; gptr[0]=&A; gptr[1]=&I;
3071: if(DEBUGMEM>1) err(warnmem,"rnfordmax");
3072: gerepilemany(av1,gptr,2);
3073: }
3074: }
3075: }
3076:
3077: static void
3078: check_pol(GEN x, long v)
3079: {
1.2 ! noro 3080: long i,tx, lx = lgef(x);
1.1 noro 3081: if (varn(x) != v)
3082: err(talker,"incorrect variable in rnf function");
3083: for (i=2; i<lx; i++)
3084: {
3085: tx = typ(x[i]);
3086: if (!is_scalar_t(tx) || tx == t_POLMOD)
3087: err(talker,"incorrect polcoeff in rnf function");
3088: }
3089: }
3090:
3091: GEN
3092: fix_relative_pol(GEN nf, GEN x, int chk_lead)
3093: {
3094: GEN xnf = (typ(nf) == t_POL)? nf: (GEN)nf[1];
1.2 ! noro 3095: long i, vnf = varn(xnf), lx = lgef(x);
1.1 noro 3096: if (typ(x) != t_POL || varn(x) >= vnf)
3097: err(talker,"incorrect polynomial in rnf function");
3098: x = dummycopy(x);
3099: for (i=2; i<lx; i++)
3100: switch(typ(x[i]))
3101: {
3102: case t_POL:
3103: check_pol((GEN)x[i], vnf);
3104: x[i] = lmodulcp((GEN)x[i], xnf); break;
3105: case t_POLMOD:
3106: if (!gegal(gmael(x,i,1), xnf)) err(consister,"rnf function");
3107: break;
3108: }
3109:
3110: if (chk_lead && !gcmp1(leading_term(x)))
3111: err(impl,"non-monic relative polynomials");
3112: return x;
3113: }
3114:
3115: static GEN
3116: rnfround2all(GEN nf, GEN pol, long all)
3117: {
1.2 ! noro 3118: gpmem_t av = avma;
! 3119: long i,j,n,N,nbidp,vpol,*ep;
! 3120: GEN A,p1,p2,nfT,list,id,I,W,pseudo,y,d,D,sym,unnf,disc;
1.1 noro 3121:
1.2 ! noro 3122: nf=checknf(nf); nfT=(GEN)nf[1]; vpol=varn(pol);
1.1 noro 3123: pol = fix_relative_pol(nf,pol,1);
1.2 ! noro 3124: N = degpol(nfT);
! 3125: n = degpol(pol);
! 3126: disc = discsr(pol); pol = lift(pol);
! 3127: list = idealfactor(nf, disc);
! 3128: ep = (long*)list[2];
! 3129: list= (GEN) list[1];
! 3130: nbidp=lg(list)-1; for(i=1;i<=nbidp;i++) ep[i] = itos((GEN)ep[i]);
1.1 noro 3131: if (DEBUGLEVEL>1)
3132: {
3133: fprintferr("Ideals to consider:\n");
3134: for (i=1; i<=nbidp; i++)
3135: if (ep[i]>1) fprintferr("%Z^%ld\n",list[i],ep[i]);
3136: flusherr();
3137: }
1.2 ! noro 3138: id = idmat(N); unnf = (GEN)id[1];
1.1 noro 3139: pseudo = NULL;
3140: for (i=1; i<=nbidp; i++)
3141: if (ep[i] > 1)
3142: {
1.2 ! noro 3143: y = rnfordmax(nf, pol, (GEN)list[i], disc);
1.1 noro 3144: pseudo = rnfjoinmodules(nf,pseudo,y);
3145: }
1.2 ! noro 3146: if (pseudo)
! 3147: {
! 3148: A = (GEN)pseudo[1];
! 3149: I = (GEN)pseudo[2];
! 3150: }
! 3151: else
1.1 noro 3152: {
1.2 ! noro 3153: I = cgetg(n+1,t_VEC); for (i=1; i<=n; i++) I[i] = (long)id;
! 3154: A = idmat_intern(n, unnf, zerocol(N));
1.1 noro 3155: }
1.2 ! noro 3156: W = mat_to_vecpol(lift_intern(basistoalg(nf,A)), vpol);
! 3157: sym = polsym_gen(pol, NULL, n-1, nfT, NULL);
! 3158: p2 = cgetg(n+1,t_MAT); for (j=1; j<=n; j++) p2[j] = lgetg(n+1,t_COL);
1.1 noro 3159: for (j=1; j<=n; j++)
3160: for (i=j; i<=n; i++)
3161: {
1.2 ! noro 3162: p1 = RXQX_mul((GEN)W[i],(GEN)W[j], nfT);
! 3163: p1 = RXQX_rem(p1, pol, nfT);
! 3164: p1 = simplify_i(quicktrace(p1,sym));
! 3165: coeff(p2,j,i)=coeff(p2,i,j) = (long)p1;
1.1 noro 3166: }
1.2 ! noro 3167: d = algtobasis_i(nf, det(p2));
1.1 noro 3168:
1.2 ! noro 3169: i=1; while (i<=n && gegal((GEN)I[i], id)) i++;
! 3170: if (i > n) D = id;
1.1 noro 3171: else
3172: {
3173: D = (GEN)I[i];
3174: for (i++; i<=n; i++)
1.2 ! noro 3175: if (!gegal((GEN)I[i], id)) D = idealmul(nf,D,(GEN)I[i]);
1.1 noro 3176: D = idealpow(nf,D,gdeux);
3177: }
1.2 ! noro 3178: p2 = core2partial(content(d));
! 3179: p2 = sqri((GEN)p2[2]);
! 3180:
1.1 noro 3181: i = all? 2: 0;
1.2 ! noro 3182: p1=cgetg(3 + i, t_VEC);
! 3183: if (i) { p1[1] = lcopy(A); p1[2] = lcopy(I); }
1.1 noro 3184: p1[1+i] = (long)idealmul(nf,D,d);
1.2 ! noro 3185: p1[2+i] = (long)gdiv(d, p2);
! 3186: return gerepileupto(av,p1);
1.1 noro 3187: }
3188:
3189: GEN
3190: rnfpseudobasis(GEN nf, GEN pol)
3191: {
3192: return rnfround2all(nf,pol,1);
3193: }
3194:
3195: GEN
3196: rnfdiscf(GEN nf, GEN pol)
3197: {
3198: return rnfround2all(nf,pol,0);
3199: }
3200:
1.2 ! noro 3201: GEN
! 3202: gen_if_principal(GEN bnf, GEN x)
! 3203: {
! 3204: gpmem_t av = avma;
! 3205: GEN z = isprincipalall(bnf,x, nf_GEN_IF_PRINCIPAL | nf_FORCE);
! 3206: if (typ(z) == t_INT) { avma = av; return NULL; }
! 3207: return z;
! 3208: }
! 3209:
! 3210: /* given bnf and a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d] it
! 3211: * does not matter), tries to simplify the HNF as much as possible. The
! 3212: * resulting matrix will be upper triangular but the diagonal coefficients
! 3213: * will not be equal to 1. The ideals are guaranteed to be integral and
! 3214: * primitive. */
1.1 noro 3215: GEN
3216: rnfsimplifybasis(GEN bnf, GEN order)
3217: {
1.2 ! noro 3218: gpmem_t av = avma;
! 3219: long j, n, l;
1.1 noro 3220: GEN p1,id,Az,Iz,nf,A,I;
3221:
1.2 ! noro 3222: bnf = checkbnf(bnf); nf = (GEN)bnf[7];
1.1 noro 3223: if (typ(order)!=t_VEC || lg(order)<3)
3224: err(talker,"not a pseudo-basis in nfsimplifybasis");
1.2 ! noro 3225: A = (GEN)order[1];
! 3226: I = (GEN)order[2]; n = lg(A)-1;
! 3227: id = idmat(degpol(nf[1]));
! 3228: Iz = cgetg(n+1,t_VEC);
! 3229: Az = cgetg(n+1,t_MAT);
1.1 noro 3230: for (j=1; j<=n; j++)
3231: {
1.2 ! noro 3232: if (gegal((GEN)I[j],id)) { Iz[j]=(long)id; Az[j]=A[j]; continue; }
! 3233:
! 3234: Iz[j] = (long)Q_primitive_part((GEN)I[j], &p1);
! 3235: Az[j] = p1? lmul((GEN)A[j],p1): A[j];
! 3236: if (p1 && gegal((GEN)Iz[j], id)) continue;
! 3237:
! 3238: p1 = gen_if_principal(bnf, (GEN)Iz[j]);
! 3239: if (p1)
1.1 noro 3240: {
1.2 ! noro 3241: Iz[j] = (long)id;
! 3242: Az[j] = (long)element_mulvec(nf,p1,(GEN)Az[j]);
1.1 noro 3243: }
3244: }
1.2 ! noro 3245: l = lg(order); p1 = cgetg(l, t_VEC);
! 3246: p1[1] = (long)Az;
! 3247: p1[2] = (long)Iz; for (j=3; j<l; j++) p1[j] = order[j];
! 3248: return gerepilecopy(av, p1);
1.1 noro 3249: }
3250:
3251: GEN
3252: rnfdet2(GEN nf, GEN A, GEN I)
3253: {
1.2 ! noro 3254: gpmem_t av,tetpil;
! 3255: long i;
1.1 noro 3256: GEN p1;
3257:
3258: nf=checknf(nf); av = tetpil = avma;
3259: p1=idealhermite(nf,det(matbasistoalg(nf,A)));
3260: for(i=1;i<lg(I);i++) { tetpil=avma; p1=idealmul(nf,p1,(GEN)I[i]); }
3261: tetpil=avma; return gerepile(av,tetpil,p1);
3262: }
3263:
3264: GEN
3265: rnfdet(GEN nf, GEN order)
3266: {
3267: if (typ(order)!=t_VEC || lg(order)<3)
3268: err(talker,"not a pseudo-matrix in rnfdet");
3269: return rnfdet2(nf,(GEN)order[1],(GEN)order[2]);
3270: }
3271:
3272: GEN
3273: rnfdet0(GEN nf, GEN x, GEN y)
3274: {
3275: return y? rnfdet2(nf,x,y): rnfdet(nf,x);
3276: }
3277:
1.2 ! noro 3278: /* Given two fractional ideals a and b, gives x in a, y in b, z in b^-1,
! 3279: t in a^-1 such that xt-yz=1. In the present version, z is in Z. */
! 3280: static GEN
! 3281: nfidealdet1(GEN nf, GEN a, GEN b)
! 3282: {
! 3283: gpmem_t av = avma;
! 3284: GEN x,p1,res,u,v,da,db;
! 3285:
! 3286: a = idealinv(nf,a);
! 3287: da = denom(a); if (!gcmp1(da)) a = gmul(da,a);
! 3288: db = denom(b); if (!gcmp1(db)) b = gmul(db,b);
! 3289: x = idealcoprime(nf,a,b);
! 3290: p1 = idealaddtoone(nf, idealmul(nf,x,a), b);
! 3291: u = (GEN)p1[1];
! 3292: v = (GEN)p1[2];
! 3293:
! 3294: res = cgetg(5,t_VEC);
! 3295: res[1] = (long)gmul(x,da);
! 3296: res[2] = (long)gdiv(v,db);
! 3297: res[3] = lnegi(db);
! 3298: res[4] = (long)element_div(nf, u, (GEN)res[1]);
! 3299: return gerepilecopy(av,res);
! 3300: }
! 3301:
! 3302: static GEN
! 3303: get_order(GEN nf, GEN O, char *s)
! 3304: {
! 3305: if (typ(O) == t_POL)
! 3306: return rnfpseudobasis(nf, O);
! 3307: if (typ(O)!=t_VEC || lg(O) < 3)
! 3308: err(talker,"not a pseudo-matrix in %s", s);
! 3309: return O;
! 3310: }
! 3311:
! 3312: /* given a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d]), gives an
! 3313: * n x n matrix (not in HNF) of a pseudo-basis and an ideal vector
! 3314: * [id,id,...,id,I] such that order = nf[7]^(n-1) x I.
! 3315: * Uses the approximation theorem ==> slow. */
1.1 noro 3316: GEN
3317: rnfsteinitz(GEN nf, GEN order)
3318: {
1.2 ! noro 3319: gpmem_t av = avma;
! 3320: long i,n,l;
1.1 noro 3321: GEN Id,A,I,p1,a,b;
3322:
3323: nf = checknf(nf);
3324: Id = idmat(degpol(nf[1]));
1.2 ! noro 3325: order = get_order(nf, order, "rnfsteinitz");
! 3326: A = matalgtobasis(nf, (GEN)order[1]);
! 3327: I = dummycopy((GEN)order[2]); n=lg(A)-1;
1.1 noro 3328: if (typ(A) != t_MAT || typ(I) != t_VEC || lg(I) != n+1)
3329: err(typeer,"rnfsteinitz");
3330: for (i=1; i<n; i++)
3331: {
1.2 ! noro 3332: GEN c1,c2;
! 3333: a = (GEN)I[i]; if (gegal(a,Id)) continue;
! 3334:
! 3335: c1 = (GEN)A[i];
! 3336: c2 = (GEN)A[i+1];
! 3337: b = (GEN)I[i+1];
! 3338: if (gegal(b,Id))
! 3339: {
! 3340: A[i] = (long)c2;
! 3341: A[i+1]= lneg(c1);
! 3342: I[i] = (long)b;
! 3343: I[i+1]= (long)a;
! 3344: }
! 3345: else
1.1 noro 3346: {
1.2 ! noro 3347: p1 = nfidealdet1(nf,a,b);
! 3348: A[i] = ladd(element_mulvec(nf, (GEN)p1[1], c1),
! 3349: element_mulvec(nf, (GEN)p1[2], c2));
! 3350: A[i+1]= ladd(element_mulvec(nf, (GEN)p1[3], c1),
! 3351: element_mulvec(nf, (GEN)p1[4], c2));
! 3352: I[i] = (long)Id;
! 3353: I[i+1]= (long)Q_primitive_part(idealmul(nf,a,b), &p1);
! 3354: if (p1) A[i+1] = (long)element_mulvec(nf, p1,(GEN)A[i+1]);
! 3355: }
! 3356: }
! 3357: l = lg(order);
! 3358: p1 = cgetg(l,t_VEC);
! 3359: p1[1]=(long)A;
! 3360: p1[2]=(long)I; for (i=3; i<l; i++) p1[i]=order[i];
! 3361: return gerepilecopy(av, p1);
1.1 noro 3362: }
3363:
1.2 ! noro 3364: /* Given bnf and either an order as output by rnfpseudobasis or a polynomial,
! 3365: * and outputs a basis if it is free, an n+1-generating set if it is not */
1.1 noro 3366: GEN
3367: rnfbasis(GEN bnf, GEN order)
3368: {
1.2 ! noro 3369: gpmem_t av = avma;
! 3370: long j, n;
! 3371: GEN nf, A, I, cl, col, a, id;
! 3372:
! 3373: bnf = checkbnf(bnf); nf = (GEN)bnf[7];
! 3374: id = idmat(degpol(nf[1]));
! 3375: order = get_order(nf, order, "rnfbasis");
! 3376: I = (GEN)order[2]; n = lg(I)-1;
1.1 noro 3377: j=1; while (j<n && gegal((GEN)I[j],id)) j++;
1.2 ! noro 3378: if (j<n)
1.1 noro 3379: {
1.2 ! noro 3380: order = rnfsteinitz(nf,order);
! 3381: I = (GEN)order[2];
1.1 noro 3382: }
1.2 ! noro 3383: A = (GEN)order[1];
! 3384: col= (GEN)A[n]; A = vecextract_i(A, 1, n-1);
! 3385: cl = (GEN)I[n];
! 3386: a = gen_if_principal(bnf, cl);
! 3387: if (!a)
1.1 noro 3388: {
1.2 ! noro 3389: GEN p1 = ideal_two_elt(nf, cl);
! 3390: A = concatsp(A, gmul((GEN)p1[1], col));
! 3391: a = (GEN)p1[2];
1.1 noro 3392: }
1.2 ! noro 3393: A = concatsp(A, element_mulvec(nf, a, col));
! 3394: return gerepilecopy(av, A);
1.1 noro 3395: }
3396:
1.2 ! noro 3397: /* Given bnf and either an order as output by rnfpseudobasis or a polynomial,
! 3398: * and outputs a basis (not pseudo) in Hermite Normal Form if it exists, zero
! 3399: * if not
1.1 noro 3400: */
3401: GEN
3402: rnfhermitebasis(GEN bnf, GEN order)
3403: {
1.2 ! noro 3404: gpmem_t av = avma;
! 3405: long j, n;
! 3406: GEN nf, A, I, a, id;
! 3407:
! 3408: bnf = checkbnf(bnf); nf = (GEN)bnf[7];
! 3409: id = idmat(degpol(nf[1]));
! 3410: order = get_order(nf, order, "rnfbasis");
! 3411: A = (GEN)order[1]; A = dummycopy(A);
! 3412: I = (GEN)order[2]; n = lg(A)-1;
1.1 noro 3413: for (j=1; j<=n; j++)
3414: {
1.2 ! noro 3415: if (gegal((GEN)I[j],id)) continue;
! 3416:
! 3417: a = gen_if_principal(bnf, (GEN)I[j]);
! 3418: if (!a) { avma = av; return gzero; }
! 3419: A[j] = (long)element_mulvec(nf, a, (GEN)A[j]);
1.1 noro 3420: }
3421: return gerepilecopy(av,A);
3422: }
3423:
1.2 ! noro 3424: static long
! 3425: _rnfisfree(GEN bnf, GEN order)
1.1 noro 3426: {
1.2 ! noro 3427: long n, j;
! 3428: GEN nf, p1, id, I;
1.1 noro 3429:
3430: bnf = checkbnf(bnf);
3431: if (gcmp1(gmael3(bnf,8,1,1))) return 1;
3432:
1.2 ! noro 3433: nf = (GEN)bnf[7]; id = idmat(degpol(nf[1]));
! 3434: order = get_order(nf, order, "rnfisfree");
! 3435: I = (GEN)order[2]; n = lg(I)-1;
1.1 noro 3436: j=1; while (j<=n && gegal((GEN)I[j],id)) j++;
1.2 ! noro 3437: if (j>n) return 1;
1.1 noro 3438:
1.2 ! noro 3439: p1 = (GEN)I[j];
1.1 noro 3440: for (j++; j<=n; j++)
1.2 ! noro 3441: if (!gegal((GEN)I[j],id)) p1 = idealmul(nf,p1,(GEN)I[j]);
! 3442: return gcmp0( isprincipal(bnf,p1) );
! 3443: }
! 3444:
! 3445: long
! 3446: rnfisfree(GEN bnf, GEN order)
! 3447: {
! 3448: gpmem_t av = avma;
! 3449: long n = _rnfisfree(bnf, order);
! 3450: avma = av; return n;
1.1 noro 3451: }
3452:
3453: /**********************************************************************/
3454: /** **/
3455: /** COMPOSITUM OF TWO NUMBER FIELDS **/
3456: /** **/
3457: /**********************************************************************/
1.2 ! noro 3458: /* modular version */
1.1 noro 3459: GEN
3460: polcompositum0(GEN A, GEN B, long flall)
3461: {
1.2 ! noro 3462: gpmem_t av = avma;
! 3463: long v, k;
1.1 noro 3464: GEN C, LPRS;
3465:
3466: if (typ(A)!=t_POL || typ(B)!=t_POL) err(typeer,"polcompositum0");
3467: if (degpol(A)<=0 || degpol(B)<=0) err(constpoler,"compositum");
3468: v = varn(A);
3469: if (varn(B) != v) err(talker,"not the same variable in compositum");
3470: C = content(A); if (!gcmp1(C)) A = gdiv(A, C);
3471: C = content(B); if (!gcmp1(C)) B = gdiv(B, C);
3472: check_pol_int(A,"compositum");
3473: check_pol_int(B,"compositum");
3474: if (!ZX_is_squarefree(A)) err(talker,"compositum: %Z not separable", A);
3475: if (!ZX_is_squarefree(B)) err(talker,"compositum: %Z not separable", B);
3476:
3477: k = 1; C = ZY_ZXY_resultant_all(A, B, &k, flall? &LPRS: NULL);
1.2 ! noro 3478: C = DDF2(C, 0); /* C = Res_Y (A, B(X + kY)) guaranteed squarefree */
! 3479: C = sort_vecpol(C);
1.1 noro 3480: if (flall)
3481: {
3482: long i,l = lg(C);
3483: GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b - k a */
3484: for (i=1; i<l; i++)
3485: { /* invmod possibly very costly */
1.2 ! noro 3486: a = gmul((GEN)LPRS[1], QX_invmod((GEN)LPRS[2], (GEN)C[i]));
1.1 noro 3487: a = gneg_i(gmod(a, (GEN)C[i]));
3488: b = gadd(polx[v], gmulsg(k,a));
3489: w = cgetg(5,t_VEC); /* [C, a, b, n ] */
1.2 ! noro 3490: w[1] = C[i];
1.1 noro 3491: w[2] = (long)to_polmod(a, (GEN)w[1]);
3492: w[3] = (long)to_polmod(b, (GEN)w[1]);
3493: w[4] = lstoi(-k); C[i] = (long)w;
3494: }
3495: }
3496: settyp(C, t_VEC); return gerepilecopy(av, C);
3497: }
3498:
3499: GEN
3500: compositum(GEN pol1,GEN pol2)
3501: {
3502: return polcompositum0(pol1,pol2,0);
3503: }
3504:
3505: GEN
3506: compositum2(GEN pol1,GEN pol2)
3507: {
3508: return polcompositum0(pol1,pol2,1);
3509: }
3510:
3511: int
3512: nfissquarefree(GEN nf, GEN x)
3513: {
1.2 ! noro 3514: gpmem_t av = avma;
1.1 noro 3515: GEN g, y = derivpol(x);
3516: if (isrational(x))
3517: g = modulargcd(x, y);
3518: else
3519: g = nfgcd(x, y, nf, NULL);
3520: avma = av; return (degpol(g) == 0);
3521: }
3522:
3523: GEN
1.2 ! noro 3524: _rnfequation(GEN A, GEN B, long *pk, GEN *pLPRS)
1.1 noro 3525: {
1.2 ! noro 3526: long k, lA, lB;
! 3527: GEN nf, C;
! 3528:
! 3529: A = get_nfpol(A, &nf);
! 3530: B = fix_relative_pol(A,B,1);
! 3531: lA = lgef(A);
! 3532: lB = lgef(B);
1.1 noro 3533: if (lA<=3 || lB<=3) err(constpoler,"rnfequation");
3534:
3535: check_pol_int(A,"rnfequation");
1.2 ! noro 3536: B = primpart(lift_intern(B));
1.1 noro 3537: for (k=2; k<lB; k++)
3538: if (lgef(B[k]) >= lA) B[k] = lres((GEN)B[k],A);
3539:
3540: if (!nfissquarefree(A,B))
3541: err(talker,"not k separable relative equation in rnfequation");
3542:
1.2 ! noro 3543: *pk = 0; C = ZY_ZXY_resultant_all(A, B, pk, pLPRS);
1.1 noro 3544: if (gsigne(leadingcoeff(C)) < 0) C = gneg_i(C);
1.2 ! noro 3545: *pk = -*pk; return primpart(C);
! 3546: }
! 3547:
! 3548: GEN
! 3549: rnfequation0(GEN A, GEN B, long flall)
! 3550: {
! 3551: gpmem_t av = avma;
! 3552: long k;
! 3553: GEN LPRS, nf, C;
! 3554:
! 3555: A = get_nfpol(A, &nf);
! 3556: C = _rnfequation(A, B, &k, flall? &LPRS: NULL);
1.1 noro 3557: if (flall)
3558: {
1.2 ! noro 3559: GEN w,a,b; /* a,b,c root of A,B,C = compositum, c = b + k a */
1.1 noro 3560: /* invmod possibly very costly */
1.2 ! noro 3561: a = gmul((GEN)LPRS[1], QX_invmod((GEN)LPRS[2], C));
1.1 noro 3562: a = gneg_i(gmod(a, C));
1.2 ! noro 3563: b = gadd(polx[varn(A)], gmulsg(k,a));
1.1 noro 3564: w = cgetg(4,t_VEC); /* [C, a, n ] */
1.2 ! noro 3565: w[1] = (long)C;
1.1 noro 3566: w[2] = (long)to_polmod(a, (GEN)w[1]);
1.2 ! noro 3567: w[3] = lstoi(k); C = w;
1.1 noro 3568: }
3569: return gerepilecopy(av, C);
3570: }
3571:
3572: GEN
3573: rnfequation(GEN nf,GEN pol2)
3574: {
3575: return rnfequation0(nf,pol2,0);
3576: }
3577:
3578: GEN
3579: rnfequation2(GEN nf,GEN pol2)
3580: {
3581: return rnfequation0(nf,pol2,1);
3582: }
3583:
3584: static GEN
3585: nftau(long r1, GEN x)
3586: {
3587: long i, ru = lg(x);
3588: GEN s;
3589:
3590: s = r1 ? (GEN)x[1] : gmul2n(greal((GEN)x[1]),1);
3591: for (i=2; i<=r1; i++) s=gadd(s,(GEN)x[i]);
3592: for ( ; i<ru; i++) s=gadd(s,gmul2n(greal((GEN)x[i]),1));
3593: return s;
3594: }
3595:
3596: static GEN
3597: nftocomplex(GEN nf, GEN x)
3598: {
3599: long ru,vnf,k;
3600: GEN p2,p3,ronf;
3601:
3602: p2 = (typ(x)==t_POLMOD)? (GEN)x[2]: gmul((GEN)nf[7],x);
3603: vnf=varn(nf[1]);
3604: ronf=(GEN)nf[6]; ru=lg(ronf); p3=cgetg(ru,t_COL);
3605: for (k=1; k<ru; k++) p3[k]=lsubst(p2,vnf,(GEN)ronf[k]);
3606: return p3;
3607: }
3608:
3609: static GEN
3610: rnfscal(GEN mth, GEN xth, GEN yth)
3611: {
3612: long n,ru,i,j,kk;
3613: GEN x,y,m,res,p1,p2;
3614:
3615: n=lg(mth)-1; ru=lg(gcoeff(mth,1,1));
3616: res=cgetg(ru,t_COL);
3617: for (kk=1; kk<ru; kk++)
3618: {
3619: m=cgetg(n+1,t_MAT);
3620: for (j=1; j<=n; j++)
3621: {
3622: p1=cgetg(n+1,t_COL); m[j]=(long)p1;
3623: for (i=1; i<=n; i++) { p2=gcoeff(mth,i,j); p1[i]=p2[kk]; }
3624: }
3625: x=cgetg(n+1,t_VEC);
3626: for (j=1; j<=n; j++) x[j]=(long)gconj((GEN)((GEN)xth[j])[kk]);
3627: y=cgetg(n+1,t_COL);
3628: for (j=1; j<=n; j++) y[j]=((GEN)yth[j])[kk];
3629: res[kk]=(long)gmul(x,gmul(m,y));
3630: }
3631: return res;
3632: }
3633:
3634: static GEN
3635: rnfdiv(GEN x, GEN y)
3636: {
3637: long i, ru = lg(x);
1.2 ! noro 3638: GEN z = cgetg(ru,t_COL);
1.1 noro 3639:
3640: for (i=1; i<ru; i++) z[i]=(long)gdiv((GEN)x[i],(GEN)y[i]);
3641: return z;
3642: }
3643:
3644: static GEN
3645: rnfmul(GEN x, GEN y)
3646: {
3647: long i, ru = lg(x);
1.2 ! noro 3648: GEN z = cgetg(ru,t_COL);
1.1 noro 3649:
3650: for (i=1; i<ru; i++) z[i]=(long)gmul((GEN)x[i],(GEN)y[i]);
3651: return z;
3652: }
3653:
3654: static GEN
3655: rnfvecmul(GEN x, GEN v)
3656: {
3657: long i, lx = lg(v);
1.2 ! noro 3658: GEN y = cgetg(lx,typ(v));
1.1 noro 3659:
3660: for (i=1; i<lx; i++) y[i]=(long)rnfmul(x,(GEN)v[i]);
3661: return y;
3662: }
3663:
3664: static GEN
1.2 ! noro 3665: allonge(GEN v, long l)
1.1 noro 3666: {
1.2 ! noro 3667: long r, r2, i;
! 3668: GEN y = cgetg(l,t_COL);
1.1 noro 3669:
1.2 ! noro 3670: r = lg(v); r2 = l-r;
! 3671: for (i=1; i < r; i++) y[i] = v[i];
! 3672: for ( ; i < l; i++) y[i] = lconj((GEN)v[i-r2]);
1.1 noro 3673: return y;
3674: }
3675:
3676: static GEN
3677: findmin(GEN nf, GEN ideal, GEN muf,long prec)
3678: {
1.2 ! noro 3679: gpmem_t av = avma;
! 3680: long i, l;
! 3681: GEN m,y, G = gmael(nf,5,2);
1.1 noro 3682:
1.2 ! noro 3683: m = gmul(G, ideal);
! 3684: m = lllintern(m,4,1,prec);
1.1 noro 3685: if (!m)
3686: {
1.2 ! noro 3687: ideal = lllint_ip(ideal,4);
! 3688: m = gmul(G, ideal);
! 3689: m = lllintern(m,4,1,prec);
1.1 noro 3690: if (!m) err(precer,"rnflllgram");
3691: }
1.2 ! noro 3692: ideal = gmul(ideal,m);
! 3693: l = lg(ideal); y = cgetg(l,t_MAT);
! 3694: for (i=1; i<l; i++)
! 3695: y[i] = (long)allonge(nftocomplex(nf,(GEN)ideal[i]),l);
! 3696: m = ground(greal(gauss(y, allonge(muf,l))));
! 3697: return gerepileupto(av, gmul(ideal,m));
1.1 noro 3698: }
3699:
3700: #define swap(x,y) { long _t=x; x=y; y=_t; }
3701:
3702: /* given a base field nf (e.g main variable y), a polynomial pol with
3703: * coefficients in nf (e.g main variable x), and an order as output
3704: * by rnfpseudobasis, outputs a reduced order.
3705: */
3706: GEN
3707: rnflllgram(GEN nf, GEN pol, GEN order,long prec)
3708: {
1.2 ! noro 3709: gpmem_t av=avma,tetpil;
! 3710: long i,j,k,l,kk,kmax,r1,ru,lx,vnf;
1.1 noro 3711: GEN p1,p2,M,I,U,ronf,poll,unro,roorder,powreorder,mth,s,MC,MPOL,MCS;
3712: GEN B,mu,Bf,temp,ideal,x,xc,xpol,muf,mufc,muno,y,z,Ikk_inv;
3713:
3714: /* Initializations and verifications */
3715:
3716: nf=checknf(nf);
3717: if (typ(order)!=t_VEC || lg(order)<3)
3718: err(talker,"not a pseudo-matrix in rnflllgram");
3719: M=(GEN)order[1]; I=(GEN)order[2]; lx=lg(I);
3720: if (lx < 3) return gcopy(order);
3721: if (lx-1 != degpol(pol)) err(consister,"rnflllgram");
3722: U=idmat(lx-1); I = dummycopy(I);
3723:
3724: /* Compute the relative T2 matrix of powers of theta */
3725:
3726: vnf=varn(nf[1]); ronf=(GEN)nf[6]; ru=lg(ronf); poll=lift(pol);
3727: r1 = nf_get_r1(nf);
3728: unro=cgetg(lx,t_COL); for (i=1; i<lx; i++) unro[i]=un;
3729: roorder=cgetg(ru,t_VEC);
3730: for (i=1; i<ru; i++)
3731: roorder[i]=lroots(gsubst(poll,vnf,(GEN)ronf[i]),prec);
3732: powreorder=cgetg(lx,t_MAT);
3733: p1=cgetg(ru,t_COL); powreorder[1]=(long)p1;
3734: for (i=1; i<ru; i++) p1[i]=(long)unro;
3735: for (k=2; k<lx; k++)
3736: {
3737: p1=cgetg(ru,t_COL); powreorder[k]=(long)p1;
3738: for (i=1; i<ru; i++)
3739: {
3740: p2=cgetg(lx,t_COL); p1[i]=(long)p2;
3741: for (j=1; j<lx; j++)
3742: p2[j] = lmul(gmael(roorder,i,j),gmael3(powreorder,k-1,i,j));
3743: }
3744: }
3745: mth=cgetg(lx,t_MAT);
3746: for (l=1; l<lx; l++)
3747: {
3748: p1=cgetg(lx,t_COL); mth[l]=(long)p1;
3749: for (k=1; k<lx; k++)
3750: {
3751: p2=cgetg(ru,t_COL); p1[k]=(long)p2;
3752: for (i=1; i<ru; i++)
3753: {
3754: s=gzero;
3755: for (j=1; j<lx; j++)
3756: s = gadd(s,gmul(gconj(gmael3(powreorder,k,i,j)),
3757: gmael3(powreorder,l,i,j)));
3758: p2[i]=(long)s;
3759: }
3760: }
3761: }
3762:
3763: /* Transform the matrix M into a matrix with coefficients in K and also
3764: with coefficients polymod */
3765:
3766: MC=cgetg(lx,t_MAT); MPOL=cgetg(lx,t_MAT);
3767: for (j=1; j<lx; j++)
3768: {
3769: p1=cgetg(lx,t_COL); MC[j]=(long)p1;
3770: p2=cgetg(lx,t_COL); MPOL[j]=(long)p2;
3771: for (i=1; i<lx; i++)
3772: {
3773: p2[i]=(long)basistoalg(nf,gcoeff(M,i,j));
3774: p1[i]=(long)nftocomplex(nf,(GEN)p2[i]);
3775: }
3776: }
3777: MCS=cgetg(lx,t_MAT);
3778:
3779: /* Start LLL algorithm */
3780:
3781: mu=cgetg(lx,t_MAT); B=cgetg(lx,t_COL);
3782: for (j=1; j<lx; j++)
3783: {
3784: p1=cgetg(lx,t_COL); mu[j]=(long)p1; for (i=1; i<lx; i++) p1[i]=zero;
3785: B[j]=zero;
3786: }
3787: kk=2; if (DEBUGLEVEL) fprintferr("kk = %ld ",kk);
3788: kmax=1; B[1]=lreal(rnfscal(mth,(GEN)MC[1],(GEN)MC[1]));
3789: MCS[1]=lcopy((GEN)MC[1]);
3790: do
3791: {
3792: if (kk>kmax)
3793: {
3794: /* Incremental Gram-Schmidt */
3795: kmax=kk; MCS[kk]=lcopy((GEN)MC[kk]);
3796: for (j=1; j<kk; j++)
3797: {
3798: coeff(mu,kk,j) = (long) rnfdiv(rnfscal(mth,(GEN)MCS[j],(GEN)MC[kk]),
3799: (GEN) B[j]);
3800: MCS[kk] = lsub((GEN) MCS[kk], rnfvecmul(gcoeff(mu,kk,j),(GEN)MCS[j]));
3801: }
3802: B[kk] = lreal(rnfscal(mth,(GEN)MCS[kk],(GEN)MCS[kk]));
3803: if (gcmp0((GEN)B[kk])) err(lllger3);
3804: }
3805:
3806: /* RED(k,k-1) */
3807: l=kk-1; Ikk_inv=idealinv(nf, (GEN)I[kk]);
3808: ideal=idealmul(nf,(GEN)I[l],Ikk_inv);
3809: x=findmin(nf,ideal,gcoeff(mu,kk,l),2*prec-2);
3810: if (!gcmp0(x))
3811: {
3812: xpol=basistoalg(nf,x); xc=nftocomplex(nf,xpol);
3813: MC[kk]=lsub((GEN)MC[kk],rnfvecmul(xc,(GEN)MC[l]));
3814: U[kk]=lsub((GEN)U[kk],gmul(xpol,(GEN)U[l]));
3815: coeff(mu,kk,l)=lsub(gcoeff(mu,kk,l),xc);
3816: for (i=1; i<l; i++)
3817: coeff(mu,kk,i)=lsub(gcoeff(mu,kk,i),rnfmul(xc,gcoeff(mu,l,i)));
3818: }
3819: /* Test LLL condition */
3820: p1=nftau(r1,gadd((GEN) B[kk],
3821: gmul(gnorml2(gcoeff(mu,kk,kk-1)),(GEN)B[kk-1])));
3822: p2=gdivgs(gmulsg(9,nftau(r1,(GEN)B[kk-1])),10);
3823: if (gcmp(p1,p2)<=0)
3824: {
3825: /* Execute SWAP(k) */
3826: k=kk;
3827: swap(MC[k-1],MC[k]);
3828: swap(U[k-1],U[k]);
3829: swap(I[k-1],I[k]);
3830: for (j=1; j<=k-2; j++) swap(coeff(mu,k-1,j),coeff(mu,k,j));
3831: muf=gcoeff(mu,k,k-1);
3832: mufc=gconj(muf); muno=greal(rnfmul(muf,mufc));
3833: Bf=gadd((GEN)B[k],rnfmul(muno,(GEN)B[k-1]));
3834: p1=rnfdiv((GEN)B[k-1],Bf);
3835: coeff(mu,k,k-1)=(long)rnfmul(mufc,p1);
3836: temp=(GEN)MCS[k-1];
3837: MCS[k-1]=ladd((GEN)MCS[k],rnfvecmul(muf,(GEN)MCS[k-1]));
3838: MCS[k]=lsub(rnfvecmul(rnfdiv((GEN)B[k],Bf),temp),
3839: rnfvecmul(gcoeff(mu,k,k-1),(GEN)MCS[k]));
3840: B[k]=(long)rnfmul((GEN)B[k],p1); B[k-1]=(long)Bf;
3841: for (i=k+1; i<=kmax; i++)
3842: {
3843: temp=gcoeff(mu,i,k);
3844: coeff(mu,i,k)=lsub(gcoeff(mu,i,k-1),rnfmul(muf,gcoeff(mu,i,k)));
3845: coeff(mu,i,k-1) = ladd(temp, rnfmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k)));
3846: }
3847: if (kk>2) { kk--; if (DEBUGLEVEL) fprintferr("%ld ",kk); }
3848: }
3849: else
3850: {
3851: for (l=kk-2; l; l--)
3852: {
3853: /* RED(k,l) */
3854: ideal=idealmul(nf,(GEN)I[l],Ikk_inv);
3855: x=findmin(nf,ideal,gcoeff(mu,kk,l),2*prec-2);
3856: if (!gcmp0(x))
3857: {
3858: xpol=basistoalg(nf,x); xc=nftocomplex(nf,xpol);
3859: MC[kk]=(long)gsub((GEN)MC[kk],rnfvecmul(xc,(GEN)MC[l]));
3860: U[kk]=(long)gsub((GEN)U[kk],gmul(xpol,(GEN)U[l]));
3861: coeff(mu,kk,l)=lsub(gcoeff(mu,kk,l),xc);
3862: for (i=1; i<l; i++)
3863: coeff(mu,kk,i) = lsub(gcoeff(mu,kk,i), rnfmul(xc,gcoeff(mu,l,i)));
3864: }
3865: }
3866: kk++; if (DEBUGLEVEL) fprintferr("%ld ",kk);
3867: }
3868: }
3869: while (kk<lx);
3870: if (DEBUGLEVEL) fprintferr("\n");
3871: p1=gmul(MPOL,U); tetpil=avma;
3872: y=cgetg(3,t_VEC); z=cgetg(3,t_VEC); y[1]=(long)z;
3873: z[2]=lcopy(I); z[1]=(long)algtobasis(nf,p1);
3874: y[2]=(long)algtobasis(nf,U);
3875: return gerepile(av,tetpil,y);
3876: }
3877:
3878: GEN
3879: rnfpolred(GEN nf, GEN pol, long prec)
3880: {
1.2 ! noro 3881: gpmem_t av = avma;
! 3882: long i, j, n, v = varn(pol);
! 3883: GEN id, al, w, I, O, bnf;
1.1 noro 3884:
3885: if (typ(pol)!=t_POL) err(typeer,"rnfpolred");
1.2 ! noro 3886: bnf = nf; nf = checknf(bnf);
! 3887: bnf = (nf == bnf)? NULL: checkbnf(bnf);
! 3888: if (degpol(pol) <= 1) return _vec(polx[v]);
! 3889:
! 3890: id = rnfpseudobasis(nf,pol);
1.1 noro 3891: if (bnf && gcmp1(gmael3(bnf,8,1,1))) /* if bnf is principal */
3892: {
1.2 ! noro 3893: GEN newI, newO, zk = idmat(degpol(nf[1]));
! 3894: O = (GEN)id[1];
! 3895: I = (GEN)id[2]; n = lg(I)-1;
! 3896: newI = cgetg(n+1,t_VEC);
! 3897: newO = cgetg(n+1,t_MAT);
1.1 noro 3898: for (j=1; j<=n; j++)
3899: {
1.2 ! noro 3900: newI[j] = (long)zk; al = gen_if_principal(bnf,(GEN)I[j]);
! 3901: newO[j] = (long)element_mulvec(nf, al, (GEN)O[j]);
1.1 noro 3902: }
1.2 ! noro 3903: id = cgetg(3,t_VEC);
! 3904: id[1] = (long)newO;
! 3905: id[2] = (long)newI;
1.1 noro 3906: }
1.2 ! noro 3907:
! 3908: id = (GEN)rnflllgram(nf,pol,id,prec)[1];
! 3909: O = (GEN)id[1];
! 3910: I = (GEN)id[2]; n = lg(I)-1;
! 3911: w = cgetg(n+1,t_VEC);
! 3912: pol = lift(pol);
1.1 noro 3913: for (j=1; j<=n; j++)
3914: {
1.2 ! noro 3915: GEN p1, newpol;
! 3916:
! 3917: p1 = (GEN)I[j]; al = gmul(gcoeff(p1,1,1),(GEN)O[j]);
! 3918: p1 = basistoalg(nf,(GEN)al[n]);
1.1 noro 3919: for (i=n-1; i; i--)
1.2 ! noro 3920: p1 = gadd(basistoalg(nf,(GEN)al[i]), gmul(polx[v],p1));
! 3921: p1 = gmodulcp(gtovec(caract2(pol,lift(p1),v)), (GEN)nf[1]);
! 3922: newpol = gtopoly(p1, v);
! 3923:
1.1 noro 3924: p1 = ggcd(newpol, derivpol(newpol));
1.2 ! noro 3925: if (degpol(p1) > 0)
1.1 noro 3926: {
1.2 ! noro 3927: newpol = gdiv(newpol, p1);
! 3928: newpol = gdiv(newpol, leading_term(newpol));
1.1 noro 3929: }
1.2 ! noro 3930: w[j] = (long)newpol;
1.1 noro 3931: }
3932: return gerepilecopy(av,w);
3933: }
3934:
3935: /* given a relative polynomial pol over nf, compute a pseudo-basis for the
3936: * extension, then an absolute basis */
1.2 ! noro 3937: static GEN
! 3938: makebasis(GEN nf, GEN pol, GEN rnfeq)
1.1 noro 3939: {
1.2 ! noro 3940: GEN elts,ids,polabs,plg,plg0,B,bs,p1,den,vbs,d,vpro;
! 3941: gpmem_t av = avma;
! 3942: long n,N,m,i,j,k, v = varn(pol);
! 3943:
! 3944: polabs= (GEN)rnfeq[1];
! 3945: plg = (GEN)rnfeq[2]; plg = lift_intern(plg);
! 3946: p1 = rnfpseudobasis(nf,pol);
1.1 noro 3947: elts= (GEN)p1[1];
3948: ids = (GEN)p1[2];
1.2 ! noro 3949: if (DEBUGLEVEL>1) fprintferr("relative basis computed\n");
! 3950: N = degpol(pol);
! 3951: n = degpol(nf[1]); m = n*N;
! 3952:
! 3953: plg0= Q_remove_denom(plg, &den); /* plg = plg0/den */
! 3954: /* nf = K = Q(a), vbs[i+1] = a^i as an elt of L = Q[X] / polabs */
! 3955: vbs = RXQ_powers(plg0, polabs, n-1);
! 3956: if (den)
! 3957: { /* restore denominators */
! 3958: vbs[2] = (long)plg; d = den;
! 3959: for (i=3; i<=n; i++) { d = mulii(d,den); vbs[i] = ldiv((GEN)vbs[i], d); }
! 3960: }
! 3961:
! 3962: /* bs = integer basis of K, as elements of L */
1.1 noro 3963: bs = gmul(vbs, vecpol_to_mat((GEN)nf[7],n));
3964:
1.2 ! noro 3965: vpro = cgetg(N+1,t_VEC);
! 3966: for (i=1; i<=N; i++) vpro[i] = lpowgs(polx[v], i-1);
! 3967: vpro = gmul(vpro, elts);
! 3968: B = cgetg(m+1, t_MAT);
! 3969: for(i=k=1; i<=N; i++)
! 3970: {
! 3971: GEN w = element_mulvec(nf, (GEN)vpro[i], (GEN)ids[i]);
! 3972: for(j=1; j<=n; j++)
! 3973: {
! 3974: p1 = gres(gmul(bs, (GEN)w[j]), polabs);
! 3975: B[k++] = (long)pol_to_vec(p1, m);
! 3976: }
! 3977: }
! 3978: B = Q_remove_denom(B, &den);
! 3979: if (den) { B = hnfmodid(B, den); B = gdiv(B, den); } else B = idmat(m);
! 3980: p1 = cgetg(3,t_VEC);
! 3981: p1[1] = (long)polabs;
! 3982: p1[2] = (long)B; return gerepilecopy(av, p1);
! 3983: }
! 3984:
! 3985: /* relative polredabs. Returns relative polynomial by default (flag = 0)
! 3986: * flag & nf_ORIG: + element (base change)
! 3987: * flag & nf_ADDZK: + integer basis
! 3988: * flag & nf_ABSOLUTE: absolute polynomial */
! 3989: GEN
! 3990: rnfpolredabs(GEN nf, GEN relpol, long flag)
! 3991: {
! 3992: GEN red, bas, z, elt, POL, pol, T, a;
! 3993: long v, fl;
! 3994: gpmem_t av = avma;
! 3995:
! 3996: if (typ(relpol)!=t_POL) err(typeer,"rnfpolredabs");
! 3997: nf = checknf(nf); v = varn(relpol);
! 3998: if (DEBUGLEVEL>1) (void)timer2();
! 3999: relpol = unifpol(nf,relpol,1);
! 4000: T = (GEN)nf[1];
! 4001: if ((flag & nf_ADDZK) && (flag != (nf_ADDZK|nf_ABSOLUTE)))
! 4002: err(impl,"this combination of flags in rnfpolredabs");
! 4003: fl = (flag & nf_ADDZK)? nf_ORIG: nf_RAW;
! 4004: if (flag & nf_PARTIALFACT)
! 4005: {
! 4006: long sa;
! 4007: bas = NULL; /* -Wall */
! 4008: POL = _rnfequation(nf, relpol, &sa, NULL);
! 4009: red = polredabs0(POL, fl | nf_PARTIALFACT);
! 4010: a = stoi(sa);
! 4011: }
! 4012: else
! 4013: {
! 4014: GEN rel, eq = rnfequation2(nf,relpol);
! 4015: POL = (GEN)eq[1];
! 4016: a = (GEN)eq[3];
! 4017: rel = poleval(relpol, gsub(polx[v],
! 4018: gmul(a, gmodulcp(polx[varn(T)],T))));
! 4019: bas = makebasis(nf, rel, eq);
! 4020: if (DEBUGLEVEL>1)
! 4021: {
! 4022: msgtimer("absolute basis");
! 4023: fprintferr("original absolute generator: %Z\n", POL);
! 4024: }
! 4025: red = polredabs0(bas, fl);
! 4026: }
! 4027: pol = (GEN)red[1];
! 4028: if (DEBUGLEVEL>1) fprintferr("reduced absolute generator: %Z\n",pol);
! 4029: if (flag & nf_ABSOLUTE)
1.1 noro 4030: {
1.2 ! noro 4031: if (flag & nf_ADDZK)
! 4032: {
! 4033: GEN t = (GEN)red[2], B = (GEN)bas[2];
! 4034: GEN v = RXQ_powers(lift_intern(t), pol, degpol(pol)-1);
! 4035: z = cgetg(3, t_VEC);
! 4036: z[1] = (long)pol;
! 4037: z[2] = lmul(v, B);
! 4038: return gerepilecopy(av, z);
! 4039: }
! 4040: return gerepilecopy(av, pol);
! 4041: }
! 4042:
! 4043: elt = eltabstorel((GEN)red[2], T, relpol, a);
! 4044:
! 4045: pol = rnfcharpoly(nf,relpol,elt,v);
! 4046: if (!(flag & nf_ORIG)) return gerepileupto(av, pol);
! 4047: z = cgetg(3,t_VEC);
! 4048: z[1] = (long)pol;
! 4049: z[2] = (long)to_polmod(modreverse_i((GEN)elt[2],(GEN)elt[1]), pol);
! 4050: return gerepilecopy(av, z);
1.1 noro 4051: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to base2.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