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