Annotation of OpenXM_contrib/pari/src/basemath/buch1.c, Revision 1.1
1.1 ! maekawa 1: /*******************************************************************/
! 2: /* */
! 3: /* CLASS GROUP AND REGULATOR (McCURLEY, BUCHMANN) */
! 4: /* QUADRATIC FIELDS */
! 5: /* */
! 6: /*******************************************************************/
! 7: /* $Id: buch1.c,v 1.2 1999/09/23 17:50:56 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: /* See buch2.c:
! 11: * precision en digits decimaux=2*(#digits decimaux de Disc)+50
! 12: * on prendra les p decomposes tels que prod(p)>lim dans la subbase
! 13: * LIMC=Max(c.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))
! 14: * LIMC2=Max(6.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))
! 15: * subbase contient les p decomposes tels que prod(p)>sqrt(Disc)
! 16: * lgsub=subbase[0]=#subbase;
! 17: * subfactorbase est la table des form[p] pour p dans subbase
! 18: * nbram est le nombre de p divisant Disc elimines dans subbase
! 19: * powsubfactorbase est la table des puissances des formes dans subfactorbase
! 20: */
! 21: #define HASHT 1024
! 22: static const long CBUCH = 15; /* of the form 2^k-1 */
! 23: static const long randshift=BITS_IN_RANDOM-1 - 4; /* BITS_IN_RANDOM-1-k */
! 24:
! 25: static long sens,KC,KC2,lgsub,limhash,RELSUP,PRECREG;
! 26: static long *primfact,*primfact1, *exprimfact,*exprimfact1, *badprim;
! 27: static long *factorbase,*numfactorbase, *subbase, *vectbase, **hashtab;
! 28: static GEN **powsubfactorbase,vperm,subfactorbase,Disc,sqrtD,isqrtD;
! 29:
! 30: GEN buchquad(GEN D, double c, double c2, long RELSUP0, long flag, long prec);
! 31: GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus);
! 32:
! 33: GEN
! 34: quadclassunit0(GEN x, long flag, GEN data, long prec)
! 35: {
! 36: long lx,RELSUP0;
! 37: double cbach, cbach2;
! 38:
! 39: if (!data) lx=1;
! 40: else
! 41: {
! 42: lx = lg(data);
! 43: if (typ(data)!=t_VEC || lx > 7)
! 44: err(talker,"incorrect parameters in quadclassunit");
! 45: if (lx > 4) lx = 4;
! 46: }
! 47: cbach = cbach2 = 0.1; RELSUP0 = 5;
! 48: switch(lx)
! 49: {
! 50: case 4: RELSUP0 = itos((GEN)data[3]);
! 51: case 3: cbach2 = gtodouble((GEN)data[2]);
! 52: case 2: cbach = gtodouble((GEN)data[1]);
! 53: }
! 54: return buchquad(x,cbach,cbach2,RELSUP0,flag,prec);
! 55: }
! 56:
! 57: /*******************************************************************/
! 58: /*******************************************************************/
! 59: /* */
! 60: /* Corps de classe de Hilbert et de rayon avec CM (Schertz) */
! 61: /* */
! 62: /*******************************************************************/
! 63: /*******************************************************************/
! 64:
! 65: int
! 66: isoforder2(GEN form)
! 67: {
! 68: GEN a=(GEN)form[1],b=(GEN)form[2],c=(GEN)form[3];
! 69: return !signe(b) || absi_equal(a,b) || egalii(a,c);
! 70: }
! 71:
! 72: /* returns an equation for the Hilbert class field of the imaginary
! 73: * quadratic field of discriminant D if flag=0, a vector of
! 74: * two-component vectors [form,g(form)] where g() is the root of the equation
! 75: * if flag is non-zero.
! 76: */
! 77: static GEN
! 78: quadhilbertimag(GEN D, GEN flag)
! 79: {
! 80: long av=avma,a,b,c,d,dover3,b2,t,h,h2,ell,l,i,i1,i2,ex,prec;
! 81: GEN z,form,L,LG,y,res,ga1,ga2,ga3,ga4,wp,court,p1,p2,qf1,qf2;
! 82: GEN u1,u2,u,w,ag,bg,al,ag2,wlf;
! 83: byteptr p = diffptr;
! 84: int raw = ((typ(flag)==t_INT && signe(flag)));
! 85:
! 86: if (DEBUGLEVEL>=2) timer2();
! 87: if (gcmpgs(D,-11)>=0) return polx[0];
! 88: d=itos(D); L=cgetg(1,t_VEC);
! 89: b2 = b = (d&1)?1:0; h=h2=0; z=gun; dover3=labs(d)/3;
! 90: while (b2 <= dover3)
! 91: {
! 92: t=(b2-d)/4;
! 93: for (a=b?b:1; a*a<=t; a++)
! 94: if (t%a==0)
! 95: {
! 96: h++; c = t/a; z = mulsi(a,z);
! 97: L = concatsp(L, qfi(stoi(a),stoi(b),stoi(c)));
! 98: if (b && a != b && a*a != t)
! 99: {
! 100: h++;
! 101: L = concatsp(L, qfi(stoi(a),stoi(-b),stoi(c)));
! 102: }
! 103: else h2++;
! 104: }
! 105: b+=2; b2=b*b;
! 106: }
! 107: if (h==1) {avma=av; return polx[0];}
! 108: if (DEBUGLEVEL>=2) msgtimer("class number = %ld",h);
! 109: wp=cgetg(1,t_VEC); wlf=cgetg(1,t_VEC); court=stoi(5);
! 110: if (typ(flag)==t_VEC)
! 111: {
! 112: for (i=1; i<lg(flag); i++)
! 113: {
! 114: ell=itos((GEN)flag[i]);
! 115: if (smodis(z,ell) && kross(d,ell) > 0)
! 116: {
! 117: court[2]=ell; form=redimag(primeform(D,court,0));
! 118: if (!gcmp1((GEN)form[1]))
! 119: {
! 120: wp = concat(wp,court); wlf = concat(wlf,form);
! 121: }
! 122: }
! 123: }
! 124: }
! 125: else
! 126: {
! 127: ell=0; ell += *p++; ell+= *p++;
! 128: while (lg(wp)<=2 || ell<=300)
! 129: {
! 130: ell += *p++; if (!*p) err(primer1);
! 131: if (smodis(z,ell) && kross(d,ell) > 0)
! 132: {
! 133: court[2]=ell; form=redimag(primeform(D,court,0));
! 134: if (!gcmp1((GEN)form[1]))
! 135: {
! 136: wp = concat(wp,court); wlf = concat(wlf,form);
! 137: }
! 138: }
! 139: }
! 140: }
! 141: l = lg(wp)-1;
! 142: if (l<2) { avma=av; return gzero; }
! 143: if (typ(flag)==t_VEC) { i1=1; i2=2; p1=(GEN)wp[1]; }
! 144: else
! 145: {
! 146: for(i=1; i<=l; i++)
! 147: if (smodis((GEN)wp[i],3) == 1) break;
! 148: i1=(i>l)?1:i; p1=(GEN)wp[i1]; form=(GEN)wlf[i1];
! 149: if (isoforder2(form))
! 150: {
! 151: if (smodis(p1,4)==3)
! 152: {
! 153: for (i=1; i<=l && (smodis((GEN)wp[i],4) == 3 ||
! 154: (isoforder2((GEN)wlf[i]) && !gegal((GEN)wlf[i],form))); i++);
! 155: if (i>l)
! 156: {
! 157: for (i=1; i<=l && isoforder2((GEN)wlf[i]) && !gegal((GEN)wlf[i],form) ;i++);
! 158: if (i>l) { avma=av; return gzero; }
! 159: }
! 160: }
! 161: else
! 162: {
! 163: for (i=1; i<=l && isoforder2((GEN)wlf[i]) && !gegal((GEN)wlf[i],form); i++);
! 164: if (i>l) { avma=av; return gzero; }
! 165: }
! 166: }
! 167: else
! 168: {
! 169: if (smodis(p1,4)==3)
! 170: {
! 171: for(i=1; i<=l; i++)
! 172: if (smodis((GEN)wp[i],4) == 1) break;
! 173: if (i>l) i=1;
! 174: }
! 175: else i=1;
! 176: }
! 177: i2=i;
! 178: }
! 179: qf1 = primeform(D,p1,0); u1 = gmodulcp((GEN)qf1[2],shifti(p1,1));
! 180: p2 = (GEN)wp[i2];
! 181: qf2 = primeform(D,p2,0); u2 = gmodulcp((GEN)qf2[2],shifti(p2,1));
! 182: ex=24/itos(ggcd(mulii(subis(p1,1),subis(p2,1)),stoi(24)));
! 183: if(DEBUGLEVEL>=2)
! 184: fprintferr("p1 = %Z, p2 = %Z, ex = %ld\n",p1,p2,ex);
! 185: if (!egalii(p1,p2)) u=lift(chinois(u1,u2));
! 186: else
! 187: {
! 188: if (!gegal(qf1,qf2)) err(bugparier,"quadhilbertimag (p1=p1, qf1!=qf2)");
! 189: u=(GEN)compimagraw(qf1,qf2)[2];
! 190: }
! 191: u = gmodulcp(u, shifti(mulii(p1,p2),1));
! 192: LG = cgetg(h+1,t_VEC);
! 193: prec = raw? DEFAULTPREC: 3;
! 194: for(;;)
! 195: {
! 196: long av0 = avma, e, emax = 0;
! 197: GEN lead, sqd = gsqrt(negi(D),prec);
! 198: for (i=1; i<=h; i++)
! 199: {
! 200: form=(GEN)L[i];
! 201: ag=(GEN)form[1]; ag2=shifti(ag,1);
! 202: bg=(GEN)form[2];
! 203: w = lift(chinois(gmodulcp(negi(bg),ag2), u));
! 204: al=cgetg(3,t_COMPLEX);
! 205: al[1]=lneg(gdiv(w,ag2));
! 206: al[2]=ldiv(sqd,ag2);
! 207: ga1 = trueeta(gdiv(al,p1),prec);
! 208: ga2 = trueeta(gdiv(al,p2),prec);
! 209: ga3 = trueeta(gdiv(al,mulii(p1,p2)),prec);
! 210: ga4 = trueeta(al,prec);
! 211: LG[i] = lpowgs(gdiv(gmul(ga1,ga2),gmul(ga3,ga4)), ex);
! 212: e = gexpo((GEN)LG[i]); if (e > 0) emax += e;
! 213: if (DEBUGLEVEL > 2) fprintferr("%ld ",i);
! 214: }
! 215: if (DEBUGLEVEL > 2) fprintferr("\n");
! 216: if (raw)
! 217: {
! 218: y=cgetg(h+1,t_VEC);
! 219: for(i=1; i<=h; i++)
! 220: {
! 221: res=cgetg(3,t_VEC); y[i]=(long)res;
! 222: res[1]=lcopy((GEN)L[i]);
! 223: res[2]=lcopy((GEN)LG[i]);
! 224: }
! 225: break;
! 226: }
! 227: if (DEBUGLEVEL>=2) msgtimer("roots");
! 228: /* to avoid integer overflow (1 + 0.) */
! 229: lead = (emax < bit_accuracy(prec))? gun: realun(prec);
! 230:
! 231: y = greal(roots_to_pol_intern(lead,LG,0,0));
! 232: y = grndtoi(y,&emax);
! 233: if (DEBUGLEVEL>=2) msgtimer("product, error bits = %ld",emax);
! 234: if (emax <= -10)
! 235: {
! 236: if (typ(flag)==t_VEC)
! 237: {
! 238: long av1 = avma;
! 239: e = degree(modulargcd(y,derivpol(y)));
! 240: avma = av1; if (e > 0) { avma=av; return gzero; }
! 241: }
! 242: break;
! 243: }
! 244: avma = av0; prec += (DEFAULTPREC-2) + (1 + (emax >> TWOPOTBITS_IN_LONG));
! 245: if (DEBUGLEVEL) err(warnprec,"quadhilbertimag",prec);
! 246: }
! 247: return gerepileupto(av,y);
! 248: }
! 249:
! 250: GEN quadhilbertreal(GEN D, long prec);
! 251:
! 252: GEN
! 253: quadhilbert(GEN D, GEN flag, long prec)
! 254: {
! 255: if (typ(D)!=t_INT)
! 256: {
! 257: D = checkbnf(D);
! 258: if (degree(gmael(D,7,1))!=2)
! 259: err(talker,"not a polynomial of degree 2 in quadhilbert");
! 260: D=gmael(D,7,3);
! 261: }
! 262: else
! 263: {
! 264: if (!isfundamental(D))
! 265: err(talker,"quadhilbert needs a fundamental discriminant");
! 266: }
! 267: if (signe(D)>0) return quadhilbertreal(D,prec);
! 268: if (!flag) flag = gzero;
! 269: return quadhilbertimag(D,flag);
! 270: }
! 271:
! 272: /* AUXILLIARY ROUTINES FOR QUADRAYIMAGWEI */
! 273:
! 274: static GEN
! 275: getal(GEN nf, GEN b, GEN a, long prec)
! 276: {
! 277: GEN p1,D,os;
! 278:
! 279: p1=idealcoprime(nf,idealdiv(nf,b,a),b);
! 280: D=(GEN)nf[3];
! 281: os=mpodd(D) ? gun : gzero; os=gmul2n(gadd(os,gsqrt(D,prec)),-1);
! 282: return gadd((GEN)p1[1],gmul((GEN)p1[2],os));
! 283: }
! 284:
! 285: static GEN
! 286: epseta(GEN D, long p, long q, GEN tau, long prec)
! 287: {
! 288: GEN p1,p2;
! 289:
! 290: if (gcmpgs(D,-4)<0)
! 291: {
! 292: p1=trueeta(gdivgs(tau,p),prec);
! 293: p2=(p==q) ? p1 : trueeta(gdivgs(tau,q),prec);
! 294: p1=gdiv(gmul(p1,p2),trueeta(gdivgs(tau,p*q),prec));
! 295: return gmul(p1,gpuigs(trueeta(tau,prec),3));
! 296: }
! 297: else return gpuigs(trueeta(tau,prec),4);
! 298: }
! 299:
! 300: static GEN
! 301: pppfun(GEN D, GEN z, GEN a, GEN den, long prec)
! 302: {
! 303: GEN om, res, y, os;
! 304:
! 305: os=mpodd(D) ? gun : gzero; os=gmul2n(gadd(os,gsqrt(D,prec)),-1);
! 306: om=cgetg(3,t_VEC);
! 307: om[1]=ladd(gcoeff(a,1,2),gmul(gcoeff(a,2,2),os));
! 308: om[2]=coeff(a,1,1);
! 309: res=gdiv(ellwp0(om,z,0,prec,0),den); y=res;
! 310: if (gcmpgs(D,-4)>=0) y=gmul(y,res);
! 311: if (gcmpgs(D,-3)==0) y=gmul(y,res);
! 312: return y;
! 313: }
! 314:
! 315: static GEN
! 316: schertzc(GEN nf, GEN a, GEN den, long prec)
! 317: {
! 318: GEN D,al,id,p1;
! 319: long k2,k3,ell,j;
! 320: byteptr p = diffptr;
! 321:
! 322: D=(GEN)nf[3];
! 323:
! 324: if (gcmpgs(D,-3)==0)
! 325: {
! 326: id=idealmul(nf,gdeux,(GEN)(primedec(nf,stoi(3))[1]));
! 327: al=getal(nf,id,a,prec);
! 328: return pppfun(D,al,a,den,prec);
! 329: }
! 330: if (gcmpgs(D,-4)==0)
! 331: {
! 332: id=idealmul(nf,(GEN)(primedec(nf,gdeux)[1]),(GEN)(primedec(nf,stoi(5))[1]));
! 333: al=getal(nf,id,a,prec);
! 334: return pppfun(D,al,a,den,prec);
! 335: }
! 336: k2=krogs(D,2); k3=krogs(D,3);
! 337: if (k2==-1)
! 338: {
! 339: if (k3==-1) return gzero;
! 340: else
! 341: {
! 342: ell=0; ell += *p++; ell += *p++;
! 343: do
! 344: {
! 345: ell += *p++;
! 346: if (!*p) err(primer1);
! 347: }
! 348: while (ell%3!=2 || krogs(D,ell)==1);
! 349: al=getal(nf,idealhermite(nf,(GEN)(primedec(nf,stoi(ell))[1])),a,prec);
! 350: p1=gzero;
! 351: for (j=1; j<=((ell-1)>>1); j++)
! 352: p1=gadd(p1,pppfun(D,gmulsg(j,al),a,den,prec));
! 353: return gneg(p1);
! 354: }
! 355: }
! 356: else
! 357: {
! 358: if (k3!=-1)
! 359: {
! 360: id=idealmul(nf,(GEN)(primedec(nf,gdeux)[1]),(GEN)(primedec(nf,stoi(3))[1]));
! 361: al=getal(nf,id,a,prec);
! 362: return pppfun(D,al,a,den,prec);
! 363: }
! 364: else
! 365: {
! 366: if (k2==1)
! 367: {
! 368: al=getal(nf,idealhermite(nf,gdeux),a,prec);
! 369: return pppfun(D,al,a,den,prec);
! 370: }
! 371: else
! 372: {
! 373: ell=0; ell += *p++; ell += *p++;
! 374: do
! 375: {
! 376: ell += *p++;
! 377: if (!*p) err(primer1);
! 378: }
! 379: while (ell%4!=3 || krogs(D,ell)==1);
! 380: al=getal(nf,idealhermite(nf,(GEN)(primedec(nf,stoi(ell))[1])),a,prec);
! 381: p1=gzero;
! 382: for (j=1; j<=((ell-1)>>1); j++)
! 383: p1=gadd(p1,pppfun(D,gmulsg(j,al),a,den,prec));
! 384: return gmulsg(-krogs(gdeux,ell),p1);
! 385: }
! 386: }
! 387: }
! 388: }
! 389:
! 390: static GEN
! 391: getallelts(GEN nf, GEN clgp)
! 392: {
! 393: GEN cyc,gen,listl,res;
! 394: long lc,i,j,k,no,k1,pro;
! 395:
! 396: cyc=(GEN)clgp[2]; gen=(GEN)clgp[3]; lc=lg(cyc)-1;
! 397: no=itos((GEN)clgp[1]);
! 398: listl=cgetg(no+1,t_VEC);
! 399: listl[1] = (long)idealhermite(nf,gun);
! 400: for (j=1; j<no; j++)
! 401: {
! 402: k = j; res = gun;
! 403: for (i=lc; k; i--)
! 404: {
! 405: pro=((GEN)cyc[i])[2]; /* attention: 1er et seul mot no code de cyc[i] */
! 406: k1 = k%pro;
! 407: if (k1) res=idealmul(nf,res,idealpows(nf,(GEN)gen[i],k1));
! 408: k /= pro;
! 409: }
! 410: listl[j+1] = (long)res;
! 411: }
! 412: return listl;
! 413: }
! 414:
! 415: /* If error and flag = 0 return error message, otherwise return empty vector */
! 416: static GEN
! 417: findbezk(GEN nf, GEN be, long flag, long prec)
! 418: {
! 419: GEN a0,b0,a,b,D,pol,y;
! 420: GEN eps=gpuigs(stoi(10),-8);
! 421: long d0,ea,eb;
! 422:
! 423: D=(GEN)nf[3]; pol=(GEN)nf[1];
! 424: a0=gmul2n(greal(be),1); a=grndtoi(a0,&ea);
! 425: b0=gdiv(gmul2n(gimag(be),1),gsqrt(negi(D),prec)); b=grndtoi(b0,&eb);
! 426: if (ea>-10 || eb>-10)
! 427: {
! 428: if (flag) return cgetg(1,t_VEC);
! 429: else err(talker,"insufficient precision in findbezk");
! 430: }
! 431: if (gcmp(gadd(gabs(gsub(a,a0),prec),gabs(gsub(b,b0),prec)),eps)>0 || mpodd(addii(a,mulii(b,D))))
! 432: {
! 433: if (flag) return cgetg(1,t_VEC);
! 434: else {outerr(be); err(talker," is not in ZK");}
! 435: }
! 436: y=cgetg(3,t_POLMOD);
! 437: y[1]=(long)pol;
! 438: d0=mpodd(D) ? -1 : 0;
! 439: y[2]=ladd(gmul(b,polx[varn(pol)]),gmul2n(gadd(a,gmulgs(b,d0)),-1));
! 440: return y;
! 441: }
! 442:
! 443: /* returns an equation for the ray class field of modulus f of the imaginary
! 444: * quadratic field bnf if flag=0, a vector of
! 445: * two-component vectors [id,g(id)] where g() is the root of the equation
! 446: * if flag is non-zero.
! 447: */
! 448: static GEN
! 449: quadrayimagwei(GEN bnr, GEN flag, long prec)
! 450: {
! 451: long av=avma,tetpil,vpol,clno,clrayno,lc,i,j,res,ell,inda,fl;
! 452: byteptr p = diffptr;
! 453: GEN allf,f,clray,bnf,nf,D,pol,fa,P2,P2new,pp,pi,pial,os,clgp,cyc,gen,listl;
! 454: GEN listray,lista,listla,pp1,la,z,pr2,listden,listc,p1,pii2,ida,ap2;
! 455: GEN om1,om2,tau,d,al,s,vpro;
! 456:
! 457: allf=conductor(bnr,gzero,1,prec);
! 458: f=gmael(allf,1,1); clray=(GEN)allf[2];
! 459: bnf=(GEN)bnr[1]; nf=(GEN)bnf[7]; D=(GEN)nf[3];
! 460: pol=(GEN)nf[1]; vpol=varn(pol);
! 461: fa=(GEN)idealfactor(nf,f)[1];
! 462: fl=itos(flag);
! 463: if (lg(fa)==1)
! 464: {
! 465: P2=quadhilbertimag(D,flag);
! 466: if (fl)
! 467: {
! 468: /* convertir les formes en ideaux */
! 469: }
! 470: tetpil=avma; return gerepile(av,tetpil,gcopy(P2));
! 471: }
! 472: os=mpodd(D) ? gun : gzero; os=gmul2n(gadd(os,gsqrt(D,prec)),-1);
! 473: if (lg(fa)==2)
! 474: {
! 475: pp=(GEN)fa[1]; pi=(GEN)pp[1];
! 476: if (fl) pial=gadd(gmul(gmael(pp,2,2),os),gmael(pp,2,1));
! 477: else pial=basistoalg(nf,(GEN)pp[2]);
! 478: }
! 479: else { pi=gun; pial=gun; }
! 480: clgp=gmael(bnf,8,1);
! 481: clno=itos((GEN)clgp[1]); cyc=(GEN)clgp[2]; gen=(GEN)clgp[3];
! 482: lc=lg(gen);
! 483: listl = getallelts(nf,clgp);
! 484: listray = getallelts(nf,clray);
! 485: clrayno=itos((GEN)clray[1]);
! 486: lista = cgetg(clrayno+1,t_VECSMALL);
! 487: listla = cgetg(clrayno+1,t_VEC);
! 488: for (i=1; i<=clrayno; i++)
! 489: {
! 490: pp = isprincipalgenforce(bnf,idealinv(nf,(GEN)listray[i]));
! 491: pp1 = (GEN)pp[1];
! 492: for (res=0,j=1; j<lc; j++)
! 493: res = res*itos((GEN)cyc[j]) + itos((GEN)pp1[j]);
! 494: lista[i] = res+1;
! 495: la = gmul((GEN)nf[7],(GEN)pp[2]);
! 496: listla[i] = lsubst(la,vpol,os);
! 497: }
! 498: z = dethnf(gmael3(bnr,2,1,1));
! 499: for (i=1; i<=clno; i++) z=mulii(z,dethnf((GEN)listl[i]));
! 500: if (gcmpgs(D,-4)<0)
! 501: {
! 502: GEN court=cgeti(3);
! 503:
! 504: court[1]=evallgefint(3) | evalsigne(1);
! 505: ell=0;
! 506: do
! 507: {
! 508: ell += *p++; court[2]=ell;
! 509: if (!*p) err(primer1);
! 510: }
! 511: while (ell%12!=11 || !gcmp1(ggcd(court,z)) || krogs(D,ell)!=1);
! 512: pr2=idealpows(nf,(GEN)(primedec(nf,stoi(ell))[1]),2);
! 513: }
! 514: else { pr2=idmat(2); ell=1; }
! 515: listden=cgetg(clno+1,t_VEC); listc=cgetg(clno+1,t_VEC);
! 516: p1=mppi(prec); setexpo(p1,2);
! 517: pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
! 518: for (i=1; i<=clno; i++)
! 519: {
! 520: ida=(GEN)listl[i]; ap2=idealmul(nf,ida,pr2);
! 521: om2=gcoeff(ida,1,1);
! 522: om1=gadd(gcoeff(ap2,1,2),gmul(gcoeff(ap2,2,2),os));
! 523: tau=gdiv(om1,om2);
! 524: d=gmul(gsqr(gdiv(pii2,om2)),epseta(D,ell,ell,tau,prec));
! 525: listden[i]=(long)d;
! 526: listc[i]=(long)schertzc(nf,(GEN)listl[i],d,prec);
! 527: }
! 528: al = gsubst(gmul((GEN)nf[7],idealcoprime(nf,f,f)),vpol,os);
! 529: P2 = fl? cgetg(clrayno+1,t_VEC): gun;
! 530: for (j=1; j<=clrayno; j++)
! 531: {
! 532: inda=lista[j];
! 533: s=pppfun(D,gdiv(al,(GEN)listla[j]),(GEN)listl[inda],(GEN)listden[inda],prec);
! 534: s=gsub(s,(GEN)listc[inda]);
! 535:
! 536: if (fl)
! 537: {
! 538: s=gmul(pial,s);
! 539: vpro=cgetg(3,t_VEC);
! 540: vpro[1]=(long)listray[j];
! 541: vpro[2]=(long)s;
! 542: P2[j]=(long)vpro;
! 543: }
! 544: else P2=gmul(P2,gsub(polx[0],gmul(pi,s)));
! 545: }
! 546: if (DEBUGLEVEL)
! 547: {
! 548: fprintferr("P2 = "); outerr(P2);
! 549: }
! 550: if (!fl)
! 551: {
! 552: P2new=gzero;
! 553: for (i=clrayno; i>=0; i--)
! 554: {
! 555: p1=findbezk(nf,truecoeff(P2,i),1,prec);
! 556: if (typ(p1)==t_VEC) {avma=av; return cgetg(1,t_VEC);}
! 557: else P2new=gadd(p1,gmul(polx[0],P2new));
! 558: }
! 559: P2=gsubst(P2new,0,gmul(gdiv(pi,pial),polx[0]));
! 560: P2=gmul(P2,gpuigs(gdiv(pial,pi),clrayno));
! 561: }
! 562: tetpil=avma;
! 563: return gerepile(av,tetpil,gcopy(P2));
! 564: }
! 565:
! 566: /* AUXILLIARY ROUTINES FOR QUADRAYSIGMA */
! 567:
! 568: /* Computes values 2*I*Pi, (om1_*om2-om1*om2_)/(2*I) and
! 569: om1_*eta2-om2_*eta1 necessary for ellphist */
! 570:
! 571: static GEN
! 572: ellphistinit(GEN om, long prec)
! 573: {
! 574: GEN p1,p2,et,om1,om2,ar,pii2,res;
! 575:
! 576: p1=mppi(prec); setexpo(p1,2);
! 577: pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
! 578: om1=(GEN)om[1]; om2=(GEN)om[2];
! 579: if (gsigne(gimag(gdiv(om1,om2)))<0)
! 580: {
! 581: p1=om1; om1=om2; om2=p1;
! 582: p1=cgetg(3,t_VEC); p1[1]=(long)om1; p1[2]=(long)om2;
! 583: }
! 584: et=elleta(om,prec);
! 585: ar=gimag(gmul(p2=gconj(om1),om2));
! 586: p1=gsub(gmul(p2,(GEN)et[2]),gmul(gconj(om2),(GEN)et[1]));
! 587: res=cgetg(4,t_VEC);
! 588: res[1]=(long)pii2; res[2]=(long)ar; res[3]=(long)p1;
! 589: return res;
! 590: }
! 591:
! 592: /* Computes log(phi^*(z,om)), using res computed by ellphistinit */
! 593:
! 594: static GEN
! 595: ellphist(GEN om, GEN res, GEN z, long prec)
! 596: {
! 597: GEN zst;
! 598:
! 599: zst=gdiv(gsub(gmul(z,(GEN)res[3]),gmul(gconj(z),(GEN)res[1])),gmul2n(gmul(gi,(GEN)res[2]),1));
! 600: return gsub(ellsigma(om,z,1,prec),gmul2n(gmul(z,zst),-1));
! 601: }
! 602:
! 603: /* Computes phi^*(la,om)/phi^*(1,om) where om is an oriented basis of the
! 604: ideal gf*gc^{-1} */
! 605:
! 606: static GEN
! 607: computeth2(GEN nf, GEN gf, GEN gc, GEN la, long prec)
! 608: {
! 609: GEN D,os,p1,p2,fdiv,omdiv,lanum,res;
! 610:
! 611: D=(GEN)nf[3];
! 612: os=mpodd(D) ? gun : gzero; os=gmul2n(gadd(os,gsqrt(D,prec)),-1);
! 613: fdiv=idealdiv(nf,gf,gc);
! 614: omdiv=cgetg(3,t_VEC); omdiv[2]=coeff(fdiv,1,1);
! 615: omdiv[1]=ladd(gmul(gcoeff(fdiv,2,2),os),gcoeff(fdiv,1,2));
! 616: la=lift(la);
! 617: lanum=gadd(gmul(truecoeff(la,1),os),truecoeff(la,0));
! 618: res=ellphistinit(omdiv,prec);
! 619: p1=gsub(ellphist(omdiv,res,lanum,prec),ellphist(omdiv,res,gun,prec));
! 620: p2=gimag(p1);
! 621: if (gexpo(greal(p1))>20 || gexpo(p2)> bit_accuracy(min(prec,lg(p2)))-10)
! 622: return cgetg(1,t_VEC);
! 623: else return gexp(p1,prec);
! 624: }
! 625:
! 626: /* Computes P_2(X)=polynomial in Z_K[X] closest to prod_gc(X-th2(gc)) where
! 627: the product is over the ray class group bnr.*/
! 628:
! 629: static GEN
! 630: computeP2(GEN bnr, GEN la, GEN flag, long prec)
! 631: {
! 632: long av=avma,tetpil,clrayno,j,fl;
! 633: GEN bnf,listray,nf,P,s,Pnew,gc,vpro,p1,gf;
! 634:
! 635: bnf=(GEN)bnr[1]; nf=(GEN)bnf[7]; gf=gmael3(bnr,2,1,1);
! 636: listray=getallelts(nf,(GEN)bnr[5]);
! 637: clrayno=lg(listray)-1;
! 638: fl=itos(flag);
! 639: if (fl) P=cgetg(clrayno+1,t_VEC);
! 640: else P=gun;
! 641: for (j=1; j<=clrayno; j++)
! 642: {
! 643: gc=(GEN)listray[j];
! 644: s=computeth2(nf,gf,gc,la,prec);
! 645: if (typ(s)==t_VEC) {avma=av; return cgetg(1,t_VEC);}
! 646: if (fl)
! 647: {
! 648: vpro=cgetg(3,t_VEC);
! 649: vpro[1]=(long)listray[j];
! 650: vpro[2]=(long)s;
! 651: P[j]=(long)vpro;
! 652: }
! 653: else P=gmul(P,gsub(polx[0],s));
! 654: }
! 655: if (!fl)
! 656: {
! 657: Pnew=gzero;
! 658: for (j=clrayno; j>=0; j--)
! 659: {
! 660: p1=findbezk(nf,truecoeff(P,j),1,prec);
! 661: if (typ(p1)==t_VEC)
! 662: {
! 663: prec=(prec<<1)-2; avma=av;
! 664: if (DEBUGLEVEL) err(warnprec,"computeP2",prec);
! 665: return computeP2(bnr,la,flag,prec);
! 666: }
! 667: Pnew=gadd(p1,gmul(polx[0],Pnew));
! 668: }
! 669: P=Pnew;
! 670: }
! 671: tetpil=avma; return gerepile(av,tetpil,gcopy(P));
! 672: }
! 673:
! 674: static GEN
! 675: compocyclo(GEN nf, long m, long d, long prec)
! 676: {
! 677: GEN p1,p2,p3,D,res,pol4,nf4;
! 678: long ell,vx,ph;
! 679:
! 680: D=(GEN)nf[3];
! 681: p1=quadhilbertimag(D, gzero);
! 682: p2=cyclo(m,0);
! 683: if (d==1)
! 684: {
! 685: ph=degree(p2);
! 686: p2=gmul(gpuigs(polx[0],ph),gsubst(p2,0,gdiv(polx[MAXVARN],polx[0])));
! 687: return gsubst(subres(p1,p2),MAXVARN,polx[0]);
! 688: }
! 689: ell = m%2 ? m : (m>>2);
! 690: if (!signe(addsi(ell,D)))
! 691: {
! 692: p2=gcoeff(nffactor(nf,p2),1,1);
! 693: ph=degree(p2);
! 694: p2=gmul(gpuigs(polx[0],ph),gsubst(p2,0,gdiv(polx[MAXVARN],polx[0])));
! 695: return gsubst(subres(p1,p2),MAXVARN,polx[0]);
! 696: }
! 697: if (ell%4==3) ell= -ell;
! 698: p3=cgetg(5,t_POL);
! 699: p3[1]=evalsigne(1)|evallgef(5)|evalvarn(0);
! 700: p3[2]=lstoi((1-ell)>>2);
! 701: p3[3]=p3[4]=un;
! 702: res=rnfequation2(nf,p3);
! 703: vx=varn((GEN)nf[1]);
! 704: pol4=gsubst((GEN)res[1],0,polx[vx]);
! 705: nf4=initalg(pol4,prec);
! 706: p1=gcoeff(nffactor(nf4,p1),1,1);
! 707: p2=gcoeff(nffactor(nf4,p2),1,1);
! 708: ph=degree(p2);
! 709: p2=gmul(gpuigs(polx[0],ph),gsubst(p2,0,gdiv(polx[MAXVARN],polx[0])));
! 710: p3=gsubst(subres(p1,p2),MAXVARN,polx[0]);
! 711: p1=gmodulcp(gsubst(lift((GEN)res[2]),0,polx[vx]),pol4);
! 712: return gsubst(lift0(p3,vx),vx,p1);
! 713: }
! 714:
! 715: static GEN
! 716: retflag(GEN x, GEN flag)
! 717: {
! 718: if (itos(flag)) err(impl,"some special cases in quadray (flag=1)");
! 719: /* to be done */
! 720: return x;
! 721: }
! 722:
! 723: /* Treat special cases directly. Exit with 0 if not special case. Internal,
! 724: no stack treatment. */
! 725: static GEN
! 726: treatspecialsigma(GEN nf, GEN gf, GEN flag, long prec)
! 727: {
! 728: GEN D,p1,p2,tryf,fa;
! 729: long Ds,flf,lfa,i;
! 730:
! 731: D=(GEN)nf[3];
! 732: if (cmpis(D,-3)==0)
! 733: {
! 734: p1=idmat(2); p2=gcoeff(gf,1,1);
! 735: if (gegal(gf,gmul(p1,p2)) && (cmpis(p2,4)==0 || cmpis(p2,5)==0 || cmpis(p2,7)==0))
! 736: return retflag(cyclo(itos(p2),0),flag);
! 737: p1=idealpows(nf,(GEN)primedec(nf,stoi(3))[1],3);
! 738: if (gegal(gf,p1))
! 739: {
! 740: p1=gcoeff(nffactor(nf,cyclo(3,0)),1,1);
! 741: p1=gneg(polcoeff0(p1,0,0)); /* should be zeta_3 */
! 742: p2=cgetg(6,t_POL);
! 743: p2[1]=evalsigne(1)|evallgef(6)|evalvarn(0);
! 744: p2[2]=(long)p1; p2[3]=p2[4]=zero; p2[5]=un;
! 745: return retflag(p2,flag);
! 746: }
! 747: return gzero;
! 748: }
! 749: if (cmpis(D,-4)==0)
! 750: {
! 751: p1=idmat(2); p2=gcoeff(gf,1,1);
! 752: if (gegal(gf,gmul(p1,p2)))
! 753: {
! 754: if (cmpis(p2,3)==0 || cmpis(p2,5)==0)
! 755: return retflag(cyclo(itos(p2),0),flag);
! 756: if (cmpis(p2,4)==0)
! 757: {
! 758: p1=gcoeff(nffactor(nf,cyclo(4,0)),1,1);
! 759: p1=gneg(polcoeff0(p1,0,0)); /* should be zeta_4=I */
! 760: p2=cgetg(5,t_POL);
! 761: p2[1]=evalsigne(1)|evallgef(5)|evalvarn(0);
! 762: p2[2]=(long)p1; p2[3]=zero; p2[4]=un;
! 763: return retflag(p2,flag);
! 764: }
! 765: }
! 766: return gzero;
! 767: }
! 768: p1=idmat(2); p2=gcoeff(gf,1,1); Ds=itos(modis(D,48));
! 769: if (gegal(gf,gmul(p1,p2)))
! 770: {
! 771: if (cmpis(p2,2)==0 && Ds%16==8)
! 772: return retflag(compocyclo(nf,4,1,prec),flag);
! 773: if (cmpis(p2,3)==0 && Ds%3==1)
! 774: return retflag(compocyclo(nf,3,1,prec),flag);
! 775: if (cmpis(p2,4)==0 && Ds%8==1)
! 776: return retflag(compocyclo(nf,4,1,prec),flag);
! 777: if (cmpis(p2,6)==0 && Ds%48==40)
! 778: return retflag(compocyclo(nf,12,1,prec),flag);
! 779: return gzero;
! 780: }
! 781: p1=gcoeff(gf,2,2);
! 782: if (gcmp1(p1)) {flf=1; tryf=p2;}
! 783: else
! 784: {
! 785: if (cmpis(p1,2) || mpodd(p2) || mpodd(gcoeff(gf,1,2))) return gzero;
! 786: flf=2; tryf=gmul2n(p2,-1);
! 787: }
! 788: fa=(GEN)factor(D)[1]; lfa=lg(fa);
! 789: for (i=1; i<lfa; i++)
! 790: if (cmpis((GEN)fa[i],3)>0 && gegal((GEN)fa[i],tryf))
! 791: {
! 792: if (flf==1) return retflag(compocyclo(nf,itos(tryf),2,prec),flag);
! 793: if (Ds%16==8) return retflag(compocyclo(nf,4*itos(tryf),2,prec),flag);
! 794: return gzero;
! 795: }
! 796: return gzero;
! 797: }
! 798:
! 799: /* Compute ray class field polynomial using sigma; if flag=1, pairs
! 800: [ideals,roots] are given instead so that congruence subgroups can be used */
! 801:
! 802: static GEN
! 803: quadrayimagsigma(GEN bnr, GEN flag, long prec)
! 804: {
! 805: long av=avma,tetpil,a,b,f2;
! 806: GEN allf,bnf,nf,pol,w,wbas,gf,la,p1,p2,y,labas,gfi;
! 807:
! 808: allf=conductor(bnr,gzero,2,prec);
! 809: bnr=(GEN)allf[2]; gf=gmael(allf,1,1);
! 810: if (gcmp1(dethnf(gf)))
! 811: {
! 812: if (typ(flag)!=t_INT) flag=(GEN)flag[2];
! 813: p1=quadhilbertimag(gmael3(bnr,1,7,3),flag);
! 814: if (itos(flag))
! 815: {
! 816: /* convertir les formes en ideaux */
! 817: }
! 818: tetpil=avma; return gerepile(av,tetpil,gcopy(p1));
! 819: }
! 820: bnf=(GEN)bnr[1]; nf=(GEN)bnf[7]; pol=(GEN)nf[1];
! 821: p1=treatspecialsigma(nf,gf,flag,prec);
! 822: if (!gcmp0(p1)) {tetpil=avma; return gerepile(av,tetpil,gcopy(p1));}
! 823: w=gmodulcp(polx[varn(pol)],pol);
! 824: wbas=algtobasis(nf,w);
! 825: f2=itos(gmul2n(gcoeff(gf,1,1),1));
! 826: gfi=invmat(gf);
! 827: for (a=0; a<f2; a++)
! 828: {
! 829: for (b=0; b<f2; b++)
! 830: {
! 831: if (DEBUGLEVEL>=2) fprintferr("[%ld,%ld] ",a,b);
! 832: la=gaddgs(gmulsg(a,w),b);
! 833: p1=gnorm(la);
! 834: if (gcmp1(modis(p1,f2)))
! 835: {
! 836: labas=gadd(gmulsg(a,wbas),algtobasis(nf,stoi(b)));
! 837: if (gcmp1(denom(gmul(gfi,gadd(labas,algtobasis(nf,stoi(-1))))))) continue;
! 838: if (gcmp1(denom(gmul(gfi,gadd(labas,algtobasis(nf,gun)))))) continue;
! 839: if (!cmpis((GEN)nf[3],-4))
! 840: {
! 841: p1=gcoeff(nffactor(nf,cyclo(4,0)),1,1);
! 842: p1=algtobasis(nf,polcoeff0(p1,0,0)); /* should be -I */
! 843: if (gcmp1(denom(gmul(gfi,gadd(labas,p1))))) continue;
! 844: if (gcmp1(denom(gmul(gfi,gsub(labas,p1))))) continue;
! 845: }
! 846: if (!cmpis((GEN)nf[3],-3))
! 847: {
! 848: p1=(GEN)nffactor(nf,cyclo(3,0))[1];
! 849: p2=algtobasis(nf,polcoeff0((GEN)p1[1],0,0)); /* -zeta_3^2 */
! 850: p1=algtobasis(nf,polcoeff0((GEN)p1[2],0,0)); /* should be -zeta_3 */
! 851: if (gcmp1(denom(gmul(gfi,gadd(labas,p1))))) continue;
! 852: if (gcmp1(denom(gmul(gfi,gsub(labas,p1))))) continue;
! 853: if (gcmp1(denom(gmul(gfi,gadd(labas,p2))))) continue;
! 854: if (gcmp1(denom(gmul(gfi,gsub(labas,p2))))) continue;
! 855: }
! 856: if (DEBUGLEVEL)
! 857: {
! 858: if (DEBUGLEVEL>=2) fprintferr("\n");
! 859: fprintferr("lambda = ");
! 860: outerr(la);
! 861: }
! 862: tetpil=avma;
! 863: y=computeP2(bnr,la,flag,prec);
! 864: return gerepile(av,tetpil,y);
! 865: }
! 866: }
! 867: }
! 868: err(talker,"bug in quadrayimagsigma, please report");
! 869: return gzero;
! 870: }
! 871:
! 872: GEN
! 873: quadray(GEN D, GEN f, GEN flag, long prec)
! 874: {
! 875: long av=avma,tetpil;
! 876: GEN bnr,y,p1,pol,bnf,flagnew;
! 877:
! 878: if (typ(D)!=t_INT)
! 879: {
! 880: bnf = checkbnf(D);
! 881: if (degree(gmael(bnf,7,1))!=2)
! 882: err(talker,"not a polynomial of degree 2 in quadray");
! 883: D=gmael(bnf,7,3);
! 884: }
! 885: else
! 886: {
! 887: if (!isfundamental(D))
! 888: err(talker,"quadray needs a fundamental discriminant");
! 889: pol=quadpoly(D); setvarn(pol, fetch_user_var("y"));
! 890: bnf=bnfinit0(pol,0,NULL,prec);
! 891: }
! 892: bnr=bnrinit0(bnf,f,1,prec);
! 893: if (gcmp1(gmael(bnr,5,1)))
! 894: {
! 895: avma=av; if (!flag || gcmp0(flag)) return polx[0];
! 896: y=cgetg(2,t_VEC); p1=cgetg(3,t_VEC); y[1]=(long)p1;
! 897: p1[1]=(long)idmat(2); p1[2]=(long)polx[0];
! 898: return y;
! 899: }
! 900: tetpil=avma;
! 901: if (signe(D)>0)
! 902: {
! 903: if (flag && !gcmp0(flag)) err(warner,"ignoring flag in quadray");
! 904: y=bnrstark(bnr,gzero,1,prec);
! 905: }
! 906: else
! 907: {
! 908: if (!flag) flag = gzero;
! 909: flagnew=flag;
! 910: if (typ(flagnew)==t_INT)
! 911: {
! 912: flagnew=absi(flagnew);
! 913: if (cmpis(flagnew,1)<=0) y=quadrayimagsigma(bnr,flagnew,prec);
! 914: else y=quadrayimagwei(bnr,mpodd(flagnew) ? gun : gzero,prec);
! 915: }
! 916: else
! 917: {
! 918: if (typ(flagnew)!=t_VEC || lg(flagnew)<=2) err(flagerr,"quadray");
! 919: y=computeP2(bnr,(GEN)flagnew[1],(GEN)flagnew[2],prec);
! 920: }
! 921: if (typ(y)==t_VEC && lg(y)==1)
! 922: {
! 923: prec=(prec<<1)-2; avma=av;
! 924: if (DEBUGLEVEL) err(warnprec,"quadray",prec);
! 925: return quadray(D,f,flag,prec);
! 926: }
! 927: }
! 928: return gerepile(av,tetpil,y);
! 929: }
! 930:
! 931: /*******************************************************************/
! 932: /* */
! 933: /* Routines related to binary quadratic forms (for internal use) */
! 934: /* */
! 935: /*******************************************************************/
! 936:
! 937: static void
! 938: rhoreal_aux2(GEN x, GEN y)
! 939: {
! 940: GEN p1,p2;
! 941: long s = signe(x[3]);
! 942:
! 943: y[1]=x[3]; setsigne(y[1],1);
! 944: p2 = (cmpii(isqrtD,(GEN)y[1]) >= 0)? isqrtD: (GEN) y[1];
! 945: p1 = shifti((GEN)y[1],1);
! 946: p2 = divii(addii(p2,(GEN)x[2]), p1);
! 947: y[2] = lsubii(mulii(p2,p1),(GEN)x[2]);
! 948:
! 949: setsigne(y[1],s);
! 950: p1 = shifti(subii(sqri((GEN)y[2]),Disc),-2);
! 951: y[3] = ldivii(p1,(GEN)y[1]);
! 952: }
! 953:
! 954: static GEN
! 955: rhoreal_aux(GEN x)
! 956: {
! 957: GEN y = cgetg(6,t_VEC);
! 958: long e;
! 959:
! 960: rhoreal_aux2(x,y);
! 961: switch(lg(x))
! 962: {
! 963: case 4: case 5: setlg(y,4); break;
! 964: case 6:
! 965: y[5]=lmulrr(divrr(addir((GEN)x[2],sqrtD),subir((GEN)x[2],sqrtD)),
! 966: (GEN)x[5]);
! 967: e = expo(y[5]);
! 968: if (e < EXP220) y[4]=x[4];
! 969: else
! 970: {
! 971: y[4]=laddsi(1,(GEN)x[4]);
! 972: setexpo(y[5], e - EXP220);
! 973: }
! 974: }
! 975: return y;
! 976: }
! 977:
! 978: static GEN
! 979: rhorealform(GEN x)
! 980: {
! 981: long av=avma,tetpil;
! 982: x = rhoreal_aux(x); tetpil=avma;
! 983: return gerepile(av,tetpil,gcopy(x));
! 984: }
! 985:
! 986: static GEN
! 987: redrealform(GEN x)
! 988: {
! 989: long l;
! 990: GEN p1;
! 991:
! 992: for(;;)
! 993: {
! 994: if (signe(x[2]) > 0 && cmpii((GEN)x[2],isqrtD) <= 0)
! 995: {
! 996: p1 = subii(isqrtD, shifti(absi((GEN)x[1]),1));
! 997: l = absi_cmp((GEN)x[2],p1);
! 998: if (l>0 || (l==0 && signe(p1)<0)) break;
! 999: }
! 1000: x = rhoreal_aux(x);
! 1001: }
! 1002: if (signe(x[1]) < 0)
! 1003: {
! 1004: if (sens || (signe(x[3])>0 && !absi_cmp((GEN)x[1],(GEN)x[3])))
! 1005: return rhoreal_aux(x); /* narrow class group, or a = -c */
! 1006: setsigne(x[1],1); setsigne(x[3],-1);
! 1007: }
! 1008: return x;
! 1009: }
! 1010:
! 1011: static GEN
! 1012: redrealform_init(GEN x)
! 1013: {
! 1014: long av=avma, tetpil;
! 1015: GEN y = cgetg(6,t_VEC);
! 1016:
! 1017: y[1]=x[1]; y[2]=x[2]; y[3]=x[3]; y[4]=zero;
! 1018: y[5]=(long)realun(PRECREG);
! 1019: y = redrealform(y); tetpil=avma;
! 1020: return gerepile(av,tetpil,gcopy(y));
! 1021: }
! 1022:
! 1023: static void
! 1024: compreal_aux(GEN x, GEN y, GEN z)
! 1025: {
! 1026: GEN s,n,d,d1,x1,x2,y1,y2,v1,v2,b3,c3,m,p1,r;
! 1027:
! 1028: s=shifti(addii((GEN)x[2],(GEN)y[2]),-1);
! 1029: n=subii((GEN)y[2],s);
! 1030: d=bezout((GEN)y[1],(GEN)x[1],&y1,&x1);
! 1031: d1=bezout(s,d,&x2,&y2);
! 1032: v1=divii((GEN)x[1],d1);
! 1033: v2=divii((GEN)y[1],d1);
! 1034: m=addii(mulii(mulii(y1,y2),n),mulii((GEN)y[3],x2));
! 1035: setsigne(m,-signe(m));
! 1036: r=modii(m,v1); p1=mulii(v2,r); b3=shifti(p1,1);
! 1037: c3=addii(mulii((GEN)y[3],d1),mulii(r,addii((GEN)y[2],p1)));
! 1038:
! 1039: z[1]=lmulii(v1,v2);
! 1040: z[2]=laddii((GEN)y[2],b3);
! 1041: z[3]=ldivii(c3,v1);
! 1042: }
! 1043:
! 1044: static GEN
! 1045: comprealform3(GEN x, GEN y)
! 1046: {
! 1047: long av = avma, tetpil;
! 1048: GEN z = cgetg(4,t_VEC);
! 1049: compreal_aux(x,y,z); z=redrealform(z); tetpil=avma;
! 1050: return gerepile(av,tetpil,gcopy(z));
! 1051: }
! 1052:
! 1053: static GEN
! 1054: comprealform5(GEN x, GEN y)
! 1055: {
! 1056: long av = avma,tetpil,e;
! 1057: GEN p1, z = cgetg(6,t_VEC);
! 1058:
! 1059: compreal_aux(x,y,z);
! 1060: z[5]=lmulrr((GEN)x[5],(GEN)y[5]);
! 1061: e=expo(z[5]); p1 = addii((GEN)x[4],(GEN)y[4]);
! 1062: if (e < EXP220) z[4] = (long)p1;
! 1063: else
! 1064: {
! 1065: z[4] = laddsi(1,p1);
! 1066: setexpo(z[5], e-EXP220);
! 1067: }
! 1068: z=redrealform(z); tetpil=avma;
! 1069: return gerepile(av,tetpil,gcopy(z));
! 1070: }
! 1071:
! 1072: static GEN
! 1073: initializeform5(long *ex)
! 1074: {
! 1075: long av = avma, i;
! 1076: GEN form = powsubfactorbase[1][ex[1]];
! 1077:
! 1078: for (i=2; i<=lgsub; i++)
! 1079: form = comprealform5(form, powsubfactorbase[i][ex[i]]);
! 1080: i=avma; return gerepile(av,i,gcopy(form));
! 1081: }
! 1082:
! 1083: /*******************************************************************/
! 1084: /* */
! 1085: /* Common subroutines */
! 1086: /* */
! 1087: /*******************************************************************/
! 1088: static void
! 1089: buch_init(void)
! 1090: {
! 1091: if (DEBUGLEVEL) timer2();
! 1092: primfact = new_chunk(100);
! 1093: primfact1 = new_chunk(100);
! 1094: exprimfact = new_chunk(100);
! 1095: exprimfact1 = new_chunk(100);
! 1096: badprim = new_chunk(100);
! 1097: hashtab = (long**) new_chunk(HASHT);
! 1098: }
! 1099:
! 1100: double
! 1101: check_bach(double cbach, double B)
! 1102: {
! 1103: if (cbach > B)
! 1104: err(talker,"sorry, buchxxx couldn't deal with this field PLEASE REPORT!");
! 1105: cbach *= 2; if (cbach > B) cbach = B;
! 1106: if (DEBUGLEVEL) fprintferr("\n*** Bach constant: %f\n", cbach);
! 1107: return cbach;
! 1108: }
! 1109:
! 1110: static long
! 1111: factorisequad(GEN f, long kcz, long limp)
! 1112: {
! 1113: long i,p,k,av,lo;
! 1114: GEN q,r, x = (GEN)f[1];
! 1115:
! 1116: if (is_pm1(x)) { primfact[0]=0; return 1; }
! 1117: av=avma; lo=0;
! 1118: if (signe(x) < 0) x = absi(x);
! 1119: for (i=1; ; i++)
! 1120: {
! 1121: p=factorbase[i]; q=dvmdis(x,p,&r);
! 1122: if (!signe(r))
! 1123: {
! 1124: k=0; while (!signe(r)) { x=q; k++; q=dvmdis(x,p,&r); }
! 1125: lo++; primfact[lo]=p; exprimfact[lo]=k;
! 1126: }
! 1127: if (cmpis(q,p)<=0) break;
! 1128: if (i==kcz) { avma=av; return 0; }
! 1129: }
! 1130: p = x[2]; avma=av;
! 1131: /* p = itos(x) if lgefint(x)=3 */
! 1132: if (lgefint(x)!=3 || p > limhash) return 0;
! 1133:
! 1134: if (p != 1 && p <= limp)
! 1135: {
! 1136: for (i=1; i<=badprim[0]; i++)
! 1137: if (p % badprim[i] == 0) return 0;
! 1138: lo++; primfact[lo]=p; exprimfact[lo]=1;
! 1139: p = 1;
! 1140: }
! 1141: primfact[0]=lo; return p;
! 1142: }
! 1143:
! 1144: static long *
! 1145: largeprime(long q, long *ex, long np, long nrho)
! 1146: {
! 1147: const long hashv = ((q&2047)-1)>>1;
! 1148: long *pt, i;
! 1149:
! 1150: for (pt = hashtab[hashv]; ; pt = (long*) pt[0])
! 1151: {
! 1152: if (!pt)
! 1153: {
! 1154: pt = (long*) gpmalloc((lgsub+4)<<TWOPOTBYTES_IN_LONG);
! 1155: *pt++ = nrho; /* nrho = pt[-3] */
! 1156: *pt++ = np; /* np = pt[-2] */
! 1157: *pt++ = q; /* q = pt[-1] */
! 1158: pt[0] = (long)hashtab[hashv];
! 1159: for (i=1; i<=lgsub; i++) pt[i]=ex[i];
! 1160: hashtab[hashv]=pt; return NULL;
! 1161: }
! 1162: if (pt[-1] == q) break;
! 1163: }
! 1164: for(i=1; i<=lgsub; i++)
! 1165: if (pt[i] != ex[i]) return pt;
! 1166: return (pt[-2]==np)? (GEN)NULL: pt;
! 1167: }
! 1168:
! 1169: static long
! 1170: badmod8(GEN x)
! 1171: {
! 1172: long r = mod8(x);
! 1173: if (!r) return 1;
! 1174: if (signe(Disc) < 0) r = 8-r;
! 1175: return (r < 4);
! 1176: }
! 1177:
! 1178: /* cree factorbase, numfactorbase, vectbase; affecte badprim */
! 1179: static void
! 1180: factorbasequad(GEN Disc, long n2, long n)
! 1181: {
! 1182: long i,p,bad, av = avma;
! 1183: byteptr d=diffptr;
! 1184:
! 1185: numfactorbase = (long*) gpmalloc(sizeof(long)*(n2+1));
! 1186: factorbase = (long*) gpmalloc(sizeof(long)*(n2+1));
! 1187: KC=0; bad=0; i=0; p = *d++;
! 1188: while (p<=n2)
! 1189: {
! 1190: switch(krogs(Disc,p))
! 1191: {
! 1192: case -1: break; /* inert */
! 1193: case 0: /* ramified */
! 1194: {
! 1195: GEN p1 = divis(Disc,p);
! 1196: if (smodis(p1,p) == 0)
! 1197: if (p!=2 || badmod8(p1)) { badprim[++bad]=p; break; }
! 1198: i++; numfactorbase[p]=i; factorbase[i] = -p; break;
! 1199: }
! 1200: default: /* split */
! 1201: i++; numfactorbase[p]=i; factorbase[i] = p;
! 1202: }
! 1203: p += *d++; if (!*d) err(primer1);
! 1204: if (KC == 0 && p>n) KC = i;
! 1205: }
! 1206: if (!KC) { free(factorbase); free(numfactorbase); return; }
! 1207: KC2 = i;
! 1208: vectbase = (long*) gpmalloc(sizeof(long)*(KC2+1));
! 1209: for (i=1; i<=KC2; i++)
! 1210: {
! 1211: p = factorbase[i];
! 1212: vectbase[i]=p; factorbase[i]=labs(p);
! 1213: }
! 1214: if (DEBUGLEVEL)
! 1215: {
! 1216: msgtimer("factor base");
! 1217: if (DEBUGLEVEL>7)
! 1218: {
! 1219: fprintferr("factorbase:\n");
! 1220: for (i=1; i<=KC; i++) fprintferr("%ld ",factorbase[i]);
! 1221: fprintferr("\n"); flusherr();
! 1222: }
! 1223: }
! 1224: avma=av; badprim[0] = bad;
! 1225: }
! 1226:
! 1227: /* cree vectbase and subfactorbase. Affecte lgsub */
! 1228: static long
! 1229: subfactorbasequad(double ll, long KC)
! 1230: {
! 1231: long i,j,k,nbidp,p,pro[100], ss;
! 1232: GEN p1,y;
! 1233: double prod;
! 1234:
! 1235: i=0; ss=0; prod=1;
! 1236: for (j=1; j<=KC && prod<=ll; j++)
! 1237: {
! 1238: p = vectbase[j];
! 1239: if (p>0) { pro[++i]=p; prod*=p; vperm[i]=j; } else ss++;
! 1240: }
! 1241: if (prod<=ll) return -1;
! 1242: nbidp=i;
! 1243: for (k=1; k<j; k++)
! 1244: if (vectbase[k]<=0) vperm[++i]=k;
! 1245:
! 1246: y=cgetg(nbidp+1,t_COL);
! 1247: if (PRECREG) /* real */
! 1248: for (j=1; j<=nbidp; j++)
! 1249: {
! 1250: p1=primeform(Disc,stoi(pro[j]),PRECREG);
! 1251: y[j] = (long) redrealform_init(p1);
! 1252: }
! 1253: else
! 1254: for (j=1; j<=nbidp; j++) /* imaginary */
! 1255: {
! 1256: p1=primeform(Disc,stoi(pro[j]),0);
! 1257: y[j] = (long)p1;
! 1258: }
! 1259: subbase = (long*) gpmalloc(sizeof(long)*(nbidp+1));
! 1260: lgsub = nbidp; for (j=1; j<=lgsub; j++) subbase[j]=pro[j];
! 1261: if (DEBUGLEVEL>7)
! 1262: {
! 1263: fprintferr("subfactorbase: ");
! 1264: for (i=1; i<=lgsub; i++)
! 1265: { fprintferr("%ld: ",subbase[i]); outerr((GEN)y[i]); }
! 1266: fprintferr("\n"); flusherr();
! 1267: }
! 1268: subfactorbase = y; return ss;
! 1269: }
! 1270:
! 1271: static void
! 1272: powsubfact(long n, long a)
! 1273: {
! 1274: GEN unform, **x = (GEN**) gpmalloc(sizeof(GEN*)*(n+1));
! 1275: long i,j;
! 1276:
! 1277: for (i=1; i<=n; i++)
! 1278: x[i] = (GEN*) gpmalloc(sizeof(GEN)*(a+1));
! 1279: if (PRECREG) /* real */
! 1280: {
! 1281: unform=cgetg(6,t_VEC);
! 1282: unform[1]=un;
! 1283: unform[2]=(mod2(Disc) == mod2(isqrtD))? (long)isqrtD: laddsi(-1,isqrtD);
! 1284: unform[3]=lshifti(subii(sqri((GEN)unform[2]),Disc),-2);
! 1285: unform[4]=zero;
! 1286: unform[5]=(long)realun(PRECREG);
! 1287: for (i=1; i<=n; i++)
! 1288: {
! 1289: x[i][0] = unform;
! 1290: for (j=1; j<=a; j++)
! 1291: x[i][j]=comprealform5(x[i][j-1],(GEN)subfactorbase[i]);
! 1292: }
! 1293: }
! 1294: else /* imaginary */
! 1295: {
! 1296: unform=cgetg(4,t_QFI);
! 1297: unform[1]=un;
! 1298: unform[2]=mod2(Disc)? un: zero;
! 1299: unform[3]=lshifti(absi(Disc),-2);
! 1300: for (i=1; i<=n; i++)
! 1301: {
! 1302: x[i][0] = unform;
! 1303: for (j=1; j<=a; j++)
! 1304: x[i][j]=compimag(x[i][j-1],(GEN)subfactorbase[i]);
! 1305: }
! 1306: }
! 1307: if (DEBUGLEVEL) msgtimer("powsubfact");
! 1308: powsubfactorbase = x;
! 1309: }
! 1310:
! 1311: static void
! 1312: desalloc(long **mat)
! 1313: {
! 1314: long i,*p,*q;
! 1315:
! 1316: free(vectbase); free(factorbase); free(numfactorbase);
! 1317: if (mat)
! 1318: {
! 1319: free(subbase);
! 1320: for (i=1; i<lg(subfactorbase); i++) free(powsubfactorbase[i]);
! 1321: for (i=1; i<lg(mat); i++) free(mat[i]);
! 1322: free(mat); free(powsubfactorbase);
! 1323: for (i=1; i<HASHT; i++)
! 1324: for (p = hashtab[i]; p; p = q) { q=(long*)p[0]; free(p-3); }
! 1325: }
! 1326: }
! 1327:
! 1328: /* L-function */
! 1329: static GEN
! 1330: lfunc(GEN Disc)
! 1331: {
! 1332: long av=avma, p;
! 1333: GEN y=realun(DEFAULTPREC);
! 1334: byteptr d=diffptr;
! 1335:
! 1336: for(p = *d++; p<=30000; p += *d++)
! 1337: {
! 1338: if (!*d) err(primer1);
! 1339: y = mulsr(p, divrs(y, p-krogs(Disc,p)));
! 1340: }
! 1341: return gerepileupto(av,y);
! 1342: }
! 1343:
! 1344: #define comp(x,y) x? (PRECREG? compreal(x,y): compimag(x,y)): y
! 1345: static GEN
! 1346: get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec)
! 1347: {
! 1348: GEN p1,p2,res,*init, u1u2=smith2(W), u1=(GEN)u1u2[1], met=(GEN)u1u2[3];
! 1349: long c,i,j, l = lg(met);
! 1350:
! 1351: u1 = reducemodmatrix(ginv(u1), W);
! 1352: for (c=1; c<l; c++)
! 1353: if (gcmp1(gcoeff(met,c,c))) break;
! 1354: if (DEBUGLEVEL) msgtimer("smith/class group");
! 1355: res=cgetg(c,t_VEC); init = (GEN*)cgetg(l,t_VEC);
! 1356: for (i=1; i<l; i++)
! 1357: init[i] = primeform(Disc,stoi(labs(vectbase[vperm[i]])),prec);
! 1358: for (j=1; j<c; j++)
! 1359: {
! 1360: p1 = NULL;
! 1361: for (i=1; i<l; i++)
! 1362: {
! 1363: p2 = gpui(init[i], gcoeff(u1,i,j), prec);
! 1364: p1 = comp(p1,p2);
! 1365: }
! 1366: res[j] = (long)p1;
! 1367: }
! 1368: if (DEBUGLEVEL) msgtimer("generators");
! 1369: *ptmet = met; return res;
! 1370: }
! 1371:
! 1372: static GEN
! 1373: extra_relations(long LIMC, long *ex, long nlze, GEN extramatc)
! 1374: {
! 1375: long av,fpc,p,ep,i,j,k,nlze2, *col, *colg, s = 0, extrarel = nlze+2;
! 1376: long MAXRELSUP = min(50,4*KC);
! 1377: GEN p1,form, extramat = cgetg(extrarel+1,t_MAT);
! 1378:
! 1379: if (DEBUGLEVEL)
! 1380: {
! 1381: fprintferr("looking for %ld extra relations\n",extrarel);
! 1382: flusherr();
! 1383: }
! 1384: for (j=1; j<=extrarel; j++) extramat[j]=lgetg(KC+1,t_COL);
! 1385: nlze2 = PRECREG? max(nlze,lgsub): min(nlze+1,KC);
! 1386: if (nlze2 < 3 && KC > 2) nlze2 = 3;
! 1387: av = avma;
! 1388: while (s<extrarel)
! 1389: {
! 1390: form = NULL;
! 1391: for (i=1; i<=nlze2; i++)
! 1392: {
! 1393: ex[i]=mymyrand()>>randshift;
! 1394: if (ex[i])
! 1395: {
! 1396: p1 = primeform(Disc,stoi(factorbase[vperm[i]]),PRECREG);
! 1397: p1 = gpuigs(p1,ex[i]); form = comp(form,p1);
! 1398: }
! 1399: }
! 1400: if (!form) continue;
! 1401:
! 1402: fpc = factorisequad(form,KC,LIMC);
! 1403: if (fpc==1)
! 1404: {
! 1405: s++; col = (GEN)extramat[s];
! 1406: for (i=1; i<=nlze2; i++) col[vperm[i]] = -ex[i];
! 1407: for ( ; i<=KC; i++) col[vperm[i]]= 0;
! 1408: for (j=1; j<=primfact[0]; j++)
! 1409: {
! 1410: p=primfact[j]; ep=exprimfact[j];
! 1411: if (smodis((GEN)form[2], p<<1) > p) ep = -ep;
! 1412: col[numfactorbase[p]] += ep;
! 1413: }
! 1414: for (i=1; i<=KC; i++)
! 1415: if (col[i]) break;
! 1416: if (i>KC)
! 1417: {
! 1418: s--; avma = av;
! 1419: if (--MAXRELSUP == 0) return NULL;
! 1420: }
! 1421: else { av = avma; if (PRECREG) coeff(extramatc,1,s) = form[4]; }
! 1422: }
! 1423: else avma = av;
! 1424: if (DEBUGLEVEL)
! 1425: {
! 1426: if (fpc == 1) fprintferr(" %ld",s);
! 1427: else if (DEBUGLEVEL>1) fprintferr(".");
! 1428: flusherr();
! 1429: }
! 1430: }
! 1431: for (j=1; j<=extrarel; j++)
! 1432: {
! 1433: colg = cgetg(KC+1,t_COL);
! 1434: col = (GEN)extramat[j]; extramat[j] = (long) colg;
! 1435: for (k=1; k<=KC; k++)
! 1436: colg[k] = lstoi(col[vperm[k]]);
! 1437: }
! 1438: if (DEBUGLEVEL)
! 1439: {
! 1440: fprintferr("\n");
! 1441: msgtimer("extra relations");
! 1442: }
! 1443: return extramat;
! 1444: }
! 1445: #undef comp
! 1446:
! 1447: /*******************************************************************/
! 1448: /* */
! 1449: /* Imaginary Quadratic fields */
! 1450: /* */
! 1451: /*******************************************************************/
! 1452:
! 1453: static GEN
! 1454: imag_random_form(long current, long *ex)
! 1455: {
! 1456: long av = avma,i;
! 1457: GEN form,pc;
! 1458:
! 1459: for(;;)
! 1460: {
! 1461: form = pc = primeform(Disc,stoi(factorbase[current]),PRECREG);
! 1462: for (i=1; i<=lgsub; i++)
! 1463: {
! 1464: ex[i] = mymyrand()>>randshift;
! 1465: if (ex[i])
! 1466: form = compimag(form,powsubfactorbase[i][ex[i]]);
! 1467: }
! 1468: if (form != pc) return form;
! 1469: avma = av; /* ex = 0, try again */
! 1470: }
! 1471: }
! 1472:
! 1473: static void
! 1474: imag_relations(long lim, long s, long LIMC, long *ex, long **mat)
! 1475: {
! 1476: static long nbtest;
! 1477: long av = avma, i,j,pp,fpc,b1,b2,ep,current, first = (s==0);
! 1478: long *col,*fpd,*oldfact,*oldexp;
! 1479: GEN pc,form,form1;
! 1480:
! 1481: if (first) nbtest = 0 ;
! 1482: while (s<lim)
! 1483: {
! 1484: avma=av; nbtest++; current = first? 1+(s%KC): 1+s-RELSUP;
! 1485: form = imag_random_form(current,ex);
! 1486: fpc = factorisequad(form,KC,LIMC);
! 1487: if (!fpc)
! 1488: {
! 1489: if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
! 1490: continue;
! 1491: }
! 1492: if (fpc > 1)
! 1493: {
! 1494: fpd = largeprime(fpc,ex,current,0);
! 1495: if (!fpd)
! 1496: {
! 1497: if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
! 1498: continue;
! 1499: }
! 1500: form1 = powsubfactorbase[1][fpd[1]];
! 1501: for (i=2; i<=lgsub; i++)
! 1502: form1 = compimag(form1,powsubfactorbase[i][fpd[i]]);
! 1503: pc=primeform(Disc,stoi(factorbase[fpd[-2]]),0);
! 1504: form1=compimag(form1,pc);
! 1505: pp = fpc << 1;
! 1506: b1=smodis((GEN)form1[2], pp);
! 1507: b2=smodis((GEN)form[2], pp);
! 1508: if (b1 != b2 && b1+b2 != pp) continue;
! 1509:
! 1510: s++; col = mat[s];
! 1511: if (DEBUGLEVEL) { fprintferr(" %ld",s); flusherr(); }
! 1512: oldfact = primfact; oldexp = exprimfact;
! 1513: primfact = primfact1; exprimfact = exprimfact1;
! 1514: factorisequad(form1,KC,LIMC);
! 1515:
! 1516: if (b1==b2)
! 1517: {
! 1518: for (i=1; i<=lgsub; i++)
! 1519: col[numfactorbase[subbase[i]]] = fpd[i]-ex[i];
! 1520: col[fpd[-2]]++;
! 1521: for (j=1; j<=primfact[0]; j++)
! 1522: {
! 1523: pp=primfact[j]; ep=exprimfact[j];
! 1524: if (smodis((GEN)form1[2], pp<<1) > pp) ep = -ep;
! 1525: col[numfactorbase[pp]] -= ep;
! 1526: }
! 1527: }
! 1528: else
! 1529: {
! 1530: for (i=1; i<=lgsub; i++)
! 1531: col[numfactorbase[subbase[i]]] = -fpd[i]-ex[i];
! 1532: col[fpd[-2]]--;
! 1533: for (j=1; j<=primfact[0]; j++)
! 1534: {
! 1535: pp=primfact[j]; ep=exprimfact[j];
! 1536: if (smodis((GEN)form1[2], pp<<1) > pp) ep = -ep;
! 1537: col[numfactorbase[pp]] += ep;
! 1538: }
! 1539: }
! 1540: primfact = oldfact; exprimfact = oldexp;
! 1541: }
! 1542: else
! 1543: {
! 1544: s++; col = mat[s];
! 1545: if (DEBUGLEVEL) { fprintferr(" %ld",s); flusherr(); }
! 1546: for (i=1; i<=lgsub; i++)
! 1547: col[numfactorbase[subbase[i]]] = -ex[i];
! 1548: }
! 1549: for (j=1; j<=primfact[0]; j++)
! 1550: {
! 1551: pp=primfact[j]; ep=exprimfact[j];
! 1552: if (smodis((GEN)form[2], pp<<1) > pp) ep = -ep;
! 1553: col[numfactorbase[pp]] += ep;
! 1554: }
! 1555: col[current]--;
! 1556: if (!first && fpc == 1 && col[current] == 0)
! 1557: {
! 1558: s--; for (i=1; i<=KC; i++) col[i]=0;
! 1559: }
! 1560: }
! 1561: if (DEBUGLEVEL)
! 1562: {
! 1563: fprintferr("\nnbrelations/nbtest = %ld/%ld\n",s,nbtest);
! 1564: msgtimer("%s relations", first? "initial": "random");
! 1565: }
! 1566: }
! 1567:
! 1568: static int
! 1569: imag_be_honest(long *ex)
! 1570: {
! 1571: long p,fpc, s = KC, nbtest = 0, av = avma;
! 1572: GEN form;
! 1573:
! 1574: while (s<KC2)
! 1575: {
! 1576: p = factorbase[s+1];
! 1577: if (DEBUGLEVEL) { fprintferr(" %ld",p); flusherr(); }
! 1578: form = imag_random_form(s+1,ex);
! 1579: fpc = factorisequad(form,s,p-1);
! 1580: if (fpc == 1) { nbtest=0; s++; }
! 1581: else
! 1582: {
! 1583: nbtest++; if (nbtest>20) return 0;
! 1584: }
! 1585: avma = av;
! 1586: }
! 1587: return 1;
! 1588: }
! 1589:
! 1590: /*******************************************************************/
! 1591: /* */
! 1592: /* Real Quadratic fields */
! 1593: /* */
! 1594: /*******************************************************************/
! 1595:
! 1596: static GEN
! 1597: real_random_form(long *ex)
! 1598: {
! 1599: long av = avma, i;
! 1600: GEN p1,form = NULL;
! 1601:
! 1602: for(;;)
! 1603: {
! 1604: for (i=1; i<=lgsub; i++)
! 1605: {
! 1606: ex[i] = mymyrand()>>randshift;
! 1607: /* if (ex[i]) KB: BUG if I put this in. Why ??? */
! 1608: {
! 1609: p1 = powsubfactorbase[i][ex[i]];
! 1610: form = form? comprealform3(form,p1): p1;
! 1611: }
! 1612: }
! 1613: if (form) return form;
! 1614: avma = av;
! 1615: }
! 1616: }
! 1617:
! 1618: static void
! 1619: real_relations(long lim, long s, long LIMC, long *ex, long **mat, GEN glog2,
! 1620: GEN vecexpo)
! 1621: {
! 1622: static long nbtest;
! 1623: long av = avma,av1,av2,tetpil,i,j,p,fpc,b1,b2,ep,current, first = (s==0);
! 1624: long *col,*fpd,*oldfact,*oldexp,limstack;
! 1625: long findecycle,nbrhocumule,nbrho;
! 1626: GEN p1,p2,form,form0,form1,form2;
! 1627:
! 1628: limstack=stack_lim(av,1);
! 1629: if (first) { nbtest = 0; current = 0; }
! 1630: while (s<lim)
! 1631: {
! 1632: form = real_random_form(ex);
! 1633: if (!first)
! 1634: {
! 1635: current = 1+s-RELSUP;
! 1636: p1=redrealform(primeform(Disc,stoi(factorbase[current]),PRECREG));
! 1637: form = comprealform3(form,p1);
! 1638: }
! 1639: form0 = form; form1 = NULL;
! 1640: findecycle = nbrhocumule = 0;
! 1641: nbrho = -1; av1 = avma;
! 1642: while (s<lim)
! 1643: {
! 1644: if (low_stack(limstack, stack_lim(av,1)))
! 1645: {
! 1646: tetpil=avma;
! 1647: if(DEBUGMEM>1) err(warnmem,"real_relations [1]");
! 1648: if (!form1) form=gerepile(av1,tetpil,gcopy(form));
! 1649: else
! 1650: {
! 1651: GEN *gptr[2]; gptr[0]=&form1; gptr[1]=&form;
! 1652: gerepilemany(av1,gptr,2);
! 1653: }
! 1654: }
! 1655: if (nbrho < 0) nbrho = 0; /* first time in */
! 1656: else
! 1657: {
! 1658: if (findecycle) break;
! 1659: form = rhorealform(form);
! 1660: nbrho++; nbrhocumule++;
! 1661: if (first)
! 1662: {
! 1663: if (absi_equal((GEN)form[1],(GEN)form0[1])
! 1664: && egalii((GEN)form[2],(GEN)form0[2])
! 1665: && (!sens || signe(form0[1])==signe(form[1]))) findecycle=1;
! 1666: }
! 1667: else
! 1668: {
! 1669: if (sens || !signe(addii((GEN)form[1],(GEN)form[3])))
! 1670: { form=rhorealform(form); nbrho++; }
! 1671: else
! 1672: { setsigne(form[1],1); setsigne(form[3],-1); }
! 1673: if (egalii((GEN)form[1],(GEN)form0[1]) &&
! 1674: egalii((GEN)form[2],(GEN)form0[2])) break;
! 1675: }
! 1676: }
! 1677: nbtest++; fpc = factorisequad(form,KC,LIMC);
! 1678: if (!fpc)
! 1679: {
! 1680: if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
! 1681: continue;
! 1682: }
! 1683: if (fpc > 1)
! 1684: {
! 1685: fpd = largeprime(fpc,ex,current,nbrhocumule);
! 1686: if (!fpd)
! 1687: {
! 1688: if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
! 1689: continue;
! 1690: }
! 1691: if (!form1) form1 = initializeform5(ex);
! 1692: if (!first)
! 1693: {
! 1694: p1 = primeform(Disc,stoi(factorbase[current]),PRECREG);
! 1695: p1 = redrealform_init(p1); form1=comprealform5(form1,p1);
! 1696: }
! 1697: av2=avma;
! 1698: for (i=1; i<=nbrho; i++)
! 1699: {
! 1700: form1 = rhorealform(form1);
! 1701: if (low_stack(limstack, stack_lim(av,1)))
! 1702: {
! 1703: if(DEBUGMEM>1) err(warnmem,"real_relations [2]");
! 1704: tetpil=avma; form1=gerepile(av2,tetpil,gcopy(form1));
! 1705: }
! 1706: }
! 1707: nbrho = 0;
! 1708:
! 1709: form2=powsubfactorbase[1][fpd[1]];
! 1710: for (i=2; i<=lgsub; i++)
! 1711: form2 = comprealform5(form2,powsubfactorbase[i][fpd[i]]);
! 1712: if (fpd[-2])
! 1713: {
! 1714: p1 = primeform(Disc,stoi(factorbase[fpd[-2]]), PRECREG);
! 1715: p1 = redrealform_init(p1); form2=comprealform5(form2,p1);
! 1716: }
! 1717: av2=avma;
! 1718: for (i=1; i<=fpd[-3]; i++)
! 1719: {
! 1720: form2 = rhorealform(form2);
! 1721: if (low_stack(limstack, stack_lim(av,1)))
! 1722: {
! 1723: if(DEBUGMEM>1) err(warnmem,"real_relations [3]");
! 1724: tetpil=avma; form2=gerepile(av2,tetpil,gcopy(form2));
! 1725: }
! 1726: }
! 1727: if (!sens && signe(addii((GEN)form2[1],(GEN)form2[3])))
! 1728: {
! 1729: setsigne(form2[1],1);
! 1730: setsigne(form2[3],-1);
! 1731: }
! 1732: p = fpc << 1;
! 1733: b1=smodis((GEN)form2[2], p);
! 1734: b2=smodis((GEN)form1[2], p);
! 1735: if (b1 != b2 && b1+b2 != p) continue;
! 1736:
! 1737: s++; col = mat[s]; if (DEBUGLEVEL) fprintferr(" %ld",s);
! 1738: oldfact = primfact; oldexp = exprimfact;
! 1739: primfact = primfact1; exprimfact = exprimfact1;
! 1740: factorisequad(form2,KC,LIMC);
! 1741: if (b1==b2)
! 1742: {
! 1743: for (i=1; i<=lgsub; i++)
! 1744: col[numfactorbase[subbase[i]]] = fpd[i]-ex[i];
! 1745: for (j=1; j<=primfact[0]; j++)
! 1746: {
! 1747: p=primfact[j]; ep=exprimfact[j];
! 1748: if (smodis((GEN)form2[2], p<<1) > p) ep = -ep;
! 1749: col[numfactorbase[p]] -= ep;
! 1750: }
! 1751: if (fpd[-2]) col[fpd[-2]]++; /* implies !first */
! 1752: p1 = subii((GEN)form1[4],(GEN)form2[4]);
! 1753: p2 = divrr((GEN)form1[5],(GEN)form2[5]);
! 1754: }
! 1755: else
! 1756: {
! 1757: for (i=1; i<=lgsub; i++)
! 1758: col[numfactorbase[subbase[i]]] = -fpd[i]-ex[i];
! 1759: for (j=1; j<=primfact[0]; j++)
! 1760: {
! 1761: p=primfact[j]; ep=exprimfact[j];
! 1762: if (smodis((GEN)form2[2], p<<1) > p) ep = -ep;
! 1763: col[numfactorbase[p]] += ep;
! 1764: }
! 1765: if (fpd[-2]) col[fpd[-2]]--;
! 1766: p1 = addii((GEN)form1[4],(GEN)form2[4]);
! 1767: p2 = mulrr((GEN)form1[5],(GEN)form2[5]);
! 1768: }
! 1769: primfact = oldfact; exprimfact = oldexp;
! 1770: }
! 1771: else
! 1772: {
! 1773: if (!form1) form1 = initializeform5(ex);
! 1774: if (!first)
! 1775: {
! 1776: p1 = primeform(Disc,stoi(factorbase[current]),PRECREG);
! 1777: p1 = redrealform_init(p1); form1=comprealform5(form1,p1);
! 1778: }
! 1779: av2=avma;
! 1780: for (i=1; i<=nbrho; i++)
! 1781: {
! 1782: form1 = rhorealform(form1);
! 1783: if (low_stack(limstack, stack_lim(av,1)))
! 1784: {
! 1785: if(DEBUGMEM>1) err(warnmem,"real_relations [4]");
! 1786: tetpil=avma; form1=gerepile(av2,tetpil,gcopy(form1));
! 1787: }
! 1788: }
! 1789: nbrho = 0;
! 1790:
! 1791: s++; col = mat[s]; if (DEBUGLEVEL) fprintferr(" %ld",s);
! 1792: for (i=1; i<=lgsub; i++)
! 1793: col[numfactorbase[subbase[i]]] = -ex[i];
! 1794: p1 = (GEN) form1[4];
! 1795: p2 = (GEN) form1[5];
! 1796: }
! 1797: for (j=1; j<=primfact[0]; j++)
! 1798: {
! 1799: p=primfact[j]; ep=exprimfact[j];
! 1800: if (smodis((GEN)form1[2], p<<1) > p) ep = -ep;
! 1801: col[numfactorbase[p]] += ep;
! 1802: }
! 1803: p1 = mpadd(mulir(mulsi(EXP220,p1), glog2), mplog(absr(p2)));
! 1804: affrr(shiftr(p1,-1), (GEN)vecexpo[s]);
! 1805: if (!first)
! 1806: {
! 1807: col[current]--;
! 1808: if (fpc == 1 && col[current] == 0)
! 1809: { s--; for (i=1; i<=KC; i++) col[i]=0; }
! 1810: break;
! 1811: }
! 1812: }
! 1813: avma = av;
! 1814: }
! 1815: if (DEBUGLEVEL)
! 1816: {
! 1817: fprintferr("\nnbrelations/nbtest = %ld/%ld\n",s,nbtest);
! 1818: msgtimer("%s relations", first? "initial": "random");
! 1819: }
! 1820: }
! 1821:
! 1822: static int
! 1823: real_be_honest(long *ex)
! 1824: {
! 1825: long p,fpc,s = KC,nbtest = 0,av = avma;
! 1826: GEN p1,form,form0;
! 1827:
! 1828: while (s<KC2)
! 1829: {
! 1830: p = factorbase[s+1];
! 1831: if (DEBUGLEVEL) { fprintferr(" %ld",p); flusherr(); }
! 1832: form = real_random_form(ex);
! 1833: p1 = redrealform(primeform(Disc,stoi(p),PRECREG));
! 1834: form = comprealform3(form,p1); form0=form;
! 1835: for(;;)
! 1836: {
! 1837: fpc = factorisequad(form,s,p-1);
! 1838: if (fpc == 1) { nbtest=0; s++; break; }
! 1839: form = rhorealform(form);
! 1840: nbtest++; if (nbtest>20) return 0;
! 1841: if (sens || !signe(addii((GEN)form[1],(GEN)form[3])))
! 1842: form = rhorealform(form);
! 1843: else
! 1844: {
! 1845: setsigne(form[1],1);
! 1846: setsigne(form[3],-1);
! 1847: }
! 1848: if (egalii((GEN)form[1],(GEN)form0[1])
! 1849: && egalii((GEN)form[2],(GEN)form0[2])) break;
! 1850: }
! 1851: avma=av;
! 1852: }
! 1853: return 1;
! 1854: }
! 1855:
! 1856: static GEN
! 1857: gcdrealnoer(GEN a,GEN b,long *pte)
! 1858: {
! 1859: long e;
! 1860: GEN k1,r;
! 1861:
! 1862: if (typ(a)==t_INT)
! 1863: {
! 1864: if (typ(b)==t_INT) return mppgcd(a,b);
! 1865: k1=cgetr(lg(b)); affir(a,k1); a=k1;
! 1866: }
! 1867: else if (typ(b)==t_INT)
! 1868: { k1=cgetr(lg(a)); affir(b,k1); b=k1; }
! 1869: if (expo(a)<-5) return absr(b);
! 1870: if (expo(b)<-5) return absr(a);
! 1871: a=absr(a); b=absr(b);
! 1872: while (expo(b) >= -5 && signe(b))
! 1873: {
! 1874: k1=gcvtoi(divrr(a,b),&e);
! 1875: if (e > 0) return NULL;
! 1876: r=subrr(a,mulir(k1,b)); a=b; b=r;
! 1877: }
! 1878: *pte=expo(b); return absr(a);
! 1879: }
! 1880:
! 1881: static GEN
! 1882: get_reg(GEN matc, long sreg)
! 1883: {
! 1884: long i,e,maxe;
! 1885: GEN reg = mpabs(gcoeff(matc,1,1));
! 1886:
! 1887: e = maxe = 0;
! 1888: for (i=2; i<=sreg; i++)
! 1889: {
! 1890: reg = gcdrealnoer(gcoeff(matc,1,i),reg,&e);
! 1891: if (!reg) return NULL;
! 1892: maxe = maxe? max(maxe,e): e;
! 1893: }
! 1894: if (DEBUGLEVEL)
! 1895: {
! 1896: if (DEBUGLEVEL>7) { fprintferr("reg = "); outerr(reg); }
! 1897: msgtimer("regulator");
! 1898: }
! 1899: return reg;
! 1900: }
! 1901:
! 1902: GEN
! 1903: buchquad(GEN D, double cbach, double cbach2, long RELSUP0, long flag, long prec)
! 1904: {
! 1905: long av0 = avma,av,tetpil,KCCO,KCCOPRO,i,j,s, *ex,**mat;
! 1906: long extrarel,nrelsup,nreldep,LIMC,LIMC2,cp,nbram,nlze;
! 1907: GEN p1,h,W,met,res,basecl,dr,c_1,pdep,C,B,extramat,extraC;
! 1908: GEN reg,vecexpo,glog2,cst;
! 1909: double drc,lim,LOGD;
! 1910:
! 1911: Disc = D; if (typ(Disc)!=t_INT) err(typeer,"buchquad");
! 1912: s=mod4(Disc);
! 1913: switch(signe(Disc))
! 1914: {
! 1915: case -1:
! 1916: if (cmpis(Disc,-4) >= 0)
! 1917: {
! 1918: p1=cgetg(6,t_VEC); p1[1]=p1[4]=p1[5]=un;
! 1919: p1[2]=p1[3]=lgetg(1,t_VEC); return p1;
! 1920: }
! 1921: if (s==2 || s==1) err(funder2,"buchquad");
! 1922: PRECREG=0; break;
! 1923:
! 1924: case 1:
! 1925: if (s==2 || s==3) err(funder2,"buchquad");
! 1926: if (flag)
! 1927: err(talker,"sorry, narrow class group not implemented. Use bnfnarrow");
! 1928: PRECREG=1; break;
! 1929:
! 1930: default: err(talker,"zero discriminant in quadclassunit");
! 1931: }
! 1932: if (carreparfait(Disc)) err(talker,"square argument in quadclassunit");
! 1933: if (!isfundamental(Disc))
! 1934: err(warner,"not a fundamental discriminant in quadclassunit");
! 1935: buch_init(); RELSUP = RELSUP0; sens = flag;
! 1936: dr=cgetr(3); affir(Disc,dr); drc=fabs(rtodbl(dr)); LOGD=log(drc);
! 1937: lim=sqrt(drc); cst = mulrr(lfunc(Disc), dbltor(lim));
! 1938: if (!PRECREG) lim /= sqrt(3.);
! 1939: cp = (long)exp(sqrt(LOGD*log(LOGD)/8.0));
! 1940: if (cp < 13) cp = 13;
! 1941: av = avma; cbach /= 2;
! 1942:
! 1943: INCREASE: avma = av; cbach = check_bach(cbach,6.);
! 1944: nreldep = nrelsup = 0;
! 1945: LIMC = (long)(cbach*LOGD*LOGD);
! 1946: if (LIMC < cp) LIMC=cp;
! 1947: LIMC2 = max(20, (long)(max(cbach,cbach2)*LOGD*LOGD));
! 1948: if (LIMC2 < LIMC) LIMC2 = LIMC;
! 1949: if (PRECREG)
! 1950: {
! 1951: PRECREG = max(prec+1, MEDDEFAULTPREC + 2*(expi(Disc)>>TWOPOTBITS_IN_LONG));
! 1952: glog2 = glog(gdeux,PRECREG);
! 1953: sqrtD = gsqrt(Disc,PRECREG);
! 1954: isqrtD = gfloor(sqrtD);
! 1955: }
! 1956: factorbasequad(Disc,LIMC2,LIMC);
! 1957: if (!KC) goto INCREASE;
! 1958:
! 1959: vperm = cgetg(KC+1,t_VECSMALL); for (i=1; i<=KC; i++) vperm[i]=i;
! 1960: nbram = subfactorbasequad(lim,KC);
! 1961: if (nbram == -1) { desalloc(NULL); goto INCREASE; }
! 1962: KCCO = KC + RELSUP;
! 1963: if (DEBUGLEVEL) { fprintferr("KC = %ld, KCCO = %ld\n",KC,KCCO); flusherr(); }
! 1964: powsubfact(lgsub,CBUCH+7);
! 1965:
! 1966: mat = (long**) gpmalloc((KCCO+1)*sizeof(long*));
! 1967: setlg(mat, KCCO+1);
! 1968: for (i=1; i<=KCCO; i++)
! 1969: {
! 1970: mat[i] = (long*) gpmalloc((KC+1)*sizeof(long));
! 1971: for (j=1; j<=KC; j++) mat[i][j]=0;
! 1972: }
! 1973: ex = new_chunk(lgsub+1);
! 1974: limhash = (LIMC<(MAXHALFULONG>>1))? LIMC*LIMC: HIGHBIT>>1;
! 1975: for (i=0; i<HASHT; i++) hashtab[i]=NULL;
! 1976:
! 1977: s = lgsub+nbram+RELSUP;
! 1978: if (PRECREG)
! 1979: {
! 1980: vecexpo=cgetg(KCCO+1,t_VEC);
! 1981: for (i=1; i<=KCCO; i++) vecexpo[i]=lgetr(PRECREG);
! 1982: real_relations(s,0,LIMC,ex,mat,glog2,vecexpo);
! 1983: real_relations(KCCO,s,LIMC,ex,mat,glog2,vecexpo);
! 1984: }
! 1985: else
! 1986: {
! 1987: imag_relations(s,0,LIMC,ex,mat);
! 1988: imag_relations(KCCO,s,LIMC,ex,mat);
! 1989: }
! 1990: if (KC2 > KC)
! 1991: {
! 1992: if (DEBUGLEVEL)
! 1993: fprintferr("be honest for primes from %ld to %ld\n",
! 1994: factorbase[KC+1],factorbase[KC2]);
! 1995: s = PRECREG? real_be_honest(ex): imag_be_honest(ex);
! 1996: if (DEBUGLEVEL)
! 1997: {
! 1998: fprintferr("\n");
! 1999: msgtimer("be honest");
! 2000: }
! 2001: if (!s) { desalloc(mat); goto INCREASE; }
! 2002: }
! 2003: C=cgetg(KCCO+1,t_MAT);
! 2004: if (PRECREG)
! 2005: {
! 2006: for (i=1; i<=KCCO; i++)
! 2007: {
! 2008: C[i]=lgetg(2,t_COL); coeff(C,1,i)=vecexpo[i];
! 2009: }
! 2010: if (DEBUGLEVEL>7) { fprintferr("C: "); outerr(C); flusherr(); }
! 2011: }
! 2012: else
! 2013: for (i=1; i<=KCCO; i++) C[i]=lgetg(1,t_COL);
! 2014: W = hnfspec(mat,vperm,&pdep,&B,&C,lgsub);
! 2015: nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;
! 2016:
! 2017: KCCOPRO=KCCO;
! 2018: if (nlze)
! 2019: {
! 2020: EXTRAREL:
! 2021: s = PRECREG? 2: 1; extrarel=nlze+2;
! 2022: extraC=cgetg(extrarel+1,t_MAT);
! 2023: for (i=1; i<=extrarel; i++) extraC[i]=lgetg(s,t_COL);
! 2024: extramat = extra_relations(LIMC,ex,nlze,extraC);
! 2025: if (!extramat) { desalloc(mat); goto INCREASE; }
! 2026: W = hnfadd(W,vperm,&pdep,&B,&C, extramat,extraC);
! 2027: nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;
! 2028: KCCOPRO += extrarel;
! 2029: if (nlze)
! 2030: {
! 2031: if (++nreldep > 5) { desalloc(mat); goto INCREASE; }
! 2032: goto EXTRAREL;
! 2033: }
! 2034: }
! 2035: /* tentative class number */
! 2036: h=gun; for (i=1; i<lg(W); i++) h=mulii(h,gcoeff(W,i,i));
! 2037: if (PRECREG)
! 2038: {
! 2039: /* tentative regulator */
! 2040: reg = get_reg(C, KCCOPRO - (lg(B)-1) - (lg(W)-1));
! 2041: if (!reg)
! 2042: {
! 2043: desalloc(mat);
! 2044: prec = (PRECREG<<1)-2; goto INCREASE;
! 2045: }
! 2046: if (gexpo(reg)<=-3)
! 2047: {
! 2048: if (++nrelsup <= 7)
! 2049: {
! 2050: if (DEBUGLEVEL) { fprintferr("regulateur nul\n"); flusherr(); }
! 2051: nlze=min(KC,nrelsup); goto EXTRAREL;
! 2052: }
! 2053: desalloc(mat); goto INCREASE;
! 2054: }
! 2055: c_1 = divrr(gmul2n(gmul(h,reg),1), cst);
! 2056: }
! 2057: else
! 2058: {
! 2059: reg = gun;
! 2060: c_1 = divrr(gmul(h,mppi(DEFAULTPREC)), cst);
! 2061: }
! 2062:
! 2063: if (gcmpgs(gmul2n(c_1,2),3)<0) { c_1=stoi(10); nrelsup=7; }
! 2064: if (gcmpgs(gmul2n(c_1,1),3)>0)
! 2065: {
! 2066: nrelsup++;
! 2067: if (nrelsup<=7)
! 2068: {
! 2069: if (DEBUGLEVEL)
! 2070: { fprintferr("***** check = %f\n\n",gtodouble(c_1)); flusherr(); }
! 2071: nlze=min(KC,nrelsup); goto EXTRAREL;
! 2072: }
! 2073: if (cbach < 5.99) { desalloc(mat); goto INCREASE; }
! 2074: err(warner,"suspicious check. Suggest increasing extra relations.");
! 2075: }
! 2076: basecl = get_clgp(Disc,W,&met,PRECREG);
! 2077: s = lg(basecl); desalloc(mat); tetpil=avma;
! 2078:
! 2079: res=cgetg(6,t_VEC);
! 2080: res[1]=lcopy(h); p1=cgetg(s,t_VEC);
! 2081: for (i=1; i<s; i++) p1[i] = (long)icopy(gcoeff(met,i,i));
! 2082: res[2]=(long)p1;
! 2083: res[3]=lcopy(basecl);
! 2084: res[4]=lcopy(reg);
! 2085: res[5]=lcopy(c_1); return gerepile(av0,tetpil,res);
! 2086: }
! 2087:
! 2088: GEN
! 2089: buchimag(GEN D, GEN c, GEN c2, GEN REL)
! 2090: {
! 2091: return buchquad(D,gtodouble(c),gtodouble(c2),itos(REL), 0,0);
! 2092: }
! 2093:
! 2094: GEN
! 2095: buchreal(GEN D, GEN sens0, GEN c, GEN c2, GEN REL, long prec)
! 2096: {
! 2097: return buchquad(D,gtodouble(c),gtodouble(c2),itos(REL), itos(sens0),prec);
! 2098: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to buch1.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