Annotation of OpenXM_contrib/pari-2.2/src/basemath/base1.c, Revision 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