Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - buch1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.1 lcov report (development 28904-c3aa21e911) Lines: 652 684 95.3 %
Date: 2023-12-04 07:51:13 Functions: 49 53 92.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : #define DEBUGLEVEL DEBUGLEVEL_quadclassunit
      18             : 
      19             : /*******************************************************************/
      20             : /*                                                                 */
      21             : /*         CLASS GROUP AND REGULATOR (McCURLEY, BUCHMANN)          */
      22             : /*                   QUADRATIC FIELDS                              */
      23             : /*                                                                 */
      24             : /*******************************************************************/
      25             : /* For largeprime() hashtable. Note that hashed pseudoprimes are odd (unless
      26             :  * 2 | index), hence the low order bit is not useful. So we hash
      27             :  * HASHBITS bits starting at bit 1, not bit 0 */
      28             : #define HASHBITS 10
      29             : static const long HASHT = 1L << HASHBITS;
      30             : 
      31             : static long
      32     2557304 : hash(long q) { return (q & ((1L << (HASHBITS+1)) - 1)) >> 1; }
      33             : #undef HASHBITS
      34             : 
      35             : /* See buch2.c:
      36             :  * B->subFB contains split p such that \prod p > sqrt(B->Disc)
      37             :  * B->powsubFB contains powers of forms in B->subFB */
      38             : #define RANDOM_BITS 4
      39             : static const long CBUCH = (1L<<RANDOM_BITS)-1;
      40             : 
      41             : struct buch_quad
      42             : {
      43             :   ulong limhash;
      44             :   long KC, KC2, PRECREG;
      45             :   long *primfact, *exprimfact, **hashtab;
      46             :   GEN FB, numFB, prodFB;
      47             :   GEN powsubFB, vperm, subFB, badprim;
      48             :   struct qfr_data *q;
      49             : };
      50             : 
      51             : /*******************************************************************/
      52             : /*                                                                 */
      53             : /*  Routines related to binary quadratic forms (for internal use)  */
      54             : /*                                                                 */
      55             : /*******************************************************************/
      56             : /* output canonical representative wrt projection Cl^+ --> Cl (a > 0) */
      57             : static GEN
      58     1164331 : qfr3_canon(GEN x, struct qfr_data *S)
      59             : {
      60     1164331 :   GEN a = gel(x,1), c = gel(x,3);
      61     1164331 :   if (signe(a) < 0) {
      62      502173 :     if (absequalii(a,c)) return qfr3_rho(x, S);
      63      502166 :     setsigne(a, 1);
      64      502166 :     setsigne(c,-1);
      65             :   }
      66     1164324 :   return x;
      67             : }
      68             : static GEN
      69        3710 : qfr3_canon_safe(GEN x, struct qfr_data *S)
      70             : {
      71        3710 :   GEN a = gel(x,1), c = gel(x,3);
      72        3710 :   if (signe(a) < 0) {
      73         224 :     if (absequalii(a,c)) return qfr3_rho(x, S);
      74         224 :     gel(x,1) = negi(a);
      75         224 :     gel(x,3) = negi(c);
      76             :   }
      77        3710 :   return x;
      78             : }
      79             : static GEN
      80     1947281 : qfr5_canon(GEN x, struct qfr_data *S)
      81             : {
      82     1947281 :   GEN a = gel(x,1), c = gel(x,3);
      83     1947281 :   if (signe(a) < 0) {
      84      868427 :     if (absequalii(a,c)) return qfr5_rho(x, S);
      85      868420 :     setsigne(a, 1);
      86      868420 :     setsigne(c,-1);
      87             :   }
      88     1947274 :   return x;
      89             : }
      90             : static GEN
      91     1730337 : QFR5_comp(GEN x,GEN y, struct qfr_data *S)
      92     1730337 : { return qfr5_canon(qfr5_comp(x,y,S), S); }
      93             : static GEN
      94     1005354 : QFR3_comp(GEN x, GEN y, struct qfr_data *S)
      95     1005354 : { return qfr3_canon(qfr3_comp(x,y,S), S); }
      96             : 
      97             : /* compute rho^n(x) */
      98             : static GEN
      99      220584 : qfr5_rho_pow(GEN x, long n, struct qfr_data *S)
     100             : {
     101             :   long i;
     102      220584 :   pari_sp av = avma;
     103     2216151 :   for (i=1; i<=n; i++)
     104             :   {
     105     1995567 :     x = qfr5_rho(x,S);
     106     1995567 :     if (gc_needed(av,1))
     107             :     {
     108           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"qfr5_rho_pow");
     109           0 :       x = gerepilecopy(av, x);
     110             :     }
     111             :   }
     112      220584 :   return gerepilecopy(av, x);
     113             : }
     114             : 
     115             : static GEN
     116      216944 : qfr5_pf(struct qfr_data *S, long p, long prec)
     117             : {
     118      216944 :   GEN y = primeform_u(S->D,p);
     119      216944 :   return qfr5_canon(qfr5_red(qfr_to_qfr5(y,prec), S), S);
     120             : }
     121             : 
     122             : static GEN
     123      158956 : qfr3_pf(struct qfr_data *S, long p)
     124             : {
     125      158956 :   GEN y = primeform_u(S->D,p);
     126      158956 :   return qfr3_canon(qfr3_red(y, S), S);
     127             : }
     128             : 
     129             : #define qfi_pf primeform_u
     130             : 
     131             : /* Warning: ex[0] not set in general */
     132             : static GEN
     133     4044133 : init_form(struct buch_quad *B, GEN ex,
     134             :           GEN (*comp)(GEN,GEN,struct qfr_data *S))
     135             : {
     136     4044133 :   long i, l = lg(B->powsubFB);
     137     4044133 :   GEN F = NULL;
     138    22848531 :   for (i=1; i<l; i++)
     139    18814018 :     if (ex[i])
     140             :     {
     141    17640220 :       GEN t = gmael(B->powsubFB,i,ex[i]);
     142    17640220 :       F = F? comp(F, t, B->q): t;
     143             :     }
     144     4034513 :   return F;
     145             : }
     146             : static GEN
     147      197442 : qfr5_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFR5_comp); }
     148             : 
     149             : static GEN
     150    11714534 : QFI_comp(GEN x, GEN y, struct qfr_data *S) { (void)S; return qfbcomp_i(x,y); }
     151             : 
     152             : static GEN
     153      158855 : qfi_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFI_comp); }
     154             : 
     155             : static GEN
     156     3686753 : random_form(struct buch_quad *B, GEN ex,
     157             :             GEN (*comp)(GEN,GEN, struct qfr_data *S))
     158             : {
     159     3686753 :   long i, l = lg(ex);
     160     3686753 :   pari_sp av = avma;
     161             :   GEN F;
     162             :   for(;;)
     163             :   {
     164    20529227 :     for (i=1; i<l; i++) ex[i] = random_bits(RANDOM_BITS);
     165     3688037 :     if ((F = init_form(B, ex, comp))) return F;
     166         448 :     set_avma(av);
     167             :   }
     168             : }
     169             : static GEN
     170      161133 : qfr3_random(struct buch_quad *B,GEN ex){ return random_form(B, ex, &QFR3_comp); }
     171             : static GEN
     172     3525750 : qfi_random(struct buch_quad *B,GEN ex) { return random_form(B, ex, &QFI_comp); }
     173             : 
     174             : /*******************************************************************/
     175             : /*                                                                 */
     176             : /*                     Common subroutines                          */
     177             : /*                                                                 */
     178             : /*******************************************************************/
     179             : /* We assume prime ideals of norm < D generate Cl(K); failure with
     180             :  * a factor base of primes with norm < LIMC <= D. Suggest new value.
     181             :  * All values of the form c * log^2 (disc K) [where D has c = 4
     182             :  * (Grenie-Molteni, under GRH)]. A value of c in [0.3, 1] should be OK for most
     183             :  * fields. No example is known where c >= 2 is necessary. */
     184             : ulong
     185        2349 : bnf_increase_LIMC(ulong LIMC, ulong D)
     186             : {
     187        2349 :   if (LIMC >= D) pari_err_BUG("Buchmann's algorithm");
     188        2349 :   if (LIMC <= D / 13.333)
     189         207 :     LIMC *= 2; /* tiny c <= 0.3 : double it */
     190             :   else
     191        2142 :     LIMC += maxuu(1, D / 20); /* large c, add 0.2 to it */
     192        2349 :   if (LIMC > D) LIMC = D;
     193        2349 :   return LIMC;
     194             : }
     195             : 
     196             : /* Is |q| <= p ? */
     197             : static int
     198       10186 : isless_iu(GEN q, ulong p) {
     199       10186 :   long l = lgefint(q);
     200       10186 :   return l==2 || (l == 3 && uel(q,2) <= p);
     201             : }
     202             : 
     203             : static GEN
     204     5015848 : Z_isquasismooth_prod(GEN N, GEN P)
     205             : {
     206     5015848 :   P = gcdii(P,N);
     207    10005705 :   while (!is_pm1(P))
     208             :   {
     209     4983147 :     N = diviiexact(N, P);
     210     4996729 :     P = gcdii(N, P);
     211             :   }
     212     5007371 :   return N;
     213             : }
     214             : 
     215             : static long
     216     5016867 : factorquad(struct buch_quad *B, GEN f, long nFB, ulong limp)
     217             : {
     218             :   ulong X;
     219     5016867 :   long i, lo = 0;
     220     5016867 :   GEN F, x = gel(f,1), FB = B->FB, P = B->primfact, E = B->exprimfact;
     221     5016867 :   if (B->badprim && !is_pm1(gcdii(x, B->badprim))) return 0;
     222     5015691 :   F =  Z_isquasismooth_prod(x, B->prodFB);
     223     5007671 :   if (cmpiu(F, B->limhash) > 0) return 0;
     224     4390725 :   for (i=1; lgefint(x) > 3; i++)
     225             :   {
     226       10186 :     ulong p = uel(FB,i), r;
     227       10186 :     GEN q = absdiviu_rem(x, p, &r);
     228       10186 :     if (!r)
     229             :     {
     230        1732 :       long k = 0;
     231        2952 :       do { k++; x = q; q = absdiviu_rem(x, p, &r); } while (!r);
     232        1732 :       lo++; P[lo] = p; E[lo] = k;
     233             :     }
     234       10186 :     if (isless_iu(q,p)) {
     235           1 :       if (lgefint(x) == 3) { X = uel(x,2); goto END; }
     236          32 :       return 0;
     237             :     }
     238       10185 :     if (i == nFB) return 0;
     239             :   }
     240     4380539 :   X = uel(x,2);
     241     4380539 :   if (X == 1) { P[0] = 0; return 1; }
     242    68516132 :   for (;; i++)
     243    68516132 :   { /* single precision affair, split for efficiency */
     244    72878052 :     ulong p = uel(FB,i);
     245    72878052 :     ulong q = X / p, r = X % p; /* gcc makes a single div */
     246    72878052 :     if (!r)
     247             :     {
     248     5619324 :       long k = 0;
     249     6784313 :       do { k++; X = q; q = X / p; r = X % p; } while (!r);
     250     5619324 :       lo++; P[lo] = p; E[lo] = k;
     251             :     }
     252    72878052 :     if (q <= p) break;
     253    68551148 :     if (i == nFB) return 0;
     254             :   }
     255     4326905 : END:
     256     4326905 :   if (X > B->limhash) return 0;
     257     4326905 :   if (X != 1 && X <= limp) { lo++; P[lo] = X; E[lo] = 1; X = 1; }
     258     4326905 :   P[0] = lo; return X;
     259             : }
     260             : 
     261             : /* Check for a "large prime relation" involving q; q may not be prime */
     262             : static long *
     263     2557351 : largeprime(struct buch_quad *B, long q, GEN ex, long np, long nrho)
     264             : {
     265     2557351 :   const long hashv = hash(q);
     266     2557290 :   long *pt, i, l = lg(B->subFB);
     267             : 
     268     2779622 :   for (pt = B->hashtab[hashv]; ; pt = (long*) pt[0])
     269             :   {
     270     2779622 :     if (!pt)
     271             :     {
     272     2358245 :       pt = (long*) pari_malloc((l+3) * sizeof(long));
     273     2358327 :       *pt++ = nrho; /* nrho = pt[-3] */
     274     2358327 :       *pt++ = np;   /* np   = pt[-2] */
     275     2358327 :       *pt++ = q;    /* q    = pt[-1] */
     276     2358327 :       pt[0] = (long)B->hashtab[hashv];
     277    14948231 :       for (i=1; i<l; i++) pt[i]=ex[i];
     278     2358327 :       B->hashtab[hashv]=pt; return NULL;
     279             :     }
     280      421377 :     if (pt[-1] == q) break;
     281             :   }
     282      242472 :   for(i=1; i<l; i++)
     283      238600 :     if (pt[i] != ex[i]) return pt;
     284        3872 :   return (pt[-2]==np)? NULL: pt;
     285             : }
     286             : 
     287             : static void
     288      169411 : clearhash(long **hash)
     289             : {
     290             :   long *pt;
     291             :   long i;
     292   173606596 :   for (i=0; i<HASHT; i++) {
     293   175795617 :     for (pt = hash[i]; pt; ) {
     294     2358432 :       void *z = (void*)(pt - 3);
     295     2358432 :       pt = (long*) pt[0]; pari_free(z);
     296             :     }
     297   173437185 :     hash[i] = NULL;
     298             :   }
     299      169770 : }
     300             : 
     301             : /* last prime stored */
     302             : ulong
     303           0 : GRH_last_prime(GRHcheck_t *S) { return (S->primes + S->nprimes-1)->p; }
     304             : /* ensure that S->primes can hold at least nb primes */
     305             : void
     306      399930 : GRH_ensure(GRHcheck_t *S, long nb)
     307             : {
     308      399930 :   if (S->maxprimes <= nb)
     309             :   {
     310           0 :     do S->maxprimes *= 2; while (S->maxprimes <= nb);
     311           0 :     pari_realloc_ip((void**)&S->primes, S->maxprimes*sizeof(*S->primes));
     312             :   }
     313      399930 : }
     314             : /* cache data for all primes up to the LIM */
     315             : static void
     316     1173412 : cache_prime_quad(GRHcheck_t *S, ulong LIM, GEN D)
     317             : {
     318             :   GRHprime_t *pr;
     319             :   long nb;
     320             : 
     321     1173412 :   if (S->limp >= LIM) return;
     322       72319 :   nb = (long)primepi_upper_bound((double)LIM); /* #{p <= LIM} <= nb */
     323       72334 :   GRH_ensure(S, nb+1); /* room for one extra prime */
     324       72334 :   for (pr = S->primes + S->nprimes;;)
     325    12647001 :   {
     326    12719335 :     ulong p = u_forprime_next(&(S->P));
     327    12719121 :     pr->p = p;
     328    12719121 :     pr->logp = log((double)p);
     329    12719121 :     pr->dec = (GEN)kroiu(D,p);
     330    12719338 :     S->nprimes++;
     331    12719338 :     pr++;
     332             :     /* store up to nextprime(LIM) included */
     333    12719338 :     if (p >= LIM) { S->limp = p; break; }
     334             :   }
     335             : }
     336             : 
     337             : static GEN
     338       70257 : compute_invresquad(GRHcheck_t *S, long LIMC)
     339             : {
     340       70257 :   pari_sp av = avma;
     341       70257 :   GEN invres = real_1(DEFAULTPREC);
     342       70254 :   double limp = log((double)LIMC) / 2;
     343             :   GRHprime_t *pr;
     344             :   long i;
     345    12718589 :   for (pr = S->primes, i = S->nprimes; i > 0; pr++, i--)
     346             :   {
     347    12648987 :     long s = (long)pr->dec;
     348    12648987 :     if (s)
     349             :     {
     350    12522664 :       ulong p = pr->p;
     351    12522664 :       if (s > 0 || pr->logp <= limp)
     352             :         /* Both p and P contribute */
     353     6369811 :         invres = mulur(p - s, divru(invres, p));
     354     6152853 :       else if (s<0)
     355             :         /* Only p contributes */
     356     6152875 :         invres = mulur(p, divru(invres, p - 1));
     357             :     }
     358             :   }
     359       69602 :   return gerepileuptoleaf(av, invres);
     360             : }
     361             : 
     362             : /* p | conductor of order of disc D ? */
     363             : static int
     364      391522 : is_bad(GEN D, ulong p)
     365             : {
     366      391522 :   pari_sp av = avma;
     367      391522 :   if (p == 2)
     368             :   {
     369       89557 :     long r = mod16(D) >> 1;
     370       89557 :     if (r && signe(D) < 0) r = 8-r;
     371       89557 :     return (r < 4);
     372             :   }
     373      301965 :   return gc_bool(av, dvdii(D, sqru(p))); /* p^2 | D ? */
     374             : }
     375             : 
     376             : /* returns the n-th suitable ideal for the factorbase */
     377             : static long
     378       70258 : nthidealquad(GEN D, long n)
     379             : {
     380       70258 :   pari_sp av = avma;
     381             :   forprime_t S;
     382             :   ulong p;
     383       70258 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
     384      357173 :   while ((p = u_forprime_next(&S)) && n > 0)
     385      286924 :     if (!is_bad(D, p) && kroiu(D, p) >= 0) n--;
     386       70250 :   return gc_long(av, p);
     387             : }
     388             : 
     389             : static int
     390     1032787 : quadGRHchk(GEN D, GRHcheck_t *S, ulong LIMC)
     391             : {
     392     1032787 :   double logC = log((double)LIMC), SA = 0, SB = 0;
     393             :   long i;
     394             : 
     395     1032787 :   cache_prime_quad(S, LIMC, D);
     396     1032784 :   for (i = 0;; i++)
     397    29332614 :   {
     398    30365398 :     GRHprime_t *pr = S->primes+i;
     399    30365398 :     ulong p = pr->p;
     400             :     long M;
     401             :     double logNP, q, A, B;
     402    30365398 :     if (p > LIMC) break;
     403    29421716 :     if ((long)pr->dec < 0)
     404             :     {
     405    14693249 :       logNP = 2 * pr->logp;
     406    14693249 :       q = 1/(double)p;
     407             :     }
     408             :     else
     409             :     {
     410    14728467 :       logNP = pr->logp;
     411    14728467 :       q = 1/sqrt((double)p);
     412             :     }
     413    29421716 :     A = logNP * q; B = logNP * A; M = (long)(logC/logNP);
     414    29421716 :     if (M > 1)
     415             :     {
     416     2496985 :       double inv1_q = 1 / (1-q);
     417     2496985 :       A *= (1 - pow(q, (double) M)) * inv1_q;
     418     2496985 :       B *= (1 - pow(q, (double) M)*(M+1 - M*q)) * inv1_q * inv1_q;
     419             :     }
     420    29421716 :     if ((long)pr->dec>0) { SA += 2*A;SB += 2*B; } else { SA += A; SB += B; }
     421    29421716 :     if (p == LIMC) break;
     422             :   }
     423     1032784 :   return GRHok(S, logC, SA, SB);
     424             : }
     425             : 
     426             : /* C2 >= C1; create B->FB, B->numFB; set B->badprim. Return L(kro_D, 1) */
     427             : static void
     428       70391 : FBquad(struct buch_quad *B, ulong C2, ulong C1, GRHcheck_t *S)
     429             : {
     430       70391 :   GEN D = B->q->D;
     431             :   long i;
     432             :   pari_sp av;
     433             :   GRHprime_t *pr;
     434             : 
     435       70391 :   cache_prime_quad(S, C2, D);
     436       70390 :   pr = S->primes;
     437       70390 :   B->numFB = cgetg(C2+1, t_VECSMALL);
     438       70390 :   B->FB    = cgetg(C2+1, t_VECSMALL);
     439       70402 :   av = avma;
     440       70402 :   B->KC = 0; i = 0;
     441       70402 :   B->badprim = gen_1;
     442     2798889 :   for (;; pr++) /* p <= C2 */
     443     2798889 :   {
     444     2869291 :     ulong p = pr->p;
     445     2869291 :     if (!B->KC && p > C1) B->KC = i;
     446     2869291 :     if (p > C2) break;
     447     2806417 :     switch ((long)pr->dec)
     448             :     {
     449     1382997 :       case -1: break; /* inert */
     450      104602 :       case  0: /* ramified */
     451      104602 :         if (is_bad(D, p)) { B->badprim = muliu(B->badprim, p); break; }
     452             :         /* fall through */
     453             :       default:  /* split */
     454     1423381 :         i++; B->numFB[p] = i; B->FB[i] = p; break;
     455             :     }
     456     2806406 :     if (p == C2)
     457             :     {
     458        7517 :       if (!B->KC) B->KC = i;
     459        7517 :       break;
     460             :     }
     461             :   }
     462       70391 :   B->KC2 = i;
     463       70391 :   setlg(B->FB, B->KC2+1);
     464       70391 :   if (B->badprim != gen_1)
     465          21 :     B->badprim = gerepileuptoint(av, B->badprim);
     466             :   else
     467             :   {
     468       70370 :     B->badprim = NULL; set_avma(av);
     469             :   }
     470       70391 :   B->prodFB = zv_prod_Z(B->FB);
     471       70390 : }
     472             : 
     473             : /* create B->vperm, return B->subFB */
     474             : static GEN
     475       70390 : subFBquad(struct buch_quad *B, GEN D, double PROD, long minSFB)
     476             : {
     477       70390 :   long i, j, lgsub = 1, ino = 1, lv = B->KC+1;
     478       70390 :   double prod = 1.;
     479             :   pari_sp av;
     480             :   GEN no;
     481             : 
     482       70390 :   B->vperm = cgetg(lv, t_VECSMALL);
     483       70389 :   av = avma;
     484       70389 :   no    = cgetg(lv, t_VECSMALL);
     485      334514 :   for (j = 1; j < lv; j++)
     486             :   {
     487      334374 :     ulong p = uel(B->FB,j);
     488      334374 :     if (!umodiu(D, p)) no[ino++] = j; /* ramified */
     489             :     else
     490             :     {
     491      254583 :       B->vperm[lgsub++] = j;
     492      254583 :       prod *= p;
     493      254583 :       if (lgsub > minSFB && prod > PROD) break;
     494             :     }
     495             :   }
     496             :   /* lgsub >= 1 otherwise quadGRHchk is false */
     497       70390 :   i = lgsub;
     498      150182 :   for (j = 1; j < ino;i++,j++) B->vperm[i] = no[j];
     499     1158913 :   for (     ; i < lv; i++)     B->vperm[i] = i;
     500       70390 :   no = gclone(vecslice(B->vperm, 1, lgsub-1));
     501       70389 :   set_avma(av); return no;
     502             : }
     503             : 
     504             : /* assume n >= 1, x[i][j] = B->subFB[i]^j, for j = 1..n */
     505             : static GEN
     506       99024 : powsubFBquad(struct buch_quad *B, long n)
     507             : {
     508       99024 :   pari_sp av = avma;
     509       99024 :   long i,j, l = lg(B->subFB);
     510       99024 :   GEN F, y, x = cgetg(l, t_VEC), D = B->q->D;
     511             : 
     512       99025 :   if (B->PRECREG) /* real */
     513             :   {
     514       39088 :     for (i=1; i<l; i++)
     515             :     {
     516       34055 :       F = qfr5_pf(B->q, B->FB[B->subFB[i]], B->PRECREG);
     517       34055 :       y = cgetg(n+1, t_VEC); gel(x,i) = y;
     518       34055 :       gel(y,1) = F;
     519      544880 :       for (j=2; j<=n; j++) gel(y,j) = QFR5_comp(gel(y,j-1), F, B->q);
     520             :     }
     521             :   }
     522             :   else /* imaginary */
     523             :   {
     524      426765 :     for (i=1; i<l; i++)
     525             :     {
     526      332787 :       F = qfi_pf(D, B->FB[B->subFB[i]]);
     527      332771 :       y = cgetg(n+1, t_VEC); gel(x,i) = y;
     528      333938 :       gel(y,1) = F;
     529     5316070 :       for (j=2; j<=n; j++) gel(y,j) = qfbcomp_i(gel(y,j-1), F);
     530             :     }
     531             :   }
     532       99011 :   x = gclone(x); set_avma(av); return x;
     533             : }
     534             : 
     535             : static void
     536       97763 : sub_fact(struct buch_quad *B, GEN col, GEN F)
     537             : {
     538       97763 :   GEN b = gel(F,2);
     539             :   long i;
     540      207856 :   for (i=1; i<=B->primfact[0]; i++)
     541             :   {
     542      110093 :     ulong p = B->primfact[i], k = B->numFB[p];
     543      110093 :     long e = B->exprimfact[i];
     544      110093 :     if (umodiu(b, p<<1) > p) e = -e;
     545      110093 :     col[k] -= e;
     546             :   }
     547       97763 : }
     548             : static void
     549     1885132 : add_fact(struct buch_quad *B, GEN col, GEN F)
     550             : {
     551     1885132 :   GEN b = gel(F,2);
     552             :   long i;
     553     5915690 :   for (i=1; i<=B->primfact[0]; i++)
     554             :   {
     555     4030444 :     ulong p = B->primfact[i], k = B->numFB[p];
     556     4030444 :     long e = B->exprimfact[i];
     557     4030444 :     if (umodiu(b, p<<1) > p) e = -e;
     558     4030558 :     col[k] += e;
     559             :   }
     560     1885246 : }
     561             : 
     562             : static GEN
     563       70257 : get_clgp(struct buch_quad *B, GEN W, GEN *ptD)
     564             : {
     565       70257 :   GEN res, init, u1, D = ZM_snf_group(W,NULL,&u1);
     566       70253 :   long i, j, l = lg(W), c = lg(D);
     567             : 
     568       70253 :   res=cgetg(c,t_VEC); init = cgetg(l,t_VEC);
     569      215966 :   for (i=1; i<l; i++) gel(init,i) = primeform_u(B->q->D, B->FB[B->vperm[i]]);
     570      195559 :   for (j=1; j<c; j++)
     571             :   {
     572      125300 :     GEN g = NULL;
     573      125300 :     if (signe(B->q->D) > 0)
     574             :     {
     575       13293 :       for (i=1; i<l; i++)
     576             :       {
     577        9583 :         GEN t, u = gcoeff(u1,i,j);
     578        9583 :         if (!signe(u)) continue;
     579        4501 :         t = qfr3_pow(gel(init,i), u, B->q);
     580        4501 :         g = g? qfr3_comp(g, t, B->q): t;
     581             :       }
     582        3710 :       g = qfr3_to_qfr(qfr3_canon_safe(qfr3_red(g, B->q), B->q), B->q->D);
     583             :     }
     584             :     else
     585             :     {
     586      422789 :       for (i=1; i<l; i++)
     587             :       {
     588      301198 :         GEN t, u = gcoeff(u1,i,j);
     589      301198 :         if (!signe(u)) continue;
     590      203219 :         t = powgi(gel(init,i), u);
     591      203225 :         g = g? qfbcomp_i(g, t): t;
     592             :       }
     593             :     }
     594      125302 :     gel(res,j) = g;
     595             :   }
     596       70259 :   *ptD = D; return res;
     597             : }
     598             : 
     599             : static long
     600       70256 : trivial_relations(struct buch_quad *B, GEN mat, GEN C)
     601             : {
     602       70256 :   long i, j = 0;
     603       70256 :   GEN col, D = B->q->D;
     604     1492631 :   for (i = 1; i <= B->KC; i++)
     605             :   { /* ramified prime ==> trivial relation */
     606     1422359 :     if (umodiu(D, B->FB[i])) continue;
     607      104203 :     col = zero_zv(B->KC);
     608      104215 :     col[i] = 2; j++;
     609      104215 :     gel(mat,j) = col;
     610      104215 :     gel(C,j) = gen_0;
     611             :   }
     612       70272 :   return j;
     613             : }
     614             : 
     615             : static void
     616           0 : dbg_all(pari_timer *T, const char *phase, long s, long n)
     617             : {
     618           0 :   err_printf("\n");
     619           0 :   timer_printf(T, "%s rel [#rel/#test = %ld/%ld]", phase,s,n);
     620           0 : }
     621             : 
     622             : /* Imaginary Quadratic fields */
     623             : 
     624             : static void
     625       98711 : imag_relations(struct buch_quad *B, long need, long *pc, ulong LIMC, GEN mat)
     626             : {
     627             :   pari_timer T;
     628       98711 :   long lgsub = lg(B->subFB), current = *pc, nbtest = 0, s = 0;
     629             :   long i, fpc;
     630             :   pari_sp av;
     631       98711 :   GEN col, form, ex = cgetg(lgsub, t_VECSMALL);
     632             : 
     633       98711 :   if (!current) current = 1;
     634       98711 :   if (DEBUGLEVEL>2) timer_start(&T);
     635       98711 :   av = avma;
     636             :   for(;;)
     637             :   {
     638     3624456 :     if (s >= need) break;
     639     3525746 :     set_avma(av);
     640     3525331 :     form = qfi_random(B,ex);
     641     3523720 :     form = qfbcomp_i(form, qfi_pf(B->q->D, B->FB[current]));
     642     3523572 :     nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
     643     3525913 :     if (!fpc)
     644             :     {
     645      287960 :       if (DEBUGLEVEL>3) err_printf(".");
     646      287960 :       if ((nbtest & 0xff) == 0 && ++current > B->KC) current = 1;
     647      287960 :       continue;
     648             :     }
     649     3237953 :     if (fpc > 1)
     650             :     {
     651     1793301 :       long *fpd = largeprime(B,fpc,ex,current,0);
     652             :       ulong b1, b2, p;
     653             :       GEN form2;
     654     1793310 :       if (!fpd)
     655             :       {
     656     1634471 :         if (DEBUGLEVEL>3) err_printf(".");
     657     1634407 :         continue;
     658             :       }
     659      158839 :       form2 = qfbcomp_i(qfi_factorback(B,fpd), qfi_pf(B->q->D, B->FB[fpd[-2]]));
     660      158850 :       p = fpc << 1;
     661      158850 :       b1 = umodiu(gel(form2,2), p);
     662      158853 :       b2 = umodiu(gel(form,2),  p);
     663      158860 :       if (b1 != b2 && b1+b2 != p) continue;
     664             : 
     665      158860 :       col = gel(mat,++s);
     666      158860 :       add_fact(B,col, form);
     667      158859 :       (void)factorquad(B,form2,B->KC,LIMC);
     668      158861 :       if (b1==b2)
     669             :       {
     670      409440 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
     671       79528 :         sub_fact(B, col, form2); col[fpd[-2]]++;
     672             :       }
     673             :       else
     674             :       {
     675      408633 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
     676       79333 :         add_fact(B, col, form2); col[fpd[-2]]--;
     677             :       }
     678      158861 :       if (DEBUGLEVEL>2) err_printf(" %ldP",s);
     679             :     }
     680             :     else
     681             :     {
     682     1444652 :       col = gel(mat,++s);
     683     6943310 :       for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
     684     1444652 :       add_fact(B, col, form);
     685     1444755 :       if (DEBUGLEVEL>2) err_printf(" %ld",s);
     686             :     }
     687     1603378 :     col[current]--;
     688     1603378 :     if (++current > B->KC) current = 1;
     689             :   }
     690       98710 :   if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
     691       98710 :   *pc = current;
     692       98710 : }
     693             : 
     694             : static int
     695           7 : imag_be_honest(struct buch_quad *B)
     696             : {
     697           7 :   long p, fpc, s = B->KC, nbtest = 0;
     698           7 :   GEN F, ex = cgetg(lg(B->subFB), t_VECSMALL);
     699           7 :   pari_sp av = avma;
     700             : 
     701         525 :   while (s<B->KC2)
     702             :   {
     703         518 :     p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
     704         518 :     F = qfbcomp_i(qfi_pf(B->q->D, p), qfi_random(B, ex));
     705         518 :     fpc = factorquad(B,F,s,p-1);
     706         518 :     if (fpc == 1) { nbtest=0; s++; }
     707             :     else
     708         392 :       if (++nbtest > 40) return 0;
     709         518 :     set_avma(av);
     710             :   }
     711           7 :   return 1;
     712             : }
     713             : 
     714             : static GEN
     715      184254 : dist(GEN e, GEN d, long prec)
     716             : {
     717      184254 :   GEN t = qfr5_dist(e, d, prec);
     718      184254 :   return signe(d) < 0 ? mkcomplex(t, gen_1): t;
     719             : }
     720             : 
     721             : /* Real Quadratic fields */
     722             : 
     723             : static void
     724        5124 : real_relations(struct buch_quad *B, long need, long *pc, long lim, ulong LIMC, GEN mat, GEN C)
     725             : {
     726             :   pari_timer T;
     727        5124 :   long lgsub = lg(B->subFB), prec = B->PRECREG, current = *pc, nbtest=0, s=0;
     728             :   long i, fpc, endcycle, rhoacc, rho;
     729             :   /* in a 2nd phase, don't include FB[current] but run along the cyle
     730             :    * ==> get more units */
     731        5124 :   int first = (current == 0);
     732             :   pari_sp av, av1;
     733        5124 :   GEN d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL);
     734             : 
     735        5124 :   if (DEBUGLEVEL>2) timer_start(&T);
     736        5124 :   if (!current) current = 1;
     737        5124 :   if (lim > need) lim = need;
     738        5124 :   av = avma;
     739             :   for(;;)
     740             :   {
     741      166236 :     if (s >= need) break;
     742      161112 :     if (first && s >= lim) {
     743        2142 :       first = 0;
     744        2142 :       if (DEBUGLEVEL>2) dbg_all(&T, "initial", s, nbtest);
     745             :     }
     746      161112 :     set_avma(av); form = qfr3_random(B, ex);
     747      161112 :     if (!first)
     748      158935 :       form = QFR3_comp(form, qfr3_pf(B->q, B->FB[current]), B->q);
     749      161112 :     av1 = avma;
     750      161112 :     form0 = form; form1 = NULL;
     751      161112 :     endcycle = rhoacc = 0;
     752      161112 :     rho = -1;
     753             : 
     754     1297569 : CYCLE:
     755     1297569 :     if (endcycle || rho > 5000)
     756             :     {
     757          21 :       if (++current > B->KC) current = 1;
     758          21 :       continue;
     759             :     }
     760     1297548 :     if (gc_needed(av,1))
     761             :     {
     762           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"real_relations");
     763           0 :       gerepileall(av1, form1? 2: 1, &form, &form1);
     764             :     }
     765     1297548 :     if (rho < 0) rho = 0; /* first time in */
     766             :     else
     767             :     {
     768     1136436 :       form = qfr3_rho(form, B->q); rho++;
     769     1136436 :       rhoacc++;
     770     1136436 :       if (first)
     771      444430 :         endcycle = (absequalii(gel(form,1),gel(form0,1))
     772      222215 :              && equalii(gel(form,2),gel(form0,2)));
     773             :       else
     774             :       {
     775      914221 :         if (absequalii(gel(form,1), gel(form,3))) /* a = -c */
     776             :         {
     777           0 :           if (absequalii(gel(form,1),gel(form0,1)) &&
     778           0 :                   equalii(gel(form,2),gel(form0,2))) continue;
     779           0 :           form = qfr3_rho(form, B->q); rho++;
     780           0 :           rhoacc++;
     781             :         }
     782             :         else
     783      914221 :           { setsigne(form[1],1); setsigne(form[3],-1); }
     784      914277 :         if (equalii(gel(form,1),gel(form0,1)) &&
     785          56 :             equalii(gel(form,2),gel(form0,2))) continue;
     786             :       }
     787             :     }
     788     1297548 :     nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
     789     1297548 :     if (!fpc)
     790             :     {
     791      385511 :       if (DEBUGLEVEL>3) err_printf(".");
     792      385511 :       goto CYCLE;
     793             :     }
     794      912037 :     if (fpc > 1)
     795             :     { /* look for Large Prime relation */
     796      764113 :       long *fpd = largeprime(B,fpc,ex,first? 0: current,rhoacc);
     797             :       ulong b1, b2, p;
     798             :       GEN form2;
     799      764113 :       if (!fpd)
     800             :       {
     801      727783 :         if (DEBUGLEVEL>3) err_printf(".");
     802      727783 :         goto CYCLE;
     803             :       }
     804       36330 :       if (!form1)
     805             :       {
     806       36330 :         form1 = qfr5_factorback(B,ex);
     807       36330 :         if (!first)
     808       36330 :           form1 = QFR5_comp(form1, qfr5_pf(B->q, B->FB[current], prec), B->q);
     809             :       }
     810       36330 :       form1 = qfr5_rho_pow(form1, rho, B->q);
     811       36330 :       rho = 0;
     812             : 
     813       36330 :       form2 = qfr5_factorback(B,fpd);
     814       36330 :       if (fpd[-2])
     815       23954 :         form2 = QFR5_comp(form2, qfr5_pf(B->q, B->FB[fpd[-2]], prec), B->q);
     816       36330 :       form2 = qfr5_rho_pow(form2, fpd[-3], B->q);
     817       36330 :       if (!absequalii(gel(form2,1),gel(form2,3)))
     818             :       {
     819       36330 :         setsigne(form2[1], 1);
     820       36330 :         setsigne(form2[3],-1);
     821             :       }
     822       36330 :       p = fpc << 1;
     823       36330 :       b1 = umodiu(gel(form2,2), p);
     824       36330 :       b2 = umodiu(gel(form1,2), p);
     825       36330 :       if (b1 != b2 && b1+b2 != p) goto CYCLE;
     826             : 
     827       36330 :       col = gel(mat,++s);
     828       36330 :       add_fact(B, col, form1);
     829       36330 :       (void)factorquad(B,form2,B->KC,LIMC);
     830       36330 :       if (b1==b2)
     831             :       {
     832      139377 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
     833       18235 :         sub_fact(B,col, form2);
     834       18235 :         if (fpd[-2]) col[fpd[-2]]++;
     835       18235 :         d = dist(subii(gel(form1,4),gel(form2,4)),
     836       18235 :                       divrr(gel(form1,5),gel(form2,5)), prec);
     837             :       }
     838             :       else
     839             :       {
     840      138670 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
     841       18095 :         add_fact(B, col, form2);
     842       18095 :         if (fpd[-2]) col[fpd[-2]]--;
     843       18095 :         d = dist(addii(gel(form1,4),gel(form2,4)),
     844       18095 :                       mulrr(gel(form1,5),gel(form2,5)), prec);
     845             :       }
     846       36330 :       if (DEBUGLEVEL>2) err_printf(" %ldP",s);
     847             :     }
     848             :     else
     849             :     { /* standard relation */
     850      147924 :       if (!form1)
     851             :       {
     852      124782 :         form1 = qfr5_factorback(B, ex);
     853      124782 :         if (!first)
     854      122605 :           form1 = QFR5_comp(form1, qfr5_pf(B->q, B->FB[current], prec), B->q);
     855             :       }
     856      147924 :       form1 = qfr5_rho_pow(form1, rho, B->q);
     857      147924 :       rho = 0;
     858             : 
     859      147924 :       col = gel(mat,++s);
     860     1138494 :       for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
     861      147924 :       add_fact(B, col, form1);
     862      147924 :       d = dist(gel(form1,4), gel(form1,5), prec);
     863      147924 :       if (DEBUGLEVEL>2) err_printf(" %ld",s);
     864             :     }
     865      184254 :     gaffect(d, gel(C,s));
     866      184254 :     if (first)
     867             :     {
     868       25319 :       if (s >= lim) continue;
     869       23163 :       goto CYCLE;
     870             :     }
     871             :     else
     872             :     {
     873      158935 :       col[current]--;
     874      158935 :       if (++current > B->KC) current = 1;
     875             :     }
     876             :   }
     877        5124 :   if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
     878        5124 :   *pc = current;
     879        5124 : }
     880             : 
     881             : static int
     882           7 : real_be_honest(struct buch_quad *B)
     883             : {
     884           7 :   long p, fpc, s = B->KC, nbtest = 0;
     885           7 :   GEN F,F0, ex = cgetg(lg(B->subFB), t_VECSMALL);
     886           7 :   pari_sp av = avma;
     887             : 
     888          28 :   while (s<B->KC2)
     889             :   {
     890          21 :     p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
     891          21 :     F = QFR3_comp(qfr3_random(B, ex), qfr3_pf(B->q, p), B->q);
     892          21 :     for (F0 = F;;)
     893             :     {
     894          42 :       fpc = factorquad(B,F,s,p-1);
     895          42 :       if (fpc == 1) { nbtest=0; s++; break; }
     896          21 :       if (++nbtest > 40) return 0;
     897          21 :       F = qfr3_canon(qfr3_rho(F, B->q), B->q);
     898          21 :       if (equalii(gel(F,1),gel(F0,1))
     899           0 :        && equalii(gel(F,2),gel(F0,2))) break;
     900             :     }
     901          21 :     set_avma(av);
     902             :   }
     903           7 :   return 1;
     904             : }
     905             : 
     906             : static GEN
     907       51814 : crabs(GEN a)
     908             : {
     909       51814 :   return signe(real_i(a)) < 0 ? gneg(a): a;
     910             : }
     911             : 
     912             : static GEN
     913       32025 : gcdreal(GEN a,GEN b)
     914             : {
     915       32025 :   if (!signe(real_i(a))) return crabs(b);
     916       31363 :   if (!signe(real_i(b))) return crabs(a);
     917       31215 :   if (expo(real_i(a))<-5) return crabs(b);
     918       11305 :   if (expo(real_i(b))<-5) return crabs(a);
     919        8456 :   a = crabs(a); b = crabs(b);
     920       18680 :   while (expo(real_i(b)) >= -5  && signe(real_i(b)))
     921             :   {
     922             :     long e;
     923       10224 :     GEN r, q = gcvtoi(divrr(real_i(a),real_i(b)),&e);
     924       10224 :     if (e > 0) return NULL;
     925       10224 :     r = gsub(a, gmul(q,b)); a = b; b = r;
     926             :   }
     927        8456 :   return crabs(a);
     928             : }
     929             : 
     930             : static int
     931       90951 : get_R(struct buch_quad *B, GEN C, long sreg, GEN z, GEN *ptR)
     932             : {
     933       90951 :   GEN R = gen_1;
     934             :   double c;
     935             :   long i;
     936       90951 :   if (B->PRECREG)
     937             :   {
     938        2877 :     R = crabs(gel(C,1));
     939       34902 :     for (i=2; i<=sreg; i++)
     940             :     {
     941       32025 :       R = gcdreal(gel(C,i), R);
     942       32025 :       if (!R) return fupb_PRECI;
     943             :     }
     944        2877 :     if (gexpo(real_i(R)) <= -3)
     945             :     {
     946           0 :       if (DEBUGLEVEL>2) err_printf("regulator is zero.\n");
     947           0 :       return fupb_RELAT;
     948             :     }
     949        2877 :     if (DEBUGLEVEL>2) err_printf("#### Tentative regulator: %Ps\n",R);
     950             :   }
     951       90951 :   c = gtodouble(gmul(z, real_i(R)));
     952       90949 :   if (c < 0.8 || c > 1.3) return fupb_RELAT;
     953       70256 :   *ptR = R; return fupb_NONE;
     954             : }
     955             : 
     956             : static int
     957       70256 : quad_be_honest(struct buch_quad *B)
     958             : {
     959             :   int r;
     960       70256 :   if (B->KC2 <= B->KC) return 1;
     961          14 :   if (DEBUGLEVEL>2)
     962           0 :     err_printf("be honest for primes from %ld to %ld\n", B->FB[B->KC+1],B->FB[B->KC2]);
     963          14 :   r = B->PRECREG? real_be_honest(B): imag_be_honest(B);
     964          14 :   if (DEBUGLEVEL>2) err_printf("\n");
     965          14 :   return r;
     966             : }
     967             : 
     968             : static GEN
     969       70429 : Buchquad_i(GEN D, double cbach, double cbach2, long prec)
     970             : {
     971       70429 :   const long MAXRELSUP = 7, SFB_MAX = 3;
     972             :   pari_timer T;
     973             :   pari_sp av, av2;
     974       70429 :   const long RELSUP = 5;
     975             :   long i, s, current, triv, sfb_trials, nrelsup, nreldep, need, nsubFB, minSFB;
     976             :   ulong low, high, LIMC0, LIMC, LIMC2, LIMCMAX, cp;
     977       70429 :   GEN W, cyc, gen, dep, mat, C, extraC, B, R, invhr, h = NULL; /*-Wall*/
     978             :   double drc, sdrc, lim, LOGD, LOGD2;
     979             :   GRHcheck_t GRHcheck;
     980             :   struct qfr_data q;
     981             :   struct buch_quad BQ;
     982       70429 :   int FIRST = 1;
     983             : 
     984       70429 :   check_quaddisc(D, &s, /*junk*/&i, "Buchquad");
     985       70432 :   R = NULL; /* -Wall */
     986       70432 :   BQ.q = &q;
     987       70432 :   q.D = D;
     988       70432 :   if (s < 0)
     989             :   {
     990       68276 :     if (abscmpiu(q.D,4) <= 0)
     991         176 :       retmkvec4(gen_1, cgetg(1,t_VEC), cgetg(1,t_VEC), gen_1);
     992       68099 :     prec = BQ.PRECREG = 0;
     993             :   } else {
     994        2156 :     BQ.PRECREG = maxss(prec+EXTRAPREC64, nbits2prec(2*expi(q.D) + 128));
     995             :   }
     996       70255 :   if (DEBUGLEVEL>2) timer_start(&T);
     997       70255 :   BQ.primfact   = new_chunk(100);
     998       70255 :   BQ.exprimfact = new_chunk(100);
     999       70256 :   BQ.hashtab = (long**) new_chunk(HASHT);
    1000    72004425 :   for (i=0; i<HASHT; i++) BQ.hashtab[i] = NULL;
    1001             : 
    1002       70258 :   drc = fabs(gtodouble(q.D));
    1003       70256 :   LOGD = log(drc);
    1004       70256 :   LOGD2 = LOGD * LOGD;
    1005             : 
    1006       70256 :   sdrc = lim = sqrt(drc);
    1007       70256 :   if (!BQ.PRECREG) lim /= sqrt(3.);
    1008       70256 :   cp = (ulong)exp(sqrt(LOGD*log(LOGD)/8.0));
    1009       70256 :   if (cp < 20) cp = 20;
    1010       70256 :   if (cbach > 6.) {
    1011           0 :     if (cbach2 < cbach) cbach2 = cbach;
    1012           0 :     cbach = 6.;
    1013             :   }
    1014       70256 :   if (cbach < 0.)
    1015           0 :     pari_err_DOMAIN("Buchquad","Bach constant","<",gen_0,dbltor(cbach));
    1016       70256 :   av = avma;
    1017       70256 :   BQ.powsubFB = BQ.subFB = NULL;
    1018       70256 :   minSFB = (expi(D) > 15)? 3: 2;
    1019       70257 :   init_GRHcheck(&GRHcheck, 2, BQ.PRECREG? 2: 0, LOGD);
    1020       70256 :   high = low = LIMC0 = maxss((long)(cbach2*LOGD2), 1);
    1021       70255 :   LIMCMAX = (long)(4.*LOGD2);
    1022             :   /* 97/1223 below to ensure a good enough approximation of residue */
    1023       70255 :   cache_prime_quad(&GRHcheck, expi(D) < 16 ? 97: 1223, D);
    1024      586625 :   while (!quadGRHchk(D, &GRHcheck, high))
    1025             :   {
    1026      516369 :     low = high;
    1027      516369 :     high *= 2;
    1028             :   }
    1029      516446 :   while (high - low > 1)
    1030             :   {
    1031      446195 :     long test = (low+high)/2;
    1032      446195 :     if (quadGRHchk(D, &GRHcheck, test))
    1033      233707 :       high = test;
    1034             :     else
    1035      212486 :       low = test;
    1036             :   }
    1037       70251 :   if (high == LIMC0+1 && quadGRHchk(D, &GRHcheck, LIMC0))
    1038           0 :     LIMC2 = LIMC0;
    1039             :   else
    1040       70251 :     LIMC2 = high;
    1041       70251 :   if (LIMC2 > LIMCMAX) LIMC2 = LIMCMAX;
    1042       70251 :   LIMC0 = (long)(cbach*LOGD2);
    1043       70251 :   LIMC = cbach ? LIMC0 : LIMC2;
    1044       70251 :   LIMC = maxss(LIMC, nthidealquad(D, 2));
    1045             : 
    1046             : /* LIMC = Max(cbach*(log D)^2, exp(sqrt(log D loglog D) / 8)) */
    1047       70256 : START:
    1048             :   do
    1049             :   {
    1050       70389 :     if (!FIRST) LIMC = bnf_increase_LIMC(LIMC,LIMCMAX);
    1051       70389 :     if (DEBUGLEVEL>2 && LIMC > LIMC0)
    1052           0 :       err_printf("%s*** Bach constant: %f\n", FIRST?"":"\n", LIMC/LOGD2);
    1053       70389 :     FIRST = 0; set_avma(av);
    1054       70388 :     guncloneNULL(BQ.subFB);
    1055       70388 :     guncloneNULL(BQ.powsubFB);
    1056       70388 :     clearhash(BQ.hashtab);
    1057       70391 :     if (LIMC < cp) LIMC = cp;
    1058       70391 :     if (LIMC2 < LIMC) LIMC2 = LIMC;
    1059       70391 :     if (BQ.PRECREG) qfr_data_init(q.D, BQ.PRECREG, &q);
    1060             : 
    1061       70391 :     FBquad(&BQ, LIMC2, LIMC, &GRHcheck);
    1062       70390 :     if (DEBUGLEVEL>2) timer_printf(&T, "factor base");
    1063       70390 :     BQ.subFB = subFBquad(&BQ, q.D, lim + 0.5, minSFB);
    1064       70390 :     if (DEBUGLEVEL>2) timer_printf(&T, "subFBquad = %Ps",
    1065             :                                  vecpermute(BQ.FB, BQ.subFB));
    1066       70390 :     nsubFB = lg(BQ.subFB) - 1;
    1067             :   }
    1068       70390 :   while (nsubFB < (expi(D) > 15 ? 3 : 2));
    1069             :   /* invhr = 2^r1 (2pi)^r2 / sqrt(D) w ~ L(chi,1) / hR */
    1070       70257 :   invhr = gmul(dbltor((BQ.PRECREG?2.:M_PI)/sdrc),
    1071             :                compute_invresquad(&GRHcheck, LIMC));
    1072       70257 :   BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
    1073       70257 :   if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
    1074       70256 :   BQ.limhash = (LIMC & HIGHMASK)? (HIGHBIT>>1): LIMC*LIMC;
    1075             : 
    1076       70256 :   need = BQ.KC + RELSUP - 2;
    1077       70256 :   current = 0;
    1078       70256 :   W = NULL;
    1079       70256 :   sfb_trials = nreldep = nrelsup = 0;
    1080       70256 :   s = nsubFB + RELSUP;
    1081       70256 :   av2 = avma;
    1082             : 
    1083             :   do
    1084             :   {
    1085      103833 :     if ((nreldep & 3) == 1 || (nrelsup & 7) == 1) {
    1086       28768 :       if (DEBUGLEVEL>2) err_printf("*** Changing sub factor base\n");
    1087       28768 :       gunclone(BQ.subFB);
    1088       28768 :       gunclone(BQ.powsubFB);
    1089       28768 :       BQ.subFB = gclone(vecslice(BQ.vperm, 1, nsubFB));
    1090       28768 :       BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
    1091       28768 :       if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
    1092       28768 :       clearhash(BQ.hashtab);
    1093             :     }
    1094      103833 :     need += 2;
    1095      103833 :     mat    = cgetg(need+1, t_MAT);
    1096      103832 :     extraC = cgetg(need+1, t_VEC);
    1097      103835 :     if (!W) { /* first time */
    1098       70258 :       C = extraC;
    1099       70258 :       triv = trivial_relations(&BQ, mat, C);
    1100       70258 :       if (DEBUGLEVEL>2) err_printf("KC = %ld, need %ld relations\n", BQ.KC, need);
    1101             :     } else {
    1102       33577 :       triv = 0;
    1103       33577 :       if (DEBUGLEVEL>2) err_printf("...need %ld more relations\n", need);
    1104             :     }
    1105      103834 :     if (BQ.PRECREG) {
    1106      189378 :       for (i = triv+1; i<=need; i++) {
    1107      184254 :         gel(mat,i) = zero_zv(BQ.KC);
    1108      184254 :         gel(extraC,i) = mkcomplex(cgetr(BQ.PRECREG), cgeti(3));
    1109             :       }
    1110        5124 :       real_relations(&BQ, need - triv, &current, s,LIMC,mat + triv,extraC + triv);
    1111             :     } else {
    1112     1702476 :       for (i = triv+1; i<=need; i++) {
    1113     1603770 :         gel(mat,i) = zero_zv(BQ.KC);
    1114     1603766 :         gel(extraC,i) = gen_0;
    1115             :       }
    1116       98706 :       imag_relations(&BQ, need - triv, &current, LIMC,mat + triv);
    1117             :     }
    1118      103833 :     if (!W)
    1119       70256 :       W = hnfspec_i(mat,BQ.vperm,&dep,&B,&C,nsubFB);
    1120             :     else
    1121       33577 :       W = hnfadd_i(W,BQ.vperm,&dep,&B,&C, mat,extraC);
    1122      103833 :     gerepileall(av2, 4, &W,&C,&B,&dep);
    1123      103835 :     need = BQ.KC - (lg(W)-1) - (lg(B)-1);
    1124      103835 :     if (need)
    1125             :     {
    1126       12884 :       if (++nreldep > 15 && cbach < 1) goto START;
    1127       12884 :       continue;
    1128             :     }
    1129             : 
    1130       90951 :     h = ZM_det_triangular(W);
    1131       90950 :     if (DEBUGLEVEL>2) err_printf("\n#### Tentative class number: %Ps\n", h);
    1132             : 
    1133       90950 :     switch(get_R(&BQ, C, (lg(C)-1) - (lg(B)-1) - (lg(W)-1), mulir(h,invhr), &R))
    1134             :     {
    1135           0 :       case fupb_PRECI:
    1136           0 :         BQ.PRECREG = precdbl(BQ.PRECREG);
    1137           0 :         FIRST = 1; goto START;
    1138             : 
    1139       20693 :       case fupb_RELAT:
    1140       20693 :         if (++nrelsup > MAXRELSUP)
    1141             :         {
    1142          35 :           if (++sfb_trials > SFB_MAX && cbach <= 1) goto START;
    1143          35 :           if (nsubFB < minss(10,BQ.KC)) nsubFB++;
    1144             :         }
    1145       20693 :         need = minss(BQ.KC, nrelsup);
    1146             :     }
    1147             :   }
    1148      103833 :   while (need);
    1149             :   /* DONE */
    1150       70256 :   if (!quad_be_honest(&BQ)) goto START;
    1151       70256 :   if (DEBUGLEVEL>2) timer_printf(&T, "be honest");
    1152       70256 :   clearhash(BQ.hashtab);
    1153       70258 :   free_GRHcheck(&GRHcheck);
    1154             : 
    1155       70258 :   gen = get_clgp(&BQ,W,&cyc);
    1156       70256 :   gunclone(BQ.subFB);
    1157       70258 :   gunclone(BQ.powsubFB);
    1158       70258 :   if (BQ.PRECREG)
    1159        2156 :     return mkvec5(h, cyc, gen, real_i(R), mpodd(imag_i(R)) ? gen_m1:gen_1);
    1160             :   else
    1161       68102 :     return mkvec4(h, cyc, gen, real_i(R));
    1162             : }
    1163             : GEN
    1164        4396 : Buchquad(GEN D, double c, double c2, long prec)
    1165             : {
    1166        4396 :   pari_sp av = avma;
    1167        4396 :   GEN z = Buchquad_i(D, c, c2, prec);
    1168        4396 :   return gerepilecopy(av, z);
    1169             : }
    1170             : 
    1171             : GEN
    1172           0 : buchimag(GEN D, GEN c, GEN c2, GEN REL)
    1173           0 : { (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),0); }
    1174             : 
    1175             : GEN
    1176           0 : buchreal(GEN D, GEN flag, GEN c, GEN c2, GEN REL, long prec) {
    1177           0 :   if (signe(flag)) pari_err_IMPL("narrow class group");
    1178           0 :   (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),prec);
    1179             : }
    1180             : 
    1181             : GEN
    1182        4396 : quadclassunit0(GEN x, long flag, GEN data, long prec)
    1183             : {
    1184             :   long lx;
    1185        4396 :   double c1 = 0.0, c2 = 0.0;
    1186             : 
    1187        4396 :   if (!data) lx=1;
    1188             :   else
    1189             :   {
    1190          28 :     lx = lg(data);
    1191          28 :     if (typ(data)!=t_VEC) pari_err_TYPE("quadclassunit", data);
    1192          28 :     if (lx > 7) pari_err_DIM("quadclassunit [tech vector]");
    1193          28 :     if (lx > 3) lx = 3;
    1194             :   }
    1195        4396 :   switch(lx)
    1196             :   {
    1197          21 :     case 3: c2 = gtodouble(gel(data,2));
    1198          28 :     case 2: c1 = gtodouble(gel(data,1));
    1199             :   }
    1200        4396 :   if (flag) pari_err_IMPL("narrow class group");
    1201        4396 :   return Buchquad(x,c1,c2,prec);
    1202             : }
    1203             : GEN
    1204       61164 : quadclassno(GEN D)
    1205             : {
    1206       61164 :   pari_sp av = avma;
    1207       61164 :   GEN h = abgrp_get_no(Buchquad_i(D, 0, 0, 0));
    1208       61168 :   return icopy_avma(h, av);
    1209             : }
    1210             : long
    1211        4868 : quadclassnos(long D)
    1212             : {
    1213        4868 :   pari_sp av = avma;
    1214        4868 :   long h = itos(abgrp_get_no(Buchquad_i(stoi(D), 0, 0, 0)));
    1215        4868 :   return gc_long(av, h);
    1216             : }

Generated by: LCOV version 1.14