[BACK]Return to buch1.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/buch1.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:03 version 1.2, 2002/09/11 07:26:50
Line 23  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
Line 23  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
   
 const int narrow = 0; /* should set narrow = flag in buchquad, but buggy */  const int narrow = 0; /* should set narrow = flag in buchquad, but buggy */
   
   /* For largeprime() hashtable. Note that hashed pseudoprimes are odd (unless
    * 2 | index), hence the low order bit is not useful. So we hash
    * HASHBITS bits starting at bit 1, not bit 0 */
   #define HASHBITS 10
   const int HASHT = 1 << HASHBITS;
   
   static long
   hash(long q) { return (q & ((1 << (HASHBITS+1)) - 1)) >> 1; }
   #undef HASHBITS
   
 /* See buch2.c:  /* See buch2.c:
  * precision en digits decimaux=2*(#digits decimaux de Disc)+50   * precision en digits decimaux=2*(#digits decimaux de Disc)+50
  * on prendra les p decomposes tels que prod(p)>lim dans la subbase   * on prendra les p decomposes tels que prod(p)>lim dans la subFB
  * LIMC=Max(c.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))   * LIMC=Max(c.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))
  * LIMC2=Max(6.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))   * LIMC2=Max(6.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc))))
  * subbase contient les p decomposes tels que prod(p)>sqrt(Disc)   * subFB contient les indices des p decomposes tels que prod(p) > sqrt(Disc)
  * lgsub=subbase[0]=#subbase;   * powsubFB est la table des puissances des formes dans subFB
  * subfactorbase est la table des form[p] pour p dans subbase  
  * nbram est le nombre de p divisant Disc elimines dans subbase  
  * powsubfactorbase est la table des puissances des formes dans subfactorbase  
  */   */
 #define HASHT 1024 /* power of 2 */  
 static const long CBUCH = 15; /* of the form 2^k-1 */  
 static const long randshift=BITS_IN_RANDOM-1 - 4; /* BITS_IN_RANDOM-1-k */  
   
 static long KC,KC2,lgsub,limhash,RELSUP,PRECREG;  #define RANDOM_BITS 4
 static long *primfact,*primfact1, *exprimfact,*exprimfact1, *badprim;  static const long CBUCH = (1<<RANDOM_BITS)-1;
 static long *factorbase,*numfactorbase, *subbase, *vectbase, **hashtab;  static const long randshift=BITS_IN_RANDOM-1 - RANDOM_BITS;
 static GEN  **powsubfactorbase,vperm,subfactorbase,Disc,sqrtD,isqrtD;  #undef RANDOM_BITS
   
   static long KC,KC2,limhash,RELSUP,PRECREG;
   static long *primfact,*exprimfact,*badprim;
   static long *FB,*numFB, **hashtab;
   static GEN  powsubFB,vperm,subFB,Disc,sqrtD,isqrtD;
   
 GEN buchquad(GEN D, double c, double c2, long RELSUP0, long flag, long prec);  GEN buchquad(GEN D, double c, double c2, long RELSUP0, long flag, long prec);
 extern GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus);  extern GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus);
 extern GEN roots_to_pol(GEN L, long v);  extern GEN colreducemodHNF(GEN x, GEN y, GEN *Q);
 extern GEN colreducemodmat(GEN x, GEN y, GEN *Q);  extern GEN quadhilbertreal(GEN D, GEN flag, long prec);
   extern void comp_gen(GEN z,GEN x,GEN y);
   extern GEN codeform5(GEN x, long prec);
   extern GEN comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD);
   extern GEN redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD);
   extern GEN rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD);
   extern GEN cgetalloc(GEN x, size_t l, long t);
   
 GEN  GEN
 quadclassunit0(GEN x, long flag, GEN data, long prec)  quadclassunit0(GEN x, long flag, GEN data, long prec)
Line 109  getallforms(GEN D, long *pth, GEN *ptz)
Line 124  getallforms(GEN D, long *pth, GEN *ptz)
   
 #define MOD4(x) ((((GEN)x)[2])&3)  #define MOD4(x) ((((GEN)x)[2])&3)
 #define MAXL 300  #define MAXL 300
 /* find P and Q two non principal prime ideals (above p,q) such that  /* find P and Q two non principal prime ideals (above p,q) such that
  *   (pq, z) = 1  [coprime to all class group representatives]   *   (pq, z) = 1  [coprime to all class group representatives]
  *   cl(P) = cl(Q) if P has order 2 in Cl(K)   *   cl(P) = cl(Q) if P has order 2 in Cl(K)
  * Try to have e = 24 / gcd(24, (p-1)(q-1)) as small as possible */   * Try to have e = 24 / gcd(24, (p-1)(q-1)) as small as possible */
Line 143  get_pq(GEN D, GEN z, GEN flag, GEN *ptp, GEN *ptq)
Line 158  get_pq(GEN D, GEN z, GEN flag, GEN *ptp, GEN *ptq)
   ell = 3;    ell = 3;
   while (l < 3 || ell<=MAXL)    while (l < 3 || ell<=MAXL)
   {    {
     ell += *diffell++; if (!*diffell) err(primer1);      NEXT_PRIME_VIADIFF_CHECK(ell, diffell);
     if (smodis(z,ell) && kross(d,ell) > 0)      if (smodis(z,ell) && kross(d,ell) > 0)
     {      {
       court[2]=ell; form = redimag(primeform(D,court,0));        court[2]=ell; form = redimag(primeform(D,court,0));
Line 210  gpq(GEN form, GEN p, GEN q, long e, GEN sqd, GEN u, lo
Line 225  gpq(GEN form, GEN p, GEN q, long e, GEN sqd, GEN u, lo
 static GEN  static GEN
 quadhilbertimag(GEN D, GEN flag)  quadhilbertimag(GEN D, GEN flag)
 {  {
   long av=avma,h,i,e,prec;    long h, i, e, prec;
     gpmem_t av=avma;
   GEN z,L,P,p,q,qfp,qfq,up,uq,u;    GEN z,L,P,p,q,qfp,qfq,up,uq,u;
   int raw = ((typ(flag)==t_INT && signe(flag)));    int raw = ((typ(flag)==t_INT && signe(flag)));
   
   if (DEBUGLEVEL>=2) timer2();    if (DEBUGLEVEL>=2) (void)timer2();
   if (gcmpgs(D,-11) >= 0) return polx[0];    if (gcmpgs(D,-11) >= 0) return polx[0];
   L = getallforms(D,&h,&z);    L = getallforms(D,&h,&z);
   if (DEBUGLEVEL>=2) msgtimer("class number = %ld",h);    if (DEBUGLEVEL>=2) msgtimer("class number = %ld",h);
Line 237  quadhilbertimag(GEN D, GEN flag)
Line 253  quadhilbertimag(GEN D, GEN flag)
   prec = raw? DEFAULTPREC: 3;    prec = raw? DEFAULTPREC: 3;
   for(;;)    for(;;)
   {    {
     long av0 = avma, ex, exmax = 0;      long ex, exmax = 0;
       gpmem_t av0 = avma;
     GEN lead, sqd = gsqrt(negi(D),prec);      GEN lead, sqd = gsqrt(negi(D),prec);
     P = cgetg(h+1,t_VEC);      P = cgetg(h+1,t_VEC);
     for (i=1; i<=h; i++)      for (i=1; i<=h; i++)
     {      {
       GEN v, s = gpq((GEN)L[i], p, q, e, sqd, u, prec);        GEN v, s = gpq((GEN)L[i], p, q, e, sqd, u, prec);
Line 250  quadhilbertimag(GEN D, GEN flag)
Line 267  quadhilbertimag(GEN D, GEN flag)
         v[2] = (long)s;          v[2] = (long)s;
       }        }
       else P[i] = (long)s;        else P[i] = (long)s;
       ex = gexpo(s); if (ex > 0) exmax += ex;        ex = gexpo(s); if (ex > 0) exmax += ex;
     }      }
     if (DEBUGLEVEL>=2) msgtimer("roots");      if (DEBUGLEVEL>=2) msgtimer("roots");
     if (raw) { P = gcopy(P); break; }      if (raw) { P = gcopy(P); break; }
Line 271  quadhilbertimag(GEN D, GEN flag)
Line 288  quadhilbertimag(GEN D, GEN flag)
   return gerepileupto(av,P);    return gerepileupto(av,P);
 }  }
   
 GEN quadhilbertreal(GEN D, GEN flag, long prec);  
   
 GEN  GEN
 quadhilbert(GEN D, GEN flag, long prec)  quadhilbert(GEN D, GEN flag, long prec)
 {  {
Line 310  get_om(GEN nf, GEN a)
Line 325  get_om(GEN nf, GEN a)
  * Set list[j + 1] = g1^e1...gk^ek where j is the integer   * Set list[j + 1] = g1^e1...gk^ek where j is the integer
  *   ek + ck [ e(k-1) + c(k-1) [... + c2 [e1]]...] */   *   ek + ck [ e(k-1) + c(k-1) [... + c2 [e1]]...] */
 static GEN  static GEN
 getallelts(GEN nf, GEN G)  getallelts(GEN bnr)
 {  {
   GEN C,c,g, *list, **pows, *gk;    GEN nf,G,C,c,g, *list, **pows, *gk;
   long lc,i,j,k,no;    long lc,i,j,k,no;
   
     nf = checknf(bnr);
     G = (GEN)bnr[5];
   
   no = itos((GEN)G[1]);    no = itos((GEN)G[1]);
   c = (GEN)G[2];    c = (GEN)G[2];
   g = (GEN)G[3]; lc = lg(c)-1;    g = (GEN)G[3]; lc = lg(c)-1;
Line 330  getallelts(GEN nf, GEN G)
Line 348  getallelts(GEN nf, GEN G)
   {    {
     c[i] = k = itos((GEN)c[i]);      c[i] = k = itos((GEN)c[i]);
     gk = (GEN*)cgetg(k, t_VEC); gk[1] = (GEN)g[i];      gk = (GEN*)cgetg(k, t_VEC); gk[1] = (GEN)g[i];
     for (j=2; j<k; j++) gk[j] = idealmul(nf, gk[j-1], gk[1]);      for (j=2; j<k; j++)
         gk[j] = idealmodidele(bnr, idealmul(nf, gk[j-1], gk[1]));
     pows[i] = gk; /* powers of g[i] */      pows[i] = gk; /* powers of g[i] */
   }    }
   
   C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc];    C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc];
   for (i=2; i<=lc; i++) C[i] = C[i-1] * c[lc-i+1];    for (i=2; i<=lc; i++) C[i] = C[i-1] * c[lc-i+1];
   /* C[i] = c(k-i+1) * ... * ck */    /* C[i] = c(k-i+1) * ... * ck */
   /* j < C[i+1] <==> j only involves g(k-i)...gk */    /* j < C[i+1] <==> j only involves g(k-i)...gk */
Line 345  getallelts(GEN nf, GEN G)
Line 364  getallelts(GEN nf, GEN G)
   {    {
     GEN p1,p2;      GEN p1,p2;
     if (j == C[i+1]) i++;      if (j == C[i+1]) i++;
     p2 = pows[lc-i][j/C[i]];      p2 = pows[lc-i][j/C[i]];
     p1 = list[j%C[i] + 1]; if (p1) p2 = idealmul(nf,p2,p1);      p1 = list[j%C[i] + 1];
       if (p1) p2 = idealmodidele(bnr, idealmul(nf,p2,p1));
     list[j + 1] = p2;      list[j + 1] = p2;
   }    }
   list[1] = idealhermite(nf,gun);    list[1] = idealhermite(nf,gun);
Line 419  get_prec(GEN P, long prec)
Line 439  get_prec(GEN P, long prec)
 }  }
   
 /* AUXILLIARY ROUTINES FOR QUADRAYSIGMA */  /* AUXILLIARY ROUTINES FOR QUADRAYSIGMA */
 GEN  
 PiI2(long prec)  
 {  
   GEN z = cgetg(3,t_COMPLEX);  
   GEN x = mppi(prec); setexpo(x,2);  
   z[1] = zero;  
   z[2] = (long)x; return z; /* = 2I Pi */  
 }  
   
 /* Compute data for ellphist */  /* Compute data for ellphist */
 static GEN  static GEN
Line 436  ellphistinit(GEN om, long prec)
Line 448  ellphistinit(GEN om, long prec)
   
   if (gsigne(gimag(gdiv(om1,om2))) < 0)    if (gsigne(gimag(gdiv(om1,om2))) < 0)
   {    {
     p1=om1; om1=om2; om2=p1;      p1 = om1; om1 = om2; om2 = p1;
     p1=cgetg(3,t_VEC); p1[1]=(long)om1; p1[2]=(long)om2;      om = cgetg(3,t_VEC);
     om = p1;      om[1] = (long)om1;
       om[2] = (long)om2;
   }    }
   om1b = gconj(om1);    om1b = gconj(om1);
   om2b = gconj(om2); res = cgetg(4,t_VEC);    om2b = gconj(om2); res = cgetg(4,t_VEC);
Line 475  computeth2(GEN om, GEN la, long prec)
Line 488  computeth2(GEN om, GEN la, long prec)
 static GEN  static GEN
 computeP2(GEN bnr, GEN la, int raw, long prec)  computeP2(GEN bnr, GEN la, int raw, long prec)
 {  {
   long av=avma, av2, clrayno,i, first = 1;    long clrayno, i, first = 1;
     gpmem_t av=avma, av2;
   GEN listray,P0,P,f,lanum, nf = checknf(bnr);    GEN listray,P0,P,f,lanum, nf = checknf(bnr);
   
   f = gmael3(bnr,2,1,1);    f = gmael3(bnr,2,1,1);
   if (typ(la) != t_COL) la = algtobasis(nf,la);    if (typ(la) != t_COL) la = algtobasis(nf,la);
   listray = getallelts(nf,(GEN)bnr[5]);    listray = getallelts(bnr);
   clrayno = lg(listray)-1; av2 = avma;    clrayno = lg(listray)-1; av2 = avma;
 PRECPB:  PRECPB:
   if (!first)    if (!first)
Line 518  do_compo(GEN x, GEN y)
Line 532  do_compo(GEN x, GEN y)
 {  {
   long a, ph = degpol(y);    long a, ph = degpol(y);
   GEN z;    GEN z;
   y = gmul(gpuigs(polx[0],ph), gsubst(y,0,gdiv(polx[MAXVARN],polx[0])));    y = gmul(gpowgs(polx[0],ph), gsubst(y,0,gdiv(polx[MAXVARN],polx[0])));
   for  (a = 0;; a = nexta(a))    for  (a = 0;; a = nexta(a))
   {    {
     if (a) x = gsubst(x, 0, gaddsg(a, polx[0]));      if (a) x = gsubst(x, 0, gaddsg(a, polx[0]));
Line 543  galoisapplypol(GEN nf, GEN s, GEN x)
Line 557  galoisapplypol(GEN nf, GEN s, GEN x)
 static GEN  static GEN
 findquad(GEN a, GEN x, GEN p)  findquad(GEN a, GEN x, GEN p)
 {  {
   long tu,tv, av = avma;    long tu, tv;
     gpmem_t av = avma;
   GEN u,v;    GEN u,v;
   if (typ(x) == t_POLMOD) x = (GEN)x[2];    if (typ(x) == t_POLMOD) x = (GEN)x[2];
   if (typ(a) == t_POLMOD) a = (GEN)a[2];    if (typ(a) == t_POLMOD) a = (GEN)a[2];
Line 645  treatspecialsigma(GEN nf, GEN gf, int raw, long prec)
Line 660  treatspecialsigma(GEN nf, GEN gf, int raw, long prec)
     return NULL;      return NULL;
   }    }
   
   p1 = gcoeff(gf,1,1);    p1 = gcoeff(gf,1,1); /* integer > 0 */
   p2 = gcoeff(gf,2,2);    p2 = gcoeff(gf,2,2);
   if (gcmp1(p2)) { fl = 0; tryf = p1; }    if (gcmp1(p2)) { fl = 0; tryf = p1; }
   else {    else {
     if (Ds % 16 != 8 || !egalii(content(gf),gdeux)) return NULL;      if (Ds % 16 != 8 || !egalii(content(gf),gdeux)) return NULL;
     fl = 1; tryf = shifti(p1,-1);      fl = 1; tryf = shifti(p1,-1);
   }    }
   if (cmpis(tryf, 3) <= 0 || !gcmp0(resii(D, tryf)) || !isprime(tryf))    /* tryf integer > 0 */
     if (cmpis(tryf, 3) <= 0 || signe(resii(D, tryf)) || !isprime(tryf))
     return NULL;      return NULL;
   
   i = itos(tryf); if (fl) i *= 4;    i = itos(tryf); if (fl) i *= 4;
Line 701  quadrayimagsigma(GEN bnr, int raw, long prec)
Line 717  quadrayimagsigma(GEN bnr, int raw, long prec)
   f2 = 2 * itos(gcoeff(f,1,1));    f2 = 2 * itos(gcoeff(f,1,1));
   u = getallrootsof1(bnf); lu = lg(u);    u = getallrootsof1(bnf); lu = lg(u);
   for (i=1; i<lu; i++)    for (i=1; i<lu; i++)
     u[i] = (long)colreducemodmat((GEN)u[i], f, NULL); /* roots of 1, mod f */      u[i] = (long)colreducemodHNF((GEN)u[i], f, NULL); /* roots of 1, mod f */
   if (DEBUGLEVEL>1)    if (DEBUGLEVEL>1)
     fprintferr("quadray: looking for [a,b] != unit mod 2f\n[a,b] = ");      fprintferr("quadray: looking for [a,b] != unit mod 2f\n[a,b] = ");
   for (a=0; a<f2; a++)    for (a=0; a<f2; a++)
Line 712  quadrayimagsigma(GEN bnr, int raw, long prec)
Line 728  quadrayimagsigma(GEN bnr, int raw, long prec)
       if (DEBUGLEVEL>1) fprintferr("[%ld,%ld] ",a,b);        if (DEBUGLEVEL>1) fprintferr("[%ld,%ld] ",a,b);
   
       labas = algtobasis(nf, la);        labas = algtobasis(nf, la);
       lamodf = colreducemodmat(labas, f, NULL);        lamodf = colreducemodHNF(labas, f, NULL);
       for (i=1; i<lu; i++)        for (i=1; i<lu; i++)
         if (gegal(lamodf, (GEN)u[i])) break;          if (gegal(lamodf, (GEN)u[i])) break;
       if (i < lu) continue; /* la = unit mod f */        if (i < lu) continue; /* la = unit mod f */
Line 731  GEN
Line 747  GEN
 quadray(GEN D, GEN f, GEN flag, long prec)  quadray(GEN D, GEN f, GEN flag, long prec)
 {  {
   GEN bnr,y,p1,pol,bnf,lambda;    GEN bnr,y,p1,pol,bnf,lambda;
   long av = avma, raw;    long raw;
     gpmem_t av = avma;
   
   if (!flag) flag = gzero;    if (!flag) flag = gzero;
   if (typ(flag)==t_INT) lambda = NULL;    if (typ(flag)==t_INT) lambda = NULL;
Line 781  quadray(GEN D, GEN f, GEN flag, long prec)
Line 798  quadray(GEN D, GEN f, GEN flag, long prec)
 /*  Routines related to binary quadratic forms (for internal use)  */  /*  Routines related to binary quadratic forms (for internal use)  */
 /*                                                                 */  /*                                                                 */
 /*******************************************************************/  /*******************************************************************/
 extern void comp_gen(GEN z,GEN x,GEN y);  
 extern GEN codeform5(GEN x, long prec);  
 extern GEN comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD);  
 extern GEN redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD);  
 extern GEN rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD);  
   
 #define rhorealform(x) rhoreal_aux(x,Disc,sqrtD,isqrtD)  #define rhorealform(x) rhoreal_aux(x,Disc,sqrtD,isqrtD)
 #define redrealform(x) gcopy(fix_signs(redrealform5(x,Disc,sqrtD,isqrtD)))  #define redrealform(x) gcopy(fix_signs(redrealform5(x,Disc,sqrtD,isqrtD)))
 #define comprealform(x,y) fix_signs(comprealform5(x,y,Disc,sqrtD,isqrtD))  
   
   static GEN
   fix_signs(GEN x)
   {
     GEN a = (GEN)x[1];
     GEN c = (GEN)x[3];
     if (signe(a) < 0)
     {
       if (narrow || absi_equal(a,c)) return rhorealform(x);
       setsigne(a,1); setsigne(c,-1);
     }
     return x;
   }
   
   static GEN
   comprealform(GEN x,GEN y) {
     return fix_signs( comprealform5(x,y,Disc,sqrtD,isqrtD) );
   }
   
 /* compute rho^n(x) */  /* compute rho^n(x) */
 static GEN  static GEN
 rhoreal_pow(GEN x, long n)  rhoreal_pow(GEN x, long n)
 {  {
   long i, av = avma, lim = stack_lim(av,1);    long i;
     gpmem_t av = avma, lim = stack_lim(av, 1);
   for (i=1; i<=n; i++)    for (i=1; i<=n; i++)
   {    {
     x = rhorealform(x);      x = rhorealform(x);
     if (low_stack(lim, stack_lim(av,1)))      if (low_stack(lim, stack_lim(av,1)))
     {      {
       if(DEBUGMEM>1) err(warnmem,"rhoreal_pow");        if(DEBUGMEM>1) err(warnmem,"rhoreal_pow");
       x = gerepilecopy(av, x);        x = gerepilecopy(av, x);
     }      }
   }    }
Line 809  rhoreal_pow(GEN x, long n)
Line 838  rhoreal_pow(GEN x, long n)
 }  }
   
 static GEN  static GEN
 fix_signs(GEN x)  realpf5(GEN D, long p)
 {  {
   GEN a = (GEN)x[1];    gpmem_t av = avma;
   GEN c = (GEN)x[3];    GEN y = primeform(D,stoi(p),PRECREG);
   if (signe(a) < 0)  
   {  
     if (narrow || absi_equal(a,c)) return rhorealform(x);  
     setsigne(a,1); setsigne(c,-1);  
   }  
   return x;  
 }  
   
 static GEN  
 redrealprimeform5(GEN Disc, long p)  
 {  
   long av = avma;  
   GEN y = primeform(Disc,stoi(p),PRECREG);  
   y = codeform5(y,PRECREG);    y = codeform5(y,PRECREG);
   return gerepileupto(av, redrealform(y));    return gerepileupto(av, redrealform(y));
 }  }
   
 static GEN  static GEN
 redrealprimeform(GEN Disc, long p)  realpf(GEN D, long p)
 {  {
   long av = avma;    gpmem_t av = avma;
   GEN y = primeform(Disc,stoi(p),PRECREG);    GEN y = primeform(D,stoi(p),PRECREG);
   return gerepileupto(av, redrealform(y));    return gerepileupto(av, redrealform(y));
 }  }
   
 static GEN  static GEN
   imagpf(GEN D, long p) { return primeform(D,stoi(p),0); }
   
   static GEN
 comprealform3(GEN x, GEN y)  comprealform3(GEN x, GEN y)
 {  {
   long av = avma;    gpmem_t av = avma;
   GEN z = cgetg(4,t_VEC); comp_gen(z,x,y);    GEN z = cgetg(4,t_VEC); comp_gen(z,x,y);
   return gerepileupto(av, redrealform(z));    return gerepileupto(av, redrealform(z));
 }  }
   
   /* Warning: ex[0] not set */
 static GEN  static GEN
 initrealform5(long *ex)  init_form(long *ex, GEN (*comp)(GEN,GEN))
 {  {
   GEN form = powsubfactorbase[1][ex[1]];    long i, l = lg(powsubFB);
   long i;    GEN F = NULL;
     for (i=1; i<l; i++)
       if (ex[i])
       {
         GEN t = gmael(powsubFB,i,ex[i]);
         F = F? comp(F,t): t;
       }
     return F;
   }
   static GEN
   initrealform5(long *ex) { return init_form(ex, &comprealform); }
   static GEN
   initimagform(long *ex) { return init_form(ex, &compimag); }
   
   for (i=2; i<=lgsub; i++)  static GEN
     form = comprealform(form, powsubfactorbase[i][ex[i]]);  random_form(GEN ex, GEN (*comp)(GEN,GEN))
   return form;  {
     long i, l = lg(ex);
     gpmem_t av = avma;
     GEN F;
     for(;;)
     {
       for (i=1; i<l; i++) ex[i] = mymyrand()>>randshift;
       if ((F = init_form(ex, comp))) return F;
       avma = av;
     }
 }  }
   
   static GEN
   real_random_form(GEN ex) { return random_form(ex, &comprealform3); }
   static GEN
   imag_random_form(GEN ex) { return random_form(ex, &compimag); }
   
 /*******************************************************************/  /*******************************************************************/
 /*                                                                 */  /*                                                                 */
 /*                     Common subroutines                          */  /*                     Common subroutines                          */
Line 865  initrealform5(long *ex)
Line 911  initrealform5(long *ex)
 static void  static void
 buch_init(void)  buch_init(void)
 {  {
   if (DEBUGLEVEL) timer2();    if (DEBUGLEVEL) (void)timer2();
   primfact  = new_chunk(100);    primfact  = new_chunk(100);
   primfact1 = new_chunk(100);  
   exprimfact  = new_chunk(100);    exprimfact  = new_chunk(100);
   exprimfact1 = new_chunk(100);  
   badprim = new_chunk(100);    badprim = new_chunk(100);
   hashtab = (long**) new_chunk(HASHT);    hashtab = (long**) new_chunk(HASHT);
 }  }
Line 885  check_bach(double cbach, double B)
Line 929  check_bach(double cbach, double B)
 }  }
   
 static long  static long
 factorisequad(GEN f, long kcz, long limp)  factorquad(GEN f, long kcz, long limp)
 {  {
   long i,p,k,av,lo;    long i, p, k, lo;
     gpmem_t av;
   GEN q,r, x = (GEN)f[1];    GEN q,r, x = (GEN)f[1];
   
   if (is_pm1(x)) { primfact[0]=0; return 1; }    if (is_pm1(x)) { primfact[0]=0; return 1; }
Line 895  factorisequad(GEN f, long kcz, long limp)
Line 940  factorisequad(GEN f, long kcz, long limp)
   if (signe(x) < 0) x = absi(x);    if (signe(x) < 0) x = absi(x);
   for (i=1; ; i++)    for (i=1; ; i++)
   {    {
     p=factorbase[i]; q=dvmdis(x,p,&r);      p=FB[i]; q=dvmdis(x,p,&r);
     if (!signe(r))      if (!signe(r))
     {      {
       k=0; while (!signe(r)) { x=q; k++; q=dvmdis(x,p,&r); }        for (k=0; !signe(r); k++) { x=q; q=dvmdis(x,p,&r); }
       lo++; primfact[lo]=p; exprimfact[lo]=k;        primfact[++lo]=p; exprimfact[lo]=k;
     }      }
     if (cmpis(q,p)<=0) break;      if (cmpis(q,p)<=0) break;
     if (i==kcz) { avma=av; return 0; }      if (i==kcz) { avma=av; return 0; }
   }    }
   p = x[2]; avma=av;    avma=av;
   /* p = itos(x) if lgefint(x)=3 */    if (lgefint(x)!=3 || (p=x[2]) > limhash) return 0;
   if (lgefint(x)!=3 || p > limhash) return 0;  
   
   if (p != 1 && p <= limp)    if (p != 1 && p <= limp)
   {    {
     for (i=1; i<=badprim[0]; i++)      for (i=1; i<=badprim[0]; i++)
       if (p % badprim[i] == 0) return 0;        if (p % badprim[i] == 0) return 0;
     lo++; primfact[lo]=p; exprimfact[lo]=1;      primfact[++lo]=p; exprimfact[lo]=1;
     p = 1;      p = 1;
   }    }
   primfact[0]=lo; return p;    primfact[0]=lo; return p;
Line 922  factorisequad(GEN f, long kcz, long limp)
Line 966  factorisequad(GEN f, long kcz, long limp)
 static long *  static long *
 largeprime(long q, long *ex, long np, long nrho)  largeprime(long q, long *ex, long np, long nrho)
 {  {
   const long hashv = ((q & (2 * HASHT - 1)) - 1) >> 1;    const long hashv = hash(q);
   long *pt, i;    long *pt, i, l = lg(subFB);
   
  /* If q = 0 (2048), very slim chance of getting a relation.  
   * And hashtab[-1] is undefined anyway */  
   if (hashv < 0) return NULL;  
   for (pt = hashtab[hashv]; ; pt = (long*) pt[0])    for (pt = hashtab[hashv]; ; pt = (long*) pt[0])
   {    {
     if (!pt)      if (!pt)
     {      {
       pt = (long*) gpmalloc((lgsub+4)<<TWOPOTBYTES_IN_LONG);        pt = (long*) gpmalloc((l+3)<<TWOPOTBYTES_IN_LONG);
       *pt++ = nrho; /* nrho = pt[-3] */        *pt++ = nrho; /* nrho = pt[-3] */
       *pt++ = np;   /* np   = pt[-2] */        *pt++ = np;   /* np   = pt[-2] */
       *pt++ = q;    /* q    = pt[-1] */        *pt++ = q;    /* q    = pt[-1] */
       pt[0] = (long)hashtab[hashv];        pt[0] = (long)hashtab[hashv];
       for (i=1; i<=lgsub; i++) pt[i]=ex[i];        for (i=1; i<l; i++) pt[i]=ex[i];
       hashtab[hashv]=pt; return NULL;        hashtab[hashv]=pt; return NULL;
     }      }
     if (pt[-1] == q) break;      if (pt[-1] == q) break;
   }    }
   for(i=1; i<=lgsub; i++)    for(i=1; i<l; i++)
     if (pt[i] != ex[i]) return pt;      if (pt[i] != ex[i]) return pt;
   return (pt[-2]==np)? (GEN)NULL: pt;    return (pt[-2]==np)? (GEN)NULL: pt;
 }  }
   
 static long  /* p | conductor of order of disc D ? */
 badmod8(GEN x)  static int
   is_bad(GEN D, ulong p)
 {  {
   long r = mod8(x);    gpmem_t av = avma;
   if (!r) return 1;    int r;
   if (signe(Disc) < 0) r = 8-r;    if (p == 2)
   return (r < 4);    {
       r = mod16(D) >> 1;
       if (r && signe(D) < 0) r = 8-r;
       return (r < 4);
     }
     r = (resii(D, muluu(p,p)) == gzero); /* p^2 | D ? */
     avma = av; return r;
 }  }
   
 /* cree factorbase, numfactorbase, vectbase; affecte badprim */  /* create FB, numFB; fill badprim. Return L(kro_D, 1) */
 static void  static GEN
 factorbasequad(GEN Disc, long n2, long n)  FBquad(GEN Disc, long n2, long n)
 {  {
   long i,p,bad, av = avma;    GEN Res = realun(DEFAULTPREC);
   byteptr d=diffptr;    long i, p, bad, s;
     gpmem_t av;
     byteptr d = diffptr;
   
   numfactorbase = (long*) gpmalloc(sizeof(long)*(n2+1));    numFB = cgetg(n2+1, t_VECSMALL);
   factorbase    = (long*) gpmalloc(sizeof(long)*(n2+1));    FB    = cgetg(n2+1, t_VECSMALL);
   KC=0; bad=0; i=0; p = *d++;    av = avma;
   while (p<=n2)    KC = 0; bad = 0; i = 0;
     if (maxprime() < n2) err(primer1);
     for (p = 0;;) /* p <= n2 */
   {    {
     switch(krogs(Disc,p))      NEXT_PRIME_VIADIFF(p, d);
       if (!KC && p > n) KC = i;
       if (p > n2) break;
       s = krogs(Disc,p);
       switch (s)
     {      {
       case -1: break; /* inert */        case -1: break; /* inert */
       case  0: /* ramified */        case  0: /* ramified */
       {          if (is_bad(Disc, (ulong)p)) { badprim[++bad]=p; break; }
         GEN p1 = divis(Disc,p);          /* fall through */
         if (smodis(p1,p) == 0)  
           if (p!=2 || badmod8(p1)) { badprim[++bad]=p; break; }  
         i++; numfactorbase[p]=i; factorbase[i] = -p; break;  
       }  
       default:  /* split */        default:  /* split */
         i++; numfactorbase[p]=i; factorbase[i] = p;          i++; numFB[p]=i; FB[i] = p; break;
     }      }
     p += *d++; if (!*d) err(primer1);      Res = mulsr(p, divrs(Res, p - s));
     if (KC == 0 && p>n) KC = i;  
   }    }
   if (!KC) { free(factorbase); free(numfactorbase); return; }    if (!KC) return NULL;
     Res = gerepileupto(av, Res);
   KC2 = i;    KC2 = i;
   vectbase = (long*) gpmalloc(sizeof(long)*(KC2+1));    setlg(FB, KC2+1);
   for (i=1; i<=KC2; i++)  
   {  
     p = factorbase[i];  
     vectbase[i]=p; factorbase[i]=labs(p);  
   }  
   if (DEBUGLEVEL)    if (DEBUGLEVEL)
   {    {
     msgtimer("factor base");      msgtimer("factor base");
     if (DEBUGLEVEL>7)      if (DEBUGLEVEL>7)
     {      {
       fprintferr("factorbase:\n");        fprintferr("FB:\n");
       for (i=1; i<=KC; i++) fprintferr("%ld ",factorbase[i]);        for (i=1; i<=KC; i++) fprintferr("%ld ", FB[i]);
       fprintferr("\n"); flusherr();        fprintferr("\n");
     }      }
   }    }
   avma=av; badprim[0] = bad;    badprim[0] = bad; return Res;
 }  }
   
 /* cree vectbase and subfactorbase. Affecte lgsub */  /* create vperm, return subFB */
 static long  static GEN
 subfactorbasequad(double ll, long KC)  subFBquad(GEN D, double PROD, long KC)
 {  {
   long i,j,k,nbidp,p,pro[100], ss;    long i, j, lgsub, ino, lv = KC+1;
   GEN y;    double prod = 1.;
   double prod;    gpmem_t av;
     GEN no;
   
   i=0; ss=0; prod=1;    vperm = cgetg(lv, t_VECSMALL);
   for (j=1; j<=KC && prod<=ll; j++)    av = avma;
     no    = cgetg(lv, t_VECSMALL); ino = 1;
     for (i=j=1; j < lv; j++)
   {    {
     p = vectbase[j];      long p = FB[j];
     if (p>0) { pro[++i]=p; prod*=p; vperm[i]=j; } else ss++;      if (smodis(D, p) == 0) no[ino++] = j; /* ramified */
       else
       {
         vperm[i] = j; i++;
         prod *= p;
         if (prod > PROD) break;
       }
   }    }
   if (prod<=ll) return -1;    if (j == lv) return NULL;
   nbidp=i;    lgsub = i;
   for (k=1; k<j; k++)    for (j = 1; j <ino; i++,j++) vperm[i] = no[j];
     if (vectbase[k]<=0) vperm[++i]=k;    for (     ; i < lv; i++)     vperm[i] = i;
     avma = av;
   y=cgetg(nbidp+1,t_COL);    if (DEBUGLEVEL) msgtimer("subFBquad (%ld elt.)", lgsub-1);
   if (PRECREG) /* real */    return vecextract_i(vperm, 1, lgsub-1);
     for (j=1; j<=nbidp; j++)  
       y[j] = (long)redrealprimeform5(Disc, pro[j]);  
   else  
     for (j=1; j<=nbidp; j++) /* imaginary */  
       y[j] = (long)primeform(Disc,stoi(pro[j]),0);  
   subbase = (long*) gpmalloc(sizeof(long)*(nbidp+1));  
   lgsub = nbidp; for (j=1; j<=lgsub; j++) subbase[j]=pro[j];  
   if (DEBUGLEVEL>7)  
   {  
     fprintferr("subfactorbase: ");  
     for (i=1; i<=lgsub; i++)  
       { fprintferr("%ld: ",subbase[i]); outerr((GEN)y[i]); }  
     fprintferr("\n"); flusherr();  
   }  
   subfactorbase = y; return ss;  
 }  }
   
 static void  /* assume n >= 1, x[i][j] = subFB[i]^j, for j = 1..n */
 powsubfact(long n, long a)  static GEN
   powsubFBquad(long n)
 {  {
   GEN unform, *y, **x = (GEN**) gpmalloc(sizeof(GEN*)*(n+1));    long i,j, l = lg(subFB);
   long i,j;    GEN F, y, x = cgetg(l, t_VEC);
   
   for (i=1; i<=n; i++)  
     x[i] = (GEN*) gpmalloc(sizeof(GEN)*(a+1));  
   if (PRECREG) /* real */    if (PRECREG) /* real */
   {    {
     unform=cgetg(6,t_VEC);      for (i=1; i<l; i++)
     unform[1]=un;  
     unform[2]=(mod2(Disc) == mod2(isqrtD))? (long)isqrtD: laddsi(-1,isqrtD);  
     unform[3]=lshifti(subii(sqri((GEN)unform[2]),Disc),-2);  
     unform[4]=zero;  
     unform[5]=(long)realun(PRECREG);  
     for (i=1; i<=n; i++)  
     {      {
       y = x[i]; y[0] = unform;        F = realpf5(Disc, FB[subFB[i]]);
       for (j=1; j<=a; j++)        y = cgetg(n+1, t_VEC); x[i] = (long)y;
         y[j] = comprealform(y[j-1],(GEN)subfactorbase[i]);        y[1] = (long)F;
         for (j=2; j<=n; j++) y[j] = (long)comprealform((GEN)y[j-1], F);
     }      }
   }    }
   else /* imaginary */    else /* imaginary */
   {    {
     unform = primeform(Disc, gun, 0);      for (i=1; i<l; i++)
     for (i=1; i<=n; i++)  
     {      {
       y = x[i]; y[0] = unform;        F = imagpf(Disc, FB[subFB[i]]);
       for (j=1; j<=a; j++)        y = cgetg(n+1, t_VEC); x[i] = (long)y;
         y[j] = compimag(y[j-1],(GEN)subfactorbase[i]);        y[1] = (long)F;
         for (j=2; j<=n; j++) y[j] = (long)compimag((GEN)y[j-1], F);
     }      }
   }    }
   if (DEBUGLEVEL) msgtimer("powsubfact");    if (DEBUGLEVEL) msgtimer("powsubFBquad");
   powsubfactorbase = x;    return x;
 }  }
   
 static void  static void
 desalloc(long **mat)  desalloc(long **mat)
 {  {
   long i,*p,*q;    long i,*p,*q;
     for (i=1; i<lg(mat); i++) free(mat[i]);
     free(mat);
     for (i=1; i<HASHT; i++)
       for (p = hashtab[i]; p; p = q) { q=(long*)p[0]; free(p-3); }
   }
   
   free(vectbase); free(factorbase); free(numfactorbase);  static void
   if (mat)  sub_fact(GEN col, GEN F)
   {
     GEN b = (GEN)F[2];
     long i, p, e;
     for (i=1; i<=primfact[0]; i++)
   {    {
     free(subbase);      p = primfact[i]; e = exprimfact[i];
     for (i=1; i<lg(subfactorbase); i++) free(powsubfactorbase[i]);      if (smodis(b, p<<1) > p) e = -e;
     for (i=1; i<lg(mat); i++) free(mat[i]);      col[numFB[p]] -= e;
     free(mat); free(powsubfactorbase);  
     for (i=1; i<HASHT; i++)  
       for (p = hashtab[i]; p; p = q) { q=(long*)p[0]; free(p-3); }  
   }    }
 }  }
   static void
 /* L-function */  add_fact(GEN col, GEN F)
 static GEN  
 lfunc(GEN Disc)  
 {  {
   long av=avma, p;    GEN b = (GEN)F[2];
   GEN y=realun(DEFAULTPREC);    long i, p, e;
   byteptr d=diffptr;    for (i=1; i<=primfact[0]; i++)
   
   for(p = *d++; p<=30000; p += *d++)  
   {    {
     if (!*d) err(primer1);      p = primfact[i]; e = exprimfact[i];
     y = mulsr(p, divrs(y, p-krogs(Disc,p)));      if (smodis(b, p<<1) > p) e = -e;
       col[numFB[p]] += e;
   }    }
   return gerepileupto(av,y);  
 }  }
   
 #define comp(x,y) x? (PRECREG? compreal(x,y): compimag(x,y)): y  #define comp(x,y) x? (PRECREG? compreal(x,y): compimag(x,y)): y
 static GEN  static GEN
 get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec)  get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec)
 {  {
   GEN p1,p2,res,*init, u1u2=smith2(W), u1=(GEN)u1u2[1], met=(GEN)u1u2[3];    GEN res,*init, u1, met = smithrel(W,NULL,&u1);
   long c,i,j, l = lg(met);    long c,i,j, l = lg(W);
   
   u1 = reducemodmatrix(ginv(u1), W);    c = lg(met);
   for (c=1; c<l; c++)  
     if (gcmp1(gcoeff(met,c,c))) break;  
   if (DEBUGLEVEL) msgtimer("smith/class group");    if (DEBUGLEVEL) msgtimer("smith/class group");
   res=cgetg(c,t_VEC); init = (GEN*)cgetg(l,t_VEC);    res=cgetg(c,t_VEC); init = (GEN*)cgetg(l,t_VEC);
   for (i=1; i<l; i++)    for (i=1; i<l; i++)
     init[i] = primeform(Disc,stoi(labs(vectbase[vperm[i]])),prec);      init[i] = primeform(Disc,stoi(FB[vperm[i]]),prec);
   for (j=1; j<c; j++)    for (j=1; j<c; j++)
   {    {
     p1 = NULL;      GEN p1 = NULL;
     for (i=1; i<l; i++)      for (i=1; i<l; i++)
     {        p1 = comp(p1, powgi(init[i], gcoeff(u1,i,j)));
       p2 = gpui(init[i], gcoeff(u1,i,j), prec);  
       p1 = comp(p1,p2);  
     }  
     res[j] = (long)p1;      res[j] = (long)p1;
   }    }
   if (DEBUGLEVEL) msgtimer("generators");    if (DEBUGLEVEL) msgtimer("generators");
Line 1142  get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec)
Line 1174  get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec)
 }  }
   
 static GEN  static GEN
 extra_relations(long LIMC, long *ex, long nlze, GEN extramatc)  extra_relations(long LIMC, long nlze, GEN *ptextraC)
 {  {
   long av,fpc,p,ep,i,j,k,nlze2, *col, *colg, s = 0, extrarel = nlze+2;    long fpc, i, k, nlze2;
     long s = 0, extrarel = nlze+2, lgsub = lg(subFB);
     gpmem_t av;
   long MAXRELSUP = min(50,4*KC);    long MAXRELSUP = min(50,4*KC);
   GEN p1,form, extramat = cgetg(extrarel+1,t_MAT);    GEN p1, form, ex, extramat, extraC, col;
   
   if (DEBUGLEVEL)    extramat = cgetg(extrarel+1, t_MAT);
   {    extraC   = cgetg(extrarel+1, t_VEC);
     fprintferr("looking for %ld extra relations\n",extrarel);    for (i=1; i<=extrarel; i++) extramat[i] = lgetg(KC+1,t_VECSMALL);
     flusherr();    if (!PRECREG)
   }      for (i=1; i<=extrarel; i++) extraC[i] = lgetg(1,t_COL);
   for (j=1; j<=extrarel; j++) extramat[j]=lgetg(KC+1,t_COL);  
   nlze2 = PRECREG? max(nlze,lgsub): min(nlze+1,KC);    if (DEBUGLEVEL) fprintferr("looking for %ld extra relations\n",extrarel);
     nlze2 = PRECREG? max(nlze,lgsub-1): min(nlze+1,KC);
   if (nlze2 < 3 && KC > 2) nlze2 = 3;    if (nlze2 < 3 && KC > 2) nlze2 = 3;
     ex = cgetg(nlze2+1, t_VECSMALL);
   av = avma;    av = avma;
   while (s<extrarel)    while (s<extrarel)
   {    {
     form = NULL;      form = NULL;
     for (i=1; i<=nlze2; i++)      for (i=1; i<=nlze2; i++)
     {      {
       ex[i]=mymyrand()>>randshift;        ex[i] = mymyrand()>>randshift;
       if (ex[i])        if (ex[i])
       {        {
         p1 = primeform(Disc,stoi(factorbase[vperm[i]]),PRECREG);          p1 = primeform(Disc,stoi(FB[vperm[i]]),PRECREG);
         p1 = gpuigs(p1,ex[i]); form = comp(form,p1);          p1 = gpowgs(p1,ex[i]); form = comp(form,p1);
       }        }
     }      }
     if (!form) continue;      if (!form) continue;
   
     fpc = factorisequad(form,KC,LIMC);      fpc = factorquad(form,KC,LIMC);
     if (fpc==1)      if (fpc==1)
     {      {
       s++; col = (GEN)extramat[s];        s++; col = (GEN)extramat[s];
       for (i=1; i<=nlze2; i++) col[vperm[i]] = -ex[i];        for (i=1; i<=nlze2; i++) col[vperm[i]] = -ex[i];
       for (   ; i<=KC; i++) col[vperm[i]]= 0;        for (   ; i<=KC;    i++) col[vperm[i]] = 0;
       for (j=1; j<=primfact[0]; j++)        add_fact(col, form);
       {  
         p=primfact[j]; ep=exprimfact[j];  
         if (smodis((GEN)form[2], p<<1) > p) ep = -ep;  
         col[numfactorbase[p]] += ep;  
       }  
       for (i=1; i<=KC; i++)        for (i=1; i<=KC; i++)
         if (col[i]) break;          if (col[i]) break;
       if (i>KC)        if (i>KC)
Line 1190  extra_relations(long LIMC, long *ex, long nlze, GEN ex
Line 1221  extra_relations(long LIMC, long *ex, long nlze, GEN ex
         s--; avma = av;          s--; avma = av;
         if (--MAXRELSUP == 0) return NULL;          if (--MAXRELSUP == 0) return NULL;
       }        }
       else { av = avma; if (PRECREG) coeff(extramatc,1,s) = form[4]; }        else { av = avma; if (PRECREG) extraC[s] = form[4]; }
     }      }
     else avma = av;      else avma = av;
     if (DEBUGLEVEL)      if (DEBUGLEVEL)
     {      {
       if (fpc == 1) fprintferr(" %ld",s);        if (fpc == 1) fprintferr(" %ld",s);
       else if (DEBUGLEVEL>1) fprintferr(".");        else if (DEBUGLEVEL>1) fprintferr(".");
       flusherr();  
     }      }
   }    }
   for (j=1; j<=extrarel; j++)  
     for (i=1; i<=extrarel; i++)
   {    {
     colg = cgetg(KC+1,t_COL);      GEN colg = cgetg(KC+1,t_COL);
     col = (GEN)extramat[j]; extramat[j] = (long) colg;      col = (GEN)extramat[i];
     for (k=1; k<=KC; k++)      extramat[i] = (long)colg;
       colg[k] = lstoi(col[vperm[k]]);      for (k=1; k<=KC; k++) colg[k] = lstoi(col[vperm[k]]);
   }    }
   if (DEBUGLEVEL)    if (DEBUGLEVEL) { fprintferr("\n"); msgtimer("extra relations"); }
     *ptextraC = extraC; return extramat;
   }
   #undef comp
   
   static long
   trivial_relations(long **mat, GEN C)
   {
     long i, s = 0;
     GEN col;
     for (i=1; i<=KC; i++)
   {    {
     fprintferr("\n");      if (smodis(Disc, FB[i])) continue;
     msgtimer("extra relations");      /* ramified prime ==> trivial relation */
       col = mat[++s];
       col[i] = 2;
       if (C) affsr(0, (GEN)C[s]);
   }    }
   return extramat;    return s;
 }  }
 #undef comp  
   
   static void
   dbg_all(long phase, long s, long n)
   {
     if (DEBUGLEVEL>1) fprintferr("\n");
     msgtimer("%s rel [#rel/#test = %ld/%ld]", phase? "random": "initial", s, n);
   }
   
   void
   wr_rel(GEN col)
   {
     long i, l = lg(col);
     fprintferr("\nrel = ");
     for (i=1; i<l; i++)
       if (col[i]) fprintferr("%ld^%ld ",i,col[i]);
     fprintferr("\n");
   }
   
   void
   dbg_rel(long s, GEN col)
   {
     if (DEBUGLEVEL == 1) fprintferr("%4ld",s);
     else { fprintferr("cglob = %ld. ", s); wr_rel(col); }
     flusherr();
   }
 /*******************************************************************/  /*******************************************************************/
 /*                                                                 */  /*                                                                 */
 /*                    Imaginary Quadratic fields                   */  /*                    Imaginary Quadratic fields                   */
 /*                                                                 */  /*                                                                 */
 /*******************************************************************/  /*******************************************************************/
   
   /* Strategy 1 until lim relation found, then Strategy 2 until LIM. */
 static GEN  static GEN
 imag_random_form(long current, long *ex)  imag_relations(long LIM, long lim, long LIMC, long **mat)
 {  {
   long av = avma,i;    long lgsub = lg(subFB), i, s, fpc, current, nbtest = 0;
   GEN form,pc;    int first = 1;
     gpmem_t av;
     GEN C, col, form, ex = cgetg(lgsub, t_VECSMALL);
     GEN dummy = cgetg(1,t_COL);
   
     C = cgetg(LIM+1, t_VEC);
     for (i=1; i<=LIM; i++) C[i] = (long)dummy;
     av = avma;
     s = trivial_relations(mat, NULL);
     lim += s; if (lim > LIM) lim = LIM;
   for(;;)    for(;;)
   {    {
     form = pc = primeform(Disc,stoi(factorbase[current]),PRECREG);      if (s >= lim) {
     for (i=1; i<=lgsub; i++)        if (s >= LIM) break;
     {        lim = LIM; first = 0;
       ex[i] = mymyrand()>>randshift;        if (DEBUGLEVEL) dbg_all(0, s, nbtest);
       if (ex[i])  
         form = compimag(form,powsubfactorbase[i][ex[i]]);  
     }      }
     if (form != pc) return form;      avma = av; current = first? 1+(s%KC): 1+s-RELSUP;
     avma = av; /* ex = 0, try again */      form = imagpf(Disc, FB[current]);
   }      form = compimag(form, imag_random_form(ex));
 }      nbtest++; fpc = factorquad(form,KC,LIMC);
   
 static void  
 imag_relations(long lim, long s, long LIMC, long *ex, long **mat)  
 {  
   static long nbtest;  
   long av = avma, i,j,pp,fpc,b1,b2,ep,current, first = (s==0);  
   long *col,*fpd,*oldfact,*oldexp;  
   GEN pc,form,form1;  
   
   if (first) nbtest = 0 ;  
   while (s<lim)  
   {  
     avma=av; nbtest++; current = first? 1+(s%KC): 1+s-RELSUP;  
     form = imag_random_form(current,ex);  
     fpc = factorisequad(form,KC,LIMC);  
     if (!fpc)      if (!fpc)
     {      {
       if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }        if (DEBUGLEVEL>1) fprintferr(".");
       continue;        continue;
     }      }
     if (fpc > 1)      if (fpc > 1)
     {      {
       fpd = largeprime(fpc,ex,current,0);        long *fpd = largeprime(fpc,ex,current,0);
         long b1, b2, p;
         GEN form2;
       if (!fpd)        if (!fpd)
       {        {
         if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }          if (DEBUGLEVEL>1) fprintferr(".");
         continue;          continue;
       }        }
       form1 = powsubfactorbase[1][fpd[1]];        form2 = initimagform(fpd);
       for (i=2; i<=lgsub; i++)        form2 = compimag(form2, imagpf(Disc, FB[fpd[-2]]));
         form1 = compimag(form1,powsubfactorbase[i][fpd[i]]);        p = fpc << 1;
       pc=primeform(Disc,stoi(factorbase[fpd[-2]]),0);        b1 = smodis((GEN)form2[2], p);
       form1=compimag(form1,pc);        b2 = smodis((GEN)form[2],  p);
       pp = fpc << 1;        if (b1 != b2 && b1+b2 != p) continue;
       b1=smodis((GEN)form1[2], pp);  
       b2=smodis((GEN)form[2],  pp);  
       if (b1 != b2 && b1+b2 != pp) continue;  
   
       s++; col = mat[s];        col = mat[++s];
       if (DEBUGLEVEL) { fprintferr(" %ld",s); flusherr(); }        add_fact(col, form);
       oldfact = primfact; oldexp = exprimfact;        (void)factorquad(form2,KC,LIMC);
       primfact = primfact1; exprimfact = exprimfact1;  
       factorisequad(form1,KC,LIMC);  
   
       if (b1==b2)        if (b1==b2)
       {        {
         for (i=1; i<=lgsub; i++)          for (i=1; i<lgsub; i++) col[subFB[i]] += fpd[i]-ex[i];
           col[numfactorbase[subbase[i]]] = fpd[i]-ex[i];          sub_fact(col, form2); col[fpd[-2]]++;
         col[fpd[-2]]++;  
         for (j=1; j<=primfact[0]; j++)  
         {  
           pp=primfact[j]; ep=exprimfact[j];  
           if (smodis((GEN)form1[2], pp<<1) > pp) ep = -ep;  
           col[numfactorbase[pp]] -= ep;  
         }  
       }        }
       else        else
       {        {
         for (i=1; i<=lgsub; i++)          for (i=1; i<lgsub; i++) col[subFB[i]] += -fpd[i]-ex[i];
           col[numfactorbase[subbase[i]]] = -fpd[i]-ex[i];          add_fact(col, form2); col[fpd[-2]]--;
         col[fpd[-2]]--;  
         for (j=1; j<=primfact[0]; j++)  
         {  
           pp=primfact[j]; ep=exprimfact[j];  
           if (smodis((GEN)form1[2], pp<<1) > pp) ep = -ep;  
           col[numfactorbase[pp]] += ep;  
         }  
       }        }
       primfact = oldfact; exprimfact = oldexp;      }
     }  
     else      else
     {      {
       s++; col = mat[s];        col = mat[++s];
       if (DEBUGLEVEL) { fprintferr(" %ld",s); flusherr(); }        for (i=1; i<lgsub; i++) col[subFB[i]] = -ex[i];
       for (i=1; i<=lgsub; i++)        add_fact(col, form);
         col[numfactorbase[subbase[i]]] = -ex[i];  
     }      }
     for (j=1; j<=primfact[0]; j++)  
     {  
       pp=primfact[j]; ep=exprimfact[j];  
       if (smodis((GEN)form[2], pp<<1) > pp) ep = -ep;  
       col[numfactorbase[pp]] += ep;  
     }  
     col[current]--;      col[current]--;
       if (DEBUGLEVEL) dbg_rel(s, col);
     if (!first && fpc == 1 && col[current] == 0)      if (!first && fpc == 1 && col[current] == 0)
     {      {
       s--; for (i=1; i<=KC; i++) col[i]=0;        s--; for (i=1; i<=KC; i++) col[i]=0;
     }      }
   }    }
   if (DEBUGLEVEL)    if (DEBUGLEVEL) dbg_all(1, s, nbtest);
   {    return C;
     fprintferr("\nnbrelations/nbtest = %ld/%ld\n",s,nbtest);  
     msgtimer("%s relations", first? "initial": "random");  
   }  
 }  }
   
 static int  static int
 imag_be_honest(long *ex)  imag_be_honest()
 {  {
   long p,fpc, s = KC, nbtest = 0, av = avma;    long p, fpc, s = KC, nbtest = 0;
   GEN form;    gpmem_t av = avma;
     GEN F, ex = cgetg(lg(subFB), t_VECSMALL);
   
   while (s<KC2)    while (s<KC2)
   {    {
     p = factorbase[s+1];      p = FB[s+1]; if (DEBUGLEVEL) fprintferr(" %ld",p);
     if (DEBUGLEVEL) { fprintferr(" %ld",p); flusherr(); }      F = compimag(imagpf(Disc, p), imag_random_form(ex));
     form = imag_random_form(s+1,ex);      fpc = factorquad(F,s,p-1);
     fpc = factorisequad(form,s,p-1);  
     if (fpc == 1) { nbtest=0; s++; }      if (fpc == 1) { nbtest=0; s++; }
     else      else
       if (++nbtest > 20) return 0;        if (++nbtest > 20) return 0;
Line 1363  imag_be_honest(long *ex)
Line 1394  imag_be_honest(long *ex)
 /*                                                                 */  /*                                                                 */
 /*******************************************************************/  /*******************************************************************/
   
 static GEN  /* (1/2) log (d * 2^{e * EXP220}) */
 real_random_form(long *ex)  GEN
   get_dist(GEN e, GEN d, long prec)
 {  {
   long av = avma, i;    GEN t = mplog(absr(d));
   GEN p1,form = NULL;    if (signe(e)) t = addrr(t, mulir(mulsi(EXP220,e), mplog2(prec)));
     return shiftr(t, -1);
   for(;;)  
   {  
     for (i=1; i<=lgsub; i++)  
     {  
       ex[i] = mymyrand()>>randshift;  
 /*    if (ex[i]) KB: less efficient if I put this in. Why ??? */  
       {  
         p1 = powsubfactorbase[i][ex[i]];  
         form = form? comprealform3(form,p1): p1;  
       }  
     }  
     if (form) return form;  
     avma = av;  
   }  
 }  }
   
 static void  static GEN
 real_relations(long lim, long s, long LIMC, long *ex, long **mat, GEN glog2,  real_relations(long LIM, long lim, long LIMC, long **mat)
                GEN vecexpo)  
 {  {
   static long nbtest;    long lgsub = lg(subFB), i, s, fpc, current, nbtest = 0, endcycle, rhoacc, rho;
   long av = avma,av1,i,j,p,fpc,b1,b2,ep,current, first = (s==0);    int first = 1;
   long *col,*fpd,*oldfact,*oldexp,limstack;    gpmem_t av, av1, limstack;
   long findecycle,nbrhocumule,nbrho;    GEN C, d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL);
   GEN p1,p2,form,form0,form1,form2;  
   
   limstack=stack_lim(av,1);    C = cgetg(LIM+1, t_VEC);
   if (first) nbtest = 0;    for (i=1; i<=LIM; i++) C[i] = lgetr(PRECREG);
   
   current = 0;    current = 0;
   p1 = NULL; /* gcc -Wall */    av = avma; limstack = stack_lim(av,1);
   while (s<lim)    s = trivial_relations(mat, C);
     lim += s; if (lim > LIM) lim = LIM;
   NEW:
     for(;;)
   {    {
       if (s >= lim) {
         if (lim == LIM) break;
         lim = LIM; first = 0;
         if (DEBUGLEVEL) dbg_all(0, s, nbtest);
       }
       avma = av;
     form = real_random_form(ex);      form = real_random_form(ex);
     if (!first)      if (!first)
     {      {
       current = 1+s-RELSUP;        current = 1+s-RELSUP;
       p1 = redrealprimeform(Disc, factorbase[current]);        form = comprealform3(form, realpf(Disc, FB[current]));
       form = comprealform3(form,p1);  
     }      }
       av1 = avma;
     form0 = form; form1 = NULL;      form0 = form; form1 = NULL;
     findecycle = nbrhocumule = 0;      endcycle = rhoacc = 0;
     nbrho = -1; av1 = avma;      rho = -1;
     while (s<lim)  
   CYCLE:
       if (endcycle) goto NEW;
       if (low_stack(limstack, stack_lim(av,1)))
     {      {
       if (low_stack(limstack, stack_lim(av,1)))        if(DEBUGMEM>1) err(warnmem,"real_relations");
       {        gerepileall(av1, form1? 2: 1, &form, &form1);
         GEN *gptr[2];      }
         int c = 1;      if (rho < 0) rho = 0; /* first time in */
         if(DEBUGMEM>1) err(warnmem,"real_relations");      else
         gptr[0]=&form; if (form1) gptr[c++]=&form1;      {
         gerepilemany(av1,gptr,c);        form = rhorealform(form); rho++;
       }        rhoacc++;
       if (nbrho < 0) nbrho = 0; /* first time in */        if (first)
           endcycle = (absi_equal((GEN)form[1],(GEN)form0[1])
                && egalii((GEN)form[2],(GEN)form0[2])
                && (!narrow || signe(form0[1])==signe(form[1])));
       else        else
       {        {
         if (findecycle) break;          if (narrow)
         form = rhorealform(form);            { form = rhorealform(form); rho++; }
         nbrho++; nbrhocumule++;          else if (absi_equal((GEN)form[1], (GEN)form[3])) /* a = -c */
         if (first)  
         {          {
           if (absi_equal((GEN)form[1],(GEN)form0[1])            if (absi_equal((GEN)form[1],(GEN)form0[1]) &&
                && egalii((GEN)form[2],(GEN)form0[2])                    egalii((GEN)form[2],(GEN)form0[2])) goto NEW;
                && (!narrow || signe(form0[1])==signe(form[1]))) findecycle=1;            form = rhorealform(form); rho++;
         }          }
         else          else
         {            { setsigne(form[1],1); setsigne(form[3],-1); }
           if (narrow)          if (egalii((GEN)form[1],(GEN)form0[1]) &&
             { form=rhorealform(form); nbrho++; }              egalii((GEN)form[2],(GEN)form0[2])) goto NEW;
           else if (absi_equal((GEN)form[1], (GEN)form[3])) /* a = -c */  
           {  
             if (absi_equal((GEN)form[1],(GEN)form0[1]) &&  
                     egalii((GEN)form[2],(GEN)form0[2])) break;  
             form=rhorealform(form); nbrho++;  
           }  
           else  
             { setsigne(form[1],1); setsigne(form[3],-1); }  
           if (egalii((GEN)form[1],(GEN)form0[1]) &&  
               egalii((GEN)form[2],(GEN)form0[2])) break;  
         }  
       }        }
       nbtest++; fpc = factorisequad(form,KC,LIMC);      }
       if (!fpc)      nbtest++; fpc = factorquad(form,KC,LIMC);
       if (!fpc)
       {
         if (DEBUGLEVEL>1) fprintferr(".");
         goto CYCLE;
       }
       if (fpc > 1)
       { /* look for Large Prime relation */
         long *fpd = largeprime(fpc,ex,current,rhoacc);
         long b1, b2, p;
         GEN form2;
         if (!fpd)
       {        {
         if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }          if (DEBUGLEVEL>1) fprintferr(".");
         continue;          goto CYCLE;
       }        }
       if (fpc > 1)        if (!form1) form1 = initrealform5(ex);
         if (!first) form1 = comprealform(form1, realpf5(Disc, FB[current]));
         form1 = rhoreal_pow(form1, rho);
         rho = 0;
   
         form2 = initrealform5(fpd);
         if (fpd[-2]) form2 = comprealform(form2, realpf5(Disc, FB[fpd[-2]]));
         form2 = rhoreal_pow(form2, fpd[-3]);
         if (!narrow && !absi_equal((GEN)form2[1],(GEN)form2[3]))
       {        {
         fpd = largeprime(fpc,ex,current,nbrhocumule);          setsigne(form2[1],1);
         if (!fpd)          setsigne(form2[3],-1);
         {        }
           if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }        p = fpc << 1;
           continue;        b1 = smodis((GEN)form2[2], p);
         }        b2 = smodis((GEN)form1[2], p);
         if (!form1) form1 = initrealform5(ex);        if (b1 != b2 && b1+b2 != p) goto CYCLE;
         if (!first)  
         {  
           p1 = redrealprimeform5(Disc, factorbase[current]);  
           form1=comprealform(form1,p1);  
         }  
         form1 = rhoreal_pow(form1, nbrho); nbrho = 0;  
         form2 = initrealform5(fpd);  
         if (fpd[-2])  
         {  
           p1 = redrealprimeform5(Disc, factorbase[fpd[-2]]);  
           form2=comprealform(form2,p1);  
         }  
         form2 = rhoreal_pow(form2, fpd[-3]);  
         if (!narrow && !absi_equal((GEN)form2[1],(GEN)form2[3]))  
         {  
           setsigne(form2[1],1);  
           setsigne(form2[3],-1);  
         }  
         p = fpc << 1;  
         b1=smodis((GEN)form2[2], p);  
         b2=smodis((GEN)form1[2], p);  
         if (b1 != b2 && b1+b2 != p) continue;  
   
         s++; col = mat[s]; if (DEBUGLEVEL) fprintferr(" %ld",s);        col = mat[++s];
         oldfact = primfact; oldexp = exprimfact;        add_fact(col, form1);
         primfact = primfact1; exprimfact = exprimfact1;        (void)factorquad(form2,KC,LIMC);
         factorisequad(form2,KC,LIMC);        if (b1==b2)
         if (b1==b2)        {
         {          for (i=1; i<lgsub; i++) col[subFB[i]] += fpd[i]-ex[i];
           for (i=1; i<=lgsub; i++)          sub_fact(col, form2);
             col[numfactorbase[subbase[i]]] = fpd[i]-ex[i];          if (fpd[-2]) col[fpd[-2]]++; /* implies !first */
           for (j=1; j<=primfact[0]; j++)          d = get_dist(subii((GEN)form1[4],(GEN)form2[4]),
           {                       divrr((GEN)form1[5],(GEN)form2[5]), PRECREG);
             p=primfact[j]; ep=exprimfact[j];  
             if (smodis((GEN)form2[2], p<<1) > p) ep = -ep;  
             col[numfactorbase[p]] -= ep;  
           }  
           if (fpd[-2]) col[fpd[-2]]++; /* implies !first */  
           p1 = subii((GEN)form1[4],(GEN)form2[4]);  
           p2 = divrr((GEN)form1[5],(GEN)form2[5]);  
         }  
         else  
         {  
           for (i=1; i<=lgsub; i++)  
             col[numfactorbase[subbase[i]]] = -fpd[i]-ex[i];  
           for (j=1; j<=primfact[0]; j++)  
           {  
             p=primfact[j]; ep=exprimfact[j];  
             if (smodis((GEN)form2[2], p<<1) > p) ep = -ep;  
             col[numfactorbase[p]] += ep;  
           }  
           if (fpd[-2]) col[fpd[-2]]--;  
           p1 = addii((GEN)form1[4],(GEN)form2[4]);  
           p2 = mulrr((GEN)form1[5],(GEN)form2[5]);  
         }  
         primfact = oldfact; exprimfact = oldexp;  
       }        }
       else        else
       {        {
         if (!form1) form1 = initrealform5(ex);          for (i=1; i<lgsub; i++) col[subFB[i]] += -fpd[i]-ex[i];
         if (!first)          add_fact(col, form2);
         {          if (fpd[-2]) col[fpd[-2]]--;
           p1 = redrealprimeform5(Disc, factorbase[current]);          d = get_dist(addii((GEN)form1[4],(GEN)form2[4]),
           form1=comprealform(form1,p1);                       mulrr((GEN)form1[5],(GEN)form2[5]), PRECREG);
         }  
         form1 = rhoreal_pow(form1,nbrho); nbrho = 0;  
   
         s++; col = mat[s]; if (DEBUGLEVEL) fprintferr(" %ld",s);  
         for (i=1; i<=lgsub; i++)  
           col[numfactorbase[subbase[i]]] = -ex[i];  
         p1 = (GEN) form1[4];  
         p2 = (GEN) form1[5];  
       }        }
       for (j=1; j<=primfact[0]; j++)      }
       else
       { /* standard relation */
         if (!form1) form1 = initrealform5(ex);
         if (!first) form1 = comprealform(form1, realpf5(Disc, FB[current]));
         form1 = rhoreal_pow(form1,rho);
         rho = 0;
   
         col = mat[++s];
         for (i=1; i<lgsub; i++) col[subFB[i]] = -ex[i];
         add_fact(col, form1);
         d = get_dist((GEN)form1[4], (GEN)form1[5], PRECREG);
       }
       if (DEBUGLEVEL) fprintferr(" %ld",s);
       affrr(d, (GEN)C[s]);
       if (first)
       {
         if (s >= lim) goto NEW;
         goto CYCLE;
       }
       else
       {
         col[current]--;
         if (fpc == 1 && col[current] == 0)
       {        {
         p=primfact[j]; ep=exprimfact[j];          s--; for (i=1; i<=KC; i++) col[i]=0;
         if (smodis((GEN)form1[2], p<<1) > p) ep = -ep;  
         col[numfactorbase[p]] += ep;  
       }        }
       p1 = mpadd(mulir(mulsi(EXP220,p1), glog2), mplog(absr(p2)));  
       affrr(shiftr(p1,-1), (GEN)vecexpo[s]);  
       if (!first)  
       {  
         col[current]--;  
         if (fpc == 1 && col[current] == 0)  
           { s--; for (i=1; i<=KC; i++) col[i]=0; }  
         break;  
       }  
     }      }
     avma = av;  
   }    }
   if (DEBUGLEVEL)    if (DEBUGLEVEL) dbg_all(1, s, nbtest);
   {    return C;
     fprintferr("\nnbrelations/nbtest = %ld/%ld\n",s,nbtest);  
     msgtimer("%s relations", first? "initial": "random");  
   }  
 }  }
   
 static int  static int
 real_be_honest(long *ex)  real_be_honest()
 {  {
   long p,fpc,s = KC,nbtest = 0,av = avma;    long p, fpc, s = KC, nbtest = 0;
   GEN p1,form,form0;    gpmem_t av = avma;
     GEN F,F0, ex = cgetg(lg(subFB), t_VECSMALL);
   
   while (s<KC2)    while (s<KC2)
   {    {
     p = factorbase[s+1];      p = FB[s+1]; if (DEBUGLEVEL) fprintferr(" %ld",p);
     if (DEBUGLEVEL) { fprintferr(" %ld",p); flusherr(); }      F = comprealform3(realpf(Disc, p), real_random_form(ex));
     form = real_random_form(ex);      for (F0 = F;;)
     p1 = redrealprimeform(Disc, p);  
     form0 = form = comprealform3(form,p1);  
     for(;;)  
     {      {
       fpc = factorisequad(form,s,p-1);        fpc = factorquad(F,s,p-1);
       if (fpc == 1) { nbtest=0; s++; break; }        if (fpc == 1) { nbtest=0; s++; break; }
       form = rhorealform(form);  
       if (++nbtest > 20) return 0;        if (++nbtest > 20) return 0;
       if (narrow || absi_equal((GEN)form[1],(GEN)form[3]))        F = fix_signs( rhorealform(F) );
         form = rhorealform(form);        if (egalii((GEN)F[1],(GEN)F0[1])
       else         && egalii((GEN)F[2],(GEN)F0[2])) break;
       {  
         setsigne(form[1],1);  
         setsigne(form[3],-1);  
       }  
       if (egalii((GEN)form[1],(GEN)form0[1])  
        && egalii((GEN)form[2],(GEN)form0[2])) break;  
     }      }
     avma = av;      avma = av;
   }    }
Line 1597  real_be_honest(long *ex)
Line 1582  real_be_honest(long *ex)
 }  }
   
 static GEN  static GEN
 gcdrealnoer(GEN a,GEN b,long *pte)  gcdreal(GEN a,GEN b)
 {  {
   long e;    long e;
   GEN k1,r;    GEN k1,r;
   
     if (!signe(a)) return mpabs(b);
     if (!signe(b)) return mpabs(a);
   
   if (typ(a)==t_INT)    if (typ(a)==t_INT)
   {    {
     if (typ(b)==t_INT) return mppgcd(a,b);      if (typ(b)==t_INT) return mppgcd(a,b);
     k1=cgetr(lg(b)); affir(a,k1); a=k1;      a = itor(a, lg(b));
   }    }
   else if (typ(b)==t_INT)    else if (typ(b)==t_INT)
     { k1=cgetr(lg(a)); affir(b,k1); b=k1; }    {
       b = itor(b, lg(a));
     }
   if (expo(a)<-5) return absr(b);    if (expo(a)<-5) return absr(b);
   if (expo(b)<-5) return absr(a);    if (expo(b)<-5) return absr(a);
   a=absr(a); b=absr(b);    a=absr(a); b=absr(b);
   while (expo(b) >= -5  && signe(b))    while (expo(b) >= -5  && signe(b))
   {    {
     k1=gcvtoi(divrr(a,b),&e);      k1 = gcvtoi(divrr(a,b),&e);
     if (e > 0) return NULL;      if (e > 0) return NULL;
     r=subrr(a,mulir(k1,b)); a=b; b=r;      r=subrr(a,mulir(k1,b)); a=b; b=r;
   }    }
   *pte=expo(b); return absr(a);    return absr(a);
 }  }
   
 static GEN  enum { RELAT, LARGE, PRECI };
 get_reg(GEN matc, long sreg)  
   static int
   get_R(GEN C, long sreg, GEN z, GEN *ptR)
 {  {
   long i,e,maxe;    GEN R = gun;
   GEN reg = mpabs(gcoeff(matc,1,1));    double c;
     long i;
   
   e = maxe = 0;    if (PRECREG)
   for (i=2; i<=sreg; i++)  
   {    {
     reg = gcdrealnoer(gcoeff(matc,1,i),reg,&e);      R = mpabs((GEN)C[1]);
     if (!reg) return NULL;      for (i=2; i<=sreg; i++)
     maxe = maxe? max(maxe,e): e;      {
         R = gcdreal((GEN)C[i], R);
         if (!R) return PRECI;
       }
       if (DEBUGLEVEL)
       {
         if (DEBUGLEVEL>7) fprintferr("R = %Z",R);
         msgtimer("regulator");
       }
       if (gexpo(R) <= -3)
       {
         if (DEBUGLEVEL) fprintferr("regulator is zero.\n");
         return RELAT;
       }
   }    }
     c = gtodouble(gmul(z, R));
     if (c < 0.75 || c > 1.5) return RELAT;
     *ptR = R; return LARGE;
   }
   
   static int
   quad_be_honest()
   {
     int r;
     if (KC2 <= KC) return 1;
   if (DEBUGLEVEL)    if (DEBUGLEVEL)
   {      fprintferr("be honest for primes from %ld to %ld\n", FB[KC+1],FB[KC2]);
     if (DEBUGLEVEL>7) { fprintferr("reg = "); outerr(reg); }    r = PRECREG? real_be_honest(): imag_be_honest();
     msgtimer("regulator");    if (DEBUGLEVEL) { fprintferr("\n"); msgtimer("be honest"); }
   }    return r;
   return reg;  
 }  }
   
 GEN  GEN
 buchquad(GEN D, double cbach, double cbach2, long RELSUP0, long flag, long prec)  buchquad(GEN D, double cbach, double cbach2, long RELSUP0, long flag, long prec)
 {  {
   long av0 = avma,av,tetpil,KCCO,KCCOPRO,i,j,s, *ex,**mat;    gpmem_t av0 = avma, av;
   long extrarel,nrelsup,nreldep,LIMC,LIMC2,cp,nbram,nlze;    long KCCO, i, j, s, **mat;
   GEN p1,h,W,met,res,basecl,dr,c_1,pdep,C,B,extramat,extraC;    long nrelsup, nreldep, LIMC, LIMC2, cp, nlze;
   GEN reg,vecexpo,glog2,cst;    GEN h, W, cyc, res, gen, dep, C, B, extramat, extraC;
   double drc,lim,LOGD;    GEN R, resc, Res, z;
     double drc, lim, LOGD, LOGD2;
     const long MAXRELSUP = 7;
   
   Disc = D; if (typ(Disc)!=t_INT) err(typeer,"buchquad");    Disc = D; if (typ(Disc)!=t_INT) err(typeer,"buchquad");
   s=mod4(Disc);    s = mod4(Disc);
   glog2 = vecexpo = NULL; /* gcc -Wall */  
   switch(signe(Disc))    switch(signe(Disc))
   {    {
     case -1:      case -1:
       if (cmpis(Disc,-4) >= 0)        if (cmpis(Disc,-4) >= 0)
       {        {
         p1=cgetg(6,t_VEC); p1[1]=p1[4]=p1[5]=un;          GEN p1=cgetg(6,t_VEC); p1[1]=p1[4]=p1[5]=un;
         p1[2]=p1[3]=lgetg(1,t_VEC); return p1;          p1[2]=p1[3]=lgetg(1,t_VEC); return p1;
       }        }
       if (s==2 || s==1) err(funder2,"buchquad");        if (s==2 || s==1) err(funder2,"buchquad");
Line 1677  buchquad(GEN D, double cbach, double cbach2, long RELS
Line 1692  buchquad(GEN D, double cbach, double cbach2, long RELS
   if (!isfundamental(Disc))    if (!isfundamental(Disc))
     err(warner,"not a fundamental discriminant in quadclassunit");      err(warner,"not a fundamental discriminant in quadclassunit");
   buch_init(); RELSUP = RELSUP0;    buch_init(); RELSUP = RELSUP0;
   dr=cgetr(3); affir(Disc,dr); drc=fabs(rtodbl(dr)); LOGD=log(drc);    drc = fabs(gtodouble(Disc));
   lim=sqrt(drc); cst = mulrr(lfunc(Disc), dbltor(lim));    LOGD = log(drc);
     LOGD2 = LOGD * LOGD;
   
     lim = sqrt(drc);
     /* resc = sqrt(D) w / 2^r1 (2pi)^r2 ~ hR / L(chi,1) */
     if (PRECREG) resc = dbltor(lim / 2.);
     else         resc = dbltor(lim / PI);
   if (!PRECREG) lim /= sqrt(3.);    if (!PRECREG) lim /= sqrt(3.);
     R = gun;
   cp = (long)exp(sqrt(LOGD*log(LOGD)/8.0));    cp = (long)exp(sqrt(LOGD*log(LOGD)/8.0));
   if (cp < 13) cp = 13;    if (cp < 13) cp = 13;
   av = avma; cbach /= 2;    av = avma; cbach /= 2;
     mat = NULL;
   
 INCREASE: avma = av; cbach = check_bach(cbach,6.);  START: avma = av; cbach = check_bach(cbach,6.);
     if (mat) { desalloc(mat); mat = NULL; }
   nreldep = nrelsup = 0;    nreldep = nrelsup = 0;
   LIMC = (long)(cbach*LOGD*LOGD);    LIMC = (long)(cbach*LOGD2);
   if (LIMC < cp) LIMC=cp;    if (LIMC < cp) { LIMC = cp; cbach = LIMC / LOGD2; }
   LIMC2 = max(20, (long)(max(cbach,cbach2)*LOGD*LOGD));    LIMC2 = (long)(max(cbach,cbach2)*LOGD2);
   if (LIMC2 < LIMC) LIMC2 = LIMC;    if (LIMC2 < LIMC) LIMC2 = LIMC;
   if (PRECREG)    if (PRECREG)
   {    {
     PRECREG = max(prec+1, MEDDEFAULTPREC + 2*(expi(Disc)>>TWOPOTBITS_IN_LONG));      PRECREG = max(prec+1, MEDDEFAULTPREC + 2*(expi(Disc)>>TWOPOTBITS_IN_LONG));
     glog2  = glog(gdeux,PRECREG);  
     sqrtD  = gsqrt(Disc,PRECREG);      sqrtD  = gsqrt(Disc,PRECREG);
     isqrtD = gfloor(sqrtD);      isqrtD = gfloor(sqrtD);
   }    }
   factorbasequad(Disc,LIMC2,LIMC);  
   if (!KC) goto INCREASE;  
   
   vperm = cgetg(KC+1,t_VECSMALL); for (i=1; i<=KC; i++) vperm[i]=i;    Res = FBquad(Disc,LIMC2,LIMC);
   nbram = subfactorbasequad(lim,KC);    if (!Res) goto START;
   if (nbram == -1) { desalloc(NULL); goto INCREASE; }    subFB = subFBquad(Disc, lim + 0.5, KC);
   KCCO = KC + RELSUP;    if (!subFB) goto START;
   if (DEBUGLEVEL) { fprintferr("KC = %ld, KCCO = %ld\n",KC,KCCO); flusherr(); }    powsubFB = powsubFBquad(CBUCH+1);
   powsubfact(lgsub,CBUCH+7);  
   
   mat = (long**) gpmalloc((KCCO+1)*sizeof(long*));  
   setlg(mat, KCCO+1);  
   for (i=1; i<=KCCO; i++)  
   {  
     mat[i] = (long*) gpmalloc((KC+1)*sizeof(long));  
     for (j=1; j<=KC; j++) mat[i][j]=0;  
   }  
   ex = new_chunk(lgsub+1);  
   limhash = ((ulong)LIMC < (MAXHALFULONG>>1))? LIMC*LIMC: HIGHBIT>>1;    limhash = ((ulong)LIMC < (MAXHALFULONG>>1))? LIMC*LIMC: HIGHBIT>>1;
   for (i=0; i<HASHT; i++) hashtab[i]=NULL;    for (i=0; i<HASHT; i++) hashtab[i]=NULL;
   
   s = lgsub+nbram+RELSUP;    KCCO = KC + RELSUP;
   if (PRECREG)    if (DEBUGLEVEL) fprintferr("KC = %ld, KCCO = %ld\n",KC,KCCO);
     mat = (long**)cgetalloc(NULL, KCCO+1, t_VEC);
     for (i=1; i<=KCCO; i++)
   {    {
     vecexpo=cgetg(KCCO+1,t_VEC);      GEN t = cgetalloc(NULL, KC+1, t_VECSMALL);
     for (i=1; i<=KCCO; i++) vecexpo[i]=lgetr(PRECREG);      for (j=1; j<=KC; j++) t[j]=0;
     real_relations(s,0,LIMC,ex,mat,glog2,vecexpo);      mat[i] = t;
     real_relations(KCCO,s,LIMC,ex,mat,glog2,vecexpo);  
   }    }
   else  
   {  
     imag_relations(s,0,LIMC,ex,mat);  
     imag_relations(KCCO,s,LIMC,ex,mat);  
   }  
   if (KC2 > KC)  
   {  
     if (DEBUGLEVEL)  
       fprintferr("be honest for primes from %ld to %ld\n",  
                   factorbase[KC+1],factorbase[KC2]);  
     s = PRECREG? real_be_honest(ex): imag_be_honest(ex);  
     if (DEBUGLEVEL)  
     {  
       fprintferr("\n");  
       msgtimer("be honest");  
     }  
     if (!s) { desalloc(mat); goto INCREASE; }  
   }  
   C=cgetg(KCCO+1,t_MAT);  
   if (PRECREG)  
   {  
     for (i=1; i<=KCCO; i++)  
     {  
       C[i]=lgetg(2,t_COL); coeff(C,1,i)=vecexpo[i];  
     }  
     if (DEBUGLEVEL>7) { fprintferr("C: "); outerr(C); flusherr(); }  
   }  
   else  
     for (i=1; i<=KCCO; i++) C[i]=lgetg(1,t_COL);  
   W = hnfspec(mat,vperm,&pdep,&B,&C,lgsub);  
   nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;  
   
   KCCOPRO=KCCO;    s = lg(subFB)-1 + RELSUP;
     C = PRECREG? real_relations(KCCO,s,LIMC,mat)
                : imag_relations(KCCO,s,LIMC,mat);
     W = hnfspec(mat,vperm,&dep,&B,&C,lg(subFB)-1);
     nlze = lg(dep)>1? lg(dep[1])-1: lg(B[1])-1;
   
   if (nlze)    if (nlze)
   {    {
 EXTRAREL:  MORE:
     s = PRECREG? 2: 1; extrarel=nlze+2;      extramat = extra_relations(LIMC,nlze, &extraC);
     extraC=cgetg(extrarel+1,t_MAT);      if (!extramat) { goto START; }
     for (i=1; i<=extrarel; i++) extraC[i]=lgetg(s,t_COL);      W = hnfadd(W,vperm,&dep,&B,&C, extramat,extraC);
     extramat = extra_relations(LIMC,ex,nlze,extraC);      nlze = lg(dep)>1? lg(dep[1])-1: lg(B[1])-1;
     if (!extramat) { desalloc(mat); goto INCREASE; }      KCCO += lg(extramat)-1;
     W = hnfadd(W,vperm,&pdep,&B,&C, extramat,extraC);  
     nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;  
     KCCOPRO += extrarel;  
     if (nlze)      if (nlze)
     {      {
       if (++nreldep > 5) { desalloc(mat); goto INCREASE; }        if (++nreldep > 5) goto START;
       goto EXTRAREL;        goto MORE;
     }      }
   }    }
   /* tentative class number */  
   h=gun; for (i=1; i<lg(W); i++) h=mulii(h,gcoeff(W,i,i));  
   if (PRECREG)  
   {  
     /* tentative regulator */  
     reg = get_reg(C, KCCOPRO - (lg(B)-1) - (lg(W)-1));  
     if (!reg)  
     {  
       desalloc(mat);  
       prec = (PRECREG<<1)-2; goto INCREASE;  
     }  
     if (gexpo(reg)<=-3)  
     {  
       if (++nrelsup <= 7)  
       {  
         if (DEBUGLEVEL) { fprintferr("regulateur nul\n"); flusherr(); }  
         nlze=min(KC,nrelsup); goto EXTRAREL;  
       }  
       desalloc(mat); goto INCREASE;  
     }  
     c_1 = divrr(gmul2n(gmul(h,reg),1), cst);  
   }  
   else  
   {  
     reg = gun;  
     c_1 = divrr(gmul(h,mppi(DEFAULTPREC)), cst);  
   }  
   
   if (gcmpgs(gmul2n(c_1,2),3)<0) { c_1=stoi(10); nrelsup=7; }    h = dethnf_i(W);
   if (gcmpgs(gmul2n(c_1,1),3)>0)    if (DEBUGLEVEL) fprintferr("\n#### Tentative class number: %Z\n", h);
   
     z = mulrr(Res, resc); /* ~ hR if enough relations, a multiple otherwise */
     switch(get_R(C, KCCO - (lg(B)-1) - (lg(W)-1), divir(h,z), &R))
   {    {
     nrelsup++;      case PRECI:
     if (nrelsup<=7)        prec = (PRECREG<<1)-2;
     {        goto START;
       if (DEBUGLEVEL)  
         { fprintferr("***** check = %f\n\n",gtodouble(c_1)); flusherr(); }      case RELAT:
       nlze=min(KC,nrelsup); goto EXTRAREL;        if (++nrelsup <= MAXRELSUP) { nlze = min(KC, nrelsup); goto MORE; }
     }        goto START;
     if (cbach < 5.99) { desalloc(mat); goto INCREASE; }  
     err(warner,"suspicious check. Suggest increasing extra relations.");  
   }    }
   basecl = get_clgp(Disc,W,&met,PRECREG);    /* DONE */
   s = lg(basecl); desalloc(mat); tetpil=avma;    if (!quad_be_honest()) goto START;
   
     gen = get_clgp(Disc,W,&cyc,PRECREG);
     desalloc(mat);
   
   res=cgetg(6,t_VEC);    res=cgetg(6,t_VEC);
   res[1]=lcopy(h); p1=cgetg(s,t_VEC);    res[1]=(long)h;
   for (i=1; i<s; i++) p1[i] = (long)icopy(gcoeff(met,i,i));    res[2]=(long)cyc;
   res[2]=(long)p1;    res[3]=(long)gen;
   res[3]=lcopy(basecl);    res[4]=(long)R;
   res[4]=lcopy(reg);    res[5]=ldiv(mpmul(R,h), z); return gerepilecopy(av0,res);
   res[5]=lcopy(c_1); return gerepile(av0,tetpil,res);  
 }  }
   
 GEN  GEN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>