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