Annotation of OpenXM_contrib/pari-2.2/src/basemath/base1.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: base1.c,v 1.41 2001/10/01 12:11:28 karim Exp $
2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /**************************************************************/
17: /* */
18: /* NUMBER FIELDS */
19: /* */
20: /**************************************************************/
21: #include "pari.h"
22: #include "parinf.h"
23: GEN idealhermite_aux(GEN nf, GEN x);
24:
25: void
26: checkrnf(GEN rnf)
27: {
28: if (typ(rnf)!=t_VEC) err(idealer1);
29: if (lg(rnf)!=12) err(idealer1);
30: }
31:
32: GEN
33: checkbnf(GEN bnf)
34: {
35: if (typ(bnf)!=t_VEC) err(idealer1);
36: switch (lg(bnf))
37: {
38: case 11: return bnf;
39: case 10:
40: if (typ(bnf[1])==t_POL)
41: err(talker,"please apply bnfinit first");
42: break;
43: case 7:
44: return checkbnf((GEN)bnf[1]);
45:
46: case 3:
47: if (typ(bnf[2])==t_POLMOD)
48: return checkbnf((GEN)bnf[1]);
49: }
50: err(idealer1);
51: return NULL; /* not reached */
52: }
53:
54: GEN
55: checkbnf_discard(GEN bnf)
56: {
57: GEN x = checkbnf(bnf);
58: if (x != bnf && lg(bnf) == 3)
59: err(warner,"non-monic polynomial. Change of variables discarded");
60: return x;
61: }
62:
63: GEN
64: checknf(GEN nf)
65: {
66: if (typ(nf)==t_POL) err(talker,"please apply nfinit first");
67: if (typ(nf)!=t_VEC) err(idealer1);
68: switch(lg(nf))
69: {
70: case 10: return nf;
71: case 11: return checknf((GEN)nf[7]);
72: case 7: return checknf((GEN)nf[1]);
73: case 3: if (typ(nf[2]) == t_POLMOD) return checknf((GEN)nf[1]);
74: }
75: err(idealer1);
76: return NULL; /* not reached */
77: }
78:
79: void
80: checkbnr(GEN bnr)
81: {
82: if (typ(bnr)!=t_VEC || lg(bnr)!=7)
83: err(talker,"incorrect bigray field");
84: (void)checkbnf((GEN)bnr[1]);
85: }
86:
87: void
88: checkbnrgen(GEN bnr)
89: {
90: checkbnr(bnr);
91: if (lg(bnr[5])<=3)
92: err(talker,"please apply bnrinit(,,1) and not bnrinit(,)");
93: }
94:
95: GEN
96: check_units(GEN bnf, char *f)
97: {
98: GEN nf, x;
99: bnf = checkbnf(bnf); nf = checknf(bnf);
100: x = (GEN)bnf[8];
101: if (lg(x) < 7 || (lg(x[5]) == 1 && lg(nf[6]) > 2))
102: err(talker,"missing units in %s", f);
103: return (GEN)x[5];
104: }
105:
106: void
107: checkid(GEN x, long N)
108: {
109: if (typ(x)!=t_MAT) err(idealer2);
110: if (lg(x) == 1 || lg(x[1]) != N+1)
111: err(talker,"incorrect matrix for ideal");
112: }
113:
114: void
115: checkbid(GEN bid)
116: {
117: if (typ(bid)!=t_VEC || lg(bid)!=6 || lg(bid[1])!=3)
118: err(talker,"incorrect bigideal");
119: }
120:
121: void
122: checkprimeid(GEN id)
123: {
124: if (typ(id) != t_VEC || lg(id) != 6)
125: err(talker,"incorrect prime ideal");
126: }
127:
128: void
129: checkprhall(GEN prhall)
130: {
131: if (typ(prhall) != t_VEC || lg(prhall) != 3 || typ(prhall[1]) != t_MAT)
132: err(talker,"incorrect prhall format");
133: }
134:
135: void
136: check_pol_int(GEN x, char *s)
137: {
138: long k = lgef(x)-1;
139: for ( ; k>1; k--)
140: if (typ(x[k])!=t_INT) err(talker,"polynomial not in Z[X] in %s",s);
141: }
142:
143: GEN
144: checknfelt_mod(GEN nf, GEN x, char *s)
145: {
146: if (!gegal((GEN)x[1],(GEN)nf[1]))
147: err(talker,"incompatible modulus in %s:\n mod = %Z,\n nf = %Z",
148: s, x[1], nf[1]);
149: return (GEN)x[2];
150: }
151:
152: GEN
153: get_primeid(GEN x)
154: {
155: if (typ(x) != t_VEC) return NULL;
156: if (lg(x) == 3) x = (GEN)x[1];
157: if (lg(x) != 6 || typ(x[3]) != t_INT) return NULL;
158: return x;
159: }
160:
161: GEN
162: get_bnf(GEN x, int *t)
163: {
164: switch(typ(x))
165: {
166: case t_POL: *t = typ_POL; return NULL;
167: case t_QUAD: *t = typ_Q ; return NULL;
168: case t_VEC:
169: switch(lg(x))
170: {
171: case 3:
172: if (typ(x[2]) != t_POLMOD) break;
173: return get_bnf((GEN)x[1],t);
174: case 6 : *t = typ_QUA; return NULL;
175: case 10: *t = typ_NF; return NULL;
176: case 11: *t = typ_BNF; return x;
177: case 7 : *t = typ_BNR;
178: x = (GEN)x[1]; if (typ(x)!=t_VEC || lg(x)!=11) break;
179: return x;
180: }
181: case t_MAT:
182: if (lg(x)==2)
183: switch(lg(x[1]))
184: {
185: case 8: case 11:
186: *t = typ_CLA; return NULL;
187: }
188: }
189: *t = typ_NULL; return NULL;
190: }
191:
192: GEN
193: get_nf(GEN x, int *t)
194: {
195: switch(typ(x))
196: {
197: case t_POL : *t = typ_POL; return NULL;
198: case t_QUAD: *t = typ_Q ; return NULL;
199: case t_VEC:
200: switch(lg(x))
201: {
202: case 3:
203: if (typ(x[2]) != t_POLMOD) break;
204: return get_nf((GEN)x[1],t);
205: case 10: *t = typ_NF; return x;
206: case 11: *t = typ_BNF;
207: x = (GEN)x[7]; if (typ(x)!=t_VEC || lg(x)!=10) break;
208: return x;
209: case 7 : *t = typ_BNR;
210: x = (GEN)x[1]; if (typ(x)!=t_VEC || lg(x)!=11) break;
211: x = (GEN)x[7]; if (typ(x)!=t_VEC || lg(x)!=10) break;
212: return x;
213: case 9 :
214: x = (GEN)x[2];
215: if (typ(x) == t_VEC && lg(x) == 4) *t = typ_GAL;
216: return NULL;
217: case 14: case 20:
218: *t = typ_ELL; return NULL;
219: }break;
220: case t_MAT:
221: if (lg(x)==2)
222: switch(lg(x[1]))
223: {
224: case 8: case 11:
225: *t = typ_CLA; return NULL;
226: }
227: }
228: *t = typ_NULL; return NULL;
229: }
230:
231: /*************************************************************************/
232: /** **/
233: /** GALOIS GROUP **/
234: /** **/
235: /*************************************************************************/
236:
237: /* exchange elements i and j in vector x */
238: static GEN
239: transroot(GEN x, int i, int j)
240: {
241: long k;
242: x = dummycopy(x);
243: k=x[i]; x[i]=x[j]; x[j]=k; return x;
244: }
245:
246: #define randshift (BITS_IN_RANDOM - 3)
247:
248: GEN
249: tschirnhaus(GEN x)
250: {
251: long a, v = varn(x), av = avma;
252: GEN u, p1 = cgetg(5,t_POL);
253:
254: if (typ(x)!=t_POL) err(notpoler,"tschirnhaus");
255: if (lgef(x) < 4) err(constpoler,"tschirnhaus");
256: if (v) { u=dummycopy(x); setvarn(u,0); x=u; }
257: p1[1] = evalsigne(1)|evalvarn(0)|evallgef(5);
258: do
259: {
260: a = mymyrand() >> randshift; if (a==0) a =1; p1[4]=lstoi(a);
261: a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[3]=lstoi(a);
262: a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[2]=lstoi(a);
263: u = caract2(x,p1,v); a=avma;
264: }
265: while (lgef(srgcd(u,derivpol(u))) > 3); /* while u not separable */
266: if (DEBUGLEVEL>1)
267: fprintferr("Tschirnhaus transform. New pol: %Z",u);
268: avma=a; return gerepileupto(av,u);
269: }
270: #undef randshift
271:
272: int
273: gpolcomp(GEN p1, GEN p2)
274: {
275: int s,j = lgef(p1)-2;
276:
277: if (lgef(p2)-2 != j)
278: err(bugparier,"gpolcomp (different degrees)");
279: for (; j>=2; j--)
280: {
281: s = absi_cmp((GEN)p1[j], (GEN)p2[j]);
282: if (s) return s;
283: }
284: return 0;
285: }
286:
287: /* assume pol in Z[X] */
288: GEN
289: primitive_pol_to_monic(GEN pol, GEN *ptlead)
290: {
291: long i,j,n = degpol(pol);
292: GEN lead,fa,e,a, POL = dummycopy(pol);
293:
294: a = POL + 2; lead = (GEN)a[n];
295: if (signe(lead) < 0) { POL = gneg_i(POL); a = POL+2; lead = negi(lead); }
296: if (is_pm1(lead)) { if (ptlead) *ptlead = NULL; return POL; }
297: fa = auxdecomp(lead,0); lead = gun;
298: e = (GEN)fa[2]; fa = (GEN)fa[1];
299: for (i=lg(e)-1; i>0;i--) e[i] = itos((GEN)e[i]);
300: for (i=lg(fa)-1; i>0; i--)
301: {
302: GEN p = (GEN)fa[i], p1,pk,pku;
303: long k = (long)ceil((double)e[i] / n);
304: long d = k * n - e[i], v, j0;
305: /* find d, k such that p^d pol(x / p^k) monic */
306: for (j=n-1; j>0; j--)
307: {
308: if (!signe(a[j])) continue;
309: v = pvaluation((GEN)a[j], p, &p1);
310: while (v + d < k * j) { k++; d += n; }
311: }
312: pk = gpowgs(p,k); j0 = d/k;
313: pku = gpowgs(p,d - k*j0);
314: /* a[j] *= p^(d - kj) */
315: for (j=j0; j>=0; j--)
316: {
317: if (j < j0) pku = mulii(pku, pk);
318: a[j] = lmulii((GEN)a[j], pku);
319: }
320: j0++;
321: pku = gpowgs(p,k*j0 - d);
322: /* a[j] /= p^(kj - d) */
323: for (j=j0; j<=n; j++)
324: {
325: if (j > j0) pku = mulii(pku, pk);
326: a[j] = ldivii((GEN)a[j], pku);
327: }
328: lead = mulii(lead, pk);
329: }
330: if (ptlead) *ptlead = lead; return POL;
331: }
332:
333: /* compute x1*x2^2 + x2*x3^2 + x3*x4^2 + x4*x1^2 */
334: static GEN
335: get_F4(GEN x)
336: {
337: GEN p1=gzero;
338: long i;
339:
340: for (i=1; i<=4; i++)
341: p1 = gadd(p1, gmul((GEN)x[i], gsqr((GEN)x[(i&3)+1])));
342: return p1;
343: }
344:
345: GEN galoisbig(GEN x, long prec);
346: GEN roots_to_pol(GEN a, long v);
347:
348: GEN
349: galois(GEN x, long prec)
350: {
351: long av=avma,av1,i,j,k,n,f,l,l2,e,e1,pr,ind;
352: GEN x1,p1,p2,p3,p4,p5,p6,w,y,z,ee;
353: static int ind5[20]={2,5,3,4, 1,3,4,5, 1,5,2,4, 1,2,3,5, 1,4,2,3};
354: static int ind6[60]={3,5,4,6, 2,6,4,5, 2,3,5,6, 2,4,3,6, 2,5,3,4,
355: 1,4,5,6, 1,5,3,6, 1,6,3,4, 1,3,4,5, 1,6,2,5,
356: 1,2,4,6, 1,5,2,4, 1,3,2,6, 1,2,3,5, 1,4,2,3};
357: if (typ(x)!=t_POL) err(notpoler,"galois");
358: n=degpol(x); if (n<=0) err(constpoler,"galois");
359: if (n>11) err(impl,"galois of degree higher than 11");
360: x = gdiv(x,content(x));
361: check_pol_int(x, "galois");
362: if (gisirreducible(x) != gun)
363: err(impl,"galois of reducible polynomial");
364:
365: if (n<4)
366: {
367: if (n<3)
368: {
369: avma=av; y=cgetg(4,t_VEC);
370: y[1] = (n==1)? un: deux;
371: y[2]=lnegi(gun);
372: }
373: else /* n=3 */
374: {
375: f=carreparfait(discsr(x));
376: avma=av; y=cgetg(4,t_VEC);
377: if (f) { y[1]=lstoi(3); y[2]=un; }
378: else { y[1]=lstoi(6); y[2]=lnegi(gun); }
379: }
380: y[3]=un; return y;
381: }
382: x1 = x = primitive_pol_to_monic(x,NULL); av1=avma;
383: if (n>7) return galoisbig(x,prec);
384: for(;;)
385: {
386: switch(n)
387: {
388: case 4: z = cgetg(7,t_VEC);
389: for(;;)
390: {
391: p1=roots(x,prec);
392: z[1] = (long)get_F4(p1);
393: z[2] = (long)get_F4(transroot(p1,1,2));
394: z[3] = (long)get_F4(transroot(p1,1,3));
395: z[4] = (long)get_F4(transroot(p1,1,4));
396: z[5] = (long)get_F4(transroot(p1,2,3));
397: z[6] = (long)get_F4(transroot(p1,3,4));
398: p4 = roots_to_pol(z,0);
399: p5=grndtoi(greal(p4),&e);
400: e1=gexpo(gimag(p4)); if (e1>e) e=e1;
401: if (e <= -10) break;
402: prec = (prec<<1)-2;
403: }
404: p6 = modulargcd(p5, derivpol(p5));
405: if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
406: p2 = (GEN)factor(p5)[1];
407: switch(lg(p2)-1)
408: {
409: case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC);
410: y[3]=un;
411: if (f) { y[2]=un; y[1]=lstoi(12); return y; }
412: y[2]=lnegi(gun); y[1]=lstoi(24); return y;
413:
414: case 2: avma=av; y=cgetg(4,t_VEC);
415: y[3]=un; y[2]=lnegi(gun); y[1]=lstoi(8); return y;
416:
417: case 3: avma=av; y=cgetg(4,t_VEC);
418: y[1]=lstoi(4); y[3]=un;
419: y[2] = (lgef(p2[1])==5)? un: lnegi(gun);
420: return y;
421:
422: default: err(bugparier,"galois (bug1)");
423: }
424:
425: case 5: z = cgetg(7,t_VEC);
426: ee = new_chunk(7); w = new_chunk(7);
427: for(;;)
428: {
429: for(;;)
430: {
431: p1=roots(x,prec);
432: for (l=1; l<=6; l++)
433: {
434: p2=(l==1)?p1: ((l<6)?transroot(p1,1,l): transroot(p1,2,5));
435: p3=gzero;
436: for (k=0,i=1; i<=5; i++,k+=4)
437: {
438: p5 = gadd(gmul((GEN)p2[ind5[k]],(GEN)p2[ind5[k+1]]),
439: gmul((GEN)p2[ind5[k+2]],(GEN)p2[ind5[k+3]]));
440: p3 = gadd(p3, gmul(gsqr((GEN)p2[i]),p5));
441: }
442: w[l]=lrndtoi(greal(p3),&e);
443: e1 = gexpo(gimag(p3)); if (e1>e) e=e1;
444: ee[l]=e; z[l] = (long)p3;
445: }
446: p4 = roots_to_pol(z,0);
447: p5=grndtoi(greal(p4),&e);
448: e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
449: if (e <= -10) break;
450: prec = (prec<<1)-2;
451: }
452: p6 = modulargcd(p5,derivpol(p5));
453: if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
454: p3=(GEN)factor(p5)[1];
455: f=carreparfait(discsr(x));
456: if (lg(p3)-1==1)
457: {
458: avma=av; y=cgetg(4,t_VEC); y[3]=un;
459: if (f) { y[2]=un; y[1]=lstoi(60); return y; }
460: else { y[2]=lneg(gun); y[1]=lstoi(120); return y; }
461: }
462: if (!f)
463: {
464: avma=av; y=cgetg(4,t_VEC);
465: y[3]=un; y[2]=lneg(gun); y[1]=lstoi(20); return y;
466: }
467: pr = - (bit_accuracy(prec) >> 1);
468: for (l=1; l<=6; l++)
469: if (ee[l] <= pr && gcmp0(poleval(p5,(GEN)w[l]))) break;
470: if (l>6) err(bugparier,"galois (bug4)");
471: p2=(l==6)? transroot(p1,2,5):transroot(p1,1,l);
472: p3=gzero;
473: for (i=1; i<=5; i++)
474: {
475: j=(i%5)+1;
476: p3=gadd(p3,gmul(gmul((GEN)p2[i],(GEN)p2[j]),
477: gsub((GEN)p2[j],(GEN)p2[i])));
478: }
479: p5=gsqr(p3); p4=grndtoi(greal(p5),&e);
480: e1 = gexpo(gimag(p5)); if (e1>e) e=e1;
481: if (e <= -10)
482: {
483: if (gcmp0(p4)) goto tchi;
484: f=carreparfait(p4); avma=av; y=cgetg(4,t_VEC);
485: y[3]=y[2]=un; y[1]=lstoi(f?5:10);
486: return y;
487: }
488: prec=(prec<<1)-2;
489: }
490:
491: case 6: z = cgetg(7, t_VEC);
492: for(;;)
493: {
494: for(;;)
495: {
496: p1=roots(x,prec);
497: for (l=1; l<=6; l++)
498: {
499: p2=(l==1)?p1:transroot(p1,1,l);
500: p3=gzero; k=0;
501: for (i=1; i<=5; i++) for (j=i+1; j<=6; j++)
502: {
503: p5=gadd(gmul((GEN)p2[ind6[k]],(GEN)p2[ind6[k+1]]),
504: gmul((GEN)p2[ind6[k+2]],(GEN)p2[ind6[k+3]]));
505: p3=gadd(p3,gmul(gsqr(gmul((GEN)p2[i],(GEN)p2[j])),p5)); k+=4;
506: }
507: z[l] = (long)p3;
508: }
509: p4=roots_to_pol(z,0);
510: p5=grndtoi(greal(p4),&e);
511: e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
512: if (e <= -10) break;
513: prec=(prec<<1)-2;
514: }
515: p6 = modulargcd(p5,derivpol(p5));
516: if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
517: p2=(GEN)factor(p5)[1];
518: switch(lg(p2)-1)
519: {
520: case 1:
521: z = cgetg(11,t_VEC); ind=0;
522: p3=gadd(gmul(gmul((GEN)p1[1],(GEN)p1[2]),(GEN)p1[3]),
523: gmul(gmul((GEN)p1[4],(GEN)p1[5]),(GEN)p1[6]));
524: z[++ind] = (long)p3;
525: for (i=1; i<=3; i++)
526: for (j=4; j<=6; j++)
527: {
528: p2=transroot(p1,i,j);
529: p3=gadd(gmul(gmul((GEN)p2[1],(GEN)p2[2]),(GEN)p2[3]),
530: gmul(gmul((GEN)p2[4],(GEN)p2[5]),(GEN)p2[6]));
531: z[++ind] = (long)p3;
532: }
533: p4 = roots_to_pol(z,0);
534: p5=grndtoi(greal(p4),&e);
535: e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
536: if (e <= -10)
537: {
538: p6 = modulargcd(p5,derivpol(p5));
539: if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
540: p2=(GEN)factor(p5)[1];
541: f=carreparfait(discsr(x));
542: avma=av; y=cgetg(4,t_VEC); y[3]=un;
543: if (lg(p2)-1==1)
544: {
545: if (f) { y[2]=un; y[1]=lstoi(360); }
546: else { y[2]=lnegi(gun); y[1]=lstoi(720); }
547: }
548: else
549: {
550: if (f) { y[2]=un; y[1]=lstoi(36); }
551: else { y[2]=lnegi(gun); y[1]=lstoi(72); }
552: }
553: return y;
554: }
555: prec=(prec<<1)-2; break;
556:
557: case 2: l2=degpol(p2[1]); if (l2>3) l2=6-l2;
558: switch(l2)
559: {
560: case 1: f=carreparfait(discsr(x));
561: avma=av; y=cgetg(4,t_VEC); y[3]=un;
562: if (f) { y[2]=un; y[1]=lstoi(60); }
563: else { y[2]=lneg(gun); y[1]=lstoi(120); }
564: return y;
565: case 2: f=carreparfait(discsr(x));
566: if (f)
567: {
568: avma=av; y=cgetg(4,t_VEC);
569: y[3]=y[2]=un; y[1]=lstoi(24);
570: }
571: else
572: {
573: p3=(lgef(p2[1])==5) ? (GEN)p2[2]:(GEN)p2[1];
574: f=carreparfait(discsr(p3));
575: avma=av; y=cgetg(4,t_VEC); y[2]=lneg(gun);
576: if (f) { y[1]=lstoi(24); y[3]=deux; }
577: else { y[1]=lstoi(48); y[3]=un; }
578: }
579: return y;
580: case 3: f=carreparfait(discsr((GEN)p2[1]))
581: || carreparfait(discsr((GEN)p2[2]));
582: avma=av; y=cgetg(4,t_VEC);
583: y[3]=un; y[2]=lneg(gun); y[1]=lstoi(f? 18: 36);
584: return y;
585: }
586: case 3:
587: for (l2=1; l2<=3; l2++)
588: if (lgef(p2[l2])>=6) p3=(GEN)p2[l2];
589: if (lgef(p3)==6)
590: {
591: f=carreparfait(discsr(p3)); avma=av; y=cgetg(4,t_VEC);
592: y[2]=lneg(gun); y[1]=lstoi(f? 6: 12);
593: }
594: else
595: {
596: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC);
597: if (f) { y[2]=un; y[1]=lstoi(12); }
598: else { y[2]=lneg(gun); y[1]=lstoi(24); }
599: }
600: y[3]=un; return y;
601: case 4: avma=av; y=cgetg(4,t_VEC);
602: y[1]=lstoi(6); y[2]=lneg(gun); y[3]=deux; return y;
603: default: err(bugparier,"galois (bug3)");
604: }
605: }
606:
607: case 7: z = cgetg(36,t_VEC);
608: for(;;)
609: {
610: ind = 0; p1=roots(x,prec); p4=gun;
611: for (i=1; i<=5; i++)
612: for (j=i+1; j<=6; j++)
613: {
614: p6 = gadd((GEN)p1[i],(GEN)p1[j]);
615: for (k=j+1; k<=7; k++)
616: z[++ind] = (long) gadd(p6,(GEN)p1[k]);
617: }
618: p4 = roots_to_pol(z,0);
619: p5 = grndtoi(greal(p4),&e);
620: e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
621: if (e <= -10) break;
622: prec = (prec<<1)-2;
623: }
624: p6 = modulargcd(p5,derivpol(p5));
625: if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
626: p2=(GEN)factor(p5)[1];
627: switch(lg(p2)-1)
628: {
629: case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC); y[3]=un;
630: if (f) { y[2]=un; y[1]=lstoi(2520); }
631: else { y[2]=lneg(gun); y[1]=lstoi(5040); }
632: return y;
633: case 2: f=degpol(p2[1]); avma=av; y=cgetg(4,t_VEC); y[3]=un;
634: if (f==7 || f==28) { y[2]=un; y[1]=lstoi(168); }
635: else { y[2]=lneg(gun); y[1]=lstoi(42); }
636: return y;
637: case 3: avma=av; y=cgetg(4,t_VEC);
638: y[3]=y[2]=un; y[1]=lstoi(21); return y;
639: case 4: avma=av; y=cgetg(4,t_VEC);
640: y[3]=un; y[2]=lneg(gun); y[1]=lstoi(14); return y;
641: case 5: avma=av; y=cgetg(4,t_VEC);
642: y[3]=y[2]=un; y[1]=lstoi(7); return y;
643: default: err(talker,"galois (bug2)");
644: }
645: }
646: tchi: avma=av1; x=tschirnhaus(x1);
647: }
648: }
649:
650: GEN
651: galoisapply(GEN nf, GEN aut, GEN x)
652: {
653: long av=avma,tetpil,lx,j,N;
654: GEN p,p1,y,pol;
655:
656: nf=checknf(nf); pol=(GEN)nf[1];
657: if (typ(aut)==t_POL) aut = gmodulcp(aut,pol);
658: else
659: {
660: if (typ(aut)!=t_POLMOD || !gegal((GEN)aut[1],pol))
661: err(talker,"incorrect galois automorphism in galoisapply");
662: }
663: switch(typ(x))
664: {
665: case t_INT: case t_INTMOD: case t_FRAC: case t_FRACN: case t_PADIC:
666: avma=av; return gcopy(x);
667:
668: case t_POLMOD: x = (GEN) x[2]; /* fall through */
669: case t_POL:
670: p1 = gsubst(x,varn(pol),aut);
671: if (typ(p1)!=t_POLMOD || !gegal((GEN)p1[1],pol)) p1 = gmodulcp(p1,pol);
672: return gerepileupto(av,p1);
673:
674: case t_VEC:
675: if (lg(x)==3)
676: {
677: y=cgetg(3,t_VEC);
678: y[1]=(long)galoisapply(nf,aut,(GEN)x[1]);
679: y[2]=lcopy((GEN)x[2]); return gerepileupto(av, y);
680: }
681: if (lg(x)!=6) err(typeer,"galoisapply");
682: y=cgetg(6,t_VEC); y[1]=x[1]; y[3]=x[3]; y[4]=x[4];
683: p = (GEN)x[1];
684: p1=centermod(galoisapply(nf,aut,(GEN)x[2]), p);
685: if (is_pm1(x[3]))
686: if (ggval(subres(gmul((GEN)nf[7],p1),pol), p) > itos((GEN)x[4]))
687: p1[1] = (signe(p1[1]) > 0)? lsub((GEN)p1[1], p)
688: : ladd((GEN)p1[1], p);
689: y[2]=(long)p1;
690: y[5]=(long)centermod(galoisapply(nf,aut,(GEN)x[5]), p);
691: return gerepilecopy(av,y);
692:
693: case t_COL:
694: N=degpol(pol);
695: if (lg(x)!=N+1) err(typeer,"galoisapply");
696: p1=galoisapply(nf,aut,gmul((GEN)nf[7],x)); tetpil=avma;
697: return gerepile(av,tetpil,algtobasis(nf,p1));
698:
699: case t_MAT:
700: lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
701: N=degpol(pol);
702: if (lg(x[1])!=N+1) err(typeer,"galoisapply");
703: p1=cgetg(lx,t_MAT);
704: for (j=1; j<lx; j++) p1[j]=(long)galoisapply(nf,aut,(GEN)x[j]);
705: if (lx==N+1) p1 = idealhermite(nf,p1);
706: return gerepileupto(av,p1);
707: }
708: err(typeer,"galoisapply");
709: return NULL; /* not reached */
710: }
711:
712: GEN pol_to_monic(GEN pol, GEN *lead);
713: int cmp_pol(GEN x, GEN y);
714:
715: GEN
716: get_nfpol(GEN x, GEN *nf)
717: {
718: if (typ(x) == t_POL) { *nf = NULL; return x; }
719: *nf = checknf(x); return (GEN)(*nf)[1];
720: }
721:
722: /* if fliso test for isomorphism, for inclusion otherwise. */
723: static GEN
724: nfiso0(GEN a, GEN b, long fliso)
725: {
726: ulong av = avma;
727: long n,m,i,vb,lx;
728: GEN nfa,nfb,p1,y,la,lb;
729:
730: a = get_nfpol(a, &nfa);
731: b = get_nfpol(b, &nfb);
732: if (fliso && nfa && !nfb) { p1=a; a=b; b=p1; p1=nfa; nfa=nfb; nfb=p1; }
733: m=degpol(a);
734: n=degpol(b); if (m<=0 || n<=0) err(constpoler,"nfiso or nfincl");
735: if (fliso)
736: { if (n!=m) return gzero; }
737: else
738: { if (n%m) return gzero; }
739:
740: if (nfb) lb = NULL; else b = pol_to_monic(b,&lb);
741: if (nfa) la = NULL; else a = pol_to_monic(a,&la);
742: if (nfa && nfb)
743: {
744: if (fliso)
745: {
746: if (!gegal((GEN)nfa[2],(GEN)nfb[2])
747: || !gegal((GEN)nfa[3],(GEN)nfb[3])) return gzero;
748: }
749: else
750: if (!divise((GEN)nfb[3], gpowgs((GEN)nfa[3],n/m))) return gzero;
751: }
752: else
753: {
754: GEN da = nfa? (GEN)nfa[3]: discsr(a);
755: GEN db = nfb? (GEN)nfb[3]: discsr(b);
756: if (fliso)
757: {
758: p1=gdiv(da,db);
759: if (typ(p1)==t_FRAC) p1=mulii((GEN)p1[1],(GEN)p1[2]);
760: if (!gcarreparfait(p1)) { avma=av; return gzero; }
761: }
762: else
763: {
764: long q=n/m;
765: GEN fa=factor(da), ex=(GEN)fa[2];
766: fa=(GEN)fa[1]; lx=lg(fa);
767: for (i=1; i<lx; i++)
768: if (mod2((GEN)ex[i]) && !divise(db,gpowgs((GEN)fa[i],q)))
769: { avma=av; return gzero; }
770: }
771: }
772: a = dummycopy(a); setvarn(a,0);
773: b = dummycopy(b); vb=varn(b);
774: if (nfb)
775: {
776: if (vb == 0) nfb = gsubst(nfb, 0, polx[MAXVARN]);
777: y = lift_intern(nfroots(nfb,a));
778: }
779: else
780: {
781: if (vb == 0) setvarn(b, fetch_var());
782: y = (GEN)polfnf(a,b)[1]; lx = lg(y);
783: for (i=1; i<lx; i++)
784: {
785: if (lgef(y[i]) != 4) { setlg(y,i); break; }
786: y[i] = (long)gneg_i(lift_intern(gmael(y,i,2)));
787: }
788: y = gen_sort(y, 0, cmp_pol);
789: settyp(y, t_VEC);
790: if (vb == 0) delete_var();
791: }
792: lx = lg(y); if (lx==1) { avma=av; return gzero; }
793: for (i=1; i<lx; i++)
794: {
795: p1 = (GEN)y[i];
796: if (typ(p1) == t_POL) setvarn(p1, vb); else p1 = scalarpol(p1, vb);
797: if (lb) p1 = poleval(p1, gmul(polx[vb],lb));
798: y[i] = la? ldiv(p1,la): (long)p1;
799: }
800: return gerepilecopy(av,y);
801: }
802:
803: GEN
804: nfisisom(GEN a, GEN b)
805: {
806: return nfiso0(a,b,1);
807: }
808:
809: GEN
810: nfisincl(GEN a, GEN b)
811: {
812: return nfiso0(a,b,0);
813: }
814:
815: /*************************************************************************/
816: /** **/
817: /** INITALG **/
818: /** **/
819: /*************************************************************************/
820: extern GEN LLL_nfbasis(GEN *x, GEN polr, GEN base, long prec);
821: extern GEN eleval(GEN f,GEN h,GEN a);
822: extern int canon_pol(GEN z);
823: extern GEN mat_to_vecpol(GEN x, long v);
824: extern GEN vecpol_to_mat(GEN v, long n);
825:
826: /* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */
827: GEN
828: quicktrace(GEN x, GEN sym)
829: {
830: GEN p1 = gzero;
831: long i;
832:
833: if (signe(x))
834: {
835: sym--;
836: for (i=lgef(x)-1; i>1; i--)
837: p1 = gadd(p1, gmul((GEN)x[i],(GEN)sym[i]));
838: }
839: return p1;
840: }
841:
842: /* T = trace(w[i]), x = sum x[i] w[i], x[i] integer */
843: static GEN
844: trace_col(GEN x, GEN T)
845: {
846: ulong av = avma;
847: GEN t = gzero;
848: long i, l = lg(x);
849:
850: t = mulii((GEN)x[1],(GEN)T[1]);
851: for (i=2; i<l; i++)
852: if (signe(x[i])) t = addii(t, mulii((GEN)x[i],(GEN)T[i]));
853: return gerepileuptoint(av, t);
854: }
855:
856: /* Seek a new, simpler, polynomial pol defining the same number field as
857: * x (assumed to be monic at this point).
858: * *ptx receives pol
859: * *ptdx receives disc(pol)
860: * *ptrev expresses the new root in terms of the old one.
861: * base updated in place.
862: * prec = 0 <==> field is totally real
863: */
864: static void
865: nfinit_reduce(long flag, GEN *px, GEN *pdx, GEN *prev, GEN *pbase, long prec)
866: {
867: GEN chbas,a,phimax,dxn,s,sn,p1,p2,p3,polmax,rev,polr;
868: GEN x = *px, dx = *pdx, base = *pbase;
869: long i,j,nmax,numb,flc,v=varn(x), n=lg(base)-1;
870:
871: if (n == 1)
872: {
873: *px = gsub(polx[v],gun); *pdx = gun;
874: *prev = polymodrecip(gmodulcp(gun, x)); return;
875: }
876: polr = prec? roots(x, prec): NULL;
877: if (polr)
878: for (s=gzero,i=1; i<=n; i++) s = gadd(s,gnorm((GEN)polr[i]));
879: else
880: s = subii(sqri((GEN)x[n+1]), shifti((GEN)x[n],1));
881:
882: chbas = LLL_nfbasis(&x,polr,base,prec);
883: if (DEBUGLEVEL) msgtimer("LLL basis");
884:
885: phimax=polx[v]; polmax=dummycopy(x);
886: nmax=(flag & nf_PARTIAL)?min(n,3):n;
887:
888: for (numb=0,i=2; i<=nmax || (!numb && i<=n); i++)
889: { /* cf pols_for_polred */
890: if (DEBUGLEVEL>2) { fprintferr("i = %ld\n",i); flusherr(); }
891: p1 = a = gmul(base,(GEN)chbas[i]); p3=content(p1);
892: if (gcmp1(p3)) p3 = NULL; else p1 = gdiv(p1,p3);
893: p1 = caract2(x,p1,v);
894: if (p3)
895: for (p2=gun, j=lgef(p1)-2; j>=2; j--)
896: {
897: p2 = gmul(p2,p3); p1[j] = lmul((GEN)p1[j], p2);
898: }
899: p2 = modulargcd(derivpol(p1),p1);
900: if (lgef(p2) > 3) continue;
901:
902: if (DEBUGLEVEL>3) outerr(p1);
903: dxn=discsr(p1); flc=absi_cmp(dxn,dx); numb++;
904: if (flc>0) continue;
905:
906: if (polr)
907: for (sn=gzero,j=1; j<=n; j++)
908: sn = gadd(sn,gnorm(poleval(a,(GEN)polr[j])));
909: else
910: sn = subii(sqri((GEN)p1[n+1]), shifti((GEN)p1[n],1));
911: if (flc>=0)
912: {
913: flc=gcmp(sn,s);
914: if (flc>0 || (!flc && gpolcomp(p1,polmax) >= 0)) continue;
915: }
916: dx=dxn; s=sn; polmax=p1; phimax=a;
917: }
918: if (!numb)
919: {
920: if (gisirreducible(x)!=gun) err(redpoler,"nfinit_reduce");
921: err(talker,"you have found a counter-example to a conjecture, "
922: "please send us\nthe polynomial as soon as possible");
923: }
924: if (phimax == polx[v]) /* no improvement */
925: rev = gmodulcp(phimax,x);
926: else
927: {
928: if (canon_pol(polmax) < 0) phimax = gneg_i(phimax);
929: if (DEBUGLEVEL>1) fprintferr("polmax = %Z\n",polmax);
930: p2 = gmodulcp(phimax,x); rev = polymodrecip(p2);
931: a = base; base = cgetg(n+1,t_VEC);
932: p1 = (GEN)rev[2];
933: for (i=1; i<=n; i++)
934: base[i] = (long)eleval(polmax, (GEN)a[i], p1);
935: p1 = vecpol_to_mat(base,n); p2=denom(p1); p1=gmul(p2,p1);
936: p1 = gdiv(hnfmodid(p1,p2), p2);
937: base = mat_to_vecpol(p1,v);
938: }
939: *px=polmax; *pdx=dx; *prev=rev; *pbase = base;
940: }
941:
942: /* pol belonging to Z[x], return a monic polynomial generating the same field
943: * as pol (x-> ax+b)) set lead = NULL if pol was monic (after dividing
944: * by the content), and to to leading coeff otherwise.
945: * No garbage collecting done.
946: */
947: GEN
948: pol_to_monic(GEN pol, GEN *lead)
949: {
950: long n = lgef(pol)-1;
951: GEN p1;
952:
953: if (n==1 || gcmp1((GEN)pol[n])) { *lead = NULL; return pol; }
954:
955: p1=content(pol); if (!gcmp1(p1)) pol = gdiv(pol,p1);
956: return primitive_pol_to_monic(pol,lead);
957: }
958:
959: /* NB: TI integral, det(TI) = d_K, ideal/dideal = codifferent */
960: GEN
961: make_MDI(GEN nf, GEN TI, GEN *ideal, GEN *dideal)
962: {
963: GEN c = content(TI);
964: GEN p1;
965:
966: *dideal = divii((GEN)nf[3], c);
967: *ideal = p1 = hnfmodid(gdiv(TI,c), *dideal);
968: return gmul(c, ideal_two_elt(nf, p1));
969: }
970:
971: /* [bas[i]/den[i]]= integer basis. roo = real part of the roots */
972: GEN
973: make_M(GEN basden,GEN roo)
974: {
975: GEN p1,d,M, bas=(GEN)basden[1], den=(GEN)basden[2];
976: long i,j, ru = lg(roo), n = lg(bas);
977: M = cgetg(n,t_MAT);
978: for (j=1; j<n; j++)
979: {
980: p1=cgetg(ru,t_COL); M[j]=(long)p1;
981: for (i=1; i<ru; i++)
982: p1[i]=(long)poleval((GEN)bas[j],(GEN)roo[i]);
983: }
984: if (den)
985: {
986: long prec = precision((GEN)roo[1]);
987: GEN invd,rd = cgetr(prec+1);
988: for (j=1; j<n; j++)
989: {
990: d = (GEN)den[j]; if (!d) continue;
991: p1 = (GEN)M[j]; affir(d,rd); invd = ginv(rd);
992: for (i=1; i<ru; i++) p1[i] = lmul((GEN)p1[i], invd);
993: }
994: }
995: if (DEBUGLEVEL>4) msgtimer("matrix M");
996: return M;
997: }
998:
999: GEN
1000: make_MC(long r1,GEN M)
1001: {
1002: long i,j,av,tetpil, n = lg(M), ru = lg(M[1]);
1003: GEN p1,p2,MC=cgetg(ru,t_MAT);
1004:
1005: for (j=1; j<ru; j++)
1006: {
1007: p1=cgetg(n,t_COL); MC[j]=(long)p1;
1008: for (i=1; i<n; i++)
1009: {
1010: av=avma; p2=gconj(gcoeff(M,j,i)); tetpil=avma;
1011: p1[i] = (j<=r1)? (long)p2: lpile(av,tetpil,gmul2n(p2,1));
1012: }
1013: }
1014: if (DEBUGLEVEL>4) msgtimer("matrix MC");
1015: return MC;
1016: }
1017:
1018: GEN
1019: get_roots(GEN x,long r1,long ru,long prec)
1020: {
1021: GEN roo = (typ(x)==t_VEC)? dummycopy(x): roots(x,prec);
1022: long i;
1023:
1024: for (i=1; i<=r1; i++) roo[i]=lreal((GEN)roo[i]);
1025: for ( ; i<=ru; i++) roo[i]=roo[(i<<1)-r1];
1026: roo[0]=evaltyp(t_VEC)|evallg(ru+1); return roo;
1027: }
1028:
1029: /* assume lg(nf) > 3 && typ(nf) = container [hopefully a genuine nf] */
1030: long
1031: nf_get_r1(GEN nf)
1032: {
1033: GEN x = (GEN)nf[2];
1034: if (typ(x) != t_VEC || lg(x) != 3 || typ(x[1]) != t_INT)
1035: err(talker,"false nf in nf_get_r1");
1036: return itos((GEN)x[1]);
1037: }
1038: long
1039: nf_get_r2(GEN nf)
1040: {
1041: GEN x = (GEN)nf[2];
1042: if (typ(x) != t_VEC || lg(x) != 3 || typ(x[2]) != t_INT)
1043: err(talker,"false nf in nf_get_r2");
1044: return itos((GEN)x[2]);
1045: }
1046:
1047: extern GEN mulmat_pol(GEN A, GEN x);
1048:
1049: static GEN
1050: get_T(GEN mul, GEN x, GEN bas, GEN den)
1051: {
1052: long i,j,n = lg(bas)-1;
1053: GEN tr,T,T1,sym;
1054: T = cgetg(n+1,t_MAT);
1055: T1 = cgetg(n+1,t_COL);
1056: sym = polsym(x, n-1);
1057:
1058: T1[1]=lstoi(n);
1059: for (i=2; i<=n; i++)
1060: {
1061: tr = quicktrace((GEN)bas[i], sym);
1062: if (den && den[i]) tr = diviiexact(tr,(GEN)den[i]);
1063: T1[i] = (long)tr; /* tr(w[i]) */
1064: }
1065: T[1] = (long)T1;
1066: for (i=2; i<=n; i++)
1067: {
1068: T[i] = lgetg(n+1,t_COL); coeff(T,1,i) = T1[i];
1069: for (j=2; j<=i; j++)
1070: coeff(T,i,j) = coeff(T,j,i) = (long)trace_col((GEN)mul[j+(i-1)*n], T1);
1071: }
1072: return T;
1073: }
1074:
1075: /* return [bas[i]*denom(bas[i]), denom(bas[i])], denom 1 is given as NULL */
1076: GEN
1077: get_bas_den(GEN bas)
1078: {
1079: GEN z,d,den, dbas = dummycopy(bas);
1080: long i, c = 0, l = lg(bas);
1081: den = cgetg(l,t_VEC);
1082: for (i=1; i<l; i++)
1083: {
1084: d = denom(content((GEN)dbas[i]));
1085: if (is_pm1(d)) d = NULL; else { dbas[i] = lmul((GEN)dbas[i],d); c++; }
1086: den[i] = (long)d;
1087: }
1088: if (!c) den = NULL; /* power basis */
1089: z = cgetg(3,t_VEC);
1090: z[1] = (long)dbas;
1091: z[2] = (long)den; return z;
1092: }
1093:
1094: /* allow x or y = NULL (act as 1) */
1095: static GEN
1096: _mulii(GEN x, GEN y)
1097: {
1098: if (!x) return y;
1099: if (!y) return x;
1100: return mulii(x,y);
1101: }
1102:
1103: GEN
1104: get_mul_table(GEN x,GEN basden,GEN invbas,GEN *T)
1105: {
1106: long i,j, n = degpol(x);
1107: GEN z,d, mul = cgetg(n*n+1,t_MAT), bas=(GEN)basden[1], den=(GEN)basden[2];
1108:
1109: for (j=1; j<=n*n; j++) mul[j]=lgetg(n+1,t_COL);
1110: for (i=1; i<=n; i++)
1111: for (j=i; j<=n; j++)
1112: {
1113: ulong av = avma;
1114: z = gres(gmul((GEN)bas[j],(GEN)bas[i]), x);
1115: z = mulmat_pol(invbas, z); /* integral column */
1116: if (den)
1117: {
1118: d = _mulii((GEN)den[i], (GEN)den[j]);
1119: if (d) z = gdivexact(z, d);
1120: }
1121: mul[j+(i-1)*n] = mul[i+(j-1)*n] = lpileupto(av,z);
1122: }
1123: if (T) *T = get_T(mul,x,bas,den);
1124: return mul;
1125: }
1126:
1127: /* fill mat = nf[5], as well as nf[8] and nf[9]
1128: * If (small) only compute a subset (use dummy 0s for the rest) */
1129: void
1130: get_nf_matrices(GEN nf, long small)
1131: {
1132: GEN x=(GEN)nf[1],dK=(GEN)nf[3],index=(GEN)nf[4],ro=(GEN)nf[6],bas=(GEN)nf[7];
1133: GEN basden,mul,invbas,M,MC,T,MDI,D,TI,A,dA,mat;
1134: long r1 = nf_get_r1(nf), n = lg(bas)-1;
1135:
1136: mat = cgetg(small? 4: 8,t_VEC); nf[5] = (long)mat;
1137: basden = get_bas_den(bas);
1138: M = make_M(basden,ro);
1139: MC = make_MC(r1,M);
1140: mat[1]=(long)M;
1141: mat[3]=(long)mulmat_real(MC,M);
1142: if (small) { nf[8]=nf[9]=mat[2]=zero; return; }
1143:
1144: invbas = QM_inv(vecpol_to_mat(bas,n), gun);
1145: mul = get_mul_table(x,basden,invbas,&T);
1146: if (DEBUGLEVEL) msgtimer("mult. table");
1147:
1148: nf[8]=(long)invbas;
1149: nf[9]=(long)mul;
1150:
1151: TI = ZM_inv(T, dK); /* dK T^-1 */
1152: MDI = make_MDI(nf,TI, &A, &dA);
1153: mat[6]=(long)TI;
1154: mat[7]=(long)MDI; /* needed in idealinv below */
1155: if (is_pm1(index))
1156: D = idealhermite_aux(nf, derivpol(x));
1157: else
1158: D = gmul(dA, idealinv(nf, A));
1159: mat[2]=(long)MC;
1160: mat[4]=(long)T;
1161: mat[5]=(long)D;
1162: if (DEBUGLEVEL) msgtimer("matrices");
1163: }
1164:
1165: /* Initialize the number field defined by the polynomial x (in variable v)
1166: * flag & nf_REGULAR
1167: * regular behaviour.
1168: * flag & nf_SMALL
1169: * compute only nf[1] (pol), nf[2] (signature), nf[5][3] (T2) and
1170: * nf[7] (integer basis), the other components are filled with gzero.
1171: * flag & nf_REDUCE
1172: * try a polred first.
1173: * flag & nf_PARTIAL
1174: * do a partial polred, not a polredabs
1175: * flag & nf_ORIG
1176: * do a polred and return [nfinit(x),Mod(a,red)], where
1177: * Mod(a,red)=Mod(v,x) (i.e return the base change).
1178: */
1179:
1180: /* here x can be a polynomial, an nf or a bnf */
1181: GEN
1182: initalgall0(GEN x, long flag, long prec)
1183: {
1184: GEN lead = NULL,nf,ro,bas,mat,rev,dK,dx,index,res;
1185: long av=avma,n,i,r1,r2,ru,PRECREG;
1186:
1187: if (DEBUGLEVEL) timer2();
1188: if (typ(x)==t_POL)
1189: {
1190: n=degpol(x); if (n<=0) err(constpoler,"initalgall0");
1191: check_pol_int(x,"initalg");
1192: if (gisirreducible(x) == gzero) err(redpoler,"nfinit");
1193: if (!gcmp1((GEN)x[n+2]))
1194: {
1195: x = pol_to_monic(x,&lead);
1196: if (!(flag & nf_SMALL))
1197: {
1198: if (!(flag & nf_REDUCE))
1199: err(warner,"non-monic polynomial. Result of the form [nf,c]");
1200: flag = flag | nf_REDUCE | nf_ORIG;
1201: }
1202: }
1203: bas = allbase4(x,0,&dK,NULL);
1204: if (DEBUGLEVEL) msgtimer("round4");
1205: dx = discsr(x); r1 = sturm(x);
1206: }
1207: else
1208: {
1209: i = lg(x);
1210: if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL)
1211: { /* polynomial + integer basis */
1212: bas=(GEN)x[2]; x=(GEN)x[1]; n=degpol(x);
1213: if (typ(bas) == t_MAT) { mat = bas; bas = mat_to_vecpol(mat,varn(x)); }
1214: else mat = vecpol_to_mat(bas,n);
1215: dx = discsr(x); r1 = sturm(x);
1216: dK = gmul(dx, gsqr(det2(mat)));
1217: }
1218: else
1219: {
1220: nf=checknf(x);
1221: bas=(GEN)nf[7]; x=(GEN)nf[1]; n=degpol(x);
1222: dK=(GEN)nf[3]; dx=mulii(dK, sqri((GEN)nf[4]));
1223: r1 = nf_get_r1(nf);
1224: }
1225: bas[1]=lpolun[varn(x)]; /* it may be gun => SEGV later */
1226: }
1227: r2=(n-r1)>>1; ru=r1+r2;
1228: PRECREG = prec + (expi(dK)>>(TWOPOTBITS_IN_LONG+1))
1229: + (long)((sqrt((double)n)+3) / sizeof(long) * 4);
1230: rev = NULL;
1231: if (flag & nf_REDUCE)
1232: {
1233: nfinit_reduce(flag, &x, &dx, &rev, &bas, r1==n? 0: prec);
1234: if (DEBUGLEVEL) msgtimer("polred");
1235: }
1236: if (!carrecomplet(divii(dx,dK),&index))
1237: err(bugparier,"nfinit (incorrect discriminant)");
1238:
1239: ro=get_roots(x,r1,ru,PRECREG);
1240: if (DEBUGLEVEL) msgtimer("roots");
1241:
1242: nf=cgetg(10,t_VEC);
1243: nf[1]=(long)x;
1244: nf[2]=lgetg(3,t_VEC);
1245: mael(nf,2,1)=lstoi(r1);
1246: mael(nf,2,2)=lstoi(r2);
1247: nf[3]=(long)dK;
1248: nf[4]=(long)index;
1249: nf[6]=(long)ro;
1250: nf[7]=(long)bas;
1251: get_nf_matrices(nf, flag & nf_SMALL);
1252:
1253: if (flag & nf_ORIG)
1254: {
1255: if (!rev) err(talker,"bad flag in initalgall0");
1256: res = cgetg(3,t_VEC);
1257: res[1]=(long)nf; nf = res;
1258: res[2]=lead? ldiv(rev,lead): (long)rev;
1259: }
1260: return gerepilecopy(av, nf);
1261: }
1262:
1263: GEN
1264: initalgred(GEN x, long prec)
1265: {
1266: return initalgall0(x,nf_REDUCE,prec);
1267: }
1268:
1269: GEN
1270: initalgred2(GEN x, long prec)
1271: {
1272: return initalgall0(x,nf_REDUCE|nf_ORIG,prec);
1273: }
1274:
1275: GEN
1276: nfinit0(GEN x, long flag,long prec)
1277: {
1278: switch(flag)
1279: {
1280: case 0:
1281: case 1: return initalgall0(x,nf_REGULAR,prec);
1282: case 2: return initalgall0(x,nf_REDUCE,prec);
1283: case 3: return initalgall0(x,nf_REDUCE|nf_ORIG,prec);
1284: case 4: return initalgall0(x,nf_REDUCE|nf_PARTIAL,prec);
1285: case 5: return initalgall0(x,nf_REDUCE|nf_ORIG|nf_PARTIAL,prec);
1286: case 6: return initalgall0(x,nf_SMALL,prec);
1287: default: err(flagerr,"nfinit");
1288: }
1289: return NULL; /* not reached */
1290: }
1291:
1292: GEN
1293: initalg(GEN x, long prec)
1294: {
1295: return initalgall0(x,nf_REGULAR,prec);
1296: }
1297:
1298: /* assume x a bnr/bnf/nf */
1299: long
1300: nfgetprec(GEN x)
1301: {
1302: GEN nf = checknf(x), ro = (GEN)nf[6];
1303: return (typ(ro)==t_VEC)?precision((GEN)ro[1]): DEFAULTPREC;
1304: }
1305:
1306: GEN
1307: nfnewprec(GEN nf, long prec)
1308: {
1309: long av=avma,r1,r2,ru,n;
1310: GEN y,pol,ro,basden,MC,mat,M;
1311:
1312: if (typ(nf) != t_VEC) err(talker,"incorrect nf in nfnewprec");
1313: if (lg(nf) == 11) return bnfnewprec(nf,prec);
1314: if (lg(nf) == 7) return bnrnewprec(nf,prec);
1315: (void)checknf(nf);
1316: y = dummycopy(nf);
1317: pol=(GEN)nf[1]; n=degpol(pol);
1318: r1 = nf_get_r1(nf);
1319: r2 = (n - r1) >> 1; ru = r1+r2;
1320: mat=dummycopy((GEN)nf[5]);
1321: ro=get_roots(pol,r1,ru,prec);
1322: y[5]=(long)mat;
1323: y[6]=(long)ro;
1324: basden = get_bas_den((GEN)nf[7]);
1325: M = make_M(basden,ro);
1326: MC = make_MC(r1,M);
1327: mat[1]=(long)M;
1328: if (typ(nf[8]) != t_INT) mat[2]=(long)MC; /* not a small nf */
1329: mat[3]=(long)mulmat_real(MC,M);
1330: return gerepilecopy(av, y);
1331: }
1332:
1333: static long
1334: nf_pm1(GEN y)
1335: {
1336: long i,l;
1337:
1338: if (!is_pm1(y[1])) return 0;
1339: l = lg(y);
1340: for (i=2; i<l; i++)
1341: if (signe(y[i])) return 0;
1342: return signe(y[1]);
1343:
1344: }
1345:
1346: static GEN
1347: is_primitive_root(GEN nf, GEN fa, GEN x, long w)
1348: {
1349: GEN y, exp = stoi(2), pp = (GEN)fa[1];
1350: long i,p, l = lg(pp);
1351:
1352: for (i=1; i<l; i++)
1353: {
1354: p = itos((GEN)pp[i]);
1355: exp[2] = w / p; y = element_pow(nf,x,exp);
1356: if (nf_pm1(y) > 0) /* y = 1 */
1357: {
1358: if (p!=2 || !gcmp1(gcoeff(fa,i,2))) return NULL;
1359: x = gneg_i(x);
1360: }
1361: }
1362: return x;
1363: }
1364:
1365: GEN
1366: rootsof1(GEN nf)
1367: {
1368: ulong av;
1369: long N,k,i,ws,prec;
1370: GEN algun,p1,y,R1,d,list,w;
1371:
1372: y=cgetg(3,t_VEC); av=avma; nf=checknf(nf);
1373: R1=gmael(nf,2,1); algun=gmael(nf,8,1);
1374: if (signe(R1))
1375: {
1376: y[1]=deux;
1377: y[2]=lneg(algun); return y;
1378: }
1379: N=degpol(nf[1]); prec=gprecision((GEN)nf[6]);
1380: #ifdef LONG_IS_32BIT
1381: if (prec < 10) prec = 10;
1382: #else
1383: if (prec < 6) prec = 6;
1384: #endif
1385: for (i=1; ; i++)
1386: {
1387: p1 = fincke_pohst(nf,stoi(N),1000,1,prec,NULL);
1388: if (p1) break;
1389: if (i == MAXITERPOL) err(accurer,"rootsof1");
1390: prec=(prec<<1)-2;
1391: if (DEBUGLEVEL) err(warnprec,"rootsof1",prec);
1392: nf=nfnewprec(nf,prec);
1393: }
1394: if (itos(ground((GEN)p1[2])) != N) err(bugparier,"rootsof1 (bug1)");
1395: w=(GEN)p1[1]; ws = itos(w);
1396: if (ws == 2)
1397: {
1398: y[1]=deux; avma=av;
1399: y[2]=lneg(algun); return y;
1400: }
1401:
1402: d = decomp(w); list = (GEN)p1[3]; k = lg(list);
1403: for (i=1; i<k; i++)
1404: {
1405: p1 = (GEN)list[i];
1406: p1 = is_primitive_root(nf,d,p1,ws);
1407: if (p1)
1408: {
1409: y[2]=lpilecopy(av,p1);
1410: y[1]=lstoi(ws); return y;
1411: }
1412: }
1413: err(bugparier,"rootsof1");
1414: return NULL; /* not reached */
1415: }
1416:
1417: /*******************************************************************/
1418: /* */
1419: /* DEDEKIND ZETA FUNCTION */
1420: /* */
1421: /*******************************************************************/
1422:
1423: ulong smulss(ulong x, ulong y, ulong *rem);
1424:
1425: static GEN
1426: dirzetak0(GEN nf, long N0)
1427: {
1428: GEN vect,p1,pol,disc,c,c2;
1429: long av=avma,i,j,k,limk,lx;
1430: ulong q,p,rem;
1431: byteptr d=diffptr;
1432: long court[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0};
1433:
1434: pol=(GEN)nf[1]; disc=(GEN)nf[4];
1435: c = (GEN) gpmalloc((N0+1)*sizeof(long));
1436: c2 = (GEN) gpmalloc((N0+1)*sizeof(long));
1437: c2[0]=c[0]=evaltyp(t_VEC) | evallg(N0+1);
1438: c2[1]=c[1]=1; for (i=2; i<=N0; i++) c[i]=0;
1439: court[2] = 0;
1440:
1441: while (court[2]<=N0)
1442: {
1443: court[2] += *d++; if (! *d) err(primer1);
1444: if (smodis(disc,court[2])) /* court does not divide index */
1445: { vect = (GEN) simplefactmod(pol,court)[1]; lx=lg(vect); }
1446: else
1447: {
1448: p1=primedec(nf,court); lx=lg(p1); vect=cgetg(lx,t_COL);
1449: for (i=1; i<lx; i++) vect[i]=mael(p1,i,4);
1450: }
1451: for (j=1; j<lx; j++)
1452: {
1453: p1=powgi(court, (GEN)vect[j]); /* p1 = court^f */
1454: if (cmpis(p1,N0) <= 0)
1455: {
1456: q=p=p1[2]; limk=N0/q;
1457: for (k=2; k<=N0; k++) c2[k]=c[k];
1458: while (q<=(ulong)N0)
1459: {
1460: for (k=1; k<=limk; k++) c2[k*q] += c[k];
1461: q = smulss(q,p,&rem);
1462: if (rem) break;
1463: limk /= p;
1464: }
1465: p1=c; c=c2; c2=p1;
1466: }
1467: }
1468: avma=av;
1469: if (DEBUGLEVEL>6) fprintferr(" %ld",court[2]);
1470: }
1471: if (DEBUGLEVEL>6) { fprintferr("\n"); flusherr(); }
1472: free(c2); return c;
1473: }
1474:
1475: GEN
1476: dirzetak(GEN nf, GEN b)
1477: {
1478: GEN z,c;
1479: long i;
1480:
1481: if (typ(b)!=t_INT) err(talker,"not an integer type in dirzetak");
1482: if (signe(b)<=0) return cgetg(1,t_VEC);
1483: nf = checknf(nf);
1484: if (is_bigint(b)) err(talker,"too many terms in dirzetak");
1485: c = dirzetak0(nf,itos(b));
1486: i = lg(c); z=cgetg(i,t_VEC);
1487: for (i-- ; i; i--) z[i]=lstoi(c[i]);
1488: free(c); return z;
1489: }
1490:
1491: GEN
1492: initzeta(GEN pol, long prec)
1493: {
1494: GEN nfz,nf,alpha,beta,mu,gr1,gr2,gru,p1,p2,cst,A0,c0,c1,c2,eps,coef;
1495: GEN limx,bnf,resi,zet,C,coeflog,racpi,aij,tabj,colzero, *tabcstn, *tabcstni;
1496: GEN c_even,ck_even,c_odd,ck_odd,serie_even,serie_odd,serie_exp,Pi;
1497: long N0,imin,imax,r1,r2,ru,R,N,i,j,k,n, av,av2,tetpil;
1498: long court[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0};
1499: stackzone *zone, *zone0, *zone1;
1500:
1501: /*************** Calcul du residu et des constantes ***************/
1502: eps=gmul2n(gun,-bit_accuracy(prec)-6); p1=dbltor(0.5);
1503: nfz=cgetg(10,t_VEC);
1504: bnf=buchinit(pol,p1,p1,prec+1); prec=(prec<<1)-1;
1505: bnf = checkbnf_discard(bnf);
1506: Pi = mppi(prec); racpi=gsqrt(Pi,prec);
1507:
1508: /* Nb de classes et regulateur */
1509: nf=(GEN)bnf[7]; N=degpol(nf[1]);
1510: r1 = nf_get_r1(nf); r2 = (N-r1)>>1;
1511: gr1=gmael(nf,2,1); gr2=gmael(nf,2,2);
1512: ru=r1+r2; R=ru+2;
1513: av=avma; p1=(GEN)bnf[8]; p2 = gmul(gmul2n(gmael(p1,1,1),r1), (GEN)p1[2]);
1514: tetpil = avma; resi=gerepile(av,tetpil,gdiv(p2, gmael(p1,4,1)));
1515:
1516: /* Calcul de N0 */
1517: cst = cgetr(prec); av = avma;
1518: mu = gadd(gmul2n(gr1,-1),gr2);
1519: alpha = gmul2n(stoi(ru+1),-1);
1520: beta = gpui(gdeux,gmul2n(gr1,-1),DEFAULTPREC);
1521: A0 = gmul2n(gpowgs(mu,R),r1);
1522: A0 = gmul(A0,gpowgs(gmul2n(Pi,1),1-ru));
1523: A0 = gsqrt(A0,DEFAULTPREC);
1524:
1525: c1 = gmul(mu,gpui(beta,ginv(mu),DEFAULTPREC));
1526: c0 = gdiv(gmul(A0,gpowgs(gmul2n(Pi,1),ru-1)),mu);
1527: c0 = gmul(c0,gpui(c1,gneg_i(alpha),DEFAULTPREC));
1528: c2 = gdiv(alpha,mu);
1529:
1530: p1 = glog(gdiv(c0,eps),DEFAULTPREC);
1531: limx = gdiv(gsub(glog(p1,DEFAULTPREC),glog(c1,DEFAULTPREC)),
1532: gadd(c2,gdiv(p1,mu)));
1533: limx = gmul(gpui(gdiv(c1,p1),mu,DEFAULTPREC),
1534: gadd(gun,gmul(alpha,limx)));
1535: p1 = gsqrt(absi((GEN)nf[3]),prec);
1536: p2 = gmul2n(gpowgs(racpi,N),r2);
1537: gaffect(gdiv(p1,p2), cst);
1538:
1539: av = avma; p1 = gfloor(gdiv(cst,limx)); N0 = p1[2];
1540: if (is_bigint(p1) || N0 > 10000000)
1541: err(talker,"discriminant too large for initzeta, sorry");
1542: if (DEBUGLEVEL>=2)
1543: { fprintferr("\ninitzeta:\nN0 = %ld\n",N0); flusherr(); timer2(); }
1544:
1545: /* Calcul de imax */
1546:
1547: imin=1; imax=1400;
1548: p1 = gmul(gpowgs(gmul2n(racpi,1),r2),gpowgs(stoi(5),r1));
1549: p1 = gdiv(p1,gmul(gmul(gsqrt(limx,DEFAULTPREC),gmul2n(eps,4)),
1550: gpowgs(racpi,3)));
1551: while (imax-imin >= 4)
1552: {
1553: long itest = (imax+imin)>>1;
1554: p2 = gmul(gpowgs(mpfactr(itest,DEFAULTPREC),r2),gpowgs(limx,itest));
1555: p2 = gmul(p2,gpowgs(mpfactr(itest/2,DEFAULTPREC),r1));
1556: if (gcmp(p2,p1) >= 0) imax=itest; else imin=itest;
1557: }
1558: imax -= (imax & 1); avma = av;
1559: if (DEBUGLEVEL>=2) { fprintferr("imax = %ld\n",imax); flusherr(); }
1560: /* Tableau des i/cst (i=1 a N0) */
1561:
1562: i = prec*N0;
1563: zone = switch_stack(NULL,i + 2*(N0+1) + 6*prec);
1564: zone1 = switch_stack(NULL,2*i);
1565: zone0 = switch_stack(NULL,2*i);
1566: switch_stack(zone,1);
1567: tabcstn = (GEN*) cgetg(N0+1,t_VEC);
1568: tabcstni = (GEN*) cgetg(N0+1,t_VEC);
1569: p1 = ginv(cst);
1570: for (i=1; i<=N0; i++) { tabcstni[i] = gun; tabcstn[i] = mulsr(i,p1); }
1571: switch_stack(zone,0);
1572:
1573: /********** Calcul des coefficients a(i,j) independants de s **********/
1574:
1575: zet=cgetg(R,t_VEC); zet[1] = lmpeuler(prec);
1576: for (i=2; i<R; i++) zet[i] = (long)gzeta(stoi(i),prec);
1577:
1578: aij=cgetg(imax+1,t_VEC);
1579: for (i=1; i<=imax; i++) aij[i] = lgetg(R,t_VEC);
1580:
1581: c_even = realun(prec);
1582: c_even = gmul2n(c_even,r1);
1583: c_odd = gmul(c_even,gpowgs(racpi,r1));
1584: if (ru&1) c_odd=gneg_i(c_odd);
1585: ck_even=cgetg(R,t_VEC); ck_odd=cgetg(r2+2,t_VEC);
1586: for (k=1; k<R; k++)
1587: {
1588: ck_even[k]=lmul((GEN)zet[k],gadd(gr2,gmul2n(gr1,-k)));
1589: if (k&1) ck_even[k]=lneg((GEN)ck_even[k]);
1590: }
1591: gru=stoi(ru);
1592: for (k=1; k<=r2+1; k++)
1593: {
1594: ck_odd[k]=lmul((GEN)zet[k], gsub(gru, gmul2n(gr1,-k)));
1595: if (k&1) ck_odd[k]=lneg((GEN)ck_odd[k]);
1596: ck_odd[k]=ladd(gru,(GEN)ck_odd[k]);
1597: }
1598: ck_odd[1]=lsub((GEN)ck_odd[1],gmul(gr1,glog(gdeux,prec)));
1599: serie_even =cgetg(ru+3,t_SER); serie_odd=cgetg(r2+3,t_SER);
1600: serie_even[1] = serie_odd[1] = evalsigne(1)+evalvalp(1);
1601: i=0;
1602:
1603: while (i<imax/2)
1604: {
1605: for (k=1; k<R; k++)
1606: serie_even[k+1]=ldivgs((GEN)ck_even[k],k);
1607: serie_exp=gmul(c_even,gexp(serie_even,0));
1608: p1=(GEN)aij[2*i+1];
1609: for (j=1; j<R; j++) p1[j]=serie_exp[ru+3-j];
1610:
1611: for (k=1; k<=r2+1; k++)
1612: serie_odd[k+1]=ldivgs((GEN)ck_odd[k],k);
1613: serie_exp=gmul(c_odd,gexp(serie_odd,0));
1614: p1=(GEN)aij[2*i+2];
1615: for (j=1; j<=r2+1; j++) p1[j]=serie_exp[r2+3-j];
1616: for ( ; j<R; j++) p1[j]=zero;
1617: i++;
1618:
1619: c_even = gdiv(c_even,gmul(gpowgs(stoi(i),ru),gpowgs(stoi(2*i-1),r2)));
1620: c_odd = gdiv(c_odd, gmul(gpowgs(stoi(i),r2),gpowgs(stoi(2*i+1),ru)));
1621: c_even = gmul2n(c_even,-r2);
1622: c_odd = gmul2n(c_odd,r1-r2);
1623: if (r1&1) { c_even=gneg_i(c_even); c_odd=gneg_i(c_odd); }
1624: p1 = gr2; p2 = gru;
1625: for (k=1; k<R; k++)
1626: {
1627: p1=gdivgs(p1,2*i-1); p2=gdivgs(p2,2*i);
1628: ck_even[k] = ladd((GEN)ck_even[k], gadd(p1,p2));
1629: }
1630: p1 = gru; p2 = gr2;
1631: for (k=1; k<=r2+1; k++)
1632: {
1633: p1=gdivgs(p1,2*i+1); p2=gdivgs(p2,2*i);
1634: ck_odd[k] = ladd((GEN)ck_odd[k], gadd(p1,p2));
1635: }
1636: }
1637: aij=gerepilecopy(av,aij);
1638: if (DEBUGLEVEL>=2) msgtimer("a(i,j)");
1639: p1=cgetg(5,t_VEC);
1640: p1[1]=lstoi(r1); p1[2]=lstoi(r2); p1[3]=lstoi(imax); p1[4]=(long)bnf;
1641: nfz[1]=(long)p1;
1642: nfz[2]=(long)resi;
1643: nfz[5]=(long)cst;
1644: nfz[6]=llog(cst,prec);
1645: nfz[7]=(long)aij;
1646:
1647: /************* Calcul du nombre d'ideaux de norme donnee *************/
1648: coef = dirzetak0(nf,N0); tabj = cgetg(N0+1,t_MAT);
1649: if (DEBUGLEVEL>=2) msgtimer("coef");
1650: colzero=cgetg(ru+2,t_COL); for (j=1; j<=ru+1; j++) colzero[j]=zero;
1651: for (i=1; i<=N0; i++)
1652: if (coef[i])
1653: {
1654: GEN t = cgetg(ru+2,t_COL);
1655: p1 = glog((GEN)tabcstn[i],prec); setsigne(p1, -signe(p1));
1656: t[1] = lstoi(coef[i]); t[2] = lmul((GEN)t[1],p1);
1657: for (j=2; j<=ru; j++)
1658: {
1659: long av2 = avma; p2 = gmul((GEN)t[j],p1);
1660: t[j+1] = (long)gerepileuptoleaf(av2, divrs(p2,j));
1661: }
1662: /* tabj[n,j]=coef(n)*ln(c/n)^(j-1)/(j-1)! */
1663: tabj[i] = (long)t;
1664: }
1665: else tabj[i]=(long)colzero;
1666: if (DEBUGLEVEL>=2) msgtimer("a(n)");
1667:
1668: coeflog=cgetg(N0+1,t_VEC); coeflog[1]=zero;
1669: for (i=2; i<=N0; i++)
1670: if (coef[i])
1671: {
1672: court[2]=i; p1=glog(court,prec);
1673: setsigne(p1,-1); coeflog[i]=(long)p1;
1674: }
1675: else coeflog[i] = zero;
1676: if (DEBUGLEVEL>=2) msgtimer("log(n)");
1677:
1678: nfz[3]=(long)tabj;
1679: p1 = cgetg(N0+1,t_VECSMALL);
1680: for (i=1; i<=N0; i++) p1[i] = coef[i];
1681: nfz[8]=(long)p1;
1682: nfz[9]=(long)coeflog;
1683:
1684: /******************** Calcul des coefficients Cik ********************/
1685:
1686: C = cgetg(ru+1,t_MAT);
1687: for (k=1; k<=ru; k++) C[k] = lgetg(imax+1,t_COL);
1688: av2 = avma;
1689: for (i=1; i<=imax; i++)
1690: {
1691: stackzone *z;
1692: for (k=1; k<=ru; k++)
1693: {
1694: p1 = NULL;
1695: for (n=1; n<=N0; n++)
1696: if (coef[n])
1697: for (j=1; j<=ru-k+1; j++)
1698: {
1699: p2 = gmul(tabcstni[n],
1700: gmul(gmael(aij,i,j+k), gmael(tabj,n,j)));
1701: p1 = p1? gadd(p1,p2): p2;
1702: }
1703: coeff(C,i,k) = p1? (long)gerepileupto(av2,p1): zero;
1704: av2 = avma;
1705: }
1706: /* use a parallel stack */
1707: z = i&1? zone1: zone0;
1708: switch_stack(z, 1);
1709: for (n=1; n<=N0; n++)
1710: if (coef[n]) tabcstni[n] = mpmul(tabcstni[n],tabcstn[n]);
1711: /* come back */
1712: switch_stack(z, 0);
1713: }
1714: nfz[4] = (long) C;
1715: if (DEBUGLEVEL>=2) msgtimer("Cik");
1716: gunclone(aij);
1717: free((void*)zone); free((void*)zone1); free((void*)zone0);
1718: free((void*)coef); return nfz;
1719: }
1720:
1721: GEN
1722: gzetakall(GEN nfz, GEN s, long flag, long prec2)
1723: {
1724: GEN resi,C,cst,cstlog,coeflog,cs,coef;
1725: GEN lambd,gammas,gammaunmoins,gammas2,gammaunmoins2,var1,var2;
1726: GEN p1,unmoins,gexpro,gar,val,valm,valk,valkm;
1727: long ts = typ(s), r1,r2,ru,imax,i,j,k,N0,sl,prec,bigprec, av = avma;
1728:
1729: if (typ(nfz)!=t_VEC || lg(nfz)!=10 || typ(nfz[1]) != t_VEC)
1730: err(talker,"not a zeta number field in zetakall");
1731: if (! is_intreal_t(ts) && ts != t_COMPLEX && ! is_frac_t(ts))
1732: err(typeer,"gzetakall");
1733: resi=(GEN)nfz[2]; C=(GEN)nfz[4]; cst=(GEN)nfz[5];
1734: cstlog=(GEN)nfz[6]; coef=(GEN)nfz[8]; coeflog=(GEN)nfz[9];
1735: r1 = itos(gmael(nfz,1,1));
1736: r2 = itos(gmael(nfz,1,2)); ru = r1+r2;
1737: imax = itos(gmael(nfz,1,3)); N0 = lg(coef)-1;
1738: bigprec = precision(cst); prec = prec2+1;
1739:
1740: if (ts==t_COMPLEX && gcmp0(gimag(s))) { s=greal(s); ts = typ(s); }
1741: if (ts==t_REAL && !signe(gfrac(s))) { s=mptrunc(s); ts = t_INT; }
1742: if (ts==t_INT)
1743: {
1744: sl=itos(s);
1745: if (sl==1) err(talker,"s = 1 is a pole (gzetakall)");
1746: if (sl==0)
1747: {
1748: avma = av;
1749: if (flag) err(talker,"s = 0 is a pole (gzetakall)");
1750: if (ru == 1) return gneg(r1? ghalf: resi);
1751: return gzero;
1752: }
1753: if (sl<0 && (r2 || !odd(sl)))
1754: {
1755: if (!flag) return gzero;
1756: s = subsi(1,s); sl = 1-sl;
1757: }
1758: unmoins=subsi(1,s);
1759: lambd = gdiv(resi, mulis(s,sl-1));
1760: gammas2=ggamma(gmul2n(s,-1),prec);
1761: gar=gpowgs(gammas2,r1);
1762: cs=gexp(gmul(cstlog,s),prec);
1763: val=s; valm=unmoins;
1764: if (sl < 0) /* r2 = 0 && odd(sl) */
1765: {
1766: gammaunmoins2=ggamma(gmul2n(unmoins,-1),prec);
1767: var1=var2=gun;
1768: for (i=2; i<=N0; i++)
1769: if (coef[i])
1770: {
1771: gexpro=gexp(gmul((GEN)coeflog[i],s),bigprec);
1772: var1 = gadd(var1,gmulsg(coef[i],gexpro));
1773: var2 = gadd(var2,gdivsg(coef[i],gmulsg(i,gexpro)));
1774: }
1775: lambd=gadd(lambd,gmul(gmul(var1,cs),gar));
1776: lambd=gadd(lambd,gmul(gmul(var2,gdiv(cst,cs)),
1777: gpowgs(gammaunmoins2,r1)));
1778: for (i=1; i<=imax; i+=2)
1779: {
1780: valk = val;
1781: valkm = valm;
1782: for (k=1; k<=ru; k++)
1783: {
1784: p1 = gcoeff(C,i,k);
1785: lambd = gsub(lambd,gdiv(p1,valk )); valk = mulii(val,valk);
1786: lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm);
1787: }
1788: val = addis(val, 2);
1789: valm = addis(valm,2);
1790: }
1791: }
1792: else
1793: {
1794: GEN tabj=(GEN)nfz[3], aij=(GEN)nfz[7];
1795:
1796: gar = gmul(gar,gpowgs(ggamma(s,prec),r2));
1797: var1=var2=gzero;
1798: for (i=1; i<=N0; i++)
1799: if (coef[i])
1800: {
1801: gexpro=gexp(gmul((GEN)coeflog[i],s),bigprec);
1802: var1 = gadd(var1,gmulsg(coef[i],gexpro));
1803: if (sl <= imax)
1804: {
1805: p1=gzero;
1806: for (j=1; j<=ru+1; j++)
1807: p1 = gadd(p1, gmul(gmael(aij,sl,j), gmael(tabj,i,j)));
1808: var2=gadd(var2,gdiv(p1,gmulsg(i,gexpro)));
1809: }
1810: }
1811: lambd=gadd(lambd,gmul(gmul(var1,cs),gar));
1812: lambd=gadd(lambd,gmul(var2,gdiv(cst,cs)));
1813: for (i=1; i<=imax; i++)
1814: {
1815: valk = val;
1816: valkm = valm;
1817: if (i == sl)
1818: for (k=1; k<=ru; k++)
1819: {
1820: p1 = gcoeff(C,i,k);
1821: lambd = gsub(lambd,gdiv(p1,valk)); valk = mulii(val,valk);
1822: }
1823: else
1824: for (k=1; k<=ru; k++)
1825: {
1826: p1 = gcoeff(C,i,k);
1827: lambd = gsub(lambd,gdiv(p1,valk )); valk = mulii(val,valk);
1828: lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm);
1829: }
1830: val = addis(val,1);
1831: valm = addis(valm,1);
1832: }
1833: }
1834: }
1835: else
1836: {
1837: GEN Pi = mppi(prec);
1838: if (is_frac_t(ts))
1839: s = gmul(s, realun(bigprec));
1840: else
1841: s = gprec_w(s, bigprec);
1842:
1843: unmoins = gsub(gun,s);
1844: lambd = gdiv(resi,gmul(s,gsub(s,gun)));
1845: gammas = ggamma(s,prec);
1846: gammas2= ggamma(gmul2n(s,-1),prec);
1847: gar = gmul(gpowgs(gammas,r2),gpowgs(gammas2,r1));
1848: cs = gexp(gmul(cstlog,s),prec);
1849: var1 = gmul(Pi,s);
1850: gammaunmoins = gdiv(Pi,gmul(gsin(var1,prec),gammas));
1851: gammaunmoins2= gdiv(gmul(gmul(gsqrt(Pi,prec),gpui(gdeux,gsub(s,gun),prec)),
1852: gammas2),
1853: gmul(gcos(gmul2n(var1,-1),prec),gammas));
1854: var1 = var2 = gun;
1855: for (i=2; i<=N0; i++)
1856: if (coef[i])
1857: {
1858: gexpro = gexp(gmul((GEN)coeflog[i],s),bigprec);
1859: var1 = gadd(var1,gmulsg(coef[i], gexpro));
1860: var2 = gadd(var2,gdivsg(coef[i], gmulsg(i,gexpro)));
1861: }
1862: lambd = gadd(lambd,gmul(gmul(var1,cs),gar));
1863: lambd = gadd(lambd,gmul(gmul(gmul(var2,gdiv(cst,cs)),
1864: gpowgs(gammaunmoins,r2)),
1865: gpowgs(gammaunmoins2,r1)));
1866: val = s;
1867: valm = unmoins;
1868: for (i=1; i<=imax; i++)
1869: {
1870: valk = val;
1871: valkm = valm;
1872: for (k=1; k<=ru; k++)
1873: {
1874: p1 = gcoeff(C,i,k);
1875: lambd = gsub(lambd,gdiv(p1,valk )); valk = gmul(val, valk);
1876: lambd = gsub(lambd,gdiv(p1,valkm)); valkm = gmul(valm,valkm);
1877: }
1878: if (r2)
1879: {
1880: val = gadd(val, gun);
1881: valm = gadd(valm,gun);
1882: }
1883: else
1884: {
1885: val = gadd(val, gdeux);
1886: valm = gadd(valm,gdeux); i++;
1887: }
1888: }
1889: }
1890: if (!flag) lambd = gdiv(lambd,gmul(gar,cs)); /* zetak */
1891: return gerepileupto(av, gprec_w(lambd, prec2));
1892: }
1893:
1894: GEN
1895: gzetak(GEN nfz, GEN s, long prec)
1896: {
1897: return gzetakall(nfz,s,0,prec);
1898: }
1899:
1900: GEN
1901: glambdak(GEN nfz, GEN s, long prec)
1902: {
1903: return gzetakall(nfz,s,1,prec);
1904: }
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