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