Annotation of OpenXM_contrib/pari-2.2/src/basemath/base1.c, Revision 1.2
1.2 ! noro 1: /* $Id: base1.c,v 1.103 2002/09/07 13:34:38 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: /* NUMBER FIELDS */
19: /* */
20: /**************************************************************/
21: #include "pari.h"
22: #include "parinf.h"
1.2 ! noro 23:
! 24: int new_galois_format = 0;
! 25:
! 26: extern GEN R_from_QR(GEN x, long prec);
! 27: extern GEN to_polmod(GEN x, GEN mod);
! 28: extern GEN roots_to_pol_r1r2(GEN a, long r1, long v);
! 29: extern GEN idealhermite_aux(GEN nf, GEN x);
! 30: extern GEN cauchy_bound(GEN p);
! 31: extern GEN galoisbig(GEN x, long prec);
! 32: extern GEN lllfp_marked(int M, GEN x, long D, long flag, long prec, int gram);
! 33: extern GEN lllint_marked(long M, GEN x, long D, int g, GEN *h, GEN *fl, GEN *B);
! 34: extern GEN mulmat_pol(GEN A, GEN x);
! 35: extern ulong smulss(ulong x, ulong y, ulong *rem);
1.1 noro 36:
37: void
38: checkrnf(GEN rnf)
39: {
1.2 ! noro 40: if (typ(rnf)!=t_VEC || lg(rnf)!=12) err(idealer1);
1.1 noro 41: }
42:
43: GEN
1.2 ! noro 44: _checkbnf(GEN bnf)
1.1 noro 45: {
1.2 ! noro 46: if (typ(bnf) == t_VEC)
! 47: switch (lg(bnf))
! 48: {
! 49: case 11: return bnf;
! 50: case 7: return checkbnf((GEN)bnf[1]);
! 51:
! 52: case 3:
! 53: if (typ(bnf[2])==t_POLMOD)
! 54: return checkbnf((GEN)bnf[1]);
! 55: }
! 56: return NULL;
! 57: }
! 58:
! 59: GEN
! 60: _checknf(GEN nf)
! 61: {
! 62: if (typ(nf)==t_VEC)
! 63: switch(lg(nf))
! 64: {
! 65: case 10: return nf;
! 66: case 11: return checknf((GEN)nf[7]);
! 67: case 7: return checknf((GEN)nf[1]);
! 68: case 3: if (typ(nf[2]) == t_POLMOD) return checknf((GEN)nf[1]);
! 69: }
! 70: return NULL;
! 71: }
! 72:
! 73: GEN
! 74: checkbnf(GEN x)
! 75: {
! 76: GEN bnf = _checkbnf(x);
! 77: if (!bnf)
1.1 noro 78: {
1.2 ! noro 79: if (_checknf(x)) err(talker,"please apply bnfinit first");
! 80: err(idealer1);
1.1 noro 81: }
1.2 ! noro 82: return bnf;
1.1 noro 83: }
84:
85: GEN
86: checkbnf_discard(GEN bnf)
87: {
88: GEN x = checkbnf(bnf);
89: if (x != bnf && lg(bnf) == 3)
90: err(warner,"non-monic polynomial. Change of variables discarded");
91: return x;
92: }
93:
94: GEN
1.2 ! noro 95: checknf(GEN x)
1.1 noro 96: {
1.2 ! noro 97: GEN nf = _checknf(x);
! 98: if (!nf)
! 99: {
! 100: if (typ(x)==t_POL) err(talker,"please apply nfinit first");
! 101: err(idealer1);
1.1 noro 102: }
1.2 ! noro 103: return nf;
1.1 noro 104: }
105:
106: void
107: checkbnr(GEN bnr)
108: {
109: if (typ(bnr)!=t_VEC || lg(bnr)!=7)
110: err(talker,"incorrect bigray field");
111: (void)checkbnf((GEN)bnr[1]);
112: }
113:
114: void
115: checkbnrgen(GEN bnr)
116: {
117: checkbnr(bnr);
118: if (lg(bnr[5])<=3)
119: err(talker,"please apply bnrinit(,,1) and not bnrinit(,)");
120: }
121:
122: GEN
123: check_units(GEN bnf, char *f)
124: {
125: GEN nf, x;
126: bnf = checkbnf(bnf); nf = checknf(bnf);
127: x = (GEN)bnf[8];
128: if (lg(x) < 7 || (lg(x[5]) == 1 && lg(nf[6]) > 2))
129: err(talker,"missing units in %s", f);
130: return (GEN)x[5];
131: }
132:
133: void
134: checkid(GEN x, long N)
135: {
136: if (typ(x)!=t_MAT) err(idealer2);
137: if (lg(x) == 1 || lg(x[1]) != N+1)
138: err(talker,"incorrect matrix for ideal");
139: }
140:
141: void
142: checkbid(GEN bid)
143: {
144: if (typ(bid)!=t_VEC || lg(bid)!=6 || lg(bid[1])!=3)
145: err(talker,"incorrect bigideal");
146: }
147:
148: void
149: checkprimeid(GEN id)
150: {
151: if (typ(id) != t_VEC || lg(id) != 6)
152: err(talker,"incorrect prime ideal");
153: }
154:
155: void
156: check_pol_int(GEN x, char *s)
157: {
158: long k = lgef(x)-1;
159: for ( ; k>1; k--)
160: if (typ(x[k])!=t_INT) err(talker,"polynomial not in Z[X] in %s",s);
161: }
162:
163: GEN
164: checknfelt_mod(GEN nf, GEN x, char *s)
165: {
166: if (!gegal((GEN)x[1],(GEN)nf[1]))
167: err(talker,"incompatible modulus in %s:\n mod = %Z,\n nf = %Z",
168: s, x[1], nf[1]);
169: return (GEN)x[2];
170: }
171:
172: GEN
173: get_primeid(GEN x)
174: {
175: if (typ(x) != t_VEC) return NULL;
176: if (lg(x) == 3) x = (GEN)x[1];
177: if (lg(x) != 6 || typ(x[3]) != t_INT) return NULL;
178: return x;
179: }
180:
181: GEN
182: get_bnf(GEN x, int *t)
183: {
184: switch(typ(x))
185: {
186: case t_POL: *t = typ_POL; return NULL;
187: case t_QUAD: *t = typ_Q ; return NULL;
188: case t_VEC:
189: switch(lg(x))
190: {
191: case 3:
192: if (typ(x[2]) != t_POLMOD) break;
193: return get_bnf((GEN)x[1],t);
194: case 6 : *t = typ_QUA; return NULL;
195: case 10: *t = typ_NF; return NULL;
196: case 11: *t = typ_BNF; return x;
197: case 7 : *t = typ_BNR;
198: x = (GEN)x[1]; if (typ(x)!=t_VEC || lg(x)!=11) break;
199: return x;
200: }
201: case t_MAT:
202: if (lg(x)==2)
203: switch(lg(x[1]))
204: {
205: case 8: case 11:
206: *t = typ_CLA; return NULL;
207: }
208: }
209: *t = typ_NULL; return NULL;
210: }
211:
212: GEN
213: get_nf(GEN x, int *t)
214: {
215: switch(typ(x))
216: {
217: case t_POL : *t = typ_POL; return NULL;
218: case t_QUAD: *t = typ_Q ; return NULL;
219: case t_VEC:
220: switch(lg(x))
221: {
222: case 3:
223: if (typ(x[2]) != t_POLMOD) break;
224: return get_nf((GEN)x[1],t);
225: case 10: *t = typ_NF; return x;
226: case 11: *t = typ_BNF;
227: x = (GEN)x[7]; if (typ(x)!=t_VEC || lg(x)!=10) break;
228: return x;
229: case 7 : *t = typ_BNR;
230: x = (GEN)x[1]; if (typ(x)!=t_VEC || lg(x)!=11) break;
231: x = (GEN)x[7]; if (typ(x)!=t_VEC || lg(x)!=10) break;
232: return x;
233: case 9 :
234: x = (GEN)x[2];
235: if (typ(x) == t_VEC && lg(x) == 4) *t = typ_GAL;
236: return NULL;
237: case 14: case 20:
238: *t = typ_ELL; return NULL;
239: }break;
240: case t_MAT:
241: if (lg(x)==2)
242: switch(lg(x[1]))
243: {
244: case 8: case 11:
245: *t = typ_CLA; return NULL;
246: }
247: }
248: *t = typ_NULL; return NULL;
249: }
250:
251: /*************************************************************************/
252: /** **/
253: /** GALOIS GROUP **/
254: /** **/
255: /*************************************************************************/
256:
257: /* exchange elements i and j in vector x */
258: static GEN
259: transroot(GEN x, int i, int j)
260: {
261: long k;
262: x = dummycopy(x);
263: k=x[i]; x[i]=x[j]; x[j]=k; return x;
264: }
265:
266: #define randshift (BITS_IN_RANDOM - 3)
267:
268: GEN
269: tschirnhaus(GEN x)
270: {
1.2 ! noro 271: gpmem_t av = avma, av2;
! 272: long a, v = varn(x);
1.1 noro 273: GEN u, p1 = cgetg(5,t_POL);
274:
275: if (typ(x)!=t_POL) err(notpoler,"tschirnhaus");
276: if (lgef(x) < 4) err(constpoler,"tschirnhaus");
277: if (v) { u=dummycopy(x); setvarn(u,0); x=u; }
278: p1[1] = evalsigne(1)|evalvarn(0)|evallgef(5);
279: do
280: {
281: a = mymyrand() >> randshift; if (a==0) a =1; p1[4]=lstoi(a);
282: a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[3]=lstoi(a);
283: a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[2]=lstoi(a);
1.2 ! noro 284: u = caract2(x,p1,v); av2=avma;
1.1 noro 285: }
286: while (lgef(srgcd(u,derivpol(u))) > 3); /* while u not separable */
287: if (DEBUGLEVEL>1)
288: fprintferr("Tschirnhaus transform. New pol: %Z",u);
1.2 ! noro 289: avma=av2; return gerepileupto(av,u);
1.1 noro 290: }
291: #undef randshift
292:
293: int
294: gpolcomp(GEN p1, GEN p2)
295: {
296: int s,j = lgef(p1)-2;
297:
298: if (lgef(p2)-2 != j)
299: err(bugparier,"gpolcomp (different degrees)");
300: for (; j>=2; j--)
301: {
302: s = absi_cmp((GEN)p1[j], (GEN)p2[j]);
303: if (s) return s;
304: }
305: return 0;
306: }
307:
1.2 ! noro 308: /* assume pol in Z[X]. Find C, L in Z such that POL = C pol(x / L) monic in Z[X].
! 309: * Return POL and set *ptlead = L */
1.1 noro 310: GEN
311: primitive_pol_to_monic(GEN pol, GEN *ptlead)
312: {
313: long i,j,n = degpol(pol);
314: GEN lead,fa,e,a, POL = dummycopy(pol);
315:
316: a = POL + 2; lead = (GEN)a[n];
317: if (signe(lead) < 0) { POL = gneg_i(POL); a = POL+2; lead = negi(lead); }
318: if (is_pm1(lead)) { if (ptlead) *ptlead = NULL; return POL; }
319: fa = auxdecomp(lead,0); lead = gun;
320: e = (GEN)fa[2]; fa = (GEN)fa[1];
321: for (i=lg(e)-1; i>0;i--) e[i] = itos((GEN)e[i]);
322: for (i=lg(fa)-1; i>0; i--)
323: {
324: GEN p = (GEN)fa[i], p1,pk,pku;
325: long k = (long)ceil((double)e[i] / n);
326: long d = k * n - e[i], v, j0;
327: /* find d, k such that p^d pol(x / p^k) monic */
328: for (j=n-1; j>0; j--)
329: {
330: if (!signe(a[j])) continue;
331: v = pvaluation((GEN)a[j], p, &p1);
332: while (v + d < k * j) { k++; d += n; }
333: }
334: pk = gpowgs(p,k); j0 = d/k;
335: pku = gpowgs(p,d - k*j0);
336: /* a[j] *= p^(d - kj) */
337: for (j=j0; j>=0; j--)
338: {
339: if (j < j0) pku = mulii(pku, pk);
340: a[j] = lmulii((GEN)a[j], pku);
341: }
342: j0++;
343: pku = gpowgs(p,k*j0 - d);
344: /* a[j] /= p^(kj - d) */
345: for (j=j0; j<=n; j++)
346: {
347: if (j > j0) pku = mulii(pku, pk);
348: a[j] = ldivii((GEN)a[j], pku);
349: }
350: lead = mulii(lead, pk);
351: }
352: if (ptlead) *ptlead = lead; return POL;
353: }
354:
355: /* compute x1*x2^2 + x2*x3^2 + x3*x4^2 + x4*x1^2 */
356: static GEN
357: get_F4(GEN x)
358: {
359: GEN p1=gzero;
360: long i;
361:
362: for (i=1; i<=4; i++)
363: p1 = gadd(p1, gmul((GEN)x[i], gsqr((GEN)x[(i&3)+1])));
364: return p1;
365: }
366:
1.2 ! noro 367: static GEN
! 368: roots_to_ZX(GEN z, long *e)
! 369: {
! 370: GEN a = roots_to_pol(z,0);
! 371: GEN b = grndtoi(greal(a),e);
! 372: long e1 = gexpo(gimag(a));
! 373: if (e1 > *e) *e = e1;
! 374: return b;
! 375: }
! 376:
! 377: static GEN
! 378: _res(long n, long s, long k)
! 379: {
! 380: GEN y = cgetg(4, t_VEC);
! 381: y[1] = lstoi(n);
! 382: y[2] = lstoi(s);
! 383: if (!new_galois_format) k = (n == 6 && (k == 6 || k == 2))? 2: 1;
! 384: y[3] = lstoi(k); return y;
! 385: }
1.1 noro 386:
387: GEN
388: galois(GEN x, long prec)
389: {
1.2 ! noro 390: gpmem_t av = avma, av1;
! 391: long i,j,k,n,f,l,l2,e,e1,pr,ind;
! 392: GEN x1,p1,p2,p3,p4,p5,w,z,ee;
1.1 noro 393: static int ind5[20]={2,5,3,4, 1,3,4,5, 1,5,2,4, 1,2,3,5, 1,4,2,3};
394: static int ind6[60]={3,5,4,6, 2,6,4,5, 2,3,5,6, 2,4,3,6, 2,5,3,4,
395: 1,4,5,6, 1,5,3,6, 1,6,3,4, 1,3,4,5, 1,6,2,5,
396: 1,2,4,6, 1,5,2,4, 1,3,2,6, 1,2,3,5, 1,4,2,3};
397: if (typ(x)!=t_POL) err(notpoler,"galois");
398: n=degpol(x); if (n<=0) err(constpoler,"galois");
399: if (n>11) err(impl,"galois of degree higher than 11");
1.2 ! noro 400: x = primpart(x);
1.1 noro 401: check_pol_int(x, "galois");
402: if (gisirreducible(x) != gun)
403: err(impl,"galois of reducible polynomial");
404:
405: if (n<4)
406: {
1.2 ! noro 407: if (n == 1) { avma = av; return _res(1, 1,1); }
! 408: if (n == 2) { avma = av; return _res(2,-1,1); }
! 409: /* n = 3 */
! 410: f = carreparfait(ZX_disc(x));
! 411: avma = av;
! 412: return f? _res(3,1,1): _res(6,-1,2);
1.1 noro 413: }
414: x1 = x = primitive_pol_to_monic(x,NULL); av1=avma;
1.2 ! noro 415: if (n > 7) return galoisbig(x,prec);
1.1 noro 416: for(;;)
417: {
1.2 ! noro 418: GEN cb = cauchy_bound(x);
1.1 noro 419: switch(n)
420: {
421: case 4: z = cgetg(7,t_VEC);
1.2 ! noro 422: prec = DEFAULTPREC + (long)((gexpo(cb)*18. / BITS_IN_LONG));
1.1 noro 423: for(;;)
424: {
425: p1=roots(x,prec);
426: z[1] = (long)get_F4(p1);
427: z[2] = (long)get_F4(transroot(p1,1,2));
428: z[3] = (long)get_F4(transroot(p1,1,3));
429: z[4] = (long)get_F4(transroot(p1,1,4));
430: z[5] = (long)get_F4(transroot(p1,2,3));
431: z[6] = (long)get_F4(transroot(p1,3,4));
1.2 ! noro 432: p5 = roots_to_ZX(z, &e);
1.1 noro 433: if (e <= -10) break;
434: prec = (prec<<1)-2;
435: }
1.2 ! noro 436: if (!ZX_is_squarefree(p5)) goto tchi;
1.1 noro 437: p2 = (GEN)factor(p5)[1];
438: switch(lg(p2)-1)
439: {
1.2 ! noro 440: case 1: f = carreparfait(ZX_disc(x)); avma = av;
! 441: return f? _res(12,1,4): _res(24,-1,5);
1.1 noro 442:
1.2 ! noro 443: case 2: avma = av; return _res(8,-1,3);
1.1 noro 444:
1.2 ! noro 445: case 3: avma = av;
! 446: return (degpol(p2[1])==2)? _res(4,1,2): _res(4,-1,1);
1.1 noro 447:
448: default: err(bugparier,"galois (bug1)");
449: }
450:
451: case 5: z = cgetg(7,t_VEC);
1.2 ! noro 452: ee= cgetg(7,t_VECSMALL);
! 453: w = cgetg(7,t_VECSMALL);
! 454: prec = DEFAULTPREC + (long)((gexpo(cb)*21. / BITS_IN_LONG));
1.1 noro 455: for(;;)
456: {
457: for(;;)
458: {
459: p1=roots(x,prec);
460: for (l=1; l<=6; l++)
461: {
462: p2=(l==1)?p1: ((l<6)?transroot(p1,1,l): transroot(p1,2,5));
463: p3=gzero;
464: for (k=0,i=1; i<=5; i++,k+=4)
465: {
466: p5 = gadd(gmul((GEN)p2[ind5[k]],(GEN)p2[ind5[k+1]]),
467: gmul((GEN)p2[ind5[k+2]],(GEN)p2[ind5[k+3]]));
468: p3 = gadd(p3, gmul(gsqr((GEN)p2[i]),p5));
469: }
470: w[l]=lrndtoi(greal(p3),&e);
471: e1 = gexpo(gimag(p3)); if (e1>e) e=e1;
472: ee[l]=e; z[l] = (long)p3;
473: }
1.2 ! noro 474: p5 = roots_to_ZX(z, &e);
1.1 noro 475: if (e <= -10) break;
476: prec = (prec<<1)-2;
477: }
1.2 ! noro 478: if (!ZX_is_squarefree(p5)) goto tchi;
1.1 noro 479: p3=(GEN)factor(p5)[1];
1.2 ! noro 480: f=carreparfait(ZX_disc(x));
1.1 noro 481: if (lg(p3)-1==1)
482: {
1.2 ! noro 483: avma = av;
! 484: return f? _res(60,1,4): _res(120,-1,5);
1.1 noro 485: }
1.2 ! noro 486: if (!f) { avma = av; return _res(20,-1,3); }
! 487:
1.1 noro 488: pr = - (bit_accuracy(prec) >> 1);
489: for (l=1; l<=6; l++)
490: if (ee[l] <= pr && gcmp0(poleval(p5,(GEN)w[l]))) break;
491: if (l>6) err(bugparier,"galois (bug4)");
492: p2=(l==6)? transroot(p1,2,5):transroot(p1,1,l);
493: p3=gzero;
494: for (i=1; i<=5; i++)
495: {
1.2 ! noro 496: j = i+1; if (j>5) j -= 5;
! 497: p3 = gadd(p3,gmul(gmul((GEN)p2[i],(GEN)p2[j]),
! 498: gsub((GEN)p2[j],(GEN)p2[i])));
1.1 noro 499: }
500: p5=gsqr(p3); p4=grndtoi(greal(p5),&e);
501: e1 = gexpo(gimag(p5)); if (e1>e) e=e1;
502: if (e <= -10)
503: {
504: if (gcmp0(p4)) goto tchi;
1.2 ! noro 505: f = carreparfait(p4); avma = av;
! 506: return f? _res(5,1,1): _res(10,1,2);
1.1 noro 507: }
508: prec=(prec<<1)-2;
509: }
510:
511: case 6: z = cgetg(7, t_VEC);
1.2 ! noro 512: prec = DEFAULTPREC + (long)((gexpo(cb)*42. / BITS_IN_LONG));
1.1 noro 513: for(;;)
514: {
515: for(;;)
516: {
517: p1=roots(x,prec);
518: for (l=1; l<=6; l++)
519: {
520: p2=(l==1)?p1:transroot(p1,1,l);
521: p3=gzero; k=0;
522: for (i=1; i<=5; i++) for (j=i+1; j<=6; j++)
523: {
524: p5=gadd(gmul((GEN)p2[ind6[k]],(GEN)p2[ind6[k+1]]),
525: gmul((GEN)p2[ind6[k+2]],(GEN)p2[ind6[k+3]]));
1.2 ! noro 526: p3=gadd(p3,gmul(gsqr(gmul((GEN)p2[i],(GEN)p2[j])),p5));
! 527: k += 4;
1.1 noro 528: }
529: z[l] = (long)p3;
530: }
1.2 ! noro 531: p5 = roots_to_ZX(z, &e);
1.1 noro 532: if (e <= -10) break;
533: prec=(prec<<1)-2;
534: }
1.2 ! noro 535: if (!ZX_is_squarefree(p5)) goto tchi;
1.1 noro 536: p2=(GEN)factor(p5)[1];
537: switch(lg(p2)-1)
538: {
539: case 1:
540: z = cgetg(11,t_VEC); ind=0;
541: p3=gadd(gmul(gmul((GEN)p1[1],(GEN)p1[2]),(GEN)p1[3]),
542: gmul(gmul((GEN)p1[4],(GEN)p1[5]),(GEN)p1[6]));
543: z[++ind] = (long)p3;
544: for (i=1; i<=3; i++)
545: for (j=4; j<=6; j++)
546: {
547: p2=transroot(p1,i,j);
548: p3=gadd(gmul(gmul((GEN)p2[1],(GEN)p2[2]),(GEN)p2[3]),
549: gmul(gmul((GEN)p2[4],(GEN)p2[5]),(GEN)p2[6]));
550: z[++ind] = (long)p3;
551: }
1.2 ! noro 552: p5 = roots_to_ZX(z, &e);
1.1 noro 553: if (e <= -10)
554: {
1.2 ! noro 555: if (!ZX_is_squarefree(p5)) goto tchi;
! 556: p2 = (GEN)factor(p5)[1];
! 557: f = carreparfait(ZX_disc(x));
! 558: avma = av;
1.1 noro 559: if (lg(p2)-1==1)
1.2 ! noro 560: return f? _res(360,1,15): _res(720,-1,16);
1.1 noro 561: else
1.2 ! noro 562: return f? _res(36,1,10): _res(72,-1,13);
1.1 noro 563: }
564: prec=(prec<<1)-2; break;
565:
566: case 2: l2=degpol(p2[1]); if (l2>3) l2=6-l2;
567: switch(l2)
568: {
1.2 ! noro 569: case 1: f = carreparfait(ZX_disc(x));
! 570: avma = av;
! 571: return f? _res(60,1,12): _res(120,-1,14);
! 572: case 2: f = carreparfait(ZX_disc(x));
! 573: if (f) { avma = av; return _res(24,1,7); }
! 574: p3 = (degpol(p2[1])==2)? (GEN)p2[2]: (GEN)p2[1];
! 575: f = carreparfait(ZX_disc(p3));
! 576: avma = av;
! 577: return f? _res(24,-1,6): _res(48,-1,11);
! 578: case 3: f = carreparfait(ZX_disc((GEN)p2[1]))
! 579: || carreparfait(ZX_disc((GEN)p2[2]));
! 580: avma = av;
! 581: return f? _res(18,-1,5): _res(36,-1,9);
1.1 noro 582: }
583: case 3:
584: for (l2=1; l2<=3; l2++)
1.2 ! noro 585: if (degpol(p2[l2]) >= 3) p3 = (GEN)p2[l2];
! 586: if (degpol(p3) == 3)
1.1 noro 587: {
1.2 ! noro 588: f = carreparfait(ZX_disc(p3)); avma = av;
! 589: return f? _res(6,-1,1): _res(12,-1,3);
1.1 noro 590: }
591: else
592: {
1.2 ! noro 593: f = carreparfait(ZX_disc(x)); avma = av;
! 594: return f? _res(12,1,4): _res(24,-1,8);
1.1 noro 595: }
1.2 ! noro 596: case 4: avma = av; return _res(6,-1,2);
1.1 noro 597: default: err(bugparier,"galois (bug3)");
598: }
599: }
600:
601: case 7: z = cgetg(36,t_VEC);
1.2 ! noro 602: prec = DEFAULTPREC + (long)((gexpo(cb)*7. / BITS_IN_LONG));
1.1 noro 603: for(;;)
604: {
1.2 ! noro 605: ind = 0; p1=roots(x,prec);
1.1 noro 606: for (i=1; i<=5; i++)
607: for (j=i+1; j<=6; j++)
608: {
1.2 ! noro 609: GEN t = gadd((GEN)p1[i],(GEN)p1[j]);
! 610: for (k=j+1; k<=7; k++) z[++ind] = ladd(t, (GEN)p1[k]);
1.1 noro 611: }
1.2 ! noro 612: p5 = roots_to_ZX(z, &e);
1.1 noro 613: if (e <= -10) break;
614: prec = (prec<<1)-2;
615: }
1.2 ! noro 616: if (!ZX_is_squarefree(p5)) goto tchi;
1.1 noro 617: p2=(GEN)factor(p5)[1];
618: switch(lg(p2)-1)
619: {
1.2 ! noro 620: case 1: f = carreparfait(ZX_disc(x)); avma = av;
! 621: return f? _res(2520,1,6): _res(5040,-1,7);
! 622: case 2: f = degpol(p2[1]); avma = av;
! 623: return (f==7 || f==28)? _res(168,1,5): _res(42,-1,4);
! 624: case 3: avma = av; return _res(21,1,3);
! 625: case 4: avma = av; return _res(14,-1,2);
! 626: case 5: avma = av; return _res(7,1,1);
1.1 noro 627: default: err(talker,"galois (bug2)");
628: }
629: }
1.2 ! noro 630: tchi: avma = av1; x = tschirnhaus(x1);
1.1 noro 631: }
632: }
633:
634: GEN
635: galoisapply(GEN nf, GEN aut, GEN x)
636: {
1.2 ! noro 637: gpmem_t av=avma,tetpil;
! 638: long lx,j,N;
1.1 noro 639: GEN p,p1,y,pol;
640:
641: nf=checknf(nf); pol=(GEN)nf[1];
642: if (typ(aut)==t_POL) aut = gmodulcp(aut,pol);
643: else
644: {
645: if (typ(aut)!=t_POLMOD || !gegal((GEN)aut[1],pol))
646: err(talker,"incorrect galois automorphism in galoisapply");
647: }
648: switch(typ(x))
649: {
650: case t_INT: case t_INTMOD: case t_FRAC: case t_FRACN: case t_PADIC:
651: avma=av; return gcopy(x);
652:
653: case t_POLMOD: x = (GEN) x[2]; /* fall through */
654: case t_POL:
655: p1 = gsubst(x,varn(pol),aut);
656: if (typ(p1)!=t_POLMOD || !gegal((GEN)p1[1],pol)) p1 = gmodulcp(p1,pol);
657: return gerepileupto(av,p1);
658:
659: case t_VEC:
660: if (lg(x)==3)
661: {
662: y=cgetg(3,t_VEC);
663: y[1]=(long)galoisapply(nf,aut,(GEN)x[1]);
664: y[2]=lcopy((GEN)x[2]); return gerepileupto(av, y);
665: }
666: if (lg(x)!=6) err(typeer,"galoisapply");
667: y=cgetg(6,t_VEC); y[1]=x[1]; y[3]=x[3]; y[4]=x[4];
668: p = (GEN)x[1];
669: p1=centermod(galoisapply(nf,aut,(GEN)x[2]), p);
670: if (is_pm1(x[3]))
671: if (ggval(subres(gmul((GEN)nf[7],p1),pol), p) > itos((GEN)x[4]))
672: p1[1] = (signe(p1[1]) > 0)? lsub((GEN)p1[1], p)
673: : ladd((GEN)p1[1], p);
674: y[2]=(long)p1;
675: y[5]=(long)centermod(galoisapply(nf,aut,(GEN)x[5]), p);
676: return gerepilecopy(av,y);
677:
678: case t_COL:
679: N=degpol(pol);
680: if (lg(x)!=N+1) err(typeer,"galoisapply");
681: p1=galoisapply(nf,aut,gmul((GEN)nf[7],x)); tetpil=avma;
682: return gerepile(av,tetpil,algtobasis(nf,p1));
683:
684: case t_MAT:
685: lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
686: N=degpol(pol);
687: if (lg(x[1])!=N+1) err(typeer,"galoisapply");
688: p1=cgetg(lx,t_MAT);
689: for (j=1; j<lx; j++) p1[j]=(long)galoisapply(nf,aut,(GEN)x[j]);
690: if (lx==N+1) p1 = idealhermite(nf,p1);
691: return gerepileupto(av,p1);
692: }
693: err(typeer,"galoisapply");
694: return NULL; /* not reached */
695: }
696:
697: GEN pol_to_monic(GEN pol, GEN *lead);
698: int cmp_pol(GEN x, GEN y);
699:
700: GEN
1.2 ! noro 701: get_bnfpol(GEN x, GEN *bnf, GEN *nf)
! 702: {
! 703: *bnf = _checkbnf(x);
! 704: *nf = _checknf(x);
! 705: if (*nf) return (GEN)(*nf)[1];
! 706: if (typ(x) != t_POL) err(idealer1);
! 707: return x;
! 708: }
! 709:
! 710: GEN
1.1 noro 711: get_nfpol(GEN x, GEN *nf)
712: {
713: if (typ(x) == t_POL) { *nf = NULL; return x; }
714: *nf = checknf(x); return (GEN)(*nf)[1];
715: }
716:
717: /* if fliso test for isomorphism, for inclusion otherwise. */
718: static GEN
719: nfiso0(GEN a, GEN b, long fliso)
720: {
1.2 ! noro 721: gpmem_t av = avma;
1.1 noro 722: long n,m,i,vb,lx;
723: GEN nfa,nfb,p1,y,la,lb;
724:
725: a = get_nfpol(a, &nfa);
726: b = get_nfpol(b, &nfb);
727: if (fliso && nfa && !nfb) { p1=a; a=b; b=p1; p1=nfa; nfa=nfb; nfb=p1; }
728: m=degpol(a);
729: n=degpol(b); if (m<=0 || n<=0) err(constpoler,"nfiso or nfincl");
730: if (fliso)
731: { if (n!=m) return gzero; }
732: else
733: { if (n%m) return gzero; }
734:
735: if (nfb) lb = NULL; else b = pol_to_monic(b,&lb);
736: if (nfa) la = NULL; else a = pol_to_monic(a,&la);
737: if (nfa && nfb)
738: {
739: if (fliso)
740: {
741: if (!gegal((GEN)nfa[2],(GEN)nfb[2])
742: || !gegal((GEN)nfa[3],(GEN)nfb[3])) return gzero;
743: }
744: else
745: if (!divise((GEN)nfb[3], gpowgs((GEN)nfa[3],n/m))) return gzero;
746: }
747: else
748: {
1.2 ! noro 749: GEN da = nfa? (GEN)nfa[3]: ZX_disc(a);
! 750: GEN db = nfb? (GEN)nfb[3]: ZX_disc(b);
1.1 noro 751: if (fliso)
752: {
753: p1=gdiv(da,db);
754: if (typ(p1)==t_FRAC) p1=mulii((GEN)p1[1],(GEN)p1[2]);
755: if (!gcarreparfait(p1)) { avma=av; return gzero; }
756: }
757: else
758: {
759: long q=n/m;
760: GEN fa=factor(da), ex=(GEN)fa[2];
761: fa=(GEN)fa[1]; lx=lg(fa);
762: for (i=1; i<lx; i++)
763: if (mod2((GEN)ex[i]) && !divise(db,gpowgs((GEN)fa[i],q)))
764: { avma=av; return gzero; }
765: }
766: }
767: a = dummycopy(a); setvarn(a,0);
768: b = dummycopy(b); vb=varn(b);
769: if (nfb)
770: {
771: if (vb == 0) nfb = gsubst(nfb, 0, polx[MAXVARN]);
772: y = lift_intern(nfroots(nfb,a));
773: }
774: else
775: {
776: if (vb == 0) setvarn(b, fetch_var());
777: y = (GEN)polfnf(a,b)[1]; lx = lg(y);
778: for (i=1; i<lx; i++)
779: {
780: if (lgef(y[i]) != 4) { setlg(y,i); break; }
781: y[i] = (long)gneg_i(lift_intern(gmael(y,i,2)));
782: }
783: y = gen_sort(y, 0, cmp_pol);
784: settyp(y, t_VEC);
1.2 ! noro 785: if (vb == 0) (void)delete_var();
1.1 noro 786: }
787: lx = lg(y); if (lx==1) { avma=av; return gzero; }
788: for (i=1; i<lx; i++)
789: {
790: p1 = (GEN)y[i];
791: if (typ(p1) == t_POL) setvarn(p1, vb); else p1 = scalarpol(p1, vb);
792: if (lb) p1 = poleval(p1, gmul(polx[vb],lb));
793: y[i] = la? ldiv(p1,la): (long)p1;
794: }
795: return gerepilecopy(av,y);
796: }
797:
798: GEN
1.2 ! noro 799: nfisisom(GEN a, GEN b) { return nfiso0(a,b,1); }
1.1 noro 800:
801: GEN
1.2 ! noro 802: nfisincl(GEN a, GEN b) { return nfiso0(a,b,0); }
1.1 noro 803:
804: /*************************************************************************/
805: /** **/
806: /** INITALG **/
807: /** **/
808: /*************************************************************************/
1.2 ! noro 809:
! 810: GEN
! 811: get_roots(GEN x,long r1,long prec)
! 812: {
! 813: GEN roo = (typ(x)==t_VEC)? dummycopy(x): roots(x,prec);
! 814: long i, ru = (lg(roo)-1 + r1) >> 1;
! 815:
! 816: for (i=1; i<=r1; i++) roo[i]=lreal((GEN)roo[i]);
! 817: for ( ; i<=ru; i++) roo[i]=roo[(i<<1)-r1];
! 818: roo[0]=evaltyp(t_VEC)|evallg(ru+1); return roo;
! 819: }
1.1 noro 820:
821: /* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */
822: GEN
823: quicktrace(GEN x, GEN sym)
824: {
825: GEN p1 = gzero;
826: long i;
827:
828: if (signe(x))
829: {
830: sym--;
831: for (i=lgef(x)-1; i>1; i--)
832: p1 = gadd(p1, gmul((GEN)x[i],(GEN)sym[i]));
833: }
834: return p1;
835: }
836:
837: /* T = trace(w[i]), x = sum x[i] w[i], x[i] integer */
838: static GEN
839: trace_col(GEN x, GEN T)
840: {
1.2 ! noro 841: gpmem_t av = avma;
1.1 noro 842: GEN t = gzero;
843: long i, l = lg(x);
844:
845: t = mulii((GEN)x[1],(GEN)T[1]);
846: for (i=2; i<l; i++)
847: if (signe(x[i])) t = addii(t, mulii((GEN)x[i],(GEN)T[i]));
848: return gerepileuptoint(av, t);
849: }
850:
851: /* pol belonging to Z[x], return a monic polynomial generating the same field
852: * as pol (x-> ax+b)) set lead = NULL if pol was monic (after dividing
853: * by the content), and to to leading coeff otherwise.
854: * No garbage collecting done.
855: */
856: GEN
857: pol_to_monic(GEN pol, GEN *lead)
858: {
859: long n = lgef(pol)-1;
860:
861: if (n==1 || gcmp1((GEN)pol[n])) { *lead = NULL; return pol; }
1.2 ! noro 862: return primitive_pol_to_monic(primpart(pol), lead);
1.1 noro 863: }
864:
1.2 ! noro 865: /* Let T = (Tr(wi wj)), then T^(-1) Z^n = codifferent = coD/dcoD
! 866: * NB: det(T) = dK, TI = dK T^(-1) integral */
! 867: static GEN
! 868: make_MDI(GEN nf, GEN TI, GEN *coD, GEN *dcoD)
1.1 noro 869: {
1.2 ! noro 870: GEN c, MDI, dK = (GEN)nf[3];
1.1 noro 871:
1.2 ! noro 872: TI = Q_primitive_part(TI, &c);
! 873: *dcoD = c? diviiexact(dK, c): dK;
! 874: *coD = hnfmodid(TI, *dcoD);
! 875: MDI = ideal_two_elt(nf, *coD);
! 876: return c? gmul(c, MDI): MDI;
1.1 noro 877: }
878:
879: /* assume lg(nf) > 3 && typ(nf) = container [hopefully a genuine nf] */
880: long
881: nf_get_r1(GEN nf)
882: {
883: GEN x = (GEN)nf[2];
884: if (typ(x) != t_VEC || lg(x) != 3 || typ(x[1]) != t_INT)
885: err(talker,"false nf in nf_get_r1");
886: return itos((GEN)x[1]);
887: }
888: long
889: nf_get_r2(GEN nf)
890: {
891: GEN x = (GEN)nf[2];
892: if (typ(x) != t_VEC || lg(x) != 3 || typ(x[2]) != t_INT)
893: err(talker,"false nf in nf_get_r2");
894: return itos((GEN)x[2]);
895: }
1.2 ! noro 896: void
! 897: nf_get_sign(GEN nf, long *r1, long *r2)
! 898: {
! 899: GEN x = (GEN)nf[2];
! 900: if (typ(x) != t_VEC || lg(x) != 3
! 901: || typ(x[1]) != t_INT || typ(x[2]) != t_INT)
! 902: err(talker,"false nf in nf_get_sign");
! 903: *r1 = itos((GEN)x[1]);
! 904: *r2 = (degpol(nf[1]) - *r1) >> 1;
! 905: }
1.1 noro 906:
907: static GEN
1.2 ! noro 908: get_Tr(GEN mul, GEN x, GEN basden)
1.1 noro 909: {
1.2 ! noro 910: GEN tr,T,T1,sym, bas = (GEN)basden[1], den = (GEN)basden[2];
1.1 noro 911: long i,j,n = lg(bas)-1;
912: T = cgetg(n+1,t_MAT);
913: T1 = cgetg(n+1,t_COL);
914: sym = polsym(x, n-1);
915:
916: T1[1]=lstoi(n);
917: for (i=2; i<=n; i++)
918: {
919: tr = quicktrace((GEN)bas[i], sym);
920: if (den && den[i]) tr = diviiexact(tr,(GEN)den[i]);
921: T1[i] = (long)tr; /* tr(w[i]) */
922: }
923: T[1] = (long)T1;
924: for (i=2; i<=n; i++)
925: {
926: T[i] = lgetg(n+1,t_COL); coeff(T,1,i) = T1[i];
927: for (j=2; j<=i; j++)
928: coeff(T,i,j) = coeff(T,j,i) = (long)trace_col((GEN)mul[j+(i-1)*n], T1);
929: }
930: return T;
931: }
932:
933: /* return [bas[i]*denom(bas[i]), denom(bas[i])], denom 1 is given as NULL */
934: GEN
935: get_bas_den(GEN bas)
936: {
1.2 ! noro 937: GEN z,b,d,den, dbas = dummycopy(bas);
! 938: long i, l = lg(bas);
! 939: int power = 1;
1.1 noro 940: den = cgetg(l,t_VEC);
941: for (i=1; i<l; i++)
942: {
1.2 ! noro 943: b = Q_remove_denom((GEN)bas[i], &d);
! 944: dbas[i]= (long)b;
! 945: den[i] = (long)d; if (d) power = 0;
1.1 noro 946: }
1.2 ! noro 947: if (power) den = NULL; /* power basis */
1.1 noro 948: z = cgetg(3,t_VEC);
949: z[1] = (long)dbas;
950: z[2] = (long)den; return z;
951: }
952:
953: /* allow x or y = NULL (act as 1) */
954: static GEN
955: _mulii(GEN x, GEN y)
956: {
957: if (!x) return y;
958: if (!y) return x;
959: return mulii(x,y);
960: }
961:
962: GEN
1.2 ! noro 963: get_mul_table(GEN x,GEN basden,GEN invbas)
1.1 noro 964: {
965: long i,j, n = degpol(x);
966: GEN z,d, mul = cgetg(n*n+1,t_MAT), bas=(GEN)basden[1], den=(GEN)basden[2];
1.2 ! noro 967:
1.1 noro 968: for (j=1; j<=n*n; j++) mul[j]=lgetg(n+1,t_COL);
969: for (i=1; i<=n; i++)
970: for (j=i; j<=n; j++)
971: {
1.2 ! noro 972: gpmem_t av = avma;
1.1 noro 973: z = gres(gmul((GEN)bas[j],(GEN)bas[i]), x);
974: z = mulmat_pol(invbas, z); /* integral column */
975: if (den)
976: {
977: d = _mulii((GEN)den[i], (GEN)den[j]);
978: if (d) z = gdivexact(z, d);
979: }
980: mul[j+(i-1)*n] = mul[i+(j-1)*n] = lpileupto(av,z);
981: }
982: return mul;
983: }
984:
1.2 ! noro 985: /* as get_Tr, mul_table not precomputed */
! 986: static GEN
! 987: make_Tr(GEN x, GEN w)
! 988: {
! 989: long i,j, n = degpol(x);
! 990: GEN p1,p2,t,d;
! 991: GEN sym = cgetg(n+2,t_VEC);
! 992: GEN den = cgetg(n+1,t_VEC);
! 993: GEN T = cgetg(n+1,t_MAT);
! 994: gpmem_t av;
! 995:
! 996: sym = polsym(x, n-1);
! 997: p1 = get_bas_den(w);
! 998: w = (GEN)p1[1];
! 999: den = (GEN)p1[2];
! 1000: for (i=1; i<=n; i++)
! 1001: {
! 1002: p1 = cgetg(n+1,t_COL); T[i] = (long)p1;
! 1003: for (j=1; j<i ; j++) p1[j] = coeff(T,i,j);
! 1004: for ( ; j<=n; j++)
! 1005: {
! 1006: av = avma;
! 1007: p2 = gres(gmul((GEN)w[i],(GEN)w[j]), x);
! 1008: t = quicktrace(p2, sym);
! 1009: if (den)
! 1010: {
! 1011: d = _mulii((GEN)den[i],(GEN)den[j]);
! 1012: if (d) t = diviiexact(t, d);
! 1013: }
! 1014: p1[j] = lpileuptoint(av, t);
! 1015: }
! 1016: }
! 1017: return T;
! 1018: }
! 1019:
! 1020: /* compute roots so that _absolute_ accuracy of M >= prec [also holds for G] */
! 1021: static void
! 1022: get_roots_for_M(nffp_t *F)
! 1023: {
! 1024: long n, eBD, er, PREC;
! 1025:
! 1026: if (F->extraprec < 0)
! 1027: { /* not initialized yet */
! 1028: n = degpol(F->x);
! 1029: eBD = 1 + gexpo((GEN)F->basden[1]);
! 1030: er = 1 + gexpo(F->ro? F->ro: cauchy_bound(F->x));
! 1031: if (er < 0) er = 0;
! 1032: F->extraprec = ((n*er + eBD + (long)log2(n)) >> TWOPOTBITS_IN_LONG);
! 1033: }
! 1034:
! 1035: PREC = F->prec + F->extraprec;
! 1036: if (F->ro && gprecision((GEN)F->ro[1]) >= PREC) return;
! 1037: F->ro = get_roots(F->x, F->r1, PREC);
! 1038: }
! 1039:
! 1040: /* [bas[i]/den[i]]= integer basis. roo = real part of the roots */
! 1041: static void
! 1042: make_M(nffp_t *F, int trunc)
! 1043: {
! 1044: GEN bas = (GEN)F->basden[1], den = (GEN)F->basden[2], ro = F->ro;
! 1045: GEN m, d, M;
! 1046: long i, j, l = lg(ro), n = lg(bas);
! 1047: M = cgetg(n,t_MAT);
! 1048: m = cgetg(l, t_COL); M[1] = (long)m;
! 1049: for (i=1; i<l; i++) m[i] = un; /* bas[1] = 1 */
! 1050: for (j=2; j<n; j++)
! 1051: {
! 1052: m = cgetg(l,t_COL); M[j] = (long)m;
! 1053: for (i=1; i<l; i++) m[i] = (long)poleval((GEN)bas[j], (GEN)ro[i]);
! 1054: }
! 1055: if (den)
! 1056: {
! 1057: GEN invd, rd = cgetr(F->prec + F->extraprec);
! 1058: for (j=2; j<n; j++)
! 1059: {
! 1060: d = (GEN)den[j]; if (!d) continue;
! 1061: m = (GEN)M[j]; affir(d,rd); invd = ginv(rd);
! 1062: for (i=1; i<l; i++) m[i] = lmul((GEN)m[i], invd);
! 1063: }
! 1064: }
! 1065:
! 1066: if (trunc && gprecision(M) > F->prec)
! 1067: {
! 1068: M = gprec_w(M, F->prec);
! 1069: F->ro = gprec_w(ro,F->prec);
! 1070: }
! 1071: if (DEBUGLEVEL>4) msgtimer("matrix M");
! 1072: F->M = M;
! 1073: }
! 1074:
! 1075: /* return G real such that G~ * G = T_2 */
! 1076: static void
! 1077: make_G(nffp_t *F)
! 1078: {
! 1079: GEN G, m, g, r, M = F->M, sqrt2 = gsqrt(gdeux, F->prec + F->extraprec);
! 1080: long i, j, k, r1 = F->r1, l = lg(M);
! 1081:
! 1082: G = cgetg(l, t_MAT);
! 1083: for (j=1; j<l; j++)
! 1084: {
! 1085: g = cgetg(l, t_COL);
! 1086: G[j] = (long)g; m = (GEN)M[j];
! 1087: for (k=i=1; i<=r1; i++) g[k++] = m[i];
! 1088: for ( ; k < l; i++)
! 1089: {
! 1090: r = (GEN)m[i];
! 1091: if (typ(r) == t_COMPLEX)
! 1092: {
! 1093: g[k++] = lmpmul(sqrt2, (GEN)r[1]);
! 1094: g[k++] = lmpmul(sqrt2, (GEN)r[2]);
! 1095: }
! 1096: else
! 1097: {
! 1098: g[k++] = lmpmul(sqrt2, r);
! 1099: g[k++] = zero;
! 1100: }
! 1101: }
! 1102: }
! 1103: F->G = G;
! 1104: }
! 1105:
! 1106: static void
! 1107: make_M_G(nffp_t *F, int trunc)
! 1108: {
! 1109: get_roots_for_M(F);
! 1110: make_M(F, trunc);
! 1111: make_G(F);
! 1112: }
! 1113:
1.1 noro 1114: void
1.2 ! noro 1115: remake_GM(GEN nf, nffp_t *F, long prec)
! 1116: {
! 1117: F->x = (GEN)nf[1];
! 1118: F->ro = NULL;
! 1119: F->r1 = nf_get_r1(nf);
! 1120: F->basden = get_bas_den((GEN)nf[7]);
! 1121: F->extraprec = -1;
! 1122: F->prec = prec; make_M_G(F, 1);
! 1123: }
! 1124:
! 1125: static void
! 1126: nffp_init(nffp_t *F, nfbasic_t *T, GEN ro, long prec)
! 1127: {
! 1128: F->x = T->x;
! 1129: F->ro = ro;
! 1130: F->r1 = T->r1;
! 1131: if (!T->basden) T->basden = get_bas_den(T->bas);
! 1132: F->basden = T->basden;
! 1133: F->extraprec = -1;
! 1134: F->prec = prec;
! 1135: }
! 1136:
! 1137: static void
! 1138: get_nf_fp_compo(nfbasic_t *T, nffp_t *F, GEN ro, long prec)
! 1139: {
! 1140: nffp_init(F,T,ro,prec);
! 1141: make_M_G(F, 1);
! 1142: }
! 1143:
! 1144: static GEN
! 1145: get_sign(long r1, long n)
! 1146: {
! 1147: GEN s = cgetg(3, t_VEC);
! 1148: s[1] = lstoi(r1);
! 1149: s[2] = lstoi((n - r1) >> 1); return s;
! 1150: }
! 1151:
! 1152: GEN
! 1153: nfbasic_to_nf(nfbasic_t *T, GEN ro, long prec)
1.1 noro 1154: {
1.2 ! noro 1155: GEN nf = cgetg(10,t_VEC);
! 1156: GEN x = T->x;
! 1157: GEN invbas,Tr,D,TI,A,dA, mat = cgetg(8,t_VEC);
! 1158: nffp_t F;
! 1159: get_nf_fp_compo(T, &F, ro, prec);
! 1160:
! 1161: nf[1] = (long)T->x;
! 1162: nf[2] = (long)get_sign(T->r1, degpol(T->x));
! 1163: nf[3] = (long)T->dK;
! 1164: nf[4] = (long)T->index;
! 1165: nf[6] = (long)F.ro;
! 1166: nf[5] = (long)mat;
! 1167: nf[7] = (long)T->bas;
1.1 noro 1168:
1.2 ! noro 1169: mat[1] = (long)F.M;
! 1170: mat[2] = (long)F.G;
! 1171:
! 1172: invbas = QM_inv(vecpol_to_mat(T->bas, lg(T->bas)-1), gun);
! 1173: nf[8] = (long)invbas;
! 1174: nf[9] = (long)get_mul_table(x, F.basden, invbas);
1.1 noro 1175: if (DEBUGLEVEL) msgtimer("mult. table");
1.2 ! noro 1176:
! 1177: Tr = get_Tr((GEN)nf[9], x, F.basden);
! 1178: TI = ZM_inv(Tr, T->dK); /* dK T^-1 */
! 1179: mat[6] = (long)TI;
! 1180: mat[7] = (long)make_MDI(nf,TI, &A, &dA); /* needed in idealinv below */
! 1181: if (is_pm1(T->index))
1.1 noro 1182: D = idealhermite_aux(nf, derivpol(x));
1183: else
1184: D = gmul(dA, idealinv(nf, A));
1.2 ! noro 1185: mat[3] = zero; /* FIXME: was gram matrix of current mat[2]. Useless */
! 1186: mat[4] = (long)Tr;
! 1187: mat[5] = (long)D;
1.1 noro 1188: if (DEBUGLEVEL) msgtimer("matrices");
1.2 ! noro 1189: return nf;
! 1190: }
! 1191:
! 1192: static GEN
! 1193: hnffromLLL(GEN nf)
! 1194: {
! 1195: GEN d, x;
! 1196: x = vecpol_to_mat((GEN)nf[7], degpol(nf[1]));
! 1197: x = Q_remove_denom(x, &d);
! 1198: if (!d) return x; /* power basis */
! 1199: return gauss(hnfmodid(x, d), x);
! 1200: }
! 1201:
! 1202: static GEN
! 1203: nfbasechange(GEN u, GEN x)
! 1204: {
! 1205: long i,lx;
! 1206: GEN y;
! 1207: switch(typ(x))
! 1208: {
! 1209: case t_COL: /* nfelt */
! 1210: return gmul(u, x);
! 1211:
! 1212: case t_MAT: /* ideal */
! 1213: y = dummycopy(x);
! 1214: lx = lg(x);
! 1215: for (i=1; i<lx; i++) y[i] = lmul(u, (GEN)y[i]);
! 1216: break;
! 1217:
! 1218: case t_VEC: /* pr */
! 1219: checkprimeid(x);
! 1220: y = dummycopy(x);
! 1221: y[2] = lmul(u, (GEN)y[2]);
! 1222: y[5] = lmul(u, (GEN)y[5]);
! 1223: break;
! 1224: default: y = x;
! 1225: }
! 1226: return y;
1.1 noro 1227: }
1228:
1.2 ! noro 1229: GEN
! 1230: nffromhnfbasis(GEN nf, GEN x)
! 1231: {
! 1232: long tx = typ(x);
! 1233: gpmem_t av = avma;
! 1234: GEN u;
! 1235: if (!is_vec_t(tx)) return gcopy(x);
! 1236: nf = checknf(nf);
! 1237: u = hnffromLLL(nf);
! 1238: return gerepilecopy(av, nfbasechange(u, x));
! 1239: }
1.1 noro 1240:
1241: GEN
1.2 ! noro 1242: nftohnfbasis(GEN nf, GEN x)
! 1243: {
! 1244: long tx = typ(x);
! 1245: gpmem_t av = avma;
! 1246: GEN u;
! 1247: if (!is_vec_t(tx)) return gcopy(x);
! 1248: nf = checknf(nf);
! 1249: u = ZM_inv(hnffromLLL(nf), gun);
! 1250: return gerepilecopy(av, nfbasechange(u, x));
! 1251: }
! 1252:
! 1253: static GEN
! 1254: get_red_G(nfbasic_t *T, GEN *pro)
! 1255: {
! 1256: GEN G, u, u0 = NULL;
! 1257: gpmem_t av;
! 1258: long i, prec, extraprec, n = degpol(T->x);
! 1259: nffp_t F;
! 1260:
! 1261: extraprec = (long) (0.5 * n * sizeof(long) / 8.);
! 1262: prec = DEFAULTPREC + extraprec;
! 1263: nffp_init(&F, T, *pro, prec);
! 1264: av = avma;
! 1265: for (i=1; ; i++)
! 1266: {
! 1267: F.prec = prec; make_M_G(&F, 0); G = F.G;
! 1268: if (u0) G = gmul(G, u0);
! 1269: if (DEBUGLEVEL)
! 1270: fprintferr("get_red_G: starting LLL, prec = %ld (%ld + %ld)\n",
! 1271: prec + F.extraprec, prec, F.extraprec);
! 1272: if ((u = lllfp_marked(1, G, 100, 2, prec, 0)))
! 1273: {
! 1274: if (typ(u) == t_MAT) break;
! 1275: u = (GEN)u[1];
! 1276: if (u0) u0 = gerepileupto(av, gmul(u0,u));
! 1277: else u0 = gerepilecopy(av, u);
! 1278: }
! 1279: if (i == MAXITERPOL) err(accurer,"red_T2");
! 1280: prec = (prec<<1)-2 + (gexpo(u0) >> TWOPOTBITS_IN_LONG);
! 1281: F.ro = NULL;
! 1282: if (DEBUGLEVEL) err(warnprec,"get_red_G", prec);
! 1283: }
! 1284: *pro = F.ro; return u0? gmul(u0,u): u;
! 1285: }
! 1286:
! 1287: /* Return the base change matrix giving an LLL-reduced basis for the
! 1288: * integer basis of nf(T).
! 1289: * *ro = roots of x, computed to precision prec [or NULL] */
! 1290: static GEN
! 1291: get_LLL_basis(nfbasic_t *T, GEN *pro)
! 1292: {
! 1293: GEN u;
! 1294: if (T->r1 != degpol(T->x)) u = get_red_G(T, pro);
! 1295: else
! 1296: {
! 1297: u = lllint_marked(1, make_Tr(T->x, T->bas), 100, 1, &u,NULL,NULL);
! 1298: if (!u) u = idmat(1);
! 1299: }
! 1300: return u;
! 1301: }
! 1302:
! 1303: static void
! 1304: set_LLL_basis(nfbasic_t *T, GEN *pro)
! 1305: {
! 1306: T->bas = gmul(T->bas, get_LLL_basis(T, pro));
! 1307: T->basden = NULL; /* recompute */
! 1308: }
! 1309:
! 1310: typedef struct {
! 1311: GEN xbest, dxbest;
! 1312: long ind, indmax, indbest;
! 1313: } ok_pol_t;
! 1314:
! 1315: /* is xn better than x ? */
! 1316: static int
! 1317: better_pol(GEN xn, GEN dxn, GEN x, GEN dx)
! 1318: {
! 1319: int fl;
! 1320: if (!x) return 1;
! 1321: fl = absi_cmp(dxn, dx);
! 1322: return (fl < 0 || (!fl && gpolcomp(xn, x) < 0));
! 1323: }
! 1324:
! 1325: static GEN
! 1326: ok_pol(void *TT, GEN xn)
! 1327: {
! 1328: ok_pol_t *T = (ok_pol_t*)TT;
! 1329: GEN dxn;
! 1330:
! 1331: if (++T->ind > T->indmax && T->xbest) return T->xbest;
! 1332:
! 1333: if (!ZX_is_squarefree(xn)) return (T->ind == T->indmax)? T->xbest: NULL;
! 1334: if (DEBUGLEVEL>3) outerr(xn);
! 1335: dxn = ZX_disc(xn);
! 1336: if (better_pol(xn, dxn, T->xbest, T->dxbest))
! 1337: {
! 1338: T->dxbest = dxn; T->xbest = xn; T->indbest = T->ind;
! 1339: }
! 1340: if (T->ind >= T->indmax) return T->xbest;
! 1341: return NULL;
! 1342: }
! 1343:
! 1344: /* z in Z[X] with positive leading coeff. Set z := z(-X) or z(X) so that
! 1345: * second coeff is > 0. In place. */
! 1346: static int
! 1347: canon_pol(GEN z)
1.1 noro 1348: {
1.2 ! noro 1349: long i,s;
1.1 noro 1350:
1.2 ! noro 1351: for (i=lgef(z)-2; i>=2; i-=2)
1.1 noro 1352: {
1.2 ! noro 1353: s = signe(z[i]);
! 1354: if (s > 0)
1.1 noro 1355: {
1.2 ! noro 1356: for (; i>=2; i-=2) z[i]=lnegi((GEN)z[i]);
! 1357: return -1;
1.1 noro 1358: }
1.2 ! noro 1359: if (s) return 1;
! 1360: }
! 1361: return 0;
! 1362: }
! 1363:
! 1364: static GEN _polred(GEN x, GEN a, GEN *pta, FP_chk_fun *CHECK);
! 1365:
! 1366: /* Seek a simpler, polynomial pol defining the same number field as
! 1367: * x (assumed to be monic at this point) */
! 1368: static GEN
! 1369: nfpolred(int part, nfbasic_t *T)
! 1370: {
! 1371: GEN x = T->x, dx = T->dx, a = T->bas;
! 1372: GEN phi, xbest, dxbest, mat, d, rev;
! 1373: long i, n = lg(a)-1, v = varn(x);
! 1374: ok_pol_t O;
! 1375: FP_chk_fun chk;
! 1376:
! 1377: if (degpol(x) == 1) { T->x = gsub(polx[v],gun); return gun; }
! 1378:
! 1379: if (!dx) dx = mulii(T->dK, sqri(T->index));
! 1380:
! 1381: O.ind = 0;
! 1382: O.indmax = part? min(n,3): n;
! 1383: O.xbest = NULL;
! 1384: chk.f = &ok_pol;
! 1385: chk.data = (void*)&O;
! 1386: if (!_polred(x, a, NULL, &chk))
! 1387: err(talker,"you found a counter-example to a conjecture, please report!");
! 1388: xbest = O.xbest; dxbest = O.dxbest;
! 1389: if (!better_pol(xbest, dxbest, x, dx)) return NULL; /* no improvement */
! 1390:
! 1391: /* update T */
! 1392: phi = (GEN)a[O.indbest];
! 1393: if (canon_pol(xbest) < 0) phi = gneg_i(phi);
! 1394: if (DEBUGLEVEL>1) fprintferr("xbest = %Z\n",xbest);
! 1395: rev = modreverse_i(phi, x);
! 1396: for (i=1; i<=n; i++) a[i] = (long)RX_RXQ_compo((GEN)a[i], rev, xbest);
! 1397: mat = vecpol_to_mat(Q_remove_denom(a, &d), n);
! 1398: if (!is_pm1(d)) mat = gdiv(hnfmodid(mat,d), d); else mat = idmat(n);
! 1399:
! 1400: (void)carrecomplet(diviiexact(dxbest,T->dK), &(T->index));
! 1401: T->bas= mat_to_vecpol(mat, v);
! 1402: T->dx = dxbest;
! 1403: T->x = xbest; return rev;
! 1404: }
! 1405:
! 1406: GEN
! 1407: get_nfindex(GEN bas)
! 1408: {
! 1409: gpmem_t av = avma;
! 1410: long n = lg(bas)-1;
! 1411: GEN d, mat = vecpol_to_mat(Q_remove_denom(bas, &d), n);
! 1412: if (!d) { avma = av; return gun; }
! 1413: return gerepileuptoint(av, diviiexact(gpowgs(d, n), det(mat)));
! 1414: }
! 1415:
! 1416: void
! 1417: nfbasic_init(GEN x, long flag, GEN fa, nfbasic_t *T)
! 1418: {
! 1419: GEN bas, dK, dx, index;
! 1420: long r1;
! 1421:
! 1422: T->basden = NULL;
! 1423: T->lead = NULL;
! 1424: if (DEBUGLEVEL) (void)timer2();
! 1425: if (typ(x) == t_POL)
! 1426: {
! 1427: check_pol_int(x, "nfinit");
! 1428: if (gisirreducible(x) == gzero) err(redpoler, "nfinit");
! 1429: x = pol_to_monic(x, &(T->lead));
! 1430: bas = allbase(x, flag, &dx, &dK, &index, &fa);
1.1 noro 1431: if (DEBUGLEVEL) msgtimer("round4");
1.2 ! noro 1432: r1 = sturm(x);
! 1433: }
! 1434: else if (typ(x) == t_VEC && lg(x)==3 && typ(x[1])==t_POL)
! 1435: { /* monic integral polynomial + integer basis */
! 1436: GEN mat;
! 1437: bas = (GEN)x[2]; x = (GEN)x[1];
! 1438: if (typ(bas) == t_MAT)
! 1439: { mat = bas; bas = mat_to_vecpol(mat,varn(x)); }
! 1440: else
! 1441: mat = vecpol_to_mat(bas, lg(bas)-1);
! 1442: index = get_nfindex(bas);
! 1443: dx = ZX_disc(x);
! 1444: dK = diviiexact(dx, sqri(index));
! 1445: r1 = sturm(x);
1.1 noro 1446: }
1447: else
1.2 ! noro 1448: { /* nf, bnf, bnr */
! 1449: GEN nf = checknf(x);
! 1450: x = (GEN)nf[1];
! 1451: dK = (GEN)nf[3];
! 1452: index = (GEN)nf[4];
! 1453: bas = (GEN)nf[7];
! 1454: r1 = nf_get_r1(nf); dx = NULL;
! 1455: }
! 1456:
! 1457: T->x = x;
! 1458: T->r1 = r1;
! 1459: T->dx = dx;
! 1460: T->dK = dK;
! 1461: T->bas = bas;
! 1462: T->index = index;
! 1463: }
! 1464:
! 1465: /* Initialize the number field defined by the polynomial x (in variable v)
! 1466: * flag & nf_RED: try a polred first.
! 1467: * flag & nf_PARTRED: as nf_RED but check the first two elements only.
! 1468: * flag & nf_ORIG
! 1469: * do a polred and return [nfinit(x), Mod(a,red)], where
! 1470: * Mod(a,red) = Mod(v,x) (i.e return the base change). */
! 1471: GEN
! 1472: _initalg(GEN x, long flag, long prec)
! 1473: {
! 1474: const gpmem_t av = avma;
! 1475: GEN nf, rev = NULL, ro = NULL;
! 1476: nfbasic_t T;
! 1477:
! 1478: nfbasic_init(x, flag, NULL, &T);
! 1479:
! 1480: if (T.lead && !(flag & (nf_RED|nf_PARTRED)))
! 1481: {
! 1482: err(warner,"non-monic polynomial. Result of the form [nf,c]");
! 1483: flag |= nf_PARTRED | nf_ORIG;
! 1484: }
! 1485: if (flag & (nf_RED|nf_PARTRED))
1.1 noro 1486: {
1.2 ! noro 1487: set_LLL_basis(&T, &ro);
! 1488: rev = nfpolred(flag & nf_PARTRED, &T);
! 1489: if (rev) ro = NULL; /* changed T.x */
! 1490: if (flag & nf_ORIG)
! 1491: {
! 1492: if (!rev) rev = polx[varn(T.x)]; /* no improvement */
! 1493: if (T.lead) rev = gdiv(rev, T.lead);
! 1494: rev = to_polmod(rev, T.x);
1.1 noro 1495: }
1496: if (DEBUGLEVEL) msgtimer("polred");
1497: }
1498:
1.2 ! noro 1499: set_LLL_basis(&T, &ro);
! 1500: if (DEBUGLEVEL) msgtimer("LLL basis");
1.1 noro 1501:
1.2 ! noro 1502: nf = nfbasic_to_nf(&T, ro, prec);
1.1 noro 1503: if (flag & nf_ORIG)
1504: {
1.2 ! noro 1505: GEN res = cgetg(3,t_VEC);
! 1506: res[1] = (long)nf;
! 1507: res[2] = (long)rev; nf = res;
1.1 noro 1508: }
1509: return gerepilecopy(av, nf);
1510: }
1511:
1512: GEN
1.2 ! noro 1513: initalgred(GEN x, long prec) { return _initalg(x, nf_RED, prec); }
! 1514: GEN
! 1515: initalgred2(GEN x, long prec) { return _initalg(x, nf_RED|nf_ORIG, prec); }
1.1 noro 1516: GEN
1.2 ! noro 1517: initalg(GEN x, long prec) { return _initalg(x, 0, prec); }
1.1 noro 1518:
1519: GEN
1520: nfinit0(GEN x, long flag,long prec)
1521: {
1522: switch(flag)
1523: {
1524: case 0:
1.2 ! noro 1525: case 1: return _initalg(x,0,prec);
! 1526: case 2: return _initalg(x,nf_RED,prec);
! 1527: case 3: return _initalg(x,nf_RED|nf_ORIG,prec);
! 1528: case 4: return _initalg(x,nf_PARTRED,prec);
! 1529: case 5: return _initalg(x,nf_PARTRED|nf_ORIG,prec);
1.1 noro 1530: default: err(flagerr,"nfinit");
1531: }
1532: return NULL; /* not reached */
1533: }
1534:
1535: /* assume x a bnr/bnf/nf */
1536: long
1537: nfgetprec(GEN x)
1538: {
1539: GEN nf = checknf(x), ro = (GEN)nf[6];
1.2 ! noro 1540: return (typ(ro)==t_VEC)? precision((GEN)ro[1]): DEFAULTPREC;
1.1 noro 1541: }
1542:
1543: GEN
1544: nfnewprec(GEN nf, long prec)
1545: {
1.2 ! noro 1546: gpmem_t av = avma;
! 1547: GEN NF;
! 1548: nffp_t F;
1.1 noro 1549:
1550: if (typ(nf) != t_VEC) err(talker,"incorrect nf in nfnewprec");
1551: if (lg(nf) == 11) return bnfnewprec(nf,prec);
1552: if (lg(nf) == 7) return bnrnewprec(nf,prec);
1553: (void)checknf(nf);
1.2 ! noro 1554: NF = dummycopy(nf);
! 1555: NF[5] = (long)dummycopy((GEN)NF[5]);
! 1556: remake_GM(NF, &F, prec);
! 1557: NF[6] = (long)F.ro;
! 1558: mael(NF,5,1) = (long)F.M;
! 1559: mael(NF,5,2) = (long)F.G;
! 1560: return gerepilecopy(av, NF);
! 1561: }
! 1562:
! 1563: /********************************************************************/
! 1564: /** **/
! 1565: /** POLRED **/
! 1566: /** **/
! 1567: /********************************************************************/
! 1568: /* remove duplicate polynomials in y, updating a (same indexes), in place */
! 1569: static void
! 1570: remove_duplicates(GEN y, GEN a)
! 1571: {
! 1572: long k, i, l = lg(y);
! 1573: gpmem_t av = avma;
! 1574: GEN z;
! 1575:
! 1576: if (l < 2) return;
! 1577: z = new_chunk(3);
! 1578: z[1] = (long)y;
! 1579: z[2] = (long)a; (void)sort_factor(z, gcmp);
! 1580: for (k=1, i=2; i<l; i++)
! 1581: if (!gegal((GEN)y[i], (GEN)y[k]))
! 1582: {
! 1583: k++;
! 1584: a[k] = a[i];
! 1585: y[k] = y[i];
! 1586: }
! 1587: l = k+1; setlg(a,l); setlg(y,l);
! 1588: avma = av;
! 1589: }
! 1590:
! 1591: /* if CHECK != NULL, return the first polynomial pol found such that
! 1592: * CHECK->f(CHECK->data, pol) != 0; return NULL if there are no such pol */
! 1593: static GEN
! 1594: _polred(GEN x, GEN a, GEN *pta, FP_chk_fun *CHECK)
! 1595: {
! 1596: long i, v = varn(x), l = lg(a);
! 1597: GEN ch, d, y = cgetg(l,t_VEC);
! 1598:
! 1599: for (i=1; i<l; i++)
! 1600: {
! 1601: if (DEBUGLEVEL>2) { fprintferr("i = %ld\n",i); flusherr(); }
! 1602: ch = QX_caract(x, (GEN)a[i], v);
! 1603: if (CHECK)
! 1604: {
! 1605: ch = CHECK->f(CHECK->data, ch);
! 1606: if (!ch) continue;
! 1607: return ch;
! 1608: }
! 1609: d = modulargcd(derivpol(ch), ch);
! 1610: if (degpol(d)) ch = gdivexact(ch,d);
! 1611:
! 1612: if (canon_pol(ch) < 0 && pta) a[i] = (long)gneg_i((GEN)a[i]);
! 1613: if (DEBUGLEVEL > 3) outerr(ch);
! 1614: y[i] = (long)ch;
! 1615: }
! 1616: if (CHECK) return NULL; /* no suitable polynomial found */
! 1617: remove_duplicates(y,a);
! 1618: if (pta) *pta = a;
! 1619: return y;
! 1620: }
! 1621:
! 1622: static GEN
! 1623: allpolred(GEN x, long flag, GEN fa, GEN *pta, FP_chk_fun *CHECK)
! 1624: {
! 1625: GEN ro = NULL;
! 1626: nfbasic_t T; nfbasic_init(x, flag, fa, &T);
! 1627: if (T.lead) err(impl,"polred for non-monic polynomial");
! 1628: set_LLL_basis(&T, &ro);
! 1629: return _polred(T.x, T.bas, pta, CHECK);
! 1630: }
! 1631:
! 1632: /* FIXME: backward compatibility */
! 1633: #define red_PARTIAL 1
! 1634: #define red_ORIG 2
! 1635:
! 1636: GEN
! 1637: polred0(GEN x, long flag, GEN fa)
! 1638: {
! 1639: gpmem_t av = avma;
! 1640: GEN y, a;
! 1641: int fl = 0;
! 1642:
! 1643: if (fa && gcmp0(fa)) fa = NULL; /* compatibility */
! 1644: if (flag & red_PARTIAL) fl |= nf_PARTIALFACT;
! 1645: if (flag & red_ORIG) fl |= nf_ORIG;
! 1646: y = allpolred(x, fl, fa, &a, NULL);
! 1647: if (fl & nf_ORIG) {
! 1648: GEN z = cgetg(3,t_MAT);
! 1649: z[1] = (long)a;
! 1650: z[2] = (long)y; y = z;
! 1651: }
1.1 noro 1652: return gerepilecopy(av, y);
1653: }
1654:
1.2 ! noro 1655: GEN
! 1656: ordred(GEN x)
! 1657: {
! 1658: gpmem_t av = avma;
! 1659: GEN y;
! 1660:
! 1661: if (typ(x) != t_POL) err(typeer,"ordred");
! 1662: if (!gcmp1(leading_term(x))) err(impl,"ordred");
! 1663: if (!signe(x)) return gcopy(x);
! 1664: y = cgetg(3,t_VEC);
! 1665: y[1] = (long)x;
! 1666: y[2] = (long)idmat(degpol(x));
! 1667: return gerepileupto(av, polred(y));
! 1668: }
! 1669:
! 1670: GEN
! 1671: polredfirstpol(GEN x, long flag, FP_chk_fun *CHECK)
! 1672: { return allpolred(x,flag,NULL,NULL,CHECK); }
! 1673: GEN
! 1674: polred(GEN x)
! 1675: { return polred0(x, 0, NULL); }
! 1676: GEN
! 1677: smallpolred(GEN x)
! 1678: { return polred0(x, red_PARTIAL, NULL); }
! 1679: GEN
! 1680: factoredpolred(GEN x, GEN fa)
! 1681: { return polred0(x, 0, fa); }
! 1682: GEN
! 1683: polred2(GEN x)
! 1684: { return polred0(x, red_ORIG, NULL); }
! 1685: GEN
! 1686: smallpolred2(GEN x)
! 1687: { return polred0(x, red_PARTIAL|red_ORIG, NULL); }
! 1688: GEN
! 1689: factoredpolred2(GEN x, GEN fa)
! 1690: { return polred0(x, red_PARTIAL, fa); }
! 1691:
! 1692: /********************************************************************/
! 1693: /** **/
! 1694: /** POLREDABS **/
! 1695: /** **/
! 1696: /********************************************************************/
! 1697:
! 1698: GEN
! 1699: T2_from_embed_norm(GEN x, long r1)
! 1700: {
! 1701: GEN p = sum(x, 1, r1);
! 1702: GEN q = sum(x, r1+1, lg(x)-1);
! 1703: if (q != gzero) p = gadd(p, gmul2n(q,1));
! 1704: return p;
! 1705: }
! 1706:
! 1707: GEN
! 1708: T2_from_embed(GEN x, long r1)
! 1709: {
! 1710: return T2_from_embed_norm(gnorm(x), r1);
! 1711: }
! 1712:
! 1713: typedef struct {
! 1714: long r1, v, prec;
! 1715: GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */
! 1716: GEN u; /* matrix giving fincke-pohst-reduced Zk basis */
! 1717: GEN M; /* embeddings of initial (LLL-reduced) Zk basis */
! 1718: GEN bound; /* T2 norm of the polynomial defining nf */
! 1719: } CG_data;
! 1720:
! 1721: /* characteristic pol of x */
! 1722: static GEN
! 1723: get_polchar(CG_data *d, GEN x)
! 1724: {
! 1725: GEN g = gmul(d->ZKembed,x);
! 1726: long e;
! 1727: g = grndtoi(roots_to_pol_r1r2(g, d->r1, d->v), &e);
! 1728: if (e > -5) err(precer, "get_polchar");
! 1729: return g;
! 1730: }
! 1731:
! 1732: /* return a defining polynomial for Q(x) */
! 1733: static GEN
! 1734: get_polmin(CG_data *d, GEN x)
! 1735: {
! 1736: GEN g = get_polchar(d,x);
! 1737: GEN h = modulargcd(g,derivpol(g));
! 1738: if (degpol(h)) g = gdivexact(g,h);
! 1739: return g;
! 1740: }
! 1741:
! 1742: /* does x generate the correct field ? */
! 1743: static GEN
! 1744: chk_gen(void *data, GEN x)
! 1745: {
! 1746: gpmem_t av = avma;
! 1747: GEN g = get_polchar((CG_data*)data,x);
! 1748: GEN h = modulargcd(g,derivpol(g));
! 1749: if (degpol(h)) { avma = av; return NULL; }
! 1750: if (DEBUGLEVEL>3) fprintferr(" generator: %Z\n",g);
! 1751: return g;
! 1752: }
! 1753:
! 1754: /* mat = base change matrix, gram = LLL-reduced T2 matrix */
! 1755: static GEN
! 1756: chk_gen_init(FP_chk_fun *chk, GEN gram, GEN mat)
! 1757: {
! 1758: CG_data *d = (CG_data*)chk->data;
! 1759: GEN P,bound,prev,x,B;
! 1760: long l = lg(gram), N = l-1,i,prec,prec2;
! 1761: int skipfirst = 0;
! 1762:
! 1763: d->u = mat;
! 1764: d->ZKembed = gmul(d->M, mat);
! 1765:
! 1766: bound = d->bound;
! 1767: prev = NULL;
! 1768: x = zerocol(N);
! 1769: for (i=1; i<l; i++)
! 1770: {
! 1771: B = gcoeff(gram,i,i);
! 1772: if (gcmp(B,bound) >= 0) continue; /* don't assume increasing norms */
! 1773:
! 1774: x[i] = un; P = get_polmin(d,x);
! 1775: x[i] = zero;
! 1776: if (degpol(P) == N)
! 1777: {
! 1778: if (gcmp(B,bound) < 0) bound = B ;
! 1779: continue;
! 1780: }
! 1781: if (DEBUGLEVEL>2) fprintferr("chk_gen_init: subfield %Z\n",P);
! 1782: if (skipfirst != i-1) continue;
! 1783:
! 1784: /* Q(w[1],...,w[i-1]) = Q[X]/(prev) is a strict subfield of nf */
! 1785: if (prev && degpol(prev) < N && !gegal(prev,P))
! 1786: {
! 1787: if (degpol(prev) * degpol(P) > 64) continue; /* too expensive */
! 1788: P = compositum(prev,P);
! 1789: P = (GEN)P[lg(P)-1];
! 1790: if (degpol(P) == N) continue;
! 1791: if (DEBUGLEVEL>2 && degpol(P) != degpol(prev))
! 1792: fprintferr("chk_gen_init: subfield %Z\n",P);
! 1793: }
! 1794: skipfirst++; prev = P;
! 1795: }
! 1796: /* x_1,...,x_skipfirst generate a strict subfield [unless N=skipfirst=1] */
! 1797: chk->skipfirst = skipfirst;
! 1798: if (DEBUGLEVEL>2) fprintferr("chk_gen_init: skipfirst = %ld\n",skipfirst);
! 1799:
! 1800: /* should be gexpo( max_k C^n_k (bound/k)^(k/2) ) */
! 1801: prec2 = (1 + (((gexpo(bound)*N)/2) >> TWOPOTBITS_IN_LONG));
! 1802: if (prec2 < 0) prec2 = 0;
! 1803: prec = 3 + prec2;
! 1804: if (DEBUGLEVEL)
! 1805: fprintferr("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec);
! 1806: if (prec > d->prec) err(bugparier, "precision problem in polredabs");
! 1807: if (prec < d->prec) d->ZKembed = gprec_w(d->ZKembed, prec);
! 1808: return bound;
! 1809: }
! 1810:
! 1811: /* store phi(beta mod z). Suitable for gerepileupto */
! 1812: static GEN
! 1813: storeeval(GEN beta, GEN z, GEN lead)
! 1814: {
! 1815: GEN y = cgetg(3,t_VEC);
! 1816: z = gcopy(z);
! 1817: beta = lead? gdiv(beta, lead): gcopy(beta);
! 1818: y[1] = (long)z;
! 1819: y[2] = (long)to_polmod(beta,z);
! 1820: return y;
! 1821: }
! 1822:
! 1823: static GEN
! 1824: storeraw(GEN beta, GEN z)
! 1825: {
! 1826: GEN y = cgetg(3,t_VEC);
! 1827: y[1] = lcopy(z);
! 1828: y[2] = lcopy(beta); return y;
! 1829: }
! 1830:
! 1831: static void
! 1832: findmindisc(GEN *py, GEN *pa)
! 1833: {
! 1834: GEN v, dmin, z, b, discs, y = *py, a = *pa;
! 1835: long i,k, l = lg(y);
! 1836:
! 1837: if (l == 2) { *py = (GEN)y[1]; *pa = (GEN)a[1]; return; }
! 1838:
! 1839: discs = cgetg(l,t_VEC);
! 1840: for (i=1; i<l; i++) discs[i] = labsi(ZX_disc((GEN)y[i]));
! 1841: v = sindexsort(discs);
! 1842: k = v[1];
! 1843: dmin = (GEN)discs[k];
! 1844: z = (GEN)y[k];
! 1845: b = (GEN)a[k];
! 1846: for (i=2; i<l; i++)
! 1847: {
! 1848: k = v[i];
! 1849: if (!egalii((GEN)discs[k], dmin)) break;
! 1850: if (gpolcomp((GEN)y[k],z) < 0) { z = (GEN)y[k]; b = (GEN)a[k]; }
! 1851: }
! 1852: *py = z; *pa = b;
! 1853: }
! 1854:
! 1855: static GEN
! 1856: storepol(GEN x, GEN z, GEN a, GEN lead, long flag)
! 1857: {
! 1858: GEN y, b;
! 1859: if (flag & nf_RAW)
! 1860: y = storeraw(a, z);
! 1861: else if (flag & nf_ORIG)
! 1862: {
! 1863: b = modreverse_i(a, x);
! 1864: y = storeeval(b,z, lead);
! 1865: }
! 1866: else
! 1867: y = gcopy(z);
! 1868: return y;
! 1869: }
! 1870:
! 1871: static GEN
! 1872: storeallpol(GEN x, GEN z, GEN a, GEN lead, long flag)
! 1873: {
! 1874: GEN y, b;
! 1875:
! 1876: if (flag & nf_RAW)
! 1877: {
! 1878: long i, c = lg(z);
! 1879: y = cgetg(c,t_VEC);
! 1880: for (i=1; i<c; i++) y[i] = (long)storeraw((GEN)a[i], (GEN)z[i]);
! 1881: }
! 1882: else if (flag & nf_ORIG)
! 1883: {
! 1884: long i, c = lg(z);
! 1885: b = cgetg(c, t_VEC);
! 1886: for (i=1; i<c; i++)
! 1887: b[i] = (long)modreverse_i((GEN)a[i], x);
! 1888:
! 1889: y = cgetg(c,t_VEC);
! 1890: for (i=1; i<c; i++) y[i] = (long)storeeval((GEN)b[i], (GEN)z[i], lead);
! 1891: }
! 1892: else
! 1893: y = gcopy(z);
! 1894: return y;
! 1895: }
! 1896:
! 1897: GEN
! 1898: _polredabs(nfbasic_t *T, GEN *u)
! 1899: {
! 1900: long i, prec, e, n = degpol(T->x);
! 1901: GEN v, ro = NULL;
! 1902: FP_chk_fun chk = { &chk_gen, &chk_gen_init, NULL, 0 };
! 1903: nffp_t F;
! 1904: CG_data d; chk.data = (void*)&d;
! 1905:
! 1906: set_LLL_basis(T, &ro);
! 1907: if (DEBUGLEVEL) msgtimer("LLL basis");
! 1908:
! 1909: /* || polchar ||_oo < 2^e */
! 1910: e = n * (gexpo(gmulsg(n, cauchy_bound(T->x))) + 1);
! 1911: prec = DEFAULTPREC + (e >> TWOPOTBITS_IN_LONG);
! 1912:
! 1913: get_nf_fp_compo(T, &F, ro, prec);
! 1914:
! 1915: d.v = varn(T->x);
! 1916: d.r1= T->r1;
! 1917: d.bound = gmul(T2_from_embed(F.ro, d.r1), dbltor(1.00000001));
! 1918: for (i=1; ; i++)
! 1919: {
! 1920: GEN R = R_from_QR(F.G, prec);
! 1921: d.prec = prec;
! 1922: d.M = F.M;
! 1923: if (R)
! 1924: {
! 1925: v = fincke_pohst(_vec(R),NULL,5000,1, 0, &chk);
! 1926: if (v) break;
! 1927: }
! 1928: if (i==MAXITERPOL) err(accurer,"polredabs0");
! 1929: prec = (prec<<1)-2;
! 1930: get_nf_fp_compo(T, &F, NULL, prec);
! 1931: if (DEBUGLEVEL) err(warnprec,"polredabs0",prec);
! 1932: }
! 1933: *u = d.u; return v;
! 1934: }
! 1935:
! 1936: GEN
! 1937: polredabs0(GEN x, long flag)
! 1938: {
! 1939: long i, l, vx;
! 1940: gpmem_t av = avma;
! 1941: GEN y, a, u;
! 1942: nfbasic_t T;
! 1943:
! 1944: nfbasic_init(x, flag & nf_PARTIALFACT, NULL, &T);
! 1945: x = T.x; vx = varn(x);
! 1946:
! 1947: if (degpol(x) == 1)
! 1948: {
! 1949: u = NULL;
! 1950: y = _vec(polx[vx]);
! 1951: a = _vec(gsub((GEN)y[1], (GEN)x[2]));
! 1952: }
! 1953: else
! 1954: {
! 1955: GEN v = _polredabs(&T, &u);
! 1956: y = (GEN)v[1];
! 1957: a = (GEN)v[2];
! 1958: }
! 1959: l = lg(a);
! 1960: for (i=1; i<l; i++)
! 1961: if (canon_pol((GEN)y[i]) < 0) a[i] = (long)gneg_i((GEN)a[i]);
! 1962: remove_duplicates(y,a);
! 1963: l = lg(a);
! 1964: if (l == 1)
! 1965: {
! 1966: y = _vec(x);
! 1967: a = _vec(polx[vx]);
! 1968: }
! 1969: if (DEBUGLEVEL) fprintferr("%ld minimal vectors found.\n",l-1);
! 1970: if (flag & nf_ALL)
! 1971: {
! 1972: if (u) for (i=1; i<l; i++) a[i] = lmul(T.bas, gmul(u, (GEN)a[i]));
! 1973: y = storeallpol(x,y,a,T.lead,flag);
! 1974: }
! 1975: else
! 1976: {
! 1977: findmindisc(&y, &a);
! 1978: if (u) a = gmul(T.bas, gmul(u, a));
! 1979: y = storepol(x,y,a,T.lead,flag);
! 1980: }
! 1981: return gerepileupto(av, y);
! 1982: }
! 1983:
! 1984: GEN
! 1985: polredabsall(GEN x, long flun) { return polredabs0(x, flun | nf_ALL); }
! 1986: GEN
! 1987: polredabs(GEN x) { return polredabs0(x,0); }
! 1988: GEN
! 1989: polredabs2(GEN x) { return polredabs0(x,nf_ORIG); }
! 1990:
1.1 noro 1991: static long
1992: nf_pm1(GEN y)
1993: {
1994: long i,l;
1995:
1996: if (!is_pm1(y[1])) return 0;
1997: l = lg(y);
1998: for (i=2; i<l; i++)
1999: if (signe(y[i])) return 0;
2000: return signe(y[1]);
2001: }
2002:
2003: static GEN
2004: is_primitive_root(GEN nf, GEN fa, GEN x, long w)
2005: {
2006: GEN y, exp = stoi(2), pp = (GEN)fa[1];
2007: long i,p, l = lg(pp);
2008:
2009: for (i=1; i<l; i++)
2010: {
2011: p = itos((GEN)pp[i]);
2012: exp[2] = w / p; y = element_pow(nf,x,exp);
2013: if (nf_pm1(y) > 0) /* y = 1 */
2014: {
2015: if (p!=2 || !gcmp1(gcoeff(fa,i,2))) return NULL;
2016: x = gneg_i(x);
2017: }
2018: }
2019: return x;
2020: }
2021:
1.2 ! noro 2022: /* only roots of 1 are +/- 1 */
! 2023: static GEN
! 2024: trivroots(GEN nf)
! 2025: {
! 2026: GEN y = cgetg(3, t_VEC);
! 2027: y[1] = deux;
! 2028: y[2] = (long)gscalcol_i(negi(gun), degpol(nf[1]));
! 2029: return y;
! 2030: }
! 2031:
1.1 noro 2032: GEN
2033: rootsof1(GEN nf)
2034: {
1.2 ! noro 2035: gpmem_t av = avma;
1.1 noro 2036: long N,k,i,ws,prec;
1.2 ! noro 2037: GEN p1,y,d,list,w;
! 2038:
! 2039: nf = checknf(nf);
! 2040: if ( nf_get_r1(nf) ) return trivroots(nf);
1.1 noro 2041:
1.2 ! noro 2042: N = degpol(nf[1]); prec = nfgetprec(nf);
1.1 noro 2043: for (i=1; ; i++)
2044: {
1.2 ! noro 2045: GEN R = R_from_QR(gmael(nf,5,2), prec);
! 2046: if (R)
! 2047: {
! 2048: p1 = fincke_pohst(_vec(R),stoi(N),1000,1, 0, NULL);
! 2049: if (p1) break;
! 2050: }
1.1 noro 2051: if (i == MAXITERPOL) err(accurer,"rootsof1");
1.2 ! noro 2052: prec = (prec<<1)-2;
1.1 noro 2053: if (DEBUGLEVEL) err(warnprec,"rootsof1",prec);
1.2 ! noro 2054: nf = nfnewprec(nf,prec);
1.1 noro 2055: }
2056: if (itos(ground((GEN)p1[2])) != N) err(bugparier,"rootsof1 (bug1)");
1.2 ! noro 2057: w = (GEN)p1[1]; ws = itos(w);
! 2058: if (ws == 2) { avma = av; return trivroots(nf); }
1.1 noro 2059:
2060: d = decomp(w); list = (GEN)p1[3]; k = lg(list);
2061: for (i=1; i<k; i++)
2062: {
2063: p1 = (GEN)list[i];
2064: p1 = is_primitive_root(nf,d,p1,ws);
1.2 ! noro 2065: if (p1) break;
1.1 noro 2066: }
1.2 ! noro 2067: if (!p1) err(bugparier,"rootsof1");
! 2068: y = cgetg(3, t_VEC);
! 2069: y[1] = lstoi(ws);
! 2070: y[2] = (long)p1; return gerepilecopy(av, y);
1.1 noro 2071: }
2072:
2073: /*******************************************************************/
2074: /* */
2075: /* DEDEKIND ZETA FUNCTION */
2076: /* */
2077: /*******************************************************************/
2078: static GEN
2079: dirzetak0(GEN nf, long N0)
2080: {
2081: GEN vect,p1,pol,disc,c,c2;
1.2 ! noro 2082: gpmem_t av=avma;
! 2083: long i,j,k,limk,lx;
1.1 noro 2084: ulong q,p,rem;
2085: byteptr d=diffptr;
1.2 ! noro 2086: long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
1.1 noro 2087:
2088: pol=(GEN)nf[1]; disc=(GEN)nf[4];
2089: c = (GEN) gpmalloc((N0+1)*sizeof(long));
2090: c2 = (GEN) gpmalloc((N0+1)*sizeof(long));
2091: c2[0]=c[0]=evaltyp(t_VEC) | evallg(N0+1);
2092: c2[1]=c[1]=1; for (i=2; i<=N0; i++) c[i]=0;
2093: court[2] = 0;
2094:
2095: while (court[2]<=N0)
2096: {
1.2 ! noro 2097: NEXT_PRIME_VIADIFF_CHECK(court[2], d);
1.1 noro 2098: if (smodis(disc,court[2])) /* court does not divide index */
2099: { vect = (GEN) simplefactmod(pol,court)[1]; lx=lg(vect); }
2100: else
2101: {
2102: p1=primedec(nf,court); lx=lg(p1); vect=cgetg(lx,t_COL);
2103: for (i=1; i<lx; i++) vect[i]=mael(p1,i,4);
2104: }
2105: for (j=1; j<lx; j++)
2106: {
2107: p1=powgi(court, (GEN)vect[j]); /* p1 = court^f */
2108: if (cmpis(p1,N0) <= 0)
2109: {
2110: q=p=p1[2]; limk=N0/q;
2111: for (k=2; k<=N0; k++) c2[k]=c[k];
2112: while (q<=(ulong)N0)
2113: {
2114: for (k=1; k<=limk; k++) c2[k*q] += c[k];
2115: q = smulss(q,p,&rem);
2116: if (rem) break;
2117: limk /= p;
2118: }
2119: p1=c; c=c2; c2=p1;
2120: }
2121: }
2122: avma=av;
2123: if (DEBUGLEVEL>6) fprintferr(" %ld",court[2]);
2124: }
2125: if (DEBUGLEVEL>6) { fprintferr("\n"); flusherr(); }
2126: free(c2); return c;
2127: }
2128:
2129: GEN
2130: dirzetak(GEN nf, GEN b)
2131: {
2132: GEN z,c;
2133: long i;
2134:
2135: if (typ(b)!=t_INT) err(talker,"not an integer type in dirzetak");
2136: if (signe(b)<=0) return cgetg(1,t_VEC);
2137: nf = checknf(nf);
2138: if (is_bigint(b)) err(talker,"too many terms in dirzetak");
2139: c = dirzetak0(nf,itos(b));
2140: i = lg(c); z=cgetg(i,t_VEC);
2141: for (i-- ; i; i--) z[i]=lstoi(c[i]);
2142: free(c); return z;
2143: }
2144:
2145: GEN
2146: initzeta(GEN pol, long prec)
2147: {
2148: GEN nfz,nf,alpha,beta,mu,gr1,gr2,gru,p1,p2,cst,A0,c0,c1,c2,eps,coef;
2149: GEN limx,bnf,resi,zet,C,coeflog,racpi,aij,tabj,colzero, *tabcstn, *tabcstni;
2150: GEN c_even,ck_even,c_odd,ck_odd,serie_even,serie_odd,serie_exp,Pi;
1.2 ! noro 2151: long N0,imin,imax,r1,r2,ru,R,N,i,j,k,n;
! 2152: gpmem_t av,av2,tetpil;
! 2153: long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
1.1 noro 2154: stackzone *zone, *zone0, *zone1;
2155:
2156: /*************** Calcul du residu et des constantes ***************/
2157: eps=gmul2n(gun,-bit_accuracy(prec)-6); p1=dbltor(0.5);
2158: nfz=cgetg(10,t_VEC);
1.2 ! noro 2159: bnf = bnfinit0(pol,2,NULL,prec+1); prec=(prec<<1)-1;
1.1 noro 2160: bnf = checkbnf_discard(bnf);
2161: Pi = mppi(prec); racpi=gsqrt(Pi,prec);
2162:
2163: /* Nb de classes et regulateur */
2164: nf=(GEN)bnf[7]; N=degpol(nf[1]);
2165: r1 = nf_get_r1(nf); r2 = (N-r1)>>1;
2166: gr1=gmael(nf,2,1); gr2=gmael(nf,2,2);
2167: ru=r1+r2; R=ru+2;
2168: av=avma; p1=(GEN)bnf[8]; p2 = gmul(gmul2n(gmael(p1,1,1),r1), (GEN)p1[2]);
2169: tetpil = avma; resi=gerepile(av,tetpil,gdiv(p2, gmael(p1,4,1)));
2170:
2171: /* Calcul de N0 */
2172: cst = cgetr(prec); av = avma;
2173: mu = gadd(gmul2n(gr1,-1),gr2);
2174: alpha = gmul2n(stoi(ru+1),-1);
2175: beta = gpui(gdeux,gmul2n(gr1,-1),DEFAULTPREC);
2176: A0 = gmul2n(gpowgs(mu,R),r1);
2177: A0 = gmul(A0,gpowgs(gmul2n(Pi,1),1-ru));
2178: A0 = gsqrt(A0,DEFAULTPREC);
2179:
2180: c1 = gmul(mu,gpui(beta,ginv(mu),DEFAULTPREC));
2181: c0 = gdiv(gmul(A0,gpowgs(gmul2n(Pi,1),ru-1)),mu);
2182: c0 = gmul(c0,gpui(c1,gneg_i(alpha),DEFAULTPREC));
2183: c2 = gdiv(alpha,mu);
2184:
2185: p1 = glog(gdiv(c0,eps),DEFAULTPREC);
2186: limx = gdiv(gsub(glog(p1,DEFAULTPREC),glog(c1,DEFAULTPREC)),
2187: gadd(c2,gdiv(p1,mu)));
2188: limx = gmul(gpui(gdiv(c1,p1),mu,DEFAULTPREC),
2189: gadd(gun,gmul(alpha,limx)));
2190: p1 = gsqrt(absi((GEN)nf[3]),prec);
2191: p2 = gmul2n(gpowgs(racpi,N),r2);
2192: gaffect(gdiv(p1,p2), cst);
2193:
2194: av = avma; p1 = gfloor(gdiv(cst,limx)); N0 = p1[2];
2195: if (is_bigint(p1) || N0 > 10000000)
2196: err(talker,"discriminant too large for initzeta, sorry");
2197: if (DEBUGLEVEL>=2)
1.2 ! noro 2198: { fprintferr("\ninitzeta:\nN0 = %ld\n",N0); flusherr(); (void)timer2(); }
1.1 noro 2199:
2200: /* Calcul de imax */
2201:
2202: imin=1; imax=1400;
2203: p1 = gmul(gpowgs(gmul2n(racpi,1),r2),gpowgs(stoi(5),r1));
2204: p1 = gdiv(p1,gmul(gmul(gsqrt(limx,DEFAULTPREC),gmul2n(eps,4)),
2205: gpowgs(racpi,3)));
2206: while (imax-imin >= 4)
2207: {
2208: long itest = (imax+imin)>>1;
2209: p2 = gmul(gpowgs(mpfactr(itest,DEFAULTPREC),r2),gpowgs(limx,itest));
2210: p2 = gmul(p2,gpowgs(mpfactr(itest/2,DEFAULTPREC),r1));
2211: if (gcmp(p2,p1) >= 0) imax=itest; else imin=itest;
2212: }
2213: imax -= (imax & 1); avma = av;
2214: if (DEBUGLEVEL>=2) { fprintferr("imax = %ld\n",imax); flusherr(); }
2215: /* Tableau des i/cst (i=1 a N0) */
2216:
2217: i = prec*N0;
2218: zone = switch_stack(NULL,i + 2*(N0+1) + 6*prec);
2219: zone1 = switch_stack(NULL,2*i);
2220: zone0 = switch_stack(NULL,2*i);
1.2 ! noro 2221: (void)switch_stack(zone,1);
1.1 noro 2222: tabcstn = (GEN*) cgetg(N0+1,t_VEC);
2223: tabcstni = (GEN*) cgetg(N0+1,t_VEC);
2224: p1 = ginv(cst);
2225: for (i=1; i<=N0; i++) { tabcstni[i] = gun; tabcstn[i] = mulsr(i,p1); }
1.2 ! noro 2226: (void)switch_stack(zone,0);
1.1 noro 2227:
2228: /********** Calcul des coefficients a(i,j) independants de s **********/
2229:
2230: zet=cgetg(R,t_VEC); zet[1] = lmpeuler(prec);
2231: for (i=2; i<R; i++) zet[i] = (long)gzeta(stoi(i),prec);
2232:
2233: aij=cgetg(imax+1,t_VEC);
2234: for (i=1; i<=imax; i++) aij[i] = lgetg(R,t_VEC);
2235:
2236: c_even = realun(prec);
2237: c_even = gmul2n(c_even,r1);
2238: c_odd = gmul(c_even,gpowgs(racpi,r1));
2239: if (ru&1) c_odd=gneg_i(c_odd);
2240: ck_even=cgetg(R,t_VEC); ck_odd=cgetg(r2+2,t_VEC);
2241: for (k=1; k<R; k++)
2242: {
2243: ck_even[k]=lmul((GEN)zet[k],gadd(gr2,gmul2n(gr1,-k)));
2244: if (k&1) ck_even[k]=lneg((GEN)ck_even[k]);
2245: }
2246: gru=stoi(ru);
2247: for (k=1; k<=r2+1; k++)
2248: {
2249: ck_odd[k]=lmul((GEN)zet[k], gsub(gru, gmul2n(gr1,-k)));
2250: if (k&1) ck_odd[k]=lneg((GEN)ck_odd[k]);
2251: ck_odd[k]=ladd(gru,(GEN)ck_odd[k]);
2252: }
1.2 ! noro 2253: ck_odd[1]=lsub((GEN)ck_odd[1],gmul(gr1,mplog2(prec)));
1.1 noro 2254: serie_even =cgetg(ru+3,t_SER); serie_odd=cgetg(r2+3,t_SER);
2255: serie_even[1] = serie_odd[1] = evalsigne(1)+evalvalp(1);
2256: i=0;
2257:
2258: while (i<imax/2)
2259: {
2260: for (k=1; k<R; k++)
2261: serie_even[k+1]=ldivgs((GEN)ck_even[k],k);
2262: serie_exp=gmul(c_even,gexp(serie_even,0));
2263: p1=(GEN)aij[2*i+1];
2264: for (j=1; j<R; j++) p1[j]=serie_exp[ru+3-j];
2265:
2266: for (k=1; k<=r2+1; k++)
2267: serie_odd[k+1]=ldivgs((GEN)ck_odd[k],k);
2268: serie_exp=gmul(c_odd,gexp(serie_odd,0));
2269: p1=(GEN)aij[2*i+2];
2270: for (j=1; j<=r2+1; j++) p1[j]=serie_exp[r2+3-j];
2271: for ( ; j<R; j++) p1[j]=zero;
2272: i++;
2273:
2274: c_even = gdiv(c_even,gmul(gpowgs(stoi(i),ru),gpowgs(stoi(2*i-1),r2)));
2275: c_odd = gdiv(c_odd, gmul(gpowgs(stoi(i),r2),gpowgs(stoi(2*i+1),ru)));
2276: c_even = gmul2n(c_even,-r2);
2277: c_odd = gmul2n(c_odd,r1-r2);
2278: if (r1&1) { c_even=gneg_i(c_even); c_odd=gneg_i(c_odd); }
2279: p1 = gr2; p2 = gru;
2280: for (k=1; k<R; k++)
2281: {
2282: p1=gdivgs(p1,2*i-1); p2=gdivgs(p2,2*i);
2283: ck_even[k] = ladd((GEN)ck_even[k], gadd(p1,p2));
2284: }
2285: p1 = gru; p2 = gr2;
2286: for (k=1; k<=r2+1; k++)
2287: {
2288: p1=gdivgs(p1,2*i+1); p2=gdivgs(p2,2*i);
2289: ck_odd[k] = ladd((GEN)ck_odd[k], gadd(p1,p2));
2290: }
2291: }
2292: aij=gerepilecopy(av,aij);
2293: if (DEBUGLEVEL>=2) msgtimer("a(i,j)");
2294: p1=cgetg(5,t_VEC);
2295: p1[1]=lstoi(r1); p1[2]=lstoi(r2); p1[3]=lstoi(imax); p1[4]=(long)bnf;
2296: nfz[1]=(long)p1;
2297: nfz[2]=(long)resi;
2298: nfz[5]=(long)cst;
2299: nfz[6]=llog(cst,prec);
2300: nfz[7]=(long)aij;
2301:
2302: /************* Calcul du nombre d'ideaux de norme donnee *************/
2303: coef = dirzetak0(nf,N0); tabj = cgetg(N0+1,t_MAT);
2304: if (DEBUGLEVEL>=2) msgtimer("coef");
2305: colzero=cgetg(ru+2,t_COL); for (j=1; j<=ru+1; j++) colzero[j]=zero;
2306: for (i=1; i<=N0; i++)
2307: if (coef[i])
2308: {
2309: GEN t = cgetg(ru+2,t_COL);
2310: p1 = glog((GEN)tabcstn[i],prec); setsigne(p1, -signe(p1));
2311: t[1] = lstoi(coef[i]); t[2] = lmul((GEN)t[1],p1);
2312: for (j=2; j<=ru; j++)
2313: {
1.2 ! noro 2314: gpmem_t av2 = avma;
! 2315: p2 = gmul((GEN)t[j], p1);
1.1 noro 2316: t[j+1] = (long)gerepileuptoleaf(av2, divrs(p2,j));
2317: }
2318: /* tabj[n,j]=coef(n)*ln(c/n)^(j-1)/(j-1)! */
2319: tabj[i] = (long)t;
2320: }
2321: else tabj[i]=(long)colzero;
2322: if (DEBUGLEVEL>=2) msgtimer("a(n)");
2323:
2324: coeflog=cgetg(N0+1,t_VEC); coeflog[1]=zero;
2325: for (i=2; i<=N0; i++)
2326: if (coef[i])
2327: {
2328: court[2]=i; p1=glog(court,prec);
2329: setsigne(p1,-1); coeflog[i]=(long)p1;
2330: }
2331: else coeflog[i] = zero;
2332: if (DEBUGLEVEL>=2) msgtimer("log(n)");
2333:
2334: nfz[3]=(long)tabj;
2335: p1 = cgetg(N0+1,t_VECSMALL);
2336: for (i=1; i<=N0; i++) p1[i] = coef[i];
2337: nfz[8]=(long)p1;
2338: nfz[9]=(long)coeflog;
2339:
2340: /******************** Calcul des coefficients Cik ********************/
2341:
2342: C = cgetg(ru+1,t_MAT);
2343: for (k=1; k<=ru; k++) C[k] = lgetg(imax+1,t_COL);
2344: av2 = avma;
2345: for (i=1; i<=imax; i++)
2346: {
2347: stackzone *z;
2348: for (k=1; k<=ru; k++)
2349: {
2350: p1 = NULL;
2351: for (n=1; n<=N0; n++)
2352: if (coef[n])
2353: for (j=1; j<=ru-k+1; j++)
2354: {
2355: p2 = gmul(tabcstni[n],
2356: gmul(gmael(aij,i,j+k), gmael(tabj,n,j)));
2357: p1 = p1? gadd(p1,p2): p2;
2358: }
2359: coeff(C,i,k) = p1? (long)gerepileupto(av2,p1): zero;
2360: av2 = avma;
2361: }
2362: /* use a parallel stack */
2363: z = i&1? zone1: zone0;
1.2 ! noro 2364: (void)switch_stack(z, 1);
1.1 noro 2365: for (n=1; n<=N0; n++)
2366: if (coef[n]) tabcstni[n] = mpmul(tabcstni[n],tabcstn[n]);
2367: /* come back */
1.2 ! noro 2368: (void)switch_stack(z, 0);
1.1 noro 2369: }
2370: nfz[4] = (long) C;
2371: if (DEBUGLEVEL>=2) msgtimer("Cik");
2372: gunclone(aij);
2373: free((void*)zone); free((void*)zone1); free((void*)zone0);
2374: free((void*)coef); return nfz;
2375: }
2376:
2377: GEN
2378: gzetakall(GEN nfz, GEN s, long flag, long prec2)
2379: {
2380: GEN resi,C,cst,cstlog,coeflog,cs,coef;
2381: GEN lambd,gammas,gammaunmoins,gammas2,gammaunmoins2,var1,var2;
2382: GEN p1,unmoins,gexpro,gar,val,valm,valk,valkm;
1.2 ! noro 2383: long ts = typ(s), r1,r2,ru,imax,i,j,k,N0,sl,prec,bigprec;
! 2384: gpmem_t av = avma;
1.1 noro 2385:
2386: if (typ(nfz)!=t_VEC || lg(nfz)!=10 || typ(nfz[1]) != t_VEC)
2387: err(talker,"not a zeta number field in zetakall");
2388: if (! is_intreal_t(ts) && ts != t_COMPLEX && ! is_frac_t(ts))
2389: err(typeer,"gzetakall");
2390: resi=(GEN)nfz[2]; C=(GEN)nfz[4]; cst=(GEN)nfz[5];
2391: cstlog=(GEN)nfz[6]; coef=(GEN)nfz[8]; coeflog=(GEN)nfz[9];
2392: r1 = itos(gmael(nfz,1,1));
2393: r2 = itos(gmael(nfz,1,2)); ru = r1+r2;
2394: imax = itos(gmael(nfz,1,3)); N0 = lg(coef)-1;
2395: bigprec = precision(cst); prec = prec2+1;
2396:
2397: if (ts==t_COMPLEX && gcmp0(gimag(s))) { s=greal(s); ts = typ(s); }
2398: if (ts==t_REAL && !signe(gfrac(s))) { s=mptrunc(s); ts = t_INT; }
2399: if (ts==t_INT)
2400: {
2401: sl=itos(s);
2402: if (sl==1) err(talker,"s = 1 is a pole (gzetakall)");
2403: if (sl==0)
2404: {
2405: avma = av;
2406: if (flag) err(talker,"s = 0 is a pole (gzetakall)");
2407: if (ru == 1) return gneg(r1? ghalf: resi);
2408: return gzero;
2409: }
2410: if (sl<0 && (r2 || !odd(sl)))
2411: {
2412: if (!flag) return gzero;
2413: s = subsi(1,s); sl = 1-sl;
2414: }
2415: unmoins=subsi(1,s);
2416: lambd = gdiv(resi, mulis(s,sl-1));
2417: gammas2=ggamma(gmul2n(s,-1),prec);
2418: gar=gpowgs(gammas2,r1);
2419: cs=gexp(gmul(cstlog,s),prec);
2420: val=s; valm=unmoins;
2421: if (sl < 0) /* r2 = 0 && odd(sl) */
2422: {
2423: gammaunmoins2=ggamma(gmul2n(unmoins,-1),prec);
2424: var1=var2=gun;
2425: for (i=2; i<=N0; i++)
2426: if (coef[i])
2427: {
2428: gexpro=gexp(gmul((GEN)coeflog[i],s),bigprec);
2429: var1 = gadd(var1,gmulsg(coef[i],gexpro));
2430: var2 = gadd(var2,gdivsg(coef[i],gmulsg(i,gexpro)));
2431: }
2432: lambd=gadd(lambd,gmul(gmul(var1,cs),gar));
2433: lambd=gadd(lambd,gmul(gmul(var2,gdiv(cst,cs)),
2434: gpowgs(gammaunmoins2,r1)));
2435: for (i=1; i<=imax; i+=2)
2436: {
2437: valk = val;
2438: valkm = valm;
2439: for (k=1; k<=ru; k++)
2440: {
2441: p1 = gcoeff(C,i,k);
2442: lambd = gsub(lambd,gdiv(p1,valk )); valk = mulii(val,valk);
2443: lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm);
2444: }
2445: val = addis(val, 2);
2446: valm = addis(valm,2);
2447: }
2448: }
2449: else
2450: {
2451: GEN tabj=(GEN)nfz[3], aij=(GEN)nfz[7];
2452:
2453: gar = gmul(gar,gpowgs(ggamma(s,prec),r2));
2454: var1=var2=gzero;
2455: for (i=1; i<=N0; i++)
2456: if (coef[i])
2457: {
2458: gexpro=gexp(gmul((GEN)coeflog[i],s),bigprec);
2459: var1 = gadd(var1,gmulsg(coef[i],gexpro));
2460: if (sl <= imax)
2461: {
2462: p1=gzero;
2463: for (j=1; j<=ru+1; j++)
2464: p1 = gadd(p1, gmul(gmael(aij,sl,j), gmael(tabj,i,j)));
2465: var2=gadd(var2,gdiv(p1,gmulsg(i,gexpro)));
2466: }
2467: }
2468: lambd=gadd(lambd,gmul(gmul(var1,cs),gar));
2469: lambd=gadd(lambd,gmul(var2,gdiv(cst,cs)));
2470: for (i=1; i<=imax; i++)
2471: {
2472: valk = val;
2473: valkm = valm;
2474: if (i == sl)
2475: for (k=1; k<=ru; k++)
2476: {
2477: p1 = gcoeff(C,i,k);
2478: lambd = gsub(lambd,gdiv(p1,valk)); valk = mulii(val,valk);
2479: }
2480: else
2481: for (k=1; k<=ru; k++)
2482: {
2483: p1 = gcoeff(C,i,k);
1.2 ! noro 2484: lambd = gsub(lambd,gdiv(p1,valk )); valk = mulii(val,valk);
1.1 noro 2485: lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm);
2486: }
2487: val = addis(val,1);
2488: valm = addis(valm,1);
2489: }
2490: }
2491: }
2492: else
2493: {
2494: GEN Pi = mppi(prec);
2495: if (is_frac_t(ts))
2496: s = gmul(s, realun(bigprec));
2497: else
2498: s = gprec_w(s, bigprec);
2499:
2500: unmoins = gsub(gun,s);
2501: lambd = gdiv(resi,gmul(s,gsub(s,gun)));
2502: gammas = ggamma(s,prec);
2503: gammas2= ggamma(gmul2n(s,-1),prec);
2504: gar = gmul(gpowgs(gammas,r2),gpowgs(gammas2,r1));
2505: cs = gexp(gmul(cstlog,s),prec);
2506: var1 = gmul(Pi,s);
2507: gammaunmoins = gdiv(Pi,gmul(gsin(var1,prec),gammas));
2508: gammaunmoins2= gdiv(gmul(gmul(gsqrt(Pi,prec),gpui(gdeux,gsub(s,gun),prec)),
2509: gammas2),
2510: gmul(gcos(gmul2n(var1,-1),prec),gammas));
2511: var1 = var2 = gun;
2512: for (i=2; i<=N0; i++)
2513: if (coef[i])
2514: {
2515: gexpro = gexp(gmul((GEN)coeflog[i],s),bigprec);
2516: var1 = gadd(var1,gmulsg(coef[i], gexpro));
2517: var2 = gadd(var2,gdivsg(coef[i], gmulsg(i,gexpro)));
2518: }
2519: lambd = gadd(lambd,gmul(gmul(var1,cs),gar));
2520: lambd = gadd(lambd,gmul(gmul(gmul(var2,gdiv(cst,cs)),
2521: gpowgs(gammaunmoins,r2)),
2522: gpowgs(gammaunmoins2,r1)));
2523: val = s;
2524: valm = unmoins;
2525: for (i=1; i<=imax; i++)
2526: {
2527: valk = val;
2528: valkm = valm;
2529: for (k=1; k<=ru; k++)
2530: {
2531: p1 = gcoeff(C,i,k);
2532: lambd = gsub(lambd,gdiv(p1,valk )); valk = gmul(val, valk);
2533: lambd = gsub(lambd,gdiv(p1,valkm)); valkm = gmul(valm,valkm);
2534: }
2535: if (r2)
2536: {
2537: val = gadd(val, gun);
2538: valm = gadd(valm,gun);
2539: }
2540: else
2541: {
2542: val = gadd(val, gdeux);
2543: valm = gadd(valm,gdeux); i++;
2544: }
2545: }
2546: }
2547: if (!flag) lambd = gdiv(lambd,gmul(gar,cs)); /* zetak */
2548: return gerepileupto(av, gprec_w(lambd, prec2));
2549: }
2550:
2551: GEN
2552: gzetak(GEN nfz, GEN s, long prec)
2553: {
2554: return gzetakall(nfz,s,0,prec);
2555: }
2556:
2557: GEN
2558: glambdak(GEN nfz, GEN s, long prec)
2559: {
2560: return gzetakall(nfz,s,1,prec);
2561: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to base1.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