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