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.2 lcov report (development 29115-f22e516b23) Lines: 652 684 95.3 %
Date: 2024-05-18 08:06:58 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     2566796 : 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     1166872 : qfr3_canon(GEN x, struct qfr_data *S)
      59             : {
      60     1166872 :   GEN a = gel(x,1), c = gel(x,3);
      61     1166872 :   if (signe(a) < 0) {
      62      403277 :     if (absequalii(a,c)) return qfr3_rho(x, S);
      63      403270 :     setsigne(a, 1);
      64      403270 :     setsigne(c,-1);
      65             :   }
      66     1166865 :   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     1952657 : qfr5_canon(GEN x, struct qfr_data *S)
      81             : {
      82     1952657 :   GEN a = gel(x,1), c = gel(x,3);
      83     1952657 :   if (signe(a) < 0) {
      84      723247 :     if (absequalii(a,c)) return qfr5_rho(x, S);
      85      723240 :     setsigne(a, 1);
      86      723240 :     setsigne(c,-1);
      87             :   }
      88     1952650 :   return x;
      89             : }
      90             : static GEN
      91     1735097 : QFR5_comp(GEN x,GEN y, struct qfr_data *S)
      92     1735097 : { return qfr5_canon(qfr5_comp(x,y,S), S); }
      93             : static GEN
      94     1007482 : QFR3_comp(GEN x, GEN y, struct qfr_data *S)
      95     1007482 : { return qfr3_canon(qfr3_comp(x,y,S), S); }
      96             : 
      97             : /* compute rho^n(x) */
      98             : static GEN
      99      220129 : qfr5_rho_pow(GEN x, long n, struct qfr_data *S)
     100             : {
     101             :   long i;
     102      220129 :   pari_sp av = avma;
     103     2194332 :   for (i=1; i<=n; i++)
     104             :   {
     105     1974203 :     x = qfr5_rho(x,S);
     106     1974203 :     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      220129 :   return gerepilecopy(av, x);
     113             : }
     114             : 
     115             : static GEN
     116      217560 : qfr5_pf(struct qfr_data *S, long p, long prec)
     117             : {
     118      217560 :   GEN y = primeform_u(S->D,p);
     119      217560 :   return qfr5_canon(qfr5_red(qfr_to_qfr5(y,prec), S), S);
     120             : }
     121             : 
     122             : static GEN
     123      159362 : qfr3_pf(struct qfr_data *S, long p)
     124             : {
     125      159362 :   GEN y = primeform_u(S->D,p);
     126      159362 :   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     4051534 : init_form(struct buch_quad *B, GEN ex,
     134             :           GEN (*comp)(GEN,GEN,struct qfr_data *S))
     135             : {
     136     4051534 :   long i, l = lg(B->powsubFB);
     137     4051534 :   GEN F = NULL;
     138    22904626 :   for (i=1; i<l; i++)
     139    18861280 :     if (ex[i])
     140             :     {
     141    17684989 :       GEN t = gmael(B->powsubFB,i,ex[i]);
     142    17684989 :       F = F? comp(F, t, B->q): t;
     143             :     }
     144     4043346 :   return F;
     145             : }
     146             : static GEN
     147      196987 : qfr5_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFR5_comp); }
     148             : 
     149             : static GEN
     150    11753191 : QFI_comp(GEN x, GEN y, struct qfr_data *S) { (void)S; return qfbcomp_i(x,y); }
     151             : 
     152             : static GEN
     153      159104 : qfi_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFI_comp); }
     154             : 
     155             : static GEN
     156     3695010 : random_form(struct buch_quad *B, GEN ex,
     157             :             GEN (*comp)(GEN,GEN, struct qfr_data *S))
     158             : {
     159     3695010 :   long i, l = lg(ex);
     160     3695010 :   pari_sp av = avma;
     161             :   GEN F;
     162             :   for(;;)
     163             :   {
     164    20588070 :     for (i=1; i<l; i++) ex[i] = random_bits(RANDOM_BITS);
     165     3695328 :     if ((F = init_form(B, ex, comp))) return F;
     166         439 :     set_avma(av);
     167             :   }
     168             : }
     169             : static GEN
     170      161539 : qfr3_random(struct buch_quad *B,GEN ex){ return random_form(B, ex, &QFR3_comp); }
     171             : static GEN
     172     3533480 : 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        2351 : bnf_increase_LIMC(ulong LIMC, ulong D)
     186             : {
     187        2351 :   if (LIMC >= D) pari_err_BUG("Buchmann's algorithm");
     188        2351 :   if (LIMC <= D / 13.333)
     189         209 :     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        2351 :   if (LIMC > D) LIMC = D;
     193        2351 :   return LIMC;
     194             : }
     195             : 
     196             : /* Is |q| <= p ? */
     197             : static int
     198       10416 : isless_iu(GEN q, ulong p) {
     199       10416 :   long l = lgefint(q);
     200       10416 :   return l==2 || (l == 3 && uel(q,2) <= p);
     201             : }
     202             : 
     203             : static GEN
     204     5026097 : Z_isquasismooth_prod(GEN N, GEN P)
     205             : {
     206     5026097 :   P = gcdii(P,N);
     207    10028433 :   while (!is_pm1(P))
     208             :   {
     209     4997227 :     N = diviiexact(N, P);
     210     5006028 :     P = gcdii(N, P);
     211             :   }
     212     5020058 :   return N;
     213             : }
     214             : 
     215             : static long
     216     5026896 : factorquad(struct buch_quad *B, GEN f, long nFB, ulong limp)
     217             : {
     218             :   ulong X;
     219     5026896 :   long i, lo = 0;
     220     5026896 :   GEN F, x = gel(f,1), FB = B->FB, P = B->primfact, E = B->exprimfact;
     221     5026896 :   if (B->badprim && !is_pm1(gcdii(x, B->badprim))) return 0;
     222     5025720 :   F =  Z_isquasismooth_prod(x, B->prodFB);
     223     5019949 :   if (cmpiu(F, B->limhash) > 0) return 0;
     224     4395199 :   for (i=1; lgefint(x) > 3; i++)
     225             :   {
     226       10416 :     ulong p = uel(FB,i), r;
     227       10416 :     GEN q = absdiviu_rem(x, p, &r);
     228       10416 :     if (!r)
     229             :     {
     230        1530 :       long k = 0;
     231        2633 :       do { k++; x = q; q = absdiviu_rem(x, p, &r); } while (!r);
     232        1530 :       lo++; P[lo] = p; E[lo] = k;
     233             :     }
     234       10416 :     if (isless_iu(q,p)) {
     235           1 :       if (lgefint(x) == 3) { X = uel(x,2); goto END; }
     236          34 :       return 0;
     237             :     }
     238       10415 :     if (i == nFB) return 0;
     239             :   }
     240     4384783 :   X = uel(x,2);
     241     4384783 :   if (X == 1) { P[0] = 0; return 1; }
     242    68716060 :   for (;; i++)
     243    68716060 :   { /* single precision affair, split for efficiency */
     244    73082283 :     ulong p = uel(FB,i);
     245    73082283 :     ulong q = X / p, r = X % p; /* gcc makes a single div */
     246    73082283 :     if (!r)
     247             :     {
     248     5631502 :       long k = 0;
     249     6798929 :       do { k++; X = q; q = X / p; r = X % p; } while (!r);
     250     5631502 :       lo++; P[lo] = p; E[lo] = k;
     251             :     }
     252    73082283 :     if (q <= p) break;
     253    68747664 :     if (i == nFB) return 0;
     254             :   }
     255     4334620 : END:
     256     4334620 :   if (X > B->limhash) return 0;
     257     4334620 :   if (X != 1 && X <= limp) { lo++; P[lo] = X; E[lo] = 1; X = 1; }
     258     4334620 :   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     2566779 : largeprime(struct buch_quad *B, long q, GEN ex, long np, long nrho)
     264             : {
     265     2566779 :   const long hashv = hash(q);
     266     2566794 :   long *pt, i, l = lg(B->subFB);
     267             : 
     268     2789700 :   for (pt = B->hashtab[hashv]; ; pt = (long*) pt[0])
     269             :   {
     270     2789700 :     if (!pt)
     271             :     {
     272     2368926 :       pt = (long*) pari_malloc((l+3) * sizeof(long));
     273     2369082 :       *pt++ = nrho; /* nrho = pt[-3] */
     274     2369082 :       *pt++ = np;   /* np   = pt[-2] */
     275     2369082 :       *pt++ = q;    /* q    = pt[-1] */
     276     2369082 :       pt[0] = (long)B->hashtab[hashv];
     277    15025019 :       for (i=1; i<l; i++) pt[i]=ex[i];
     278     2369082 :       B->hashtab[hashv]=pt; return NULL;
     279             :     }
     280      420774 :     if (pt[-1] == q) break;
     281             :   }
     282      236297 :   for(i=1; i<l; i++)
     283      233001 :     if (pt[i] != ex[i]) return pt;
     284        3296 :   return (pt[-2]==np)? NULL: pt;
     285             : }
     286             : 
     287             : static void
     288      169449 : clearhash(long **hash)
     289             : {
     290             :   long *pt;
     291             :   long i;
     292   173646800 :   for (i=0; i<HASHT; i++) {
     293   175846455 :     for (pt = hash[i]; pt; ) {
     294     2369104 :       void *z = (void*)(pt - 3);
     295     2369104 :       pt = (long*) pt[0]; pari_free(z);
     296             :     }
     297   173477351 :     hash[i] = NULL;
     298             :   }
     299      170278 : }
     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      399971 : GRH_ensure(GRHcheck_t *S, long nb)
     307             : {
     308      399971 :   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      399971 : }
     314             : /* cache data for all primes up to the LIM */
     315             : static void
     316     1173458 : cache_prime_quad(GRHcheck_t *S, ulong LIM, GEN D)
     317             : {
     318             :   GRHprime_t *pr;
     319             :   long nb;
     320             : 
     321     1173458 :   if (S->limp >= LIM) return;
     322       72332 :   nb = (long)primepi_upper_bound((double)LIM); /* #{p <= LIM} <= nb */
     323       72337 :   GRH_ensure(S, nb+1); /* room for one extra prime */
     324       72336 :   for (pr = S->primes + S->nprimes;;)
     325    12648096 :   {
     326    12720432 :     ulong p = u_forprime_next(&(S->P));
     327    12720062 :     pr->p = p;
     328    12720062 :     pr->logp = log((double)p);
     329    12720062 :     pr->dec = (GEN)kroiu(D,p);
     330    12720435 :     S->nprimes++;
     331    12720435 :     pr++;
     332             :     /* store up to nextprime(LIM) included */
     333    12720435 :     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       70256 :   double limp = log((double)LIMC) / 2;
     343             :   GRHprime_t *pr;
     344             :   long i;
     345    12672749 :   for (pr = S->primes, i = S->nprimes; i > 0; pr++, i--)
     346             :   {
     347    12609371 :     long s = (long)pr->dec;
     348    12609371 :     if (s)
     349             :     {
     350    12494433 :       ulong p = pr->p;
     351    12494433 :       if (s > 0 || pr->logp <= limp)
     352             :         /* Both p and P contribute */
     353     6346275 :         invres = mulur(p - s, divru(invres, p));
     354     6148158 :       else if (s<0)
     355             :         /* Only p contributes */
     356     6148164 :         invres = mulur(p, divru(invres, p - 1));
     357             :     }
     358             :   }
     359       63378 :   return gerepileuptoleaf(av, invres);
     360             : }
     361             : 
     362             : /* p | conductor of order of disc D ? */
     363             : static int
     364      391514 : is_bad(GEN D, ulong p)
     365             : {
     366      391514 :   pari_sp av = avma;
     367      391514 :   if (p == 2)
     368             :   {
     369       89558 :     long r = mod16(D) >> 1;
     370       89558 :     if (r && signe(D) < 0) r = 8-r;
     371       89558 :     return (r < 4);
     372             :   }
     373      301956 :   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       70257 : nthidealquad(GEN D, long n)
     379             : {
     380       70257 :   pari_sp av = avma;
     381             :   forprime_t S;
     382             :   ulong p;
     383       70257 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
     384      357173 :   while ((p = u_forprime_next(&S)) && n > 0)
     385      286927 :     if (!is_bad(D, p) && kroiu(D, p) >= 0) n--;
     386       70252 :   return gc_long(av, p);
     387             : }
     388             : 
     389             : static int
     390     1032821 : quadGRHchk(GEN D, GRHcheck_t *S, ulong LIMC)
     391             : {
     392     1032821 :   double logC = log((double)LIMC), SA = 0, SB = 0;
     393             :   long i;
     394             : 
     395     1032821 :   cache_prime_quad(S, LIMC, D);
     396     1032816 :   for (i = 0;; i++)
     397    29334413 :   {
     398    30367229 :     GRHprime_t *pr = S->primes+i;
     399    30367229 :     ulong p = pr->p;
     400             :     long M;
     401             :     double logNP, q, A, B;
     402    30367229 :     if (p > LIMC) break;
     403    29423531 :     if ((long)pr->dec < 0)
     404             :     {
     405    14693485 :       logNP = 2 * pr->logp;
     406    14693485 :       q = 1/(double)p;
     407             :     }
     408             :     else
     409             :     {
     410    14730046 :       logNP = pr->logp;
     411    14730046 :       q = 1/sqrt((double)p);
     412             :     }
     413    29423531 :     A = logNP * q; B = logNP * A; M = (long)(logC/logNP);
     414    29423531 :     if (M > 1)
     415             :     {
     416     2497008 :       double inv1_q = 1 / (1-q);
     417     2497008 :       A *= (1 - pow(q, (double) M)) * inv1_q;
     418     2497008 :       B *= (1 - pow(q, (double) M)*(M+1 - M*q)) * inv1_q * inv1_q;
     419             :     }
     420    29423531 :     if ((long)pr->dec>0) { SA += 2*A;SB += 2*B; } else { SA += A; SB += B; }
     421    29423531 :     if (p == LIMC) break;
     422             :   }
     423     1032816 :   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       70391 :   pr = S->primes;
     437       70391 :   B->numFB = cgetg(C2+1, t_VECSMALL);
     438       70391 :   B->FB    = cgetg(C2+1, t_VECSMALL);
     439       70277 :   av = avma;
     440       70277 :   B->KC = 0; i = 0;
     441       70277 :   B->badprim = gen_1;
     442     2798589 :   for (;; pr++) /* p <= C2 */
     443     2798589 :   {
     444     2868866 :     ulong p = pr->p;
     445     2868866 :     if (!B->KC && p > C1) B->KC = i;
     446     2868866 :     if (p > C2) break;
     447     2805993 :     switch ((long)pr->dec)
     448             :     {
     449     1382887 :       case -1: break; /* inert */
     450      104594 :       case  0: /* ramified */
     451      104594 :         if (is_bad(D, p)) { B->badprim = muliu(B->badprim, p); break; }
     452             :         /* fall through */
     453             :       default:  /* split */
     454     1423191 :         i++; B->numFB[p] = i; B->FB[i] = p; break;
     455             :     }
     456     2806106 :     if (p == C2)
     457             :     {
     458        7517 :       if (!B->KC) B->KC = i;
     459        7517 :       break;
     460             :     }
     461             :   }
     462       70390 :   B->KC2 = i;
     463       70390 :   setlg(B->FB, B->KC2+1);
     464       70390 :   if (B->badprim != gen_1)
     465          21 :     B->badprim = gerepileuptoint(av, B->badprim);
     466             :   else
     467             :   {
     468       70369 :     B->badprim = NULL; set_avma(av);
     469             :   }
     470       70390 :   B->prodFB = zv_prod_Z(B->FB);
     471       70388 : }
     472             : 
     473             : /* create B->vperm, return B->subFB */
     474             : static GEN
     475       70389 : subFBquad(struct buch_quad *B, GEN D, double PROD, long minSFB)
     476             : {
     477       70389 :   long i, j, lgsub = 1, ino = 1, lv = B->KC+1;
     478       70389 :   double prod = 1.;
     479             :   pari_sp av;
     480             :   GEN no;
     481             : 
     482       70389 :   B->vperm = cgetg(lv, t_VECSMALL);
     483       70389 :   av = avma;
     484       70389 :   no    = cgetg(lv, t_VECSMALL);
     485      334513 :   for (j = 1; j < lv; j++)
     486             :   {
     487      334373 :     ulong p = uel(B->FB,j);
     488      334373 :     if (!umodiu(D, p)) no[ino++] = j; /* ramified */
     489             :     else
     490             :     {
     491      254581 :       B->vperm[lgsub++] = j;
     492      254581 :       prod *= p;
     493      254581 :       if (lgsub > minSFB && prod > PROD) break;
     494             :     }
     495             :   }
     496             :   /* lgsub >= 1 otherwise quadGRHchk is false */
     497       70389 :   i = lgsub;
     498      150180 :   for (j = 1; j < ino;i++,j++) B->vperm[i] = no[j];
     499     1158901 :   for (     ; i < lv; i++)     B->vperm[i] = i;
     500       70389 :   no = gclone(vecslice(B->vperm, 1, lgsub-1));
     501       70390 :   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       99063 : powsubFBquad(struct buch_quad *B, long n)
     507             : {
     508       99063 :   pari_sp av = avma;
     509       99063 :   long i,j, l = lg(B->subFB);
     510       99063 :   GEN F, y, x = cgetg(l, t_VEC), D = B->q->D;
     511             : 
     512       99064 :   if (B->PRECREG) /* real */
     513             :   {
     514       39627 :     for (i=1; i<l; i++)
     515             :     {
     516       34510 :       F = qfr5_pf(B->q, B->FB[B->subFB[i]], B->PRECREG);
     517       34510 :       y = cgetg(n+1, t_VEC); gel(x,i) = y;
     518       34510 :       gel(y,1) = F;
     519      552160 :       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      426254 :     for (i=1; i<l; i++)
     525             :     {
     526      332322 :       F = qfi_pf(D, B->FB[B->subFB[i]]);
     527      332282 :       y = cgetg(n+1, t_VEC); gel(x,i) = y;
     528      333008 :       gel(y,1) = F;
     529     5309537 :       for (j=2; j<=n; j++) gel(y,j) = qfbcomp_i(gel(y,j-1), F);
     530             :     }
     531             :   }
     532       99049 :   x = gclone(x); set_avma(av); return x;
     533             : }
     534             : 
     535             : static void
     536       97570 : sub_fact(struct buch_quad *B, GEN col, GEN F)
     537             : {
     538       97570 :   GEN b = gel(F,2);
     539             :   long i;
     540      207876 :   for (i=1; i<=B->primfact[0]; i++)
     541             :   {
     542      110306 :     ulong p = B->primfact[i], k = B->numFB[p];
     543      110306 :     long e = B->exprimfact[i];
     544      110306 :     if (umodiu(b, p<<1) > p) e = -e;
     545      110306 :     col[k] -= e;
     546             :   }
     547       97570 : }
     548             : static void
     549     1885453 : add_fact(struct buch_quad *B, GEN col, GEN F)
     550             : {
     551     1885453 :   GEN b = gel(F,2);
     552             :   long i;
     553     5917911 :   for (i=1; i<=B->primfact[0]; i++)
     554             :   {
     555     4032398 :     ulong p = B->primfact[i], k = B->numFB[p];
     556     4032398 :     long e = B->exprimfact[i];
     557     4032398 :     if (umodiu(b, p<<1) > p) e = -e;
     558     4032458 :     col[k] += e;
     559             :   }
     560     1885513 : }
     561             : 
     562             : static GEN
     563       70258 : get_clgp(struct buch_quad *B, GEN W, GEN *ptD)
     564             : {
     565       70258 :   GEN res, init, u1, D = ZM_snf_group(W,NULL,&u1);
     566       70257 :   long i, j, l = lg(W), c = lg(D);
     567             : 
     568       70257 :   res=cgetg(c,t_VEC); init = cgetg(l,t_VEC);
     569      215916 :   for (i=1; i<l; i++) gel(init,i) = primeform_u(B->q->D, B->FB[B->vperm[i]]);
     570      195554 :   for (j=1; j<c; j++)
     571             :   {
     572      125298 :     GEN g = NULL;
     573      125298 :     if (signe(B->q->D) > 0)
     574             :     {
     575       13328 :       for (i=1; i<l; i++)
     576             :       {
     577        9618 :         GEN t, u = gcoeff(u1,i,j);
     578        9618 :         if (!signe(u)) continue;
     579        4543 :         t = qfr3_pow(gel(init,i), u, B->q);
     580        4543 :         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      422708 :       for (i=1; i<l; i++)
     587             :       {
     588      301123 :         GEN t, u = gcoeff(u1,i,j);
     589      301123 :         if (!signe(u)) continue;
     590      203146 :         t = powgi(gel(init,i), u);
     591      203151 :         g = g? qfbcomp_i(g, t): t;
     592             :       }
     593             :     }
     594      125299 :     gel(res,j) = g;
     595             :   }
     596       70256 :   *ptD = D; return res;
     597             : }
     598             : 
     599             : static long
     600       70255 : trivial_relations(struct buch_quad *B, GEN mat, GEN C)
     601             : {
     602       70255 :   long i, j = 0;
     603       70255 :   GEN col, D = B->q->D;
     604     1492657 :   for (i = 1; i <= B->KC; i++)
     605             :   { /* ramified prime ==> trivial relation */
     606     1422385 :     if (umodiu(D, B->FB[i])) continue;
     607      104201 :     col = zero_zv(B->KC);
     608      104217 :     col[i] = 2; j++;
     609      104217 :     gel(mat,j) = col;
     610      104217 :     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       98743 : imag_relations(struct buch_quad *B, long need, long *pc, ulong LIMC, GEN mat)
     626             : {
     627             :   pari_timer T;
     628       98743 :   long lgsub = lg(B->subFB), current = *pc, nbtest = 0, s = 0;
     629             :   long i, fpc;
     630             :   pari_sp av;
     631       98743 :   GEN col, form, ex = cgetg(lgsub, t_VECSMALL);
     632             : 
     633       98743 :   if (!current) current = 1;
     634       98743 :   if (DEBUGLEVEL>2) timer_start(&T);
     635       98743 :   av = avma;
     636             :   for(;;)
     637             :   {
     638     3631854 :     if (s >= need) break;
     639     3533110 :     set_avma(av);
     640     3533016 :     form = qfi_random(B,ex);
     641     3531381 :     form = qfbcomp_i(form, qfi_pf(B->q->D, B->FB[current]));
     642     3530894 :     nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
     643     3532520 :     if (!fpc)
     644             :     {
     645      286481 :       if (DEBUGLEVEL>3) err_printf(".");
     646      286481 :       if ((nbtest & 0xff) == 0 && ++current > B->KC) current = 1;
     647      286481 :       continue;
     648             :     }
     649     3246039 :     if (fpc > 1)
     650             :     {
     651     1801670 :       long *fpd = largeprime(B,fpc,ex,current,0);
     652             :       ulong b1, b2, p;
     653             :       GEN form2;
     654     1802008 :       if (!fpd)
     655             :       {
     656     1642935 :         if (DEBUGLEVEL>3) err_printf(".");
     657     1642931 :         continue;
     658             :       }
     659      159073 :       form2 = qfbcomp_i(qfi_factorback(B,fpd), qfi_pf(B->q->D, B->FB[fpd[-2]]));
     660      159112 :       p = fpc << 1;
     661      159112 :       b1 = umodiu(gel(form2,2), p);
     662      159116 :       b2 = umodiu(gel(form,2),  p);
     663      159116 :       if (b1 != b2 && b1+b2 != p) continue;
     664             : 
     665      159116 :       col = gel(mat,++s);
     666      159116 :       add_fact(B,col, form);
     667      159115 :       (void)factorquad(B,form2,B->KC,LIMC);
     668      159115 :       if (b1==b2)
     669             :       {
     670      412274 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
     671       79944 :         sub_fact(B, col, form2); col[fpd[-2]]++;
     672             :       }
     673             :       else
     674             :       {
     675      407600 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
     676       79171 :         add_fact(B, col, form2); col[fpd[-2]]--;
     677             :       }
     678      159116 :       if (DEBUGLEVEL>2) err_printf(" %ldP",s);
     679             :     }
     680             :     else
     681             :     {
     682     1444369 :       col = gel(mat,++s);
     683     6941543 :       for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
     684     1444369 :       add_fact(B, col, form);
     685     1444729 :       if (DEBUGLEVEL>2) err_printf(" %ld",s);
     686             :     }
     687     1603699 :     col[current]--;
     688     1603699 :     if (++current > B->KC) current = 1;
     689             :   }
     690       98744 :   if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
     691       98744 :   *pc = current;
     692       98744 : }
     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      184660 : dist(GEN e, GEN d, long prec)
     716             : {
     717      184660 :   GEN t = qfr5_dist(e, d, prec);
     718      184660 :   return signe(d) < 0 ? mkcomplex(t, gen_1): t;
     719             : }
     720             : 
     721             : /* Real Quadratic fields */
     722             : 
     723             : static void
     724        5138 : real_relations(struct buch_quad *B, long need, long *pc, long lim, ulong LIMC, GEN mat, GEN C)
     725             : {
     726             :   pari_timer T;
     727        5138 :   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        5138 :   int first = (current == 0);
     732             :   pari_sp av, av1;
     733        5138 :   GEN d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL);
     734             : 
     735        5138 :   if (DEBUGLEVEL>2) timer_start(&T);
     736        5138 :   if (!current) current = 1;
     737        5138 :   if (lim > need) lim = need;
     738        5138 :   av = avma;
     739             :   for(;;)
     740             :   {
     741      166656 :     if (s >= need) break;
     742      161518 :     if (first && s >= lim) {
     743        2142 :       first = 0;
     744        2142 :       if (DEBUGLEVEL>2) dbg_all(&T, "initial", s, nbtest);
     745             :     }
     746      161518 :     set_avma(av); form = qfr3_random(B, ex);
     747      161518 :     if (!first)
     748      159341 :       form = QFR3_comp(form, qfr3_pf(B->q, B->FB[current]), B->q);
     749      161518 :     av1 = avma;
     750      161518 :     form0 = form; form1 = NULL;
     751      161518 :     endcycle = rhoacc = 0;
     752      161518 :     rho = -1;
     753             : 
     754     1300936 : CYCLE:
     755     1300936 :     if (endcycle || rho > 5000)
     756             :     {
     757          21 :       if (++current > B->KC) current = 1;
     758          21 :       continue;
     759             :     }
     760     1300915 :     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     1300915 :     if (rho < 0) rho = 0; /* first time in */
     766             :     else
     767             :     {
     768     1139397 :       form = qfr3_rho(form, B->q); rho++;
     769     1139397 :       rhoacc++;
     770     1139397 :       if (first)
     771      442526 :         endcycle = (absequalii(gel(form,1),gel(form0,1))
     772      221263 :              && equalii(gel(form,2),gel(form0,2)));
     773             :       else
     774             :       {
     775      918134 :         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      918134 :           { setsigne(form[1],1); setsigne(form[3],-1); }
     784      918176 :         if (equalii(gel(form,1),gel(form0,1)) &&
     785          42 :             equalii(gel(form,2),gel(form0,2))) continue;
     786             :       }
     787             :     }
     788     1300915 :     nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
     789     1300915 :     if (!fpc)
     790             :     {
     791      386778 :       if (DEBUGLEVEL>3) err_printf(".");
     792      386778 :       goto CYCLE;
     793             :     }
     794      914137 :     if (fpc > 1)
     795             :     { /* look for Large Prime relation */
     796      764946 :       long *fpd = largeprime(B,fpc,ex,first? 0: current,rhoacc);
     797             :       ulong b1, b2, p;
     798             :       GEN form2;
     799      764946 :       if (!fpd)
     800             :       {
     801      729477 :         if (DEBUGLEVEL>3) err_printf(".");
     802      729477 :         goto CYCLE;
     803             :       }
     804       35469 :       if (!form1)
     805             :       {
     806       35469 :         form1 = qfr5_factorback(B,ex);
     807       35469 :         if (!first)
     808       35469 :           form1 = QFR5_comp(form1, qfr5_pf(B->q, B->FB[current], prec), B->q);
     809             :       }
     810       35469 :       form1 = qfr5_rho_pow(form1, rho, B->q);
     811       35469 :       rho = 0;
     812             : 
     813       35469 :       form2 = qfr5_factorback(B,fpd);
     814       35469 :       if (fpd[-2])
     815       23709 :         form2 = QFR5_comp(form2, qfr5_pf(B->q, B->FB[fpd[-2]], prec), B->q);
     816       35469 :       form2 = qfr5_rho_pow(form2, fpd[-3], B->q);
     817       35469 :       if (!absequalii(gel(form2,1),gel(form2,3)))
     818             :       {
     819       35469 :         setsigne(form2[1], 1);
     820       35469 :         setsigne(form2[3],-1);
     821             :       }
     822       35469 :       p = fpc << 1;
     823       35469 :       b1 = umodiu(gel(form2,2), p);
     824       35469 :       b2 = umodiu(gel(form1,2), p);
     825       35469 :       if (b1 != b2 && b1+b2 != p) goto CYCLE;
     826             : 
     827       35469 :       col = gel(mat,++s);
     828       35469 :       add_fact(B, col, form1);
     829       35469 :       (void)factorquad(B,form2,B->KC,LIMC);
     830       35469 :       if (b1==b2)
     831             :       {
     832      135100 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
     833       17626 :         sub_fact(B,col, form2);
     834       17626 :         if (fpd[-2]) col[fpd[-2]]++;
     835       17626 :         d = dist(subii(gel(form1,4),gel(form2,4)),
     836       17626 :                       divrr(gel(form1,5),gel(form2,5)), prec);
     837             :       }
     838             :       else
     839             :       {
     840      136780 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
     841       17843 :         add_fact(B, col, form2);
     842       17843 :         if (fpd[-2]) col[fpd[-2]]--;
     843       17843 :         d = dist(addii(gel(form1,4),gel(form2,4)),
     844       17843 :                       mulrr(gel(form1,5),gel(form2,5)), prec);
     845             :       }
     846       35469 :       if (DEBUGLEVEL>2) err_printf(" %ldP",s);
     847             :     }
     848             :     else
     849             :     { /* standard relation */
     850      149191 :       if (!form1)
     851             :       {
     852      126049 :         form1 = qfr5_factorback(B, ex);
     853      126049 :         if (!first)
     854      123872 :           form1 = QFR5_comp(form1, qfr5_pf(B->q, B->FB[current], prec), B->q);
     855             :       }
     856      149191 :       form1 = qfr5_rho_pow(form1, rho, B->q);
     857      149191 :       rho = 0;
     858             : 
     859      149191 :       col = gel(mat,++s);
     860     1147377 :       for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
     861      149191 :       add_fact(B, col, form1);
     862      149191 :       d = dist(gel(form1,4), gel(form1,5), prec);
     863      149191 :       if (DEBUGLEVEL>2) err_printf(" %ld",s);
     864             :     }
     865      184660 :     gaffect(d, gel(C,s));
     866      184660 :     if (first)
     867             :     {
     868       25319 :       if (s >= lim) continue;
     869       23163 :       goto CYCLE;
     870             :     }
     871             :     else
     872             :     {
     873      159341 :       col[current]--;
     874      159341 :       if (++current > B->KC) current = 1;
     875             :     }
     876             :   }
     877        5138 :   if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
     878        5138 :   *pc = current;
     879        5138 : }
     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          49 :       fpc = factorquad(B,F,s,p-1);
     895          49 :       if (fpc == 1) { nbtest=0; s++; break; }
     896          28 :       if (++nbtest > 40) return 0;
     897          28 :       F = qfr3_canon(qfr3_rho(F, B->q), B->q);
     898          28 :       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       54922 : crabs(GEN a)
     908             : {
     909       54922 :   return signe(real_i(a)) < 0 ? gneg(a): a;
     910             : }
     911             : 
     912             : static GEN
     913       33901 : gcdreal(GEN a,GEN b)
     914             : {
     915       33901 :   if (!signe(real_i(a))) return crabs(b);
     916       32977 :   if (!signe(real_i(b))) return crabs(a);
     917       32862 :   if (expo(real_i(a))<-5) return crabs(b);
     918       11970 :   if (expo(real_i(b))<-5) return crabs(a);
     919        9023 :   a = crabs(a); b = crabs(b);
     920       19893 :   while (expo(real_i(b)) >= -5  && signe(real_i(b)))
     921             :   {
     922             :     long e;
     923       10870 :     GEN r, q = gcvtoi(divrr(real_i(a),real_i(b)),&e);
     924       10870 :     if (e > 0) return NULL;
     925       10870 :     r = gsub(a, gmul(q,b)); a = b; b = r;
     926             :   }
     927        9023 :   return crabs(a);
     928             : }
     929             : 
     930             : static int
     931       91149 : get_R(struct buch_quad *B, GEN C, long sreg, GEN z, GEN *ptR)
     932             : {
     933       91149 :   GEN R = gen_1;
     934             :   double c;
     935             :   long i;
     936       91149 :   if (B->PRECREG)
     937             :   {
     938        2975 :     R = crabs(gel(C,1));
     939       36876 :     for (i=2; i<=sreg; i++)
     940             :     {
     941       33901 :       R = gcdreal(gel(C,i), R);
     942       33901 :       if (!R) return fupb_PRECI;
     943             :     }
     944        2975 :     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        2975 :     if (DEBUGLEVEL>2) err_printf("#### Tentative regulator: %Ps\n",R);
     950             :   }
     951       91149 :   c = gtodouble(gmul(z, real_i(R)));
     952       91151 :   if (c < 0.8 || c > 1.3) return fupb_RELAT;
     953       70254 :   *ptR = R; return fupb_NONE;
     954             : }
     955             : 
     956             : static int
     957       70254 : quad_be_honest(struct buch_quad *B)
     958             : {
     959             :   int r;
     960       70254 :   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       70433 : Buchquad_i(GEN D, double cbach, double cbach2, long prec)
     970             : {
     971       70433 :   const long MAXRELSUP = 7, SFB_MAX = 3;
     972             :   pari_timer T;
     973             :   pari_sp av, av2;
     974       70433 :   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       70433 :   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       70433 :   int FIRST = 1;
     983             : 
     984       70433 :   check_quaddisc(D, &s, /*junk*/&i, "Buchquad");
     985       70433 :   R = NULL; /* -Wall */
     986       70433 :   BQ.q = &q;
     987       70433 :   q.D = D;
     988       70433 :   if (s < 0)
     989             :   {
     990       68277 :     if (abscmpiu(q.D,4) <= 0)
     991         175 :       retmkvec4(gen_1, cgetg(1,t_VEC), cgetg(1,t_VEC), gen_1);
     992       68102 :     prec = BQ.PRECREG = 0;
     993             :   } else {
     994        2156 :     BQ.PRECREG = maxss(prec+EXTRAPREC64, nbits2prec(2*expi(q.D) + 128));
     995             :   }
     996       70258 :   if (DEBUGLEVEL>2) timer_start(&T);
     997       70258 :   BQ.primfact   = new_chunk(100);
     998       70258 :   BQ.exprimfact = new_chunk(100);
     999       70258 :   BQ.hashtab = (long**) new_chunk(HASHT);
    1000    72010604 :   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       70256 :   LIMCMAX = (long)(4.*LOGD2);
    1022             :   /* 97/1223 below to ensure a good enough approximation of residue */
    1023       70256 :   cache_prime_quad(&GRHcheck, expi(D) < 16 ? 97: 1223, D);
    1024      586645 :   while (!quadGRHchk(D, &GRHcheck, high))
    1025             :   {
    1026      516387 :     low = high;
    1027      516387 :     high *= 2;
    1028             :   }
    1029      516452 :   while (high - low > 1)
    1030             :   {
    1031      446199 :     long test = (low+high)/2;
    1032      446199 :     if (quadGRHchk(D, &GRHcheck, test))
    1033      233709 :       high = test;
    1034             :     else
    1035      212490 :       low = test;
    1036             :   }
    1037       70253 :   if (high == LIMC0+1 && quadGRHchk(D, &GRHcheck, LIMC0))
    1038           0 :     LIMC2 = LIMC0;
    1039             :   else
    1040       70253 :     LIMC2 = high;
    1041       70253 :   if (LIMC2 > LIMCMAX) LIMC2 = LIMCMAX;
    1042       70253 :   LIMC0 = (long)(cbach*LOGD2);
    1043       70253 :   LIMC = cbach ? LIMC0 : LIMC2;
    1044       70253 :   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       70389 :     guncloneNULL(BQ.subFB);
    1055       70389 :     guncloneNULL(BQ.powsubFB);
    1056       70389 :     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       70389 :     if (DEBUGLEVEL>2) timer_printf(&T, "factor base");
    1063       70389 :     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       70255 :   BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
    1073       70256 :   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      103880 :     if ((nreldep & 3) == 1 || (nrelsup & 7) == 1) {
    1086       28809 :       if (DEBUGLEVEL>2) err_printf("*** Changing sub factor base\n");
    1087       28809 :       gunclone(BQ.subFB);
    1088       28809 :       gunclone(BQ.powsubFB);
    1089       28809 :       BQ.subFB = gclone(vecslice(BQ.vperm, 1, nsubFB));
    1090       28809 :       BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
    1091       28809 :       if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
    1092       28809 :       clearhash(BQ.hashtab);
    1093             :     }
    1094      103880 :     need += 2;
    1095      103880 :     mat    = cgetg(need+1, t_MAT);
    1096      103880 :     extraC = cgetg(need+1, t_VEC);
    1097      103879 :     if (!W) { /* first time */
    1098       70255 :       C = extraC;
    1099       70255 :       triv = trivial_relations(&BQ, mat, C);
    1100       70257 :       if (DEBUGLEVEL>2) err_printf("KC = %ld, need %ld relations\n", BQ.KC, need);
    1101             :     } else {
    1102       33624 :       triv = 0;
    1103       33624 :       if (DEBUGLEVEL>2) err_printf("...need %ld more relations\n", need);
    1104             :     }
    1105      103882 :     if (BQ.PRECREG) {
    1106      189798 :       for (i = triv+1; i<=need; i++) {
    1107      184660 :         gel(mat,i) = zero_zv(BQ.KC);
    1108      184660 :         gel(extraC,i) = mkcomplex(cgetr(BQ.PRECREG), cgeti(3));
    1109             :       }
    1110        5138 :       real_relations(&BQ, need - triv, &current, s,LIMC,mat + triv,extraC + triv);
    1111             :     } else {
    1112     1702753 :       for (i = triv+1; i<=need; i++) {
    1113     1604007 :         gel(mat,i) = zero_zv(BQ.KC);
    1114     1604009 :         gel(extraC,i) = gen_0;
    1115             :       }
    1116       98746 :       imag_relations(&BQ, need - triv, &current, LIMC,mat + triv);
    1117             :     }
    1118      103882 :     if (!W)
    1119       70258 :       W = hnfspec_i(mat,BQ.vperm,&dep,&B,&C,nsubFB);
    1120             :     else
    1121       33624 :       W = hnfadd_i(W,BQ.vperm,&dep,&B,&C, mat,extraC);
    1122      103882 :     gerepileall(av2, 4, &W,&C,&B,&dep);
    1123      103882 :     need = BQ.KC - (lg(W)-1) - (lg(B)-1);
    1124      103882 :     if (need)
    1125             :     {
    1126       12727 :       if (++nreldep > 15 && cbach < 1) goto START;
    1127       12727 :       continue;
    1128             :     }
    1129             : 
    1130       91155 :     h = ZM_det_triangular(W);
    1131       91150 :     if (DEBUGLEVEL>2) err_printf("\n#### Tentative class number: %Ps\n", h);
    1132             : 
    1133       91150 :     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       20897 :       case fupb_RELAT:
    1140       20897 :         if (++nrelsup > MAXRELSUP)
    1141             :         {
    1142          38 :           if (++sfb_trials > SFB_MAX && cbach <= 1) goto START;
    1143          38 :           if (nsubFB < minss(10,BQ.KC)) nsubFB++;
    1144             :         }
    1145       20897 :         need = minss(BQ.KC, nrelsup);
    1146             :     }
    1147             :   }
    1148      103878 :   while (need);
    1149             :   /* DONE */
    1150       70254 :   if (!quad_be_honest(&BQ)) goto START;
    1151       70254 :   if (DEBUGLEVEL>2) timer_printf(&T, "be honest");
    1152       70254 :   clearhash(BQ.hashtab);
    1153       70257 :   free_GRHcheck(&GRHcheck);
    1154             : 
    1155       70258 :   gen = get_clgp(&BQ,W,&cyc);
    1156       70257 :   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       61166 : quadclassno(GEN D)
    1205             : {
    1206       61166 :   pari_sp av = avma;
    1207       61166 :   GEN h = abgrp_get_no(Buchquad_i(D, 0, 0, 0));
    1208       61169 :   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