Annotation of OpenXM_contrib/pari-2.2/src/basemath/buch1.c, Revision 1.2
1.2 ! noro 1: /* $Id: buch1.c,v 1.51 2002/09/05 10:55:35 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /*******************************************************************/
17: /* */
18: /* CLASS GROUP AND REGULATOR (McCURLEY, BUCHMANN) */
19: /* QUADRATIC FIELDS */
20: /* */
21: /*******************************************************************/
22: #include "pari.h"
23:
24: const int narrow = 0; /* should set narrow = flag in buchquad, but buggy */
25:
1.2 ! noro 26: /* For largeprime() hashtable. Note that hashed pseudoprimes are odd (unless
! 27: * 2 | index), hence the low order bit is not useful. So we hash
! 28: * HASHBITS bits starting at bit 1, not bit 0 */
! 29: #define HASHBITS 10
! 30: const int HASHT = 1 << HASHBITS;
! 31:
! 32: static long
! 33: hash(long q) { return (q & ((1 << (HASHBITS+1)) - 1)) >> 1; }
! 34: #undef HASHBITS
! 35:
1.1 noro 36: /* See buch2.c:
37: * precision en digits decimaux=2*(#digits decimaux de Disc)+50
1.2 ! noro 38: * on prendra les p decomposes tels que prod(p)>lim dans la subFB
1.1 noro 39: * LIMC=Max(c.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))
40: * LIMC2=Max(6.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))
1.2 ! noro 41: * subFB contient les indices des p decomposes tels que prod(p) > sqrt(Disc)
! 42: * powsubFB est la table des puissances des formes dans subFB
1.1 noro 43: */
1.2 ! noro 44:
! 45: #define RANDOM_BITS 4
! 46: static const long CBUCH = (1<<RANDOM_BITS)-1;
! 47: static const long randshift=BITS_IN_RANDOM-1 - RANDOM_BITS;
! 48: #undef RANDOM_BITS
! 49:
! 50: static long KC,KC2,limhash,RELSUP,PRECREG;
! 51: static long *primfact,*exprimfact,*badprim;
! 52: static long *FB,*numFB, **hashtab;
! 53: static GEN powsubFB,vperm,subFB,Disc,sqrtD,isqrtD;
1.1 noro 54:
55: GEN buchquad(GEN D, double c, double c2, long RELSUP0, long flag, long prec);
56: extern GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus);
1.2 ! noro 57: extern GEN colreducemodHNF(GEN x, GEN y, GEN *Q);
! 58: extern GEN quadhilbertreal(GEN D, GEN flag, long prec);
! 59: extern void comp_gen(GEN z,GEN x,GEN y);
! 60: extern GEN codeform5(GEN x, long prec);
! 61: extern GEN comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD);
! 62: extern GEN redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD);
! 63: extern GEN rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD);
! 64: extern GEN cgetalloc(GEN x, size_t l, long t);
1.1 noro 65:
66: GEN
67: quadclassunit0(GEN x, long flag, GEN data, long prec)
68: {
69: long lx,RELSUP0;
70: double cbach, cbach2;
71:
72: if (!data) lx=1;
73: else
74: {
75: lx = lg(data);
76: if (typ(data)!=t_VEC || lx > 7)
77: err(talker,"incorrect parameters in quadclassunit");
78: if (lx > 4) lx = 4;
79: }
80: cbach = cbach2 = 0.1; RELSUP0 = 5;
81: switch(lx)
82: {
83: case 4: RELSUP0 = itos((GEN)data[3]);
84: case 3: cbach2 = gtodouble((GEN)data[2]);
85: case 2: cbach = gtodouble((GEN)data[1]);
86: }
87: return buchquad(x,cbach,cbach2,RELSUP0,flag,prec);
88: }
89:
90: /*******************************************************************/
91: /* */
92: /* Hilbert and Ray Class field using CM (Schertz) */
93: /* */
94: /*******************************************************************/
95:
96: int
97: isoforder2(GEN form)
98: {
99: GEN a=(GEN)form[1],b=(GEN)form[2],c=(GEN)form[3];
100: return !signe(b) || absi_equal(a,b) || egalii(a,c);
101: }
102:
103: GEN
104: getallforms(GEN D, long *pth, GEN *ptz)
105: {
106: long d = itos(D), t, b2, a, b, c, h, dover3 = labs(d)/3;
107: GEN z, L=cgetg(labs(d), t_VEC);
108: b2 = b = (d&1); h = 0; z=gun;
109: while (b2 <= dover3)
110: {
111: t = (b2-d)/4;
112: for (a=b?b:1; a*a<=t; a++)
113: if (t%a == 0)
114: {
115: c = t/a; z = mulsi(a,z);
116: L[++h] = (long)qfi(stoi(a),stoi(b),stoi(c));
117: if (b && a != b && a*a != t)
118: L[++h] = (long)qfi(stoi(a),stoi(-b),stoi(c));
119: }
120: b+=2; b2=b*b;
121: }
122: *pth = h; *ptz = z; setlg(L,h+1); return L;
123: }
124:
125: #define MOD4(x) ((((GEN)x)[2])&3)
126: #define MAXL 300
1.2 ! noro 127: /* find P and Q two non principal prime ideals (above p,q) such that
1.1 noro 128: * (pq, z) = 1 [coprime to all class group representatives]
129: * cl(P) = cl(Q) if P has order 2 in Cl(K)
130: * Try to have e = 24 / gcd(24, (p-1)(q-1)) as small as possible */
131: void
132: get_pq(GEN D, GEN z, GEN flag, GEN *ptp, GEN *ptq)
133: {
134: GEN wp=cgetg(MAXL,t_VEC), wlf=cgetg(MAXL,t_VEC), court=icopy(gun);
135: GEN p, form;
136: long i, ell, l = 1, d = itos(D);
137: byteptr diffell = diffptr + 2;
138:
139: if (typ(flag)==t_VEC)
140: { /* assume flag = [p,q]. Check nevertheless */
141: for (i=1; i<lg(flag); i++)
142: {
143: ell = itos((GEN)flag[i]);
144: if (smodis(z,ell) && kross(d,ell) > 0)
145: {
146: form = redimag(primeform(D,(GEN)flag[i],0));
147: if (!gcmp1((GEN)form[1])) {
148: wp[l] = flag[i];
149: l++; if (l == 3) break;
150: }
151: }
152: }
153: if (l<3) err(talker,"[quadhilbert] incorrect values in flag: %Z", flag);
154: *ptp = (GEN)wp[1];
155: *ptq = (GEN)wp[2]; return;
156: }
157:
158: ell = 3;
159: while (l < 3 || ell<=MAXL)
160: {
1.2 ! noro 161: NEXT_PRIME_VIADIFF_CHECK(ell, diffell);
1.1 noro 162: if (smodis(z,ell) && kross(d,ell) > 0)
163: {
164: court[2]=ell; form = redimag(primeform(D,court,0));
165: if (!gcmp1((GEN)form[1])) {
166: wp[l] = licopy(court);
167: wlf[l] = (long)form; l++;
168: }
169: }
170: }
171: setlg(wp,l); setlg(wlf,l);
172:
173: for (i=1; i<l; i++)
174: if (smodis((GEN)wp[i],3) == 1) break;
175: if (i==l) i = 1;
176: p = (GEN)wp[i]; form = (GEN)wlf[i];
177: i = l;
178: if (isoforder2(form))
179: {
180: long oki = 0;
181: for (i=1; i<l; i++)
182: if (gegal((GEN)wlf[i],form))
183: {
184: if (MOD4(p) == 1 || MOD4(wp[i]) == 1) break;
185: if (!oki) oki = i; /* not too bad, still e = 2 */
186: }
187: if (i==l) i = oki;
188: if (!i) err(bugparier,"quadhilbertimag (can't find p,q)");
189: }
190: else
191: {
192: if (MOD4(p) == 3)
193: {
194: for (i=1; i<l; i++)
195: if (MOD4(wp[i]) == 1) break;
196: }
197: if (i==l) i = 1;
198: }
199: *ptp = p;
200: *ptq = (GEN)wp[i];
201: }
202: #undef MAXL
203:
204: static GEN
205: gpq(GEN form, GEN p, GEN q, long e, GEN sqd, GEN u, long prec)
206: {
207: GEN a2 = shifti((GEN)form[1], 1);
208: GEN b = (GEN)form[2], p1,p2,p3,p4;
209: GEN w = lift(chinois(gmodulcp(negi(b),a2), u));
210: GEN al = cgetg(3,t_COMPLEX);
211: al[1] = lneg(gdiv(w,a2));
212: al[2] = ldiv(sqd,a2);
213: p1 = trueeta(gdiv(al,p),prec);
214: p2 = egalii(p,q)? p1: trueeta(gdiv(al,q),prec);
215: p3 = trueeta(gdiv(al,mulii(p,q)),prec);
216: p4 = trueeta(al,prec);
217: return gpowgs(gdiv(gmul(p1,p2),gmul(p3,p4)), e);
218: }
219:
220: /* returns an equation for the Hilbert class field of the imaginary
221: * quadratic field of discriminant D if flag=0, a vector of
222: * two-component vectors [form,g(form)] where g() is the root of the equation
223: * if flag is non-zero.
224: */
225: static GEN
226: quadhilbertimag(GEN D, GEN flag)
227: {
1.2 ! noro 228: long h, i, e, prec;
! 229: gpmem_t av=avma;
1.1 noro 230: GEN z,L,P,p,q,qfp,qfq,up,uq,u;
231: int raw = ((typ(flag)==t_INT && signe(flag)));
232:
1.2 ! noro 233: if (DEBUGLEVEL>=2) (void)timer2();
1.1 noro 234: if (gcmpgs(D,-11) >= 0) return polx[0];
235: L = getallforms(D,&h,&z);
236: if (DEBUGLEVEL>=2) msgtimer("class number = %ld",h);
237: if (h == 1) { avma=av; return polx[0]; }
238:
239: get_pq(D, z, flag, &p, &q);
240: e = 24 / cgcd((smodis(p,24)-1) * (smodis(q,24)-1), 24);
241: if(DEBUGLEVEL>=2) fprintferr("p = %Z, q = %Z, e = %ld\n",p,q,e);
242: qfp = primeform(D,p,0); up = gmodulcp((GEN)qfp[2],shifti(p,1));
243: if (egalii(p,q))
244: {
245: u = (GEN)compimagraw(qfp,qfp)[2];
246: u = gmodulcp(u, shifti(mulii(p,q),1));
247: }
248: else
249: {
250: qfq = primeform(D,q,0); uq = gmodulcp((GEN)qfq[2],shifti(q,1));
251: u = chinois(up,uq);
252: }
253: prec = raw? DEFAULTPREC: 3;
254: for(;;)
255: {
1.2 ! noro 256: long ex, exmax = 0;
! 257: gpmem_t av0 = avma;
1.1 noro 258: GEN lead, sqd = gsqrt(negi(D),prec);
1.2 ! noro 259: P = cgetg(h+1,t_VEC);
1.1 noro 260: for (i=1; i<=h; i++)
261: {
262: GEN v, s = gpq((GEN)L[i], p, q, e, sqd, u, prec);
263: if (raw)
264: {
265: v = cgetg(3,t_VEC); P[i] = (long)v;
266: v[1] = L[i];
267: v[2] = (long)s;
268: }
269: else P[i] = (long)s;
1.2 ! noro 270: ex = gexpo(s); if (ex > 0) exmax += ex;
1.1 noro 271: }
272: if (DEBUGLEVEL>=2) msgtimer("roots");
273: if (raw) { P = gcopy(P); break; }
274: /* to avoid integer overflow (1 + 0.) */
275: lead = (exmax < bit_accuracy(prec))? gun: realun(prec);
276:
277: P = greal(roots_to_pol_intern(lead,P,0,0));
278: P = grndtoi(P,&exmax);
279: if (DEBUGLEVEL>=2) msgtimer("product, error bits = %ld",exmax);
280: if (exmax <= -10)
281: {
282: if (typ(flag)==t_VEC && !issquarefree(P)) { avma=av; return gzero; }
283: break;
284: }
285: avma = av0; prec += (DEFAULTPREC-2) + (1 + (exmax >> TWOPOTBITS_IN_LONG));
286: if (DEBUGLEVEL) err(warnprec,"quadhilbertimag",prec);
287: }
288: return gerepileupto(av,P);
289: }
290:
291: GEN
292: quadhilbert(GEN D, GEN flag, long prec)
293: {
294: if (!flag) flag = gzero;
295: if (typ(D)!=t_INT)
296: {
297: D = checkbnf(D);
298: if (degpol(gmael(D,7,1)) != 2)
299: err(talker,"not a polynomial of degree 2 in quadhilbert");
300: D=gmael(D,7,3);
301: }
302: else
303: {
304: if (!isfundamental(D))
305: err(talker,"quadhilbert needs a fundamental discriminant");
306: }
307: return (signe(D)>0)? quadhilbertreal(D,flag,prec)
308: : quadhilbertimag(D,flag);
309: }
310:
311: /* AUXILLIARY ROUTINES FOR QUADRAYIMAGWEI */
312: #define to_approx(nf,a) ((GEN)gmul(gmael((nf),5,1), (a))[1])
313:
314: /* Z-basis for a (over C) */
315: static GEN
316: get_om(GEN nf, GEN a)
317: {
318: GEN om = cgetg(3,t_VEC);
319: om[1] = (long)to_approx(nf,(GEN)a[2]);
320: om[2] = (long)to_approx(nf,(GEN)a[1]);
321: return om;
322: }
323:
324: /* Compute all elts in class group G = [|G|,c,g], c=cyclic factors, g=gens.
325: * Set list[j + 1] = g1^e1...gk^ek where j is the integer
326: * ek + ck [ e(k-1) + c(k-1) [... + c2 [e1]]...] */
327: static GEN
1.2 ! noro 328: getallelts(GEN bnr)
1.1 noro 329: {
1.2 ! noro 330: GEN nf,G,C,c,g, *list, **pows, *gk;
1.1 noro 331: long lc,i,j,k,no;
332:
1.2 ! noro 333: nf = checknf(bnr);
! 334: G = (GEN)bnr[5];
! 335:
1.1 noro 336: no = itos((GEN)G[1]);
337: c = (GEN)G[2];
338: g = (GEN)G[3]; lc = lg(c)-1;
339: list = (GEN*) cgetg(no+1,t_VEC);
340: if (!lc)
341: {
342: list[1] = idealhermite(nf,gun);
343: return (GEN)list;
344: }
345: pows = (GEN**)cgetg(lc+1,t_VEC);
346: c = dummycopy(c); settyp(c, t_VECSMALL);
347: for (i=1; i<=lc; i++)
348: {
349: c[i] = k = itos((GEN)c[i]);
350: gk = (GEN*)cgetg(k, t_VEC); gk[1] = (GEN)g[i];
1.2 ! noro 351: for (j=2; j<k; j++)
! 352: gk[j] = idealmodidele(bnr, idealmul(nf, gk[j-1], gk[1]));
1.1 noro 353: pows[i] = gk; /* powers of g[i] */
354: }
355:
1.2 ! noro 356: C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc];
1.1 noro 357: for (i=2; i<=lc; i++) C[i] = C[i-1] * c[lc-i+1];
358: /* C[i] = c(k-i+1) * ... * ck */
359: /* j < C[i+1] <==> j only involves g(k-i)...gk */
360: i = 1; list[1] = 0; /* dummy */
361: for (j=1; j < C[1]; j++)
362: list[j + 1] = pows[lc][j];
363: for ( ; j<no; j++)
364: {
365: GEN p1,p2;
366: if (j == C[i+1]) i++;
1.2 ! noro 367: p2 = pows[lc-i][j/C[i]];
! 368: p1 = list[j%C[i] + 1];
! 369: if (p1) p2 = idealmodidele(bnr, idealmul(nf,p2,p1));
1.1 noro 370: list[j + 1] = p2;
371: }
372: list[1] = idealhermite(nf,gun);
373: return (GEN)list;
374: }
375:
376: /* x quadratic integer (approximate), recognize it. If error return NULL */
377: static GEN
378: findbezk(GEN nf, GEN x)
379: {
380: GEN a,b, M = gmael(nf,5,1), y = cgetg(3, t_COL), u = gcoeff(M,1,2);
381: long ea,eb;
382:
383: b = grndtoi(gdiv(gimag(x), gimag(u)), &eb);
384: a = grndtoi(greal(gsub(x, gmul(b,u))),&ea);
385: if (ea>-20 || eb>-20) return NULL;
386: if (!signe(b)) return a;
387: y[1] = (long)a;
388: y[2] = (long)b; return basistoalg(nf,y);
389: }
390:
391: static GEN
392: findbezk_pol(GEN nf, GEN x)
393: {
394: long i, lx = lgef(x);
395: GEN y = cgetg(lx,t_POL);
396: for (i=2; i<lx; i++)
397: if (! (y[i] = (long)findbezk(nf,(GEN)x[i])) ) return NULL;
398: y[1] = x[1]; return y;
399: }
400:
401: GEN
402: form_to_ideal(GEN x)
403: {
404: long tx = typ(x);
405: GEN b,c, y = cgetg(3, t_MAT);
406: if (tx != t_QFR && tx != t_QFI) err(typeer,"form_to_ideal");
407: c = cgetg(3,t_COL); y[1] = (long)c;
408: c[1] = x[1]; c[2] = zero;
409: c = cgetg(3,t_COL); y[2] = (long)c;
410: b = negi((GEN)x[2]);
411: if (mpodd(b)) b = addis(b,1);
412: c[1] = lshifti(b,-1); c[2] = un; return y;
413: }
414:
415: /* P as output by quadhilbertimag, convert forms to ideals */
416: static void
417: convert_to_id(GEN P)
418: {
419: long i,l = lg(P);
420: for (i=1; i<l; i++)
421: {
422: GEN p1 = (GEN)P[i];
423: p1[1] = (long)form_to_ideal((GEN)p1[1]);
424: }
425: }
426:
427: /* P approximation computed at initial precision prec. Compute needed prec
428: * to know P with 1 word worth of trailing decimals */
429: static long
430: get_prec(GEN P, long prec)
431: {
432: long k = gprecision(P);
433: if (k == 3) return (prec<<1)-2; /* approximation not trustworthy */
434: k = prec - k; /* lost precision when computing P */
435: if (k < 0) k = 0;
436: k += MEDDEFAULTPREC + (gexpo(P) >> TWOPOTBITS_IN_LONG);
437: if (k <= prec) k = (prec<<1)-2; /* dubious: old prec should have worked */
438: return k;
439: }
440:
441: /* AUXILLIARY ROUTINES FOR QUADRAYSIGMA */
442:
443: /* Compute data for ellphist */
444: static GEN
445: ellphistinit(GEN om, long prec)
446: {
447: GEN p1,res,om1b,om2b, om1 = (GEN)om[1], om2 = (GEN)om[2];
448:
449: if (gsigne(gimag(gdiv(om1,om2))) < 0)
450: {
1.2 ! noro 451: p1 = om1; om1 = om2; om2 = p1;
! 452: om = cgetg(3,t_VEC);
! 453: om[1] = (long)om1;
! 454: om[2] = (long)om2;
1.1 noro 455: }
456: om1b = gconj(om1);
457: om2b = gconj(om2); res = cgetg(4,t_VEC);
458: res[1] = ldivgs(elleisnum(om,2,0,prec),12);
459: res[2] = ldiv(PiI2(prec), gmul(om2, gimag(gmul(om1b,om2))));
460: res[3] = (long)om2b; return res;
461: }
462:
463: /* Computes log(phi^*(z,om)), using res computed by ellphistinit */
464: static GEN
465: ellphist(GEN om, GEN res, GEN z, long prec)
466: {
467: GEN u = gimag(gmul(z, (GEN)res[3]));
468: GEN zst = gsub(gmul(u, (GEN)res[2]), gmul(z,(GEN)res[1]));
469: return gsub(ellsigma(om,z,1,prec),gmul2n(gmul(z,zst),-1));
470: }
471:
472: /* Computes phi^*(la,om)/phi^*(1,om) where (1,om) is an oriented basis of the
473: ideal gf*gc^{-1} */
474: static GEN
475: computeth2(GEN om, GEN la, long prec)
476: {
477: GEN p1,p2,res = ellphistinit(om,prec);
478:
479: p1 = gsub(ellphist(om,res,la,prec), ellphist(om,res,gun,prec));
480: p2 = gimag(p1);
481: if (gexpo(greal(p1))>20 || gexpo(p2)> bit_accuracy(min(prec,lg(p2)))-10)
482: return NULL;
483: return gexp(p1,prec);
484: }
485:
486: /* Computes P_2(X)=polynomial in Z_K[X] closest to prod_gc(X-th2(gc)) where
487: the product is over the ray class group bnr.*/
488: static GEN
489: computeP2(GEN bnr, GEN la, int raw, long prec)
490: {
1.2 ! noro 491: long clrayno, i, first = 1;
! 492: gpmem_t av=avma, av2;
1.1 noro 493: GEN listray,P0,P,f,lanum, nf = checknf(bnr);
1.2 ! noro 494:
1.1 noro 495: f = gmael3(bnr,2,1,1);
496: if (typ(la) != t_COL) la = algtobasis(nf,la);
1.2 ! noro 497: listray = getallelts(bnr);
1.1 noro 498: clrayno = lg(listray)-1; av2 = avma;
499: PRECPB:
500: if (!first)
501: {
502: if (DEBUGLEVEL) err(warnprec,"computeP2",prec);
503: nf = gerepileupto(av2, nfnewprec(checknf(bnr),prec));
504: }
505: first = 0; lanum = to_approx(nf,la);
506: P = cgetg(clrayno+1,t_VEC);
507: for (i=1; i<=clrayno; i++)
508: {
509: GEN om = get_om(nf, idealdiv(nf,f,(GEN)listray[i]));
510: GEN v, s = computeth2(om,lanum,prec);
511: if (!s) { prec = (prec<<1)-2; goto PRECPB; }
512: if (raw)
513: {
514: v = cgetg(3,t_VEC); P[i] = (long)v;
515: v[1] = (long)listray[i];
516: v[2] = (long)s;
517: }
518: else P[i] = (long)s;
519: }
520: if (!raw)
521: {
522: P0 = roots_to_pol(P, 0);
523: P = findbezk_pol(nf, P0);
524: if (!P) { prec = get_prec(P0, prec); goto PRECPB; }
525: }
526: return gerepilecopy(av, P);
527: }
528:
529: #define nexta(a) (a>0 ? -a : 1-a)
530: static GEN
531: do_compo(GEN x, GEN y)
532: {
533: long a, ph = degpol(y);
534: GEN z;
1.2 ! noro 535: y = gmul(gpowgs(polx[0],ph), gsubst(y,0,gdiv(polx[MAXVARN],polx[0])));
1.1 noro 536: for (a = 0;; a = nexta(a))
537: {
538: if (a) x = gsubst(x, 0, gaddsg(a, polx[0]));
539: z = subres(x,y);
540: z = gsubst(z, MAXVARN, polx[0]);
541: if (issquarefree(z)) return z;
542: }
543: }
544: #undef nexta
545:
546: static GEN
547: galoisapplypol(GEN nf, GEN s, GEN x)
548: {
549: long i, lx = lg(x);
550: GEN y = cgetg(lx,t_POL);
551:
552: for (i=2; i<lx; i++) y[i] = (long)galoisapply(nf,s,(GEN)x[i]);
553: y[1] = x[1]; return y;
554: }
555:
556: /* x quadratic, write it as ua + v, u,v rational */
557: static GEN
558: findquad(GEN a, GEN x, GEN p)
559: {
1.2 ! noro 560: long tu, tv;
! 561: gpmem_t av = avma;
1.1 noro 562: GEN u,v;
563: if (typ(x) == t_POLMOD) x = (GEN)x[2];
564: if (typ(a) == t_POLMOD) a = (GEN)a[2];
565: u = poldivres(x, a, &v);
566: u = simplify(u); tu = typ(u);
567: v = simplify(v); tv = typ(v);
568: if (!is_scalar_t(tu) || !is_scalar_t(tv))
569: err(talker, "incorrect data in findquad");
570: x = v;
571: if (!gcmp0(u)) x = gadd(gmul(u, polx[varn(a)]), x);
572: if (typ(x) == t_POL) x = gmodulcp(x,p);
573: return gerepileupto(av, x);
574: }
575:
576: static GEN
577: findquad_pol(GEN nf, GEN a, GEN x)
578: {
579: long i, lx = lg(x);
580: GEN p = (GEN)nf[1], y = cgetg(lx,t_POL);
581: for (i=2; i<lx; i++) y[i] = (long)findquad(a, (GEN)x[i], p);
582: y[1] = x[1]; return y;
583: }
584:
585: static GEN
586: compocyclo(GEN nf, long m, long d, long prec)
587: {
588: GEN sb,a,b,s,p1,p2,p3,res,polL,polLK,nfL, D = (GEN)nf[3];
589: long ell,vx;
590:
591: p1 = quadhilbertimag(D, gzero);
592: p2 = cyclo(m,0);
593: if (d==1) return do_compo(p1,p2);
594:
595: ell = m&1 ? m : (m>>2);
596: if (!cmpsi(-ell,D)) /* ell = |D| */
597: {
598: p2 = gcoeff(nffactor(nf,p2),1,1);
599: return do_compo(p1,p2);
600: }
601: if (ell%4 == 3) ell = -ell;
602: /* nf = K = Q(a), L = K(b) quadratic extension = Q(t) */
603: polLK = quadpoly(stoi(ell)); /* relative polynomial */
604: res = rnfequation2(nf, polLK);
605: vx = varn(nf[1]);
606: polL = gsubst((GEN)res[1],0,polx[vx]); /* = charpoly(t) */
607: a = gsubst(lift((GEN)res[2]), 0,polx[vx]);
608: b = gsub(polx[vx], gmul((GEN)res[3], a));
609: nfL = initalg(polL,prec);
610: p1 = gcoeff(nffactor(nfL,p1),1,1);
611: p2 = gcoeff(nffactor(nfL,p2),1,1);
612: p3 = do_compo(p1,p2); /* relative equation over L */
613: /* compute non trivial s in Gal(L / K) */
614: sb= gneg(gadd(b, truecoeff(polLK,1))); /* s(b) = Tr(b) - b */
615: s = gadd(polx[vx], gsub(sb, b)); /* s(t) = t + s(b) - b */
616: p3 = gmul(p3, galoisapplypol(nfL, s, p3));
617: return findquad_pol(nf, a, p3);
618: }
619:
620: /* I integral ideal in HNF. (x) = I, x small in Z ? */
621: static long
622: isZ(GEN I)
623: {
624: GEN x = gcoeff(I,1,1);
625: if (signe(gcoeff(I,1,2)) || !egalii(x, gcoeff(I,2,2))) return 0;
626: return is_bigint(x)? -1: itos(x);
627: }
628:
629: /* Treat special cases directly. return NULL if not special case */
630: static GEN
631: treatspecialsigma(GEN nf, GEN gf, int raw, long prec)
632: {
633: GEN p1,p2,tryf, D = (GEN)nf[3];
634: long Ds,fl,i;
635:
636: if (raw)
637: err(talker,"special case in Schertz's theorem. Odd flag meaningless");
638: i = isZ(gf);
639: if (cmpis(D,-3)==0)
640: {
641: if (i == 4 || i == 5 || i == 7) return cyclo(i,0);
642: if (cmpis(gcoeff(gf,1,1), 9) || cmpis(content(gf),3)) return NULL;
643: p1 = (GEN)nfroots(nf,cyclo(3,0))[1]; /* f = P_3^3 */
644: return gadd(gpowgs(polx[0],3), p1); /* x^3+j */
645: }
646: if (cmpis(D,-4)==0)
647: {
648: if (i == 3 || i == 5) return cyclo(i,0);
649: if (i != 4) return NULL;
650: p1 = (GEN)nfroots(nf,cyclo(4,0))[1];
651: return gadd(gpowgs(polx[0],2), p1); /* x^2+i */
652: }
653: Ds = smodis(D,48);
654: if (i)
655: {
656: if (i==2 && Ds%16== 8) return compocyclo(nf, 4,1,prec);
657: if (i==3 && Ds% 3== 1) return compocyclo(nf, 3,1,prec);
658: if (i==4 && Ds% 8== 1) return compocyclo(nf, 4,1,prec);
659: if (i==6 && Ds ==40) return compocyclo(nf,12,1,prec);
660: return NULL;
661: }
662:
1.2 ! noro 663: p1 = gcoeff(gf,1,1); /* integer > 0 */
1.1 noro 664: p2 = gcoeff(gf,2,2);
665: if (gcmp1(p2)) { fl = 0; tryf = p1; }
666: else {
667: if (Ds % 16 != 8 || !egalii(content(gf),gdeux)) return NULL;
668: fl = 1; tryf = shifti(p1,-1);
669: }
1.2 ! noro 670: /* tryf integer > 0 */
! 671: if (cmpis(tryf, 3) <= 0 || signe(resii(D, tryf)) || !isprime(tryf))
1.1 noro 672: return NULL;
673:
674: i = itos(tryf); if (fl) i *= 4;
675: return compocyclo(nf,i,2,prec);
676: }
677:
678: static GEN
679: getallrootsof1(GEN bnf)
680: {
681: GEN z, u, nf = checknf(bnf), racunit = gmael3(bnf,8,4,2);
682: long i, n = itos(gmael3(bnf,8,4,1));
683:
684: u = cgetg(n+1, t_VEC);
685: z = basistoalg(nf, racunit);
686: for (i=1; ; i++)
687: {
688: u[i] = (long)algtobasis(nf,z);
689: if (i == n) return u;
690: z = gmul(z, racunit);
691: }
692: }
693:
694: /* Compute ray class field polynomial using sigma; if raw=1, pairs
695: [ideals,roots] are given instead so that congruence subgroups can be used */
696: static GEN
697: quadrayimagsigma(GEN bnr, int raw, long prec)
698: {
699: GEN allf,bnf,nf,pol,w,f,la,P,labas,lamodf,u;
700: long a,b,f2,i,lu;
701:
702: allf = conductor(bnr,gzero,2); bnr = (GEN)allf[2];
703: f = gmael(allf,1,1);
704: bnf= (GEN)bnr[1];
705: nf = (GEN)bnf[7];
706: pol= (GEN)nf[1];
707: if (gcmp1(gcoeff(f,1,1))) /* f = 1 ? */
708: {
709: P = quadhilbertimag((GEN)nf[3], stoi(raw));
710: if (raw) convert_to_id(P);
711: return gcopy(P);
712: }
713: P = treatspecialsigma(nf,f,raw,prec);
714: if (P) return P;
715:
716: w = gmodulcp(polx[varn(pol)], pol);
717: f2 = 2 * itos(gcoeff(f,1,1));
718: u = getallrootsof1(bnf); lu = lg(u);
719: for (i=1; i<lu; i++)
1.2 ! noro 720: u[i] = (long)colreducemodHNF((GEN)u[i], f, NULL); /* roots of 1, mod f */
1.1 noro 721: if (DEBUGLEVEL>1)
722: fprintferr("quadray: looking for [a,b] != unit mod 2f\n[a,b] = ");
723: for (a=0; a<f2; a++)
724: for (b=0; b<f2; b++)
725: {
726: la = gaddgs(gmulsg(a,w),b);
727: if (smodis(gnorm(la), f2) != 1) continue;
728: if (DEBUGLEVEL>1) fprintferr("[%ld,%ld] ",a,b);
729:
730: labas = algtobasis(nf, la);
1.2 ! noro 731: lamodf = colreducemodHNF(labas, f, NULL);
1.1 noro 732: for (i=1; i<lu; i++)
733: if (gegal(lamodf, (GEN)u[i])) break;
734: if (i < lu) continue; /* la = unit mod f */
735: if (DEBUGLEVEL)
736: {
737: if (DEBUGLEVEL>1) fprintferr("\n");
738: fprintferr("lambda = %Z\n",la);
739: }
740: return computeP2(bnr,labas,raw,prec);
741: }
742: err(bugparier,"quadrayimagsigma");
743: return NULL;
744: }
745:
746: GEN
747: quadray(GEN D, GEN f, GEN flag, long prec)
748: {
749: GEN bnr,y,p1,pol,bnf,lambda;
1.2 ! noro 750: long raw;
! 751: gpmem_t av = avma;
1.1 noro 752:
753: if (!flag) flag = gzero;
754: if (typ(flag)==t_INT) lambda = NULL;
755: else
756: {
757: if (typ(flag)!=t_VEC || lg(flag)!=3) err(flagerr,"quadray");
758: lambda = (GEN)flag[1]; flag = (GEN)flag[2];
759: if (typ(flag)!=t_INT) err(flagerr,"quadray");
760: }
761: if (typ(D)!=t_INT)
762: {
763: bnf = checkbnf(D);
764: if (degpol(gmael(bnf,7,1)) != 2)
765: err(talker,"not a polynomial of degree 2 in quadray");
766: D=gmael(bnf,7,3);
767: }
768: else
769: {
770: if (!isfundamental(D))
771: err(talker,"quadray needs a fundamental discriminant");
772: pol = quadpoly(D); setvarn(pol, fetch_user_var("y"));
773: bnf = bnfinit0(pol,signe(D)>0?1:0,NULL,prec);
774: }
775: raw = (mpodd(flag) && signe(D) < 0);
776: bnr = bnrinit0(bnf,f,1);
777: if (gcmp1(gmael(bnr,5,1)))
778: {
779: avma = av; if (!raw) return polx[0];
780: y = cgetg(2,t_VEC); p1 = cgetg(3,t_VEC); y[1] = (long)p1;
781: p1[1]=(long)idmat(2);
782: p1[2]=(long)polx[0]; return y;
783: }
784: if (signe(D) > 0)
785: y = bnrstark(bnr,gzero, gcmp0(flag)?1:5, prec);
786: else
787: {
788: if (lambda)
789: y = computeP2(bnr,lambda,raw,prec);
790: else
791: y = quadrayimagsigma(bnr,raw,prec);
792: }
793: return gerepileupto(av, y);
794: }
795:
796: /*******************************************************************/
797: /* */
798: /* Routines related to binary quadratic forms (for internal use) */
799: /* */
800: /*******************************************************************/
801: #define rhorealform(x) rhoreal_aux(x,Disc,sqrtD,isqrtD)
802: #define redrealform(x) gcopy(fix_signs(redrealform5(x,Disc,sqrtD,isqrtD)))
1.2 ! noro 803:
! 804: static GEN
! 805: fix_signs(GEN x)
! 806: {
! 807: GEN a = (GEN)x[1];
! 808: GEN c = (GEN)x[3];
! 809: if (signe(a) < 0)
! 810: {
! 811: if (narrow || absi_equal(a,c)) return rhorealform(x);
! 812: setsigne(a,1); setsigne(c,-1);
! 813: }
! 814: return x;
! 815: }
! 816:
! 817: static GEN
! 818: comprealform(GEN x,GEN y) {
! 819: return fix_signs( comprealform5(x,y,Disc,sqrtD,isqrtD) );
! 820: }
1.1 noro 821:
822: /* compute rho^n(x) */
1.2 ! noro 823: static GEN
1.1 noro 824: rhoreal_pow(GEN x, long n)
825: {
1.2 ! noro 826: long i;
! 827: gpmem_t av = avma, lim = stack_lim(av, 1);
1.1 noro 828: for (i=1; i<=n; i++)
829: {
830: x = rhorealform(x);
831: if (low_stack(lim, stack_lim(av,1)))
832: {
1.2 ! noro 833: if(DEBUGMEM>1) err(warnmem,"rhoreal_pow");
1.1 noro 834: x = gerepilecopy(av, x);
835: }
836: }
837: return gerepilecopy(av, x);
838: }
839:
840: static GEN
1.2 ! noro 841: realpf5(GEN D, long p)
1.1 noro 842: {
1.2 ! noro 843: gpmem_t av = avma;
! 844: GEN y = primeform(D,stoi(p),PRECREG);
! 845: y = codeform5(y,PRECREG);
! 846: return gerepileupto(av, redrealform(y));
1.1 noro 847: }
848:
849: static GEN
1.2 ! noro 850: realpf(GEN D, long p)
1.1 noro 851: {
1.2 ! noro 852: gpmem_t av = avma;
! 853: GEN y = primeform(D,stoi(p),PRECREG);
1.1 noro 854: return gerepileupto(av, redrealform(y));
855: }
856:
857: static GEN
1.2 ! noro 858: imagpf(GEN D, long p) { return primeform(D,stoi(p),0); }
1.1 noro 859:
860: static GEN
861: comprealform3(GEN x, GEN y)
862: {
1.2 ! noro 863: gpmem_t av = avma;
1.1 noro 864: GEN z = cgetg(4,t_VEC); comp_gen(z,x,y);
865: return gerepileupto(av, redrealform(z));
866: }
867:
1.2 ! noro 868: /* Warning: ex[0] not set */
1.1 noro 869: static GEN
1.2 ! noro 870: init_form(long *ex, GEN (*comp)(GEN,GEN))
1.1 noro 871: {
1.2 ! noro 872: long i, l = lg(powsubFB);
! 873: GEN F = NULL;
! 874: for (i=1; i<l; i++)
! 875: if (ex[i])
! 876: {
! 877: GEN t = gmael(powsubFB,i,ex[i]);
! 878: F = F? comp(F,t): t;
! 879: }
! 880: return F;
! 881: }
! 882: static GEN
! 883: initrealform5(long *ex) { return init_form(ex, &comprealform); }
! 884: static GEN
! 885: initimagform(long *ex) { return init_form(ex, &compimag); }
1.1 noro 886:
1.2 ! noro 887: static GEN
! 888: random_form(GEN ex, GEN (*comp)(GEN,GEN))
! 889: {
! 890: long i, l = lg(ex);
! 891: gpmem_t av = avma;
! 892: GEN F;
! 893: for(;;)
! 894: {
! 895: for (i=1; i<l; i++) ex[i] = mymyrand()>>randshift;
! 896: if ((F = init_form(ex, comp))) return F;
! 897: avma = av;
! 898: }
1.1 noro 899: }
900:
1.2 ! noro 901: static GEN
! 902: real_random_form(GEN ex) { return random_form(ex, &comprealform3); }
! 903: static GEN
! 904: imag_random_form(GEN ex) { return random_form(ex, &compimag); }
! 905:
1.1 noro 906: /*******************************************************************/
907: /* */
908: /* Common subroutines */
909: /* */
910: /*******************************************************************/
911: static void
912: buch_init(void)
913: {
1.2 ! noro 914: if (DEBUGLEVEL) (void)timer2();
1.1 noro 915: primfact = new_chunk(100);
916: exprimfact = new_chunk(100);
917: badprim = new_chunk(100);
918: hashtab = (long**) new_chunk(HASHT);
919: }
920:
921: double
922: check_bach(double cbach, double B)
923: {
924: if (cbach >= B)
925: err(talker,"sorry, buchxxx couldn't deal with this field PLEASE REPORT!");
926: cbach *= 2; if (cbach > B) cbach = B;
927: if (DEBUGLEVEL) fprintferr("\n*** Bach constant: %f\n", cbach);
928: return cbach;
929: }
930:
931: static long
1.2 ! noro 932: factorquad(GEN f, long kcz, long limp)
1.1 noro 933: {
1.2 ! noro 934: long i, p, k, lo;
! 935: gpmem_t av;
1.1 noro 936: GEN q,r, x = (GEN)f[1];
937:
938: if (is_pm1(x)) { primfact[0]=0; return 1; }
939: av=avma; lo=0;
940: if (signe(x) < 0) x = absi(x);
941: for (i=1; ; i++)
942: {
1.2 ! noro 943: p=FB[i]; q=dvmdis(x,p,&r);
1.1 noro 944: if (!signe(r))
945: {
1.2 ! noro 946: for (k=0; !signe(r); k++) { x=q; q=dvmdis(x,p,&r); }
! 947: primfact[++lo]=p; exprimfact[lo]=k;
1.1 noro 948: }
949: if (cmpis(q,p)<=0) break;
950: if (i==kcz) { avma=av; return 0; }
951: }
1.2 ! noro 952: avma=av;
! 953: if (lgefint(x)!=3 || (p=x[2]) > limhash) return 0;
1.1 noro 954:
955: if (p != 1 && p <= limp)
956: {
957: for (i=1; i<=badprim[0]; i++)
958: if (p % badprim[i] == 0) return 0;
1.2 ! noro 959: primfact[++lo]=p; exprimfact[lo]=1;
1.1 noro 960: p = 1;
961: }
962: primfact[0]=lo; return p;
963: }
964:
965: /* q may not be prime, but check for a "large prime relation" involving q */
966: static long *
967: largeprime(long q, long *ex, long np, long nrho)
968: {
1.2 ! noro 969: const long hashv = hash(q);
! 970: long *pt, i, l = lg(subFB);
1.1 noro 971:
972: for (pt = hashtab[hashv]; ; pt = (long*) pt[0])
973: {
974: if (!pt)
975: {
1.2 ! noro 976: pt = (long*) gpmalloc((l+3)<<TWOPOTBYTES_IN_LONG);
1.1 noro 977: *pt++ = nrho; /* nrho = pt[-3] */
978: *pt++ = np; /* np = pt[-2] */
979: *pt++ = q; /* q = pt[-1] */
980: pt[0] = (long)hashtab[hashv];
1.2 ! noro 981: for (i=1; i<l; i++) pt[i]=ex[i];
1.1 noro 982: hashtab[hashv]=pt; return NULL;
983: }
984: if (pt[-1] == q) break;
985: }
1.2 ! noro 986: for(i=1; i<l; i++)
1.1 noro 987: if (pt[i] != ex[i]) return pt;
988: return (pt[-2]==np)? (GEN)NULL: pt;
989: }
990:
1.2 ! noro 991: /* p | conductor of order of disc D ? */
! 992: static int
! 993: is_bad(GEN D, ulong p)
1.1 noro 994: {
1.2 ! noro 995: gpmem_t av = avma;
! 996: int r;
! 997: if (p == 2)
! 998: {
! 999: r = mod16(D) >> 1;
! 1000: if (r && signe(D) < 0) r = 8-r;
! 1001: return (r < 4);
! 1002: }
! 1003: r = (resii(D, muluu(p,p)) == gzero); /* p^2 | D ? */
! 1004: avma = av; return r;
1.1 noro 1005: }
1006:
1.2 ! noro 1007: /* create FB, numFB; fill badprim. Return L(kro_D, 1) */
! 1008: static GEN
! 1009: FBquad(GEN Disc, long n2, long n)
1.1 noro 1010: {
1.2 ! noro 1011: GEN Res = realun(DEFAULTPREC);
! 1012: long i, p, bad, s;
! 1013: gpmem_t av;
! 1014: byteptr d = diffptr;
1.1 noro 1015:
1.2 ! noro 1016: numFB = cgetg(n2+1, t_VECSMALL);
! 1017: FB = cgetg(n2+1, t_VECSMALL);
! 1018: av = avma;
! 1019: KC = 0; bad = 0; i = 0;
! 1020: if (maxprime() < n2) err(primer1);
! 1021: for (p = 0;;) /* p <= n2 */
! 1022: {
! 1023: NEXT_PRIME_VIADIFF(p, d);
! 1024: if (!KC && p > n) KC = i;
! 1025: if (p > n2) break;
! 1026: s = krogs(Disc,p);
! 1027: switch (s)
1.1 noro 1028: {
1029: case -1: break; /* inert */
1030: case 0: /* ramified */
1.2 ! noro 1031: if (is_bad(Disc, (ulong)p)) { badprim[++bad]=p; break; }
! 1032: /* fall through */
1.1 noro 1033: default: /* split */
1.2 ! noro 1034: i++; numFB[p]=i; FB[i] = p; break;
1.1 noro 1035: }
1.2 ! noro 1036: Res = mulsr(p, divrs(Res, p - s));
1.1 noro 1037: }
1.2 ! noro 1038: if (!KC) return NULL;
! 1039: Res = gerepileupto(av, Res);
1.1 noro 1040: KC2 = i;
1.2 ! noro 1041: setlg(FB, KC2+1);
1.1 noro 1042: if (DEBUGLEVEL)
1043: {
1044: msgtimer("factor base");
1045: if (DEBUGLEVEL>7)
1046: {
1.2 ! noro 1047: fprintferr("FB:\n");
! 1048: for (i=1; i<=KC; i++) fprintferr("%ld ", FB[i]);
! 1049: fprintferr("\n");
1.1 noro 1050: }
1051: }
1.2 ! noro 1052: badprim[0] = bad; return Res;
1.1 noro 1053: }
1054:
1.2 ! noro 1055: /* create vperm, return subFB */
! 1056: static GEN
! 1057: subFBquad(GEN D, double PROD, long KC)
1.1 noro 1058: {
1.2 ! noro 1059: long i, j, lgsub, ino, lv = KC+1;
! 1060: double prod = 1.;
! 1061: gpmem_t av;
! 1062: GEN no;
1.1 noro 1063:
1.2 ! noro 1064: vperm = cgetg(lv, t_VECSMALL);
! 1065: av = avma;
! 1066: no = cgetg(lv, t_VECSMALL); ino = 1;
! 1067: for (i=j=1; j < lv; j++)
1.1 noro 1068: {
1.2 ! noro 1069: long p = FB[j];
! 1070: if (smodis(D, p) == 0) no[ino++] = j; /* ramified */
! 1071: else
! 1072: {
! 1073: vperm[i] = j; i++;
! 1074: prod *= p;
! 1075: if (prod > PROD) break;
! 1076: }
1.1 noro 1077: }
1.2 ! noro 1078: if (j == lv) return NULL;
! 1079: lgsub = i;
! 1080: for (j = 1; j <ino; i++,j++) vperm[i] = no[j];
! 1081: for ( ; i < lv; i++) vperm[i] = i;
! 1082: avma = av;
! 1083: if (DEBUGLEVEL) msgtimer("subFBquad (%ld elt.)", lgsub-1);
! 1084: return vecextract_i(vperm, 1, lgsub-1);
1.1 noro 1085: }
1086:
1.2 ! noro 1087: /* assume n >= 1, x[i][j] = subFB[i]^j, for j = 1..n */
! 1088: static GEN
! 1089: powsubFBquad(long n)
1.1 noro 1090: {
1.2 ! noro 1091: long i,j, l = lg(subFB);
! 1092: GEN F, y, x = cgetg(l, t_VEC);
1.1 noro 1093:
1094: if (PRECREG) /* real */
1095: {
1.2 ! noro 1096: for (i=1; i<l; i++)
! 1097: {
! 1098: F = realpf5(Disc, FB[subFB[i]]);
! 1099: y = cgetg(n+1, t_VEC); x[i] = (long)y;
! 1100: y[1] = (long)F;
! 1101: for (j=2; j<=n; j++) y[j] = (long)comprealform((GEN)y[j-1], F);
1.1 noro 1102: }
1103: }
1104: else /* imaginary */
1105: {
1.2 ! noro 1106: for (i=1; i<l; i++)
1.1 noro 1107: {
1.2 ! noro 1108: F = imagpf(Disc, FB[subFB[i]]);
! 1109: y = cgetg(n+1, t_VEC); x[i] = (long)y;
! 1110: y[1] = (long)F;
! 1111: for (j=2; j<=n; j++) y[j] = (long)compimag((GEN)y[j-1], F);
1.1 noro 1112: }
1113: }
1.2 ! noro 1114: if (DEBUGLEVEL) msgtimer("powsubFBquad");
! 1115: return x;
1.1 noro 1116: }
1117:
1118: static void
1119: desalloc(long **mat)
1120: {
1121: long i,*p,*q;
1.2 ! noro 1122: for (i=1; i<lg(mat); i++) free(mat[i]);
! 1123: free(mat);
! 1124: for (i=1; i<HASHT; i++)
! 1125: for (p = hashtab[i]; p; p = q) { q=(long*)p[0]; free(p-3); }
! 1126: }
1.1 noro 1127:
1.2 ! noro 1128: static void
! 1129: sub_fact(GEN col, GEN F)
! 1130: {
! 1131: GEN b = (GEN)F[2];
! 1132: long i, p, e;
! 1133: for (i=1; i<=primfact[0]; i++)
! 1134: {
! 1135: p = primfact[i]; e = exprimfact[i];
! 1136: if (smodis(b, p<<1) > p) e = -e;
! 1137: col[numFB[p]] -= e;
1.1 noro 1138: }
1139: }
1.2 ! noro 1140: static void
! 1141: add_fact(GEN col, GEN F)
1.1 noro 1142: {
1.2 ! noro 1143: GEN b = (GEN)F[2];
! 1144: long i, p, e;
! 1145: for (i=1; i<=primfact[0]; i++)
! 1146: {
! 1147: p = primfact[i]; e = exprimfact[i];
! 1148: if (smodis(b, p<<1) > p) e = -e;
! 1149: col[numFB[p]] += e;
1.1 noro 1150: }
1151: }
1152:
1153: #define comp(x,y) x? (PRECREG? compreal(x,y): compimag(x,y)): y
1154: static GEN
1155: get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec)
1156: {
1.2 ! noro 1157: GEN res,*init, u1, met = smithrel(W,NULL,&u1);
! 1158: long c,i,j, l = lg(W);
1.1 noro 1159:
1.2 ! noro 1160: c = lg(met);
1.1 noro 1161: if (DEBUGLEVEL) msgtimer("smith/class group");
1162: res=cgetg(c,t_VEC); init = (GEN*)cgetg(l,t_VEC);
1163: for (i=1; i<l; i++)
1.2 ! noro 1164: init[i] = primeform(Disc,stoi(FB[vperm[i]]),prec);
1.1 noro 1165: for (j=1; j<c; j++)
1166: {
1.2 ! noro 1167: GEN p1 = NULL;
1.1 noro 1168: for (i=1; i<l; i++)
1.2 ! noro 1169: p1 = comp(p1, powgi(init[i], gcoeff(u1,i,j)));
1.1 noro 1170: res[j] = (long)p1;
1171: }
1172: if (DEBUGLEVEL) msgtimer("generators");
1173: *ptmet = met; return res;
1174: }
1175:
1176: static GEN
1.2 ! noro 1177: extra_relations(long LIMC, long nlze, GEN *ptextraC)
1.1 noro 1178: {
1.2 ! noro 1179: long fpc, i, k, nlze2;
! 1180: long s = 0, extrarel = nlze+2, lgsub = lg(subFB);
! 1181: gpmem_t av;
1.1 noro 1182: long MAXRELSUP = min(50,4*KC);
1.2 ! noro 1183: GEN p1, form, ex, extramat, extraC, col;
1.1 noro 1184:
1.2 ! noro 1185: extramat = cgetg(extrarel+1, t_MAT);
! 1186: extraC = cgetg(extrarel+1, t_VEC);
! 1187: for (i=1; i<=extrarel; i++) extramat[i] = lgetg(KC+1,t_VECSMALL);
! 1188: if (!PRECREG)
! 1189: for (i=1; i<=extrarel; i++) extraC[i] = lgetg(1,t_COL);
! 1190:
! 1191: if (DEBUGLEVEL) fprintferr("looking for %ld extra relations\n",extrarel);
! 1192: nlze2 = PRECREG? max(nlze,lgsub-1): min(nlze+1,KC);
1.1 noro 1193: if (nlze2 < 3 && KC > 2) nlze2 = 3;
1.2 ! noro 1194: ex = cgetg(nlze2+1, t_VECSMALL);
1.1 noro 1195: av = avma;
1196: while (s<extrarel)
1197: {
1198: form = NULL;
1199: for (i=1; i<=nlze2; i++)
1200: {
1.2 ! noro 1201: ex[i] = mymyrand()>>randshift;
1.1 noro 1202: if (ex[i])
1203: {
1.2 ! noro 1204: p1 = primeform(Disc,stoi(FB[vperm[i]]),PRECREG);
! 1205: p1 = gpowgs(p1,ex[i]); form = comp(form,p1);
1.1 noro 1206: }
1207: }
1208: if (!form) continue;
1209:
1.2 ! noro 1210: fpc = factorquad(form,KC,LIMC);
1.1 noro 1211: if (fpc==1)
1212: {
1213: s++; col = (GEN)extramat[s];
1214: for (i=1; i<=nlze2; i++) col[vperm[i]] = -ex[i];
1.2 ! noro 1215: for ( ; i<=KC; i++) col[vperm[i]] = 0;
! 1216: add_fact(col, form);
1.1 noro 1217: for (i=1; i<=KC; i++)
1218: if (col[i]) break;
1219: if (i>KC)
1220: {
1221: s--; avma = av;
1222: if (--MAXRELSUP == 0) return NULL;
1223: }
1.2 ! noro 1224: else { av = avma; if (PRECREG) extraC[s] = form[4]; }
1.1 noro 1225: }
1226: else avma = av;
1227: if (DEBUGLEVEL)
1228: {
1229: if (fpc == 1) fprintferr(" %ld",s);
1230: else if (DEBUGLEVEL>1) fprintferr(".");
1231: }
1232: }
1.2 ! noro 1233:
! 1234: for (i=1; i<=extrarel; i++)
1.1 noro 1235: {
1.2 ! noro 1236: GEN colg = cgetg(KC+1,t_COL);
! 1237: col = (GEN)extramat[i];
! 1238: extramat[i] = (long)colg;
! 1239: for (k=1; k<=KC; k++) colg[k] = lstoi(col[vperm[k]]);
1.1 noro 1240: }
1.2 ! noro 1241: if (DEBUGLEVEL) { fprintferr("\n"); msgtimer("extra relations"); }
! 1242: *ptextraC = extraC; return extramat;
! 1243: }
! 1244: #undef comp
! 1245:
! 1246: static long
! 1247: trivial_relations(long **mat, GEN C)
! 1248: {
! 1249: long i, s = 0;
! 1250: GEN col;
! 1251: for (i=1; i<=KC; i++)
! 1252: {
! 1253: if (smodis(Disc, FB[i])) continue;
! 1254: /* ramified prime ==> trivial relation */
! 1255: col = mat[++s];
! 1256: col[i] = 2;
! 1257: if (C) affsr(0, (GEN)C[s]);
1.1 noro 1258: }
1.2 ! noro 1259: return s;
! 1260: }
! 1261:
! 1262: static void
! 1263: dbg_all(long phase, long s, long n)
! 1264: {
! 1265: if (DEBUGLEVEL>1) fprintferr("\n");
! 1266: msgtimer("%s rel [#rel/#test = %ld/%ld]", phase? "random": "initial", s, n);
! 1267: }
! 1268:
! 1269: void
! 1270: wr_rel(GEN col)
! 1271: {
! 1272: long i, l = lg(col);
! 1273: fprintferr("\nrel = ");
! 1274: for (i=1; i<l; i++)
! 1275: if (col[i]) fprintferr("%ld^%ld ",i,col[i]);
! 1276: fprintferr("\n");
1.1 noro 1277: }
1278:
1.2 ! noro 1279: void
! 1280: dbg_rel(long s, GEN col)
! 1281: {
! 1282: if (DEBUGLEVEL == 1) fprintferr("%4ld",s);
! 1283: else { fprintferr("cglob = %ld. ", s); wr_rel(col); }
! 1284: flusherr();
! 1285: }
1.1 noro 1286: /*******************************************************************/
1287: /* */
1288: /* Imaginary Quadratic fields */
1289: /* */
1290: /*******************************************************************/
1291:
1.2 ! noro 1292: /* Strategy 1 until lim relation found, then Strategy 2 until LIM. */
1.1 noro 1293: static GEN
1.2 ! noro 1294: imag_relations(long LIM, long lim, long LIMC, long **mat)
1.1 noro 1295: {
1.2 ! noro 1296: long lgsub = lg(subFB), i, s, fpc, current, nbtest = 0;
! 1297: int first = 1;
! 1298: gpmem_t av;
! 1299: GEN C, col, form, ex = cgetg(lgsub, t_VECSMALL);
! 1300: GEN dummy = cgetg(1,t_COL);
1.1 noro 1301:
1.2 ! noro 1302: C = cgetg(LIM+1, t_VEC);
! 1303: for (i=1; i<=LIM; i++) C[i] = (long)dummy;
! 1304: av = avma;
! 1305: s = trivial_relations(mat, NULL);
! 1306: lim += s; if (lim > LIM) lim = LIM;
1.1 noro 1307: for(;;)
1308: {
1.2 ! noro 1309: if (s >= lim) {
! 1310: if (s >= LIM) break;
! 1311: lim = LIM; first = 0;
! 1312: if (DEBUGLEVEL) dbg_all(0, s, nbtest);
! 1313: }
! 1314: avma = av; current = first? 1+(s%KC): 1+s-RELSUP;
! 1315: form = imagpf(Disc, FB[current]);
! 1316: form = compimag(form, imag_random_form(ex));
! 1317: nbtest++; fpc = factorquad(form,KC,LIMC);
1.1 noro 1318: if (!fpc)
1319: {
1.2 ! noro 1320: if (DEBUGLEVEL>1) fprintferr(".");
1.1 noro 1321: continue;
1322: }
1323: if (fpc > 1)
1324: {
1.2 ! noro 1325: long *fpd = largeprime(fpc,ex,current,0);
! 1326: long b1, b2, p;
! 1327: GEN form2;
1.1 noro 1328: if (!fpd)
1329: {
1.2 ! noro 1330: if (DEBUGLEVEL>1) fprintferr(".");
1.1 noro 1331: continue;
1332: }
1.2 ! noro 1333: form2 = initimagform(fpd);
! 1334: form2 = compimag(form2, imagpf(Disc, FB[fpd[-2]]));
! 1335: p = fpc << 1;
! 1336: b1 = smodis((GEN)form2[2], p);
! 1337: b2 = smodis((GEN)form[2], p);
! 1338: if (b1 != b2 && b1+b2 != p) continue;
! 1339:
! 1340: col = mat[++s];
! 1341: add_fact(col, form);
! 1342: (void)factorquad(form2,KC,LIMC);
1.1 noro 1343: if (b1==b2)
1344: {
1.2 ! noro 1345: for (i=1; i<lgsub; i++) col[subFB[i]] += fpd[i]-ex[i];
! 1346: sub_fact(col, form2); col[fpd[-2]]++;
1.1 noro 1347: }
1348: else
1349: {
1.2 ! noro 1350: for (i=1; i<lgsub; i++) col[subFB[i]] += -fpd[i]-ex[i];
! 1351: add_fact(col, form2); col[fpd[-2]]--;
1.1 noro 1352: }
1.2 ! noro 1353: }
1.1 noro 1354: else
1355: {
1.2 ! noro 1356: col = mat[++s];
! 1357: for (i=1; i<lgsub; i++) col[subFB[i]] = -ex[i];
! 1358: add_fact(col, form);
1.1 noro 1359: }
1360: col[current]--;
1.2 ! noro 1361: if (DEBUGLEVEL) dbg_rel(s, col);
1.1 noro 1362: if (!first && fpc == 1 && col[current] == 0)
1363: {
1364: s--; for (i=1; i<=KC; i++) col[i]=0;
1365: }
1366: }
1.2 ! noro 1367: if (DEBUGLEVEL) dbg_all(1, s, nbtest);
! 1368: return C;
1.1 noro 1369: }
1370:
1371: static int
1.2 ! noro 1372: imag_be_honest()
1.1 noro 1373: {
1.2 ! noro 1374: long p, fpc, s = KC, nbtest = 0;
! 1375: gpmem_t av = avma;
! 1376: GEN F, ex = cgetg(lg(subFB), t_VECSMALL);
1.1 noro 1377:
1378: while (s<KC2)
1379: {
1.2 ! noro 1380: p = FB[s+1]; if (DEBUGLEVEL) fprintferr(" %ld",p);
! 1381: F = compimag(imagpf(Disc, p), imag_random_form(ex));
! 1382: fpc = factorquad(F,s,p-1);
1.1 noro 1383: if (fpc == 1) { nbtest=0; s++; }
1384: else
1385: if (++nbtest > 20) return 0;
1386: avma = av;
1387: }
1388: return 1;
1389: }
1390:
1391: /*******************************************************************/
1392: /* */
1393: /* Real Quadratic fields */
1394: /* */
1395: /*******************************************************************/
1396:
1.2 ! noro 1397: /* (1/2) log (d * 2^{e * EXP220}) */
! 1398: GEN
! 1399: get_dist(GEN e, GEN d, long prec)
! 1400: {
! 1401: GEN t = mplog(absr(d));
! 1402: if (signe(e)) t = addrr(t, mulir(mulsi(EXP220,e), mplog2(prec)));
! 1403: return shiftr(t, -1);
! 1404: }
! 1405:
1.1 noro 1406: static GEN
1.2 ! noro 1407: real_relations(long LIM, long lim, long LIMC, long **mat)
1.1 noro 1408: {
1.2 ! noro 1409: long lgsub = lg(subFB), i, s, fpc, current, nbtest = 0, endcycle, rhoacc, rho;
! 1410: int first = 1;
! 1411: gpmem_t av, av1, limstack;
! 1412: GEN C, d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL);
1.1 noro 1413:
1.2 ! noro 1414: C = cgetg(LIM+1, t_VEC);
! 1415: for (i=1; i<=LIM; i++) C[i] = lgetr(PRECREG);
! 1416:
! 1417: current = 0;
! 1418: av = avma; limstack = stack_lim(av,1);
! 1419: s = trivial_relations(mat, C);
! 1420: lim += s; if (lim > LIM) lim = LIM;
! 1421: NEW:
1.1 noro 1422: for(;;)
1423: {
1.2 ! noro 1424: if (s >= lim) {
! 1425: if (lim == LIM) break;
! 1426: lim = LIM; first = 0;
! 1427: if (DEBUGLEVEL) dbg_all(0, s, nbtest);
1.1 noro 1428: }
1429: avma = av;
1430: form = real_random_form(ex);
1431: if (!first)
1432: {
1433: current = 1+s-RELSUP;
1.2 ! noro 1434: form = comprealform3(form, realpf(Disc, FB[current]));
1.1 noro 1435: }
1.2 ! noro 1436: av1 = avma;
1.1 noro 1437: form0 = form; form1 = NULL;
1.2 ! noro 1438: endcycle = rhoacc = 0;
! 1439: rho = -1;
! 1440:
! 1441: CYCLE:
! 1442: if (endcycle) goto NEW;
! 1443: if (low_stack(limstack, stack_lim(av,1)))
! 1444: {
! 1445: if(DEBUGMEM>1) err(warnmem,"real_relations");
! 1446: gerepileall(av1, form1? 2: 1, &form, &form1);
! 1447: }
! 1448: if (rho < 0) rho = 0; /* first time in */
! 1449: else
1.1 noro 1450: {
1.2 ! noro 1451: form = rhorealform(form); rho++;
! 1452: rhoacc++;
! 1453: if (first)
! 1454: endcycle = (absi_equal((GEN)form[1],(GEN)form0[1])
! 1455: && egalii((GEN)form[2],(GEN)form0[2])
! 1456: && (!narrow || signe(form0[1])==signe(form[1])));
1.1 noro 1457: else
1458: {
1.2 ! noro 1459: if (narrow)
! 1460: { form = rhorealform(form); rho++; }
! 1461: else if (absi_equal((GEN)form[1], (GEN)form[3])) /* a = -c */
1.1 noro 1462: {
1.2 ! noro 1463: if (absi_equal((GEN)form[1],(GEN)form0[1]) &&
! 1464: egalii((GEN)form[2],(GEN)form0[2])) goto NEW;
! 1465: form = rhorealform(form); rho++;
1.1 noro 1466: }
1467: else
1.2 ! noro 1468: { setsigne(form[1],1); setsigne(form[3],-1); }
! 1469: if (egalii((GEN)form[1],(GEN)form0[1]) &&
! 1470: egalii((GEN)form[2],(GEN)form0[2])) goto NEW;
1.1 noro 1471: }
1.2 ! noro 1472: }
! 1473: nbtest++; fpc = factorquad(form,KC,LIMC);
! 1474: if (!fpc)
! 1475: {
! 1476: if (DEBUGLEVEL>1) fprintferr(".");
! 1477: goto CYCLE;
! 1478: }
! 1479: if (fpc > 1)
! 1480: { /* look for Large Prime relation */
! 1481: long *fpd = largeprime(fpc,ex,current,rhoacc);
! 1482: long b1, b2, p;
! 1483: GEN form2;
! 1484: if (!fpd)
1.1 noro 1485: {
1.2 ! noro 1486: if (DEBUGLEVEL>1) fprintferr(".");
! 1487: goto CYCLE;
1.1 noro 1488: }
1.2 ! noro 1489: if (!form1) form1 = initrealform5(ex);
! 1490: if (!first) form1 = comprealform(form1, realpf5(Disc, FB[current]));
! 1491: form1 = rhoreal_pow(form1, rho);
! 1492: rho = 0;
! 1493:
! 1494: form2 = initrealform5(fpd);
! 1495: if (fpd[-2]) form2 = comprealform(form2, realpf5(Disc, FB[fpd[-2]]));
! 1496: form2 = rhoreal_pow(form2, fpd[-3]);
! 1497: if (!narrow && !absi_equal((GEN)form2[1],(GEN)form2[3]))
! 1498: {
! 1499: setsigne(form2[1],1);
! 1500: setsigne(form2[3],-1);
! 1501: }
! 1502: p = fpc << 1;
! 1503: b1 = smodis((GEN)form2[2], p);
! 1504: b2 = smodis((GEN)form1[2], p);
! 1505: if (b1 != b2 && b1+b2 != p) goto CYCLE;
! 1506:
! 1507: col = mat[++s];
! 1508: add_fact(col, form1);
! 1509: (void)factorquad(form2,KC,LIMC);
! 1510: if (b1==b2)
1.1 noro 1511: {
1.2 ! noro 1512: for (i=1; i<lgsub; i++) col[subFB[i]] += fpd[i]-ex[i];
! 1513: sub_fact(col, form2);
! 1514: if (fpd[-2]) col[fpd[-2]]++; /* implies !first */
! 1515: d = get_dist(subii((GEN)form1[4],(GEN)form2[4]),
! 1516: divrr((GEN)form1[5],(GEN)form2[5]), PRECREG);
1.1 noro 1517: }
1518: else
1519: {
1.2 ! noro 1520: for (i=1; i<lgsub; i++) col[subFB[i]] += -fpd[i]-ex[i];
! 1521: add_fact(col, form2);
! 1522: if (fpd[-2]) col[fpd[-2]]--;
! 1523: d = get_dist(addii((GEN)form1[4],(GEN)form2[4]),
! 1524: mulrr((GEN)form1[5],(GEN)form2[5]), PRECREG);
1.1 noro 1525: }
1.2 ! noro 1526: }
! 1527: else
! 1528: { /* standard relation */
! 1529: if (!form1) form1 = initrealform5(ex);
! 1530: if (!first) form1 = comprealform(form1, realpf5(Disc, FB[current]));
! 1531: form1 = rhoreal_pow(form1,rho);
! 1532: rho = 0;
! 1533:
! 1534: col = mat[++s];
! 1535: for (i=1; i<lgsub; i++) col[subFB[i]] = -ex[i];
! 1536: add_fact(col, form1);
! 1537: d = get_dist((GEN)form1[4], (GEN)form1[5], PRECREG);
! 1538: }
! 1539: if (DEBUGLEVEL) fprintferr(" %ld",s);
! 1540: affrr(d, (GEN)C[s]);
! 1541: if (first)
! 1542: {
! 1543: if (s >= lim) goto NEW;
! 1544: goto CYCLE;
! 1545: }
! 1546: else
! 1547: {
! 1548: col[current]--;
! 1549: if (fpc == 1 && col[current] == 0)
1.1 noro 1550: {
1.2 ! noro 1551: s--; for (i=1; i<=KC; i++) col[i]=0;
1.1 noro 1552: }
1553: }
1554: }
1.2 ! noro 1555: if (DEBUGLEVEL) dbg_all(1, s, nbtest);
! 1556: return C;
1.1 noro 1557: }
1558:
1559: static int
1.2 ! noro 1560: real_be_honest()
1.1 noro 1561: {
1.2 ! noro 1562: long p, fpc, s = KC, nbtest = 0;
! 1563: gpmem_t av = avma;
! 1564: GEN F,F0, ex = cgetg(lg(subFB), t_VECSMALL);
1.1 noro 1565:
1566: while (s<KC2)
1567: {
1.2 ! noro 1568: p = FB[s+1]; if (DEBUGLEVEL) fprintferr(" %ld",p);
! 1569: F = comprealform3(realpf(Disc, p), real_random_form(ex));
! 1570: for (F0 = F;;)
1.1 noro 1571: {
1.2 ! noro 1572: fpc = factorquad(F,s,p-1);
1.1 noro 1573: if (fpc == 1) { nbtest=0; s++; break; }
1574: if (++nbtest > 20) return 0;
1.2 ! noro 1575: F = fix_signs( rhorealform(F) );
! 1576: if (egalii((GEN)F[1],(GEN)F0[1])
! 1577: && egalii((GEN)F[2],(GEN)F0[2])) break;
1.1 noro 1578: }
1579: avma = av;
1580: }
1581: return 1;
1582: }
1583:
1584: static GEN
1.2 ! noro 1585: gcdreal(GEN a,GEN b)
1.1 noro 1586: {
1587: long e;
1588: GEN k1,r;
1589:
1.2 ! noro 1590: if (!signe(a)) return mpabs(b);
! 1591: if (!signe(b)) return mpabs(a);
! 1592:
1.1 noro 1593: if (typ(a)==t_INT)
1594: {
1595: if (typ(b)==t_INT) return mppgcd(a,b);
1.2 ! noro 1596: a = itor(a, lg(b));
1.1 noro 1597: }
1598: else if (typ(b)==t_INT)
1.2 ! noro 1599: {
! 1600: b = itor(b, lg(a));
! 1601: }
1.1 noro 1602: if (expo(a)<-5) return absr(b);
1603: if (expo(b)<-5) return absr(a);
1604: a=absr(a); b=absr(b);
1605: while (expo(b) >= -5 && signe(b))
1606: {
1.2 ! noro 1607: k1 = gcvtoi(divrr(a,b),&e);
1.1 noro 1608: if (e > 0) return NULL;
1609: r=subrr(a,mulir(k1,b)); a=b; b=r;
1610: }
1.2 ! noro 1611: return absr(a);
1.1 noro 1612: }
1613:
1.2 ! noro 1614: enum { RELAT, LARGE, PRECI };
! 1615:
! 1616: static int
! 1617: get_R(GEN C, long sreg, GEN z, GEN *ptR)
1.1 noro 1618: {
1.2 ! noro 1619: GEN R = gun;
! 1620: double c;
! 1621: long i;
1.1 noro 1622:
1.2 ! noro 1623: if (PRECREG)
1.1 noro 1624: {
1.2 ! noro 1625: R = mpabs((GEN)C[1]);
! 1626: for (i=2; i<=sreg; i++)
! 1627: {
! 1628: R = gcdreal((GEN)C[i], R);
! 1629: if (!R) return PRECI;
! 1630: }
! 1631: if (DEBUGLEVEL)
! 1632: {
! 1633: if (DEBUGLEVEL>7) fprintferr("R = %Z",R);
! 1634: msgtimer("regulator");
! 1635: }
! 1636: if (gexpo(R) <= -3)
! 1637: {
! 1638: if (DEBUGLEVEL) fprintferr("regulator is zero.\n");
! 1639: return RELAT;
! 1640: }
1.1 noro 1641: }
1.2 ! noro 1642: c = gtodouble(gmul(z, R));
! 1643: if (c < 0.75 || c > 1.5) return RELAT;
! 1644: *ptR = R; return LARGE;
! 1645: }
! 1646:
! 1647: static int
! 1648: quad_be_honest()
! 1649: {
! 1650: int r;
! 1651: if (KC2 <= KC) return 1;
1.1 noro 1652: if (DEBUGLEVEL)
1.2 ! noro 1653: fprintferr("be honest for primes from %ld to %ld\n", FB[KC+1],FB[KC2]);
! 1654: r = PRECREG? real_be_honest(): imag_be_honest();
! 1655: if (DEBUGLEVEL) { fprintferr("\n"); msgtimer("be honest"); }
! 1656: return r;
1.1 noro 1657: }
1658:
1659: GEN
1660: buchquad(GEN D, double cbach, double cbach2, long RELSUP0, long flag, long prec)
1661: {
1.2 ! noro 1662: gpmem_t av0 = avma, av;
! 1663: long KCCO, i, j, s, **mat;
! 1664: long nrelsup, nreldep, LIMC, LIMC2, cp, nlze;
! 1665: GEN h, W, cyc, res, gen, dep, C, B, extramat, extraC;
! 1666: GEN R, resc, Res, z;
! 1667: double drc, lim, LOGD, LOGD2;
! 1668: const long MAXRELSUP = 7;
1.1 noro 1669:
1670: Disc = D; if (typ(Disc)!=t_INT) err(typeer,"buchquad");
1.2 ! noro 1671: s = mod4(Disc);
1.1 noro 1672: switch(signe(Disc))
1673: {
1674: case -1:
1675: if (cmpis(Disc,-4) >= 0)
1676: {
1.2 ! noro 1677: GEN p1=cgetg(6,t_VEC); p1[1]=p1[4]=p1[5]=un;
1.1 noro 1678: p1[2]=p1[3]=lgetg(1,t_VEC); return p1;
1679: }
1680: if (s==2 || s==1) err(funder2,"buchquad");
1681: PRECREG=0; break;
1682:
1683: case 1:
1684: if (s==2 || s==3) err(funder2,"buchquad");
1685: if (flag)
1686: err(talker,"sorry, narrow class group not implemented. Use bnfnarrow");
1687: PRECREG=1; break;
1688:
1689: default: err(talker,"zero discriminant in quadclassunit");
1690: }
1691: if (carreparfait(Disc)) err(talker,"square argument in quadclassunit");
1692: if (!isfundamental(Disc))
1693: err(warner,"not a fundamental discriminant in quadclassunit");
1694: buch_init(); RELSUP = RELSUP0;
1.2 ! noro 1695: drc = fabs(gtodouble(Disc));
! 1696: LOGD = log(drc);
! 1697: LOGD2 = LOGD * LOGD;
! 1698:
! 1699: lim = sqrt(drc);
! 1700: /* resc = sqrt(D) w / 2^r1 (2pi)^r2 ~ hR / L(chi,1) */
! 1701: if (PRECREG) resc = dbltor(lim / 2.);
! 1702: else resc = dbltor(lim / PI);
1.1 noro 1703: if (!PRECREG) lim /= sqrt(3.);
1.2 ! noro 1704: R = gun;
1.1 noro 1705: cp = (long)exp(sqrt(LOGD*log(LOGD)/8.0));
1706: if (cp < 13) cp = 13;
1707: av = avma; cbach /= 2;
1.2 ! noro 1708: mat = NULL;
1.1 noro 1709:
1.2 ! noro 1710: START: avma = av; cbach = check_bach(cbach,6.);
! 1711: if (mat) { desalloc(mat); mat = NULL; }
1.1 noro 1712: nreldep = nrelsup = 0;
1.2 ! noro 1713: LIMC = (long)(cbach*LOGD2);
! 1714: if (LIMC < cp) { LIMC = cp; cbach = LIMC / LOGD2; }
! 1715: LIMC2 = (long)(max(cbach,cbach2)*LOGD2);
1.1 noro 1716: if (LIMC2 < LIMC) LIMC2 = LIMC;
1717: if (PRECREG)
1718: {
1719: PRECREG = max(prec+1, MEDDEFAULTPREC + 2*(expi(Disc)>>TWOPOTBITS_IN_LONG));
1720: sqrtD = gsqrt(Disc,PRECREG);
1721: isqrtD = gfloor(sqrtD);
1722: }
1723:
1.2 ! noro 1724: Res = FBquad(Disc,LIMC2,LIMC);
! 1725: if (!Res) goto START;
! 1726: subFB = subFBquad(Disc, lim + 0.5, KC);
! 1727: if (!subFB) goto START;
! 1728: powsubFB = powsubFBquad(CBUCH+1);
! 1729: limhash = ((ulong)LIMC < (MAXHALFULONG>>1))? LIMC*LIMC: HIGHBIT>>1;
! 1730: for (i=0; i<HASHT; i++) hashtab[i]=NULL;
! 1731:
1.1 noro 1732: KCCO = KC + RELSUP;
1.2 ! noro 1733: if (DEBUGLEVEL) fprintferr("KC = %ld, KCCO = %ld\n",KC,KCCO);
! 1734: mat = (long**)cgetalloc(NULL, KCCO+1, t_VEC);
1.1 noro 1735: for (i=1; i<=KCCO; i++)
1736: {
1.2 ! noro 1737: GEN t = cgetalloc(NULL, KC+1, t_VECSMALL);
! 1738: for (j=1; j<=KC; j++) t[j]=0;
! 1739: mat[i] = t;
1.1 noro 1740: }
1741:
1.2 ! noro 1742: s = lg(subFB)-1 + RELSUP;
! 1743: C = PRECREG? real_relations(KCCO,s,LIMC,mat)
! 1744: : imag_relations(KCCO,s,LIMC,mat);
! 1745: W = hnfspec(mat,vperm,&dep,&B,&C,lg(subFB)-1);
! 1746: nlze = lg(dep)>1? lg(dep[1])-1: lg(B[1])-1;
1.1 noro 1747:
1748: if (nlze)
1749: {
1.2 ! noro 1750: MORE:
! 1751: extramat = extra_relations(LIMC,nlze, &extraC);
! 1752: if (!extramat) { goto START; }
! 1753: W = hnfadd(W,vperm,&dep,&B,&C, extramat,extraC);
! 1754: nlze = lg(dep)>1? lg(dep[1])-1: lg(B[1])-1;
! 1755: KCCO += lg(extramat)-1;
1.1 noro 1756: if (nlze)
1757: {
1.2 ! noro 1758: if (++nreldep > 5) goto START;
! 1759: goto MORE;
1.1 noro 1760: }
1761: }
1.2 ! noro 1762:
! 1763: h = dethnf_i(W);
! 1764: if (DEBUGLEVEL) fprintferr("\n#### Tentative class number: %Z\n", h);
! 1765:
! 1766: z = mulrr(Res, resc); /* ~ hR if enough relations, a multiple otherwise */
! 1767: switch(get_R(C, KCCO - (lg(B)-1) - (lg(W)-1), divir(h,z), &R))
1.1 noro 1768: {
1.2 ! noro 1769: case PRECI:
! 1770: prec = (PRECREG<<1)-2;
! 1771: goto START;
! 1772:
! 1773: case RELAT:
! 1774: if (++nrelsup <= MAXRELSUP) { nlze = min(KC, nrelsup); goto MORE; }
! 1775: goto START;
1.1 noro 1776: }
1.2 ! noro 1777: /* DONE */
! 1778: if (!quad_be_honest()) goto START;
1.1 noro 1779:
1.2 ! noro 1780: gen = get_clgp(Disc,W,&cyc,PRECREG);
! 1781: desalloc(mat);
1.1 noro 1782:
1783: res=cgetg(6,t_VEC);
1.2 ! noro 1784: res[1]=(long)h;
! 1785: res[2]=(long)cyc;
! 1786: res[3]=(long)gen;
! 1787: res[4]=(long)R;
! 1788: res[5]=ldiv(mpmul(R,h), z); return gerepilecopy(av0,res);
1.1 noro 1789: }
1790:
1791: GEN
1792: buchimag(GEN D, GEN c, GEN c2, GEN REL)
1793: {
1794: return buchquad(D,gtodouble(c),gtodouble(c2),itos(REL), 0,0);
1795: }
1796:
1797: GEN
1798: buchreal(GEN D, GEN sens0, GEN c, GEN c2, GEN REL, long prec)
1799: {
1800: return buchquad(D,gtodouble(c),gtodouble(c2),itos(REL), itos(sens0),prec);
1801: }
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-2.2 / src / basemath