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.12.0 lcov report (development 23339-b1c33c51a) Lines: 634 669 94.8 %
Date: 2018-12-11 05:41:34 Functions: 44 47 93.6 %
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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : #include "pari.h"
      14             : #include "paripriv.h"
      15             : 
      16             : /*******************************************************************/
      17             : /*                                                                 */
      18             : /*         CLASS GROUP AND REGULATOR (McCURLEY, BUCHMANN)          */
      19             : /*                   QUADRATIC FIELDS                              */
      20             : /*                                                                 */
      21             : /*******************************************************************/
      22             : /* For largeprime() hashtable. Note that hashed pseudoprimes are odd (unless
      23             :  * 2 | index), hence the low order bit is not useful. So we hash
      24             :  * HASHBITS bits starting at bit 1, not bit 0 */
      25             : #define HASHBITS 10
      26             : static const long HASHT = 1L << HASHBITS;
      27             : 
      28             : static long
      29     1557367 : hash(long q) { return (q & ((1L << (HASHBITS+1)) - 1)) >> 1; }
      30             : #undef HASHBITS
      31             : 
      32             : /* See buch2.c:
      33             :  * B->subFB contains split p such that \prod p > sqrt(B->Disc)
      34             :  * B->powsubFB contains powers of forms in B->subFB */
      35             : #define RANDOM_BITS 4
      36             : static const long CBUCH = (1L<<RANDOM_BITS)-1;
      37             : 
      38             : struct buch_quad
      39             : {
      40             :   ulong limhash;
      41             :   long KC, KC2, PRECREG;
      42             :   long *primfact, *exprimfact, **hashtab;
      43             :   GEN FB, numFB;
      44             :   GEN powsubFB, vperm, subFB, badprim;
      45             :   struct qfr_data *QFR;
      46             : };
      47             : 
      48             : /*******************************************************************/
      49             : /*                                                                 */
      50             : /*  Routines related to binary quadratic forms (for internal use)  */
      51             : /*                                                                 */
      52             : /*******************************************************************/
      53             : /* output canonical representative wrt projection Cl^+ --> Cl (a > 0) */
      54             : static GEN
      55     1164331 : qfr3_canon(GEN x, struct qfr_data *S)
      56             : {
      57     1164331 :   GEN a = gel(x,1), c = gel(x,3);
      58     1164331 :   if (signe(a) < 0) {
      59      502173 :     if (absequalii(a,c)) return qfr3_rho(x, S);
      60      502166 :     setsigne(a, 1);
      61      502166 :     setsigne(c,-1);
      62             :   }
      63     1164324 :   return x;
      64             : }
      65             : static GEN
      66        3710 : qfr3_canon_safe(GEN x, struct qfr_data *S)
      67             : {
      68        3710 :   GEN a = gel(x,1), c = gel(x,3);
      69        3710 :   if (signe(a) < 0) {
      70         189 :     if (absequalii(a,c)) return qfr3_rho(x, S);
      71         189 :     gel(x,1) = negi(a);
      72         189 :     gel(x,3) = negi(c);
      73             :   }
      74        3710 :   return x;
      75             : }
      76             : static GEN
      77     1947281 : qfr5_canon(GEN x, struct qfr_data *S)
      78             : {
      79     1947281 :   GEN a = gel(x,1), c = gel(x,3);
      80     1947281 :   if (signe(a) < 0) {
      81      868427 :     if (absequalii(a,c)) return qfr5_rho(x, S);
      82      868420 :     setsigne(a, 1);
      83      868420 :     setsigne(c,-1);
      84             :   }
      85     1947274 :   return x;
      86             : }
      87             : static GEN
      88     1730337 : QFR5_comp(GEN x,GEN y, struct qfr_data *S)
      89     1730337 : { return qfr5_canon(qfr5_comp(x,y,S), S); }
      90             : static GEN
      91     1005354 : QFR3_comp(GEN x, GEN y, struct qfr_data *S)
      92     1005354 : { return qfr3_canon(qfr3_comp(x,y,S), S); }
      93             : 
      94             : /* compute rho^n(x) */
      95             : static GEN
      96      220584 : qfr5_rho_pow(GEN x, long n, struct qfr_data *S)
      97             : {
      98             :   long i;
      99      220584 :   pari_sp av = avma;
     100     2216151 :   for (i=1; i<=n; i++)
     101             :   {
     102     1995567 :     x = qfr5_rho(x,S);
     103     1995567 :     if (gc_needed(av,1))
     104             :     {
     105           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"qfr5_rho_pow");
     106           0 :       x = gerepilecopy(av, x);
     107             :     }
     108             :   }
     109      220584 :   return gerepilecopy(av, x);
     110             : }
     111             : 
     112             : static GEN
     113      216944 : qfr5_pf(struct qfr_data *S, long p, long prec)
     114             : {
     115      216944 :   GEN y = primeform_u(S->D,p);
     116      216944 :   return qfr5_canon(qfr5_red(qfr_to_qfr5(y,prec), S), S);
     117             : }
     118             : 
     119             : static GEN
     120      158956 : qfr3_pf(struct qfr_data *S, long p)
     121             : {
     122      158956 :   GEN y = primeform_u(S->D,p);
     123      158956 :   return qfr3_canon(qfr3_red(y, S), S);
     124             : }
     125             : 
     126             : #define qfi_pf primeform_u
     127             : 
     128             : /* Warning: ex[0] not set in general */
     129             : static GEN
     130     1804157 : init_form(struct buch_quad *B, GEN ex,
     131             :           GEN (*comp)(GEN,GEN,struct qfr_data *S))
     132             : {
     133     1804157 :   long i, l = lg(B->powsubFB);
     134     1804157 :   GEN F = NULL;
     135    12716951 :   for (i=1; i<l; i++)
     136    10912794 :     if (ex[i])
     137             :     {
     138    10231424 :       GEN t = gmael(B->powsubFB,i,ex[i]);
     139    10231424 :       F = F? comp(F, t, B->QFR): t;
     140             :     }
     141     1804157 :   return F;
     142             : }
     143             : static GEN
     144      197442 : qfr5_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFR5_comp); }
     145             : 
     146             : static GEN
     147     6544381 : QFI_comp(GEN x, GEN y, struct qfr_data *S) { (void)S; return qficomp(x,y); }
     148             : 
     149             : static GEN
     150       47614 : qfi_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFI_comp); }
     151             : 
     152             : static GEN
     153     1558966 : random_form(struct buch_quad *B, GEN ex,
     154             :             GEN (*comp)(GEN,GEN, struct qfr_data *S))
     155             : {
     156     1558966 :   long i, l = lg(ex);
     157     1558966 :   pari_sp av = avma;
     158             :   GEN F;
     159             :   for(;;)
     160             :   {
     161     1559236 :     for (i=1; i<l; i++) ex[i] = random_bits(RANDOM_BITS);
     162     3118067 :     if ((F = init_form(B, ex, comp))) return F;
     163         135 :     set_avma(av);
     164             :   }
     165             : }
     166             : static GEN
     167      161133 : qfr3_random(struct buch_quad *B,GEN ex){ return random_form(B, ex, &QFR3_comp); }
     168             : static GEN
     169     1397833 : qfi_random(struct buch_quad *B,GEN ex) { return random_form(B, ex, &QFI_comp); }
     170             : 
     171             : /*******************************************************************/
     172             : /*                                                                 */
     173             : /*                     Common subroutines                          */
     174             : /*                                                                 */
     175             : /*******************************************************************/
     176             : long
     177         350 : bnf_increase_LIMC(long LIMC, long LIMCMAX)
     178             : {
     179         350 :   if (LIMC >= LIMCMAX) pari_err_BUG("Buchmann's algorithm");
     180         350 :   if (LIMC <= LIMCMAX/40) /* cbach <= 0.3 */
     181          28 :     LIMC *= 2;
     182         322 :   else if (LIMCMAX < 60) /* \Delta_K <= 9 */
     183           0 :     LIMC++;
     184             :   else
     185         322 :     LIMC += LIMCMAX / 60; /* cbach += 0.2 */
     186         350 :   if (LIMC > LIMCMAX) LIMC = LIMCMAX;
     187         350 :   return LIMC;
     188             : }
     189             : 
     190             : /* Is |q| <= p ? */
     191             : static int
     192       73056 : isless_iu(GEN q, ulong p) {
     193       73056 :   long l = lgefint(q);
     194       73056 :   return l==2 || (l == 3 && uel(q,2) <= p);
     195             : }
     196             : 
     197             : static long
     198     2779332 : factorquad(struct buch_quad *B, GEN f, long nFB, ulong limp)
     199             : {
     200             :   ulong X;
     201     2779332 :   long i, lo = 0;
     202     2779332 :   GEN x = gel(f,1), FB = B->FB, P = B->primfact, E = B->exprimfact;
     203             : 
     204     5703992 :   for (i=1; lgefint(x) > 3; i++)
     205             :   {
     206       73056 :     ulong p = uel(FB,i), r;
     207       73056 :     GEN q = absdiviu_rem(x, p, &r);
     208       73056 :     if (!r)
     209             :     {
     210        3822 :       long k = 0;
     211        6154 :       do { k++; x = q; q = absdiviu_rem(x, p, &r); } while (!r);
     212        3822 :       lo++; P[lo] = p; E[lo] = k;
     213             :     }
     214       73056 :     if (isless_iu(q,p)) {
     215           1 :       if (lgefint(x) == 3) { X = uel(x,2); goto END; }
     216         391 :       return 0;
     217             :     }
     218       73055 :     if (i == nFB) return 0;
     219             :   }
     220     2778940 :   X = uel(x,2);
     221     2778940 :   if (X == 1) { P[0] = 0; return 1; }
     222   116522154 :   for (;; i++)
     223   116522154 :   { /* single precision affair, split for efficiency */
     224   119294129 :     ulong p = uel(FB,i);
     225   119294129 :     ulong q = X / p, r = X % p; /* gcc makes a single div */
     226   119294129 :     if (!r)
     227             :     {
     228     4382031 :       long k = 0;
     229     5324715 :       do { k++; X = q; q = X / p; r = X % p; } while (!r);
     230     4382031 :       lo++; P[lo] = p; E[lo] = k;
     231             :     }
     232   119294129 :     if (q <= p) break;
     233   117194022 :     if (i == nFB) return 0;
     234             :   }
     235             : END:
     236     2100108 :   if (X > B->limhash) return 0;
     237     2100080 :   if (X != 1 && X <= limp)
     238             :   {
     239      419477 :     if (B->badprim && ugcdui(X, B->badprim) > 1) return 0;
     240      418726 :     lo++; P[lo] = X; E[lo] = 1;
     241      418726 :     X = 1;
     242             :   }
     243     2099329 :   P[0] = lo; return X;
     244             : }
     245             : 
     246             : /* Check for a "large prime relation" involving q; q may not be prime */
     247             : static long *
     248     1557367 : largeprime(struct buch_quad *B, long q, GEN ex, long np, long nrho)
     249             : {
     250     1557367 :   const long hashv = hash(q);
     251     1557367 :   long *pt, i, l = lg(B->subFB);
     252             : 
     253     1777958 :   for (pt = B->hashtab[hashv]; ; pt = (long*) pt[0])
     254             :   {
     255     1998549 :     if (!pt)
     256             :     {
     257     1469748 :       pt = (long*) pari_malloc((l+3) * sizeof(long));
     258     1469748 :       *pt++ = nrho; /* nrho = pt[-3] */
     259     1469748 :       *pt++ = np;   /* np   = pt[-2] */
     260     1469748 :       *pt++ = q;    /* q    = pt[-1] */
     261     1469748 :       pt[0] = (long)B->hashtab[hashv];
     262     1469748 :       for (i=1; i<l; i++) pt[i]=ex[i];
     263     1469748 :       B->hashtab[hashv]=pt; return NULL;
     264             :     }
     265      308210 :     if (pt[-1] == q) break;
     266             :   }
     267      122707 :   for(i=1; i<l; i++)
     268      119032 :     if (pt[i] != ex[i]) return pt;
     269        3675 :   return (pt[-2]==np)? NULL: pt;
     270             : }
     271             : 
     272             : static void
     273       39638 : clearhash(long **hash)
     274             : {
     275             :   long *pt;
     276             :   long i;
     277    40628950 :   for (i=0; i<HASHT; i++) {
     278    82648372 :     for (pt = hash[i]; pt; ) {
     279     1469748 :       void *z = (void*)(pt - 3);
     280     1469748 :       pt = (long*) pt[0]; pari_free(z);
     281             :     }
     282    40589312 :     hash[i] = NULL;
     283             :   }
     284       39638 : }
     285             : 
     286             : /* last prime stored */
     287             : ulong
     288       16274 : GRH_last_prime(GRHcheck_t *S) { return (S->primes + S->nprimes-1)->p; }
     289             : /* ensure that S->primes can hold at least nb primes */
     290             : void
     291       71236 : GRH_ensure(GRHcheck_t *S, long nb)
     292             : {
     293       71236 :   if (S->maxprimes <= nb)
     294             :   {
     295           0 :     do S->maxprimes *= 2; while (S->maxprimes <= nb);
     296           0 :     S->primes = (GRHprime_t*)pari_realloc((void*)S->primes,
     297           0 :                                           S->maxprimes*sizeof(*S->primes));
     298             :   }
     299       71236 : }
     300             : /* cache data for all primes up to the LIM */
     301             : static void
     302      250311 : cache_prime_quad(GRHcheck_t *S, ulong LIM, GEN D)
     303             : {
     304             :   GRHprime_t *pr;
     305             :   long nb;
     306             : 
     307      250311 :   if (S->limp >= LIM) return;
     308       17870 :   nb = (long)primepi_upper_bound((double)LIM); /* #{p <= LIM} <= nb */
     309       17870 :   GRH_ensure(S, nb+1); /* room for one extra prime */
     310       17870 :   for (pr = S->primes + S->nprimes;;)
     311     2208196 :   {
     312     2226066 :     ulong p = u_forprime_next(&(S->P));
     313     2226066 :     pr->p = p;
     314     2226066 :     pr->logp = log((double)p);
     315     2226066 :     pr->dec = (GEN)kroiu(D,p);
     316     2226066 :     S->nprimes++;
     317     2226066 :     pr++;
     318             :     /* store up to nextprime(LIM) included */
     319     2226066 :     if (p >= LIM) { S->limp = p; break; }
     320             :   }
     321             : }
     322             : 
     323             : static GEN
     324       16274 : compute_invresquad(GRHcheck_t *S)
     325             : {
     326       16274 :   pari_sp av = avma;
     327       16274 :   GEN invres = real_1(DEFAULTPREC);
     328       16274 :   GRHprime_t *pr = S->primes;
     329       16274 :   long i = S->nprimes, LIMC = GRH_last_prime(S)+diffptr[i]-1; /* nextprime(p+1)-1*/
     330       16274 :   double limp = log((double)LIMC) / 2;
     331     2242340 :   for (; i > 0; pr++, i--)
     332             :   {
     333     2226066 :     long s = (long)pr->dec;
     334     2226066 :     if (s)
     335             :     {
     336     2197684 :       ulong p = pr->p;
     337     2197684 :       if (s>0 || pr->logp <= limp)
     338             :         /* Both p and P contribute */
     339     1150112 :         invres = mulur(p - s, divru(invres, p));
     340     1047572 :       else if (s<0)
     341             :         /* Only p contributes */
     342     1047572 :         invres = mulur(p, divru(invres, p - 1));
     343             :     }
     344             :   }
     345       16274 :   return gerepileuptoleaf(av, invres);
     346             : }
     347             : 
     348             : /* p | conductor of order of disc D ? */
     349             : static int
     350       80907 : is_bad(GEN D, ulong p)
     351             : {
     352       80907 :   pari_sp av = avma;
     353       80907 :   if (p == 2)
     354             :   {
     355       21636 :     long r = mod16(D) >> 1;
     356       21636 :     if (r && signe(D) < 0) r = 8-r;
     357       21636 :     return (r < 4);
     358             :   }
     359       59271 :   return gc_bool(av, dvdii(D, sqru(p))); /* p^2 | D ? */
     360             : }
     361             : 
     362             : /* returns the n-th suitable ideal for the factorbase */
     363             : static long
     364       16274 : nthidealquad(GEN D, long n)
     365             : {
     366       16274 :   pari_sp av = avma;
     367             :   forprime_t S;
     368             :   ulong p;
     369       16274 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
     370       88481 :   while ((p = u_forprime_next(&S)) && n > 0)
     371       55933 :     if (!is_bad(D, p) && kroiu(D, p) >= 0) n--;
     372       16274 :   return gc_long(av, p);
     373             : }
     374             : 
     375             : static int
     376      217588 : quadGRHchk(GEN D, GRHcheck_t *S, ulong LIMC)
     377             : {
     378      217588 :   double logC = log((double)LIMC), SA = 0, SB = 0;
     379             :   long i;
     380             : 
     381      217588 :   cache_prime_quad(S, LIMC, D);
     382    10531147 :   for (i = 0;; i++)
     383    10313559 :   {
     384    10531147 :     GRHprime_t *pr = S->primes+i;
     385    10531147 :     ulong p = pr->p;
     386             :     long M;
     387             :     double logNP, q, A, B;
     388    10531147 :     if (p > LIMC) break;
     389    10336049 :     if ((long)pr->dec < 0)
     390             :     {
     391     5079042 :       logNP = 2 * pr->logp;
     392     5079042 :       q = 1/(double)p;
     393             :     }
     394             :     else
     395             :     {
     396     5257007 :       logNP = pr->logp;
     397     5257007 :       q = 1/sqrt((double)p);
     398             :     }
     399    10336049 :     A = logNP * q; B = logNP * A; M = (long)(logC/logNP);
     400    10336049 :     if (M > 1)
     401             :     {
     402      614351 :       double inv1_q = 1 / (1-q);
     403      614351 :       A *= (1 - pow(q, (double) M)) * inv1_q;
     404      614351 :       B *= (1 - pow(q, (double) M)*(M+1 - M*q)) * inv1_q * inv1_q;
     405             :     }
     406    10336049 :     if ((long)pr->dec>0) { SA += 2*A;SB += 2*B; } else { SA += A; SB += B; }
     407    10336049 :     if (p == LIMC) break;
     408             :   }
     409      217588 :   return GRHok(S, logC, SA, SB);
     410             : }
     411             : 
     412             : /* C2 >= C1; create B->FB, B->numFB; set B->badprim. Return L(kro_D, 1) */
     413             : static void
     414       16449 : FBquad(struct buch_quad *B, ulong C2, ulong C1, GRHcheck_t *S)
     415             : {
     416       16449 :   GEN D = B->QFR->D;
     417             :   long i;
     418             :   pari_sp av;
     419             :   GRHprime_t *pr;
     420             : 
     421       16449 :   cache_prime_quad(S, C2, D);
     422       16449 :   pr = S->primes;
     423       16449 :   B->numFB = cgetg(C2+1, t_VECSMALL);
     424       16449 :   B->FB    = cgetg(C2+1, t_VECSMALL);
     425       16449 :   av = avma;
     426       16449 :   B->KC = 0; i = 0;
     427       16449 :   B->badprim = gen_1;
     428      896937 :   for (;; pr++) /* p <= C2 */
     429      896937 :   {
     430      913386 :     ulong p = pr->p;
     431      913386 :     if (!B->KC && p > C1) B->KC = i;
     432      913386 :     if (p > C2) break;
     433      898206 :     switch ((long)pr->dec)
     434             :     {
     435      436726 :       case -1: break; /* inert */
     436             :       case  0: /* ramified */
     437       24974 :         if (is_bad(D, p)) { B->badprim = muliu(B->badprim, p); break; }
     438             :         /* fall through */
     439             :       default:  /* split */
     440      461256 :         i++; B->numFB[p] = i; B->FB[i] = p; break;
     441             :     }
     442      898206 :     if (p == C2)
     443             :     {
     444        1269 :       if (!B->KC) B->KC = i;
     445        1269 :       break;
     446             :     }
     447             :   }
     448       16449 :   B->KC2 = i;
     449       16449 :   setlg(B->FB, B->KC2+1);
     450       16449 :   if (B->badprim != gen_1)
     451         196 :     B->badprim = gerepileuptoint(av, B->badprim);
     452             :   else
     453             :   {
     454       16253 :     B->badprim = NULL; set_avma(av);
     455             :   }
     456       16449 : }
     457             : 
     458             : /* create B->vperm, return B->subFB */
     459             : static GEN
     460       16449 : subFBquad(struct buch_quad *B, GEN D, double PROD, long minSFB)
     461             : {
     462       16449 :   long i, j, lgsub = 1, ino = 1, lv = B->KC+1;
     463       16449 :   double prod = 1.;
     464             :   pari_sp av;
     465             :   GEN no;
     466             : 
     467       16449 :   B->vperm = cgetg(lv, t_VECSMALL);
     468       16449 :   av = avma;
     469       16449 :   no    = cgetg(lv, t_VECSMALL);
     470       84037 :   for (j = 1; j < lv; j++)
     471             :   {
     472       83855 :     ulong p = uel(B->FB,j);
     473       83855 :     if (!umodiu(D, p)) no[ino++] = j; /* ramified */
     474             :     else
     475             :     {
     476       64706 :       B->vperm[lgsub++] = j;
     477       64706 :       prod *= p;
     478       64706 :       if (lgsub > minSFB && prod > PROD) break;
     479             :     }
     480             :   }
     481             :   /* lgsub >= 1 otherwise quadGRHchk is false */
     482       16449 :   i = lgsub;
     483       16449 :   for (j = 1; j < ino;i++,j++) B->vperm[i] = no[j];
     484       16449 :   for (     ; i < lv; i++)     B->vperm[i] = i;
     485       16449 :   no = gclone(vecslice(B->vperm, 1, lgsub-1));
     486       16449 :   set_avma(av); return no;
     487             : }
     488             : 
     489             : /* assume n >= 1, x[i][j] = B->subFB[i]^j, for j = 1..n */
     490             : static GEN
     491       23189 : powsubFBquad(struct buch_quad *B, long n)
     492             : {
     493       23189 :   pari_sp av = avma;
     494       23189 :   long i,j, l = lg(B->subFB);
     495       23189 :   GEN F, y, x = cgetg(l, t_VEC), D = B->QFR->D;
     496             : 
     497       23189 :   if (B->PRECREG) /* real */
     498             :   {
     499       39088 :     for (i=1; i<l; i++)
     500             :     {
     501       34055 :       F = qfr5_pf(B->QFR, B->FB[B->subFB[i]], B->PRECREG);
     502       34055 :       y = cgetg(n+1, t_VEC); gel(x,i) = y;
     503       34055 :       gel(y,1) = F;
     504       34055 :       for (j=2; j<=n; j++) gel(y,j) = QFR5_comp(gel(y,j-1), F, B->QFR);
     505             :     }
     506             :   }
     507             :   else /* imaginary */
     508             :   {
     509       85425 :     for (i=1; i<l; i++)
     510             :     {
     511       67269 :       F = qfi_pf(D, B->FB[B->subFB[i]]);
     512       67269 :       y = cgetg(n+1, t_VEC); gel(x,i) = y;
     513       67269 :       gel(y,1) = F;
     514       67269 :       for (j=2; j<=n; j++) gel(y,j) = qficomp(gel(y,j-1), F);
     515             :     }
     516             :   }
     517       23189 :   x = gclone(x); set_avma(av); return x;
     518             : }
     519             : 
     520             : static void
     521       41958 : sub_fact(struct buch_quad *B, GEN col, GEN F)
     522             : {
     523       41958 :   GEN b = gel(F,2);
     524             :   long i;
     525      128739 :   for (i=1; i<=B->primfact[0]; i++)
     526             :   {
     527       86781 :     ulong p = B->primfact[i], k = B->numFB[p];
     528       86781 :     long e = B->exprimfact[i];
     529       86781 :     if (umodiu(b, p<<1) > p) e = -e;
     530       86781 :     col[k] -= e;
     531             :   }
     532       41958 : }
     533             : static void
     534      590367 : add_fact(struct buch_quad *B, GEN col, GEN F)
     535             : {
     536      590367 :   GEN b = gel(F,2);
     537             :   long i;
     538     2296829 :   for (i=1; i<=B->primfact[0]; i++)
     539             :   {
     540     1706462 :     ulong p = B->primfact[i], k = B->numFB[p];
     541     1706462 :     long e = B->exprimfact[i];
     542     1706462 :     if (umodiu(b, p<<1) > p) e = -e;
     543     1706462 :     col[k] += e;
     544             :   }
     545      590367 : }
     546             : 
     547             : static GEN
     548       16274 : get_clgp(struct buch_quad *B, GEN W, GEN *ptD, long prec)
     549             : {
     550       16274 :   GEN res, init, u1, D = ZM_snf_group(W,NULL,&u1), Z = prec? real_0(prec): NULL;
     551       16274 :   long i, j, l = lg(W), c = lg(D);
     552             : 
     553       16274 :   res=cgetg(c,t_VEC); init = cgetg(l,t_VEC);
     554       16274 :   for (i=1; i<l; i++) gel(init,i) = primeform_u(B->QFR->D, B->FB[B->vperm[i]]);
     555       44587 :   for (j=1; j<c; j++)
     556             :   {
     557       28313 :     GEN g = NULL;
     558       28313 :     if (prec)
     559             :     {
     560       13293 :       for (i=1; i<l; i++)
     561             :       {
     562        9583 :         GEN t, u = gcoeff(u1,i,j);
     563        9583 :         if (!signe(u)) continue;
     564        4389 :         t = qfr3_pow(gel(init,i), u, B->QFR);
     565        4389 :         g = g? qfr3_comp(g, t, B->QFR): t;
     566             :       }
     567        3710 :       g = qfr3_to_qfr(qfr3_canon_safe(qfr3_red(g, B->QFR), B->QFR), Z);
     568             :     }
     569             :     else
     570             :     {
     571       86176 :       for (i=1; i<l; i++)
     572             :       {
     573       61573 :         GEN t, u = gcoeff(u1,i,j);
     574       61573 :         if (!signe(u)) continue;
     575       39264 :         t = powgi(gel(init,i), u);
     576       39264 :         g = g? qficomp(g, t): t;
     577             :       }
     578             :     }
     579       28313 :     gel(res,j) = g;
     580             :   }
     581       16274 :   *ptD = D; return res;
     582             : }
     583             : 
     584             : static long
     585       16274 : trivial_relations(struct buch_quad *B, GEN mat, GEN C)
     586             : {
     587       16274 :   long i, j = 0;
     588       16274 :   GEN col, D = B->QFR->D;
     589      476361 :   for (i = 1; i <= B->KC; i++)
     590             :   { /* ramified prime ==> trivial relation */
     591      460087 :     if (umodiu(D, B->FB[i])) continue;
     592       24260 :     col = zero_zv(B->KC);
     593       24260 :     col[i] = 2; j++;
     594       24260 :     gel(mat,j) = col;
     595       24260 :     gel(C,j) = gen_0;
     596             :   }
     597       16274 :   return j;
     598             : }
     599             : 
     600             : static void
     601           0 : dbg_all(pari_timer *T, const char *phase, long s, long n)
     602             : {
     603           0 :   err_printf("\n");
     604           0 :   timer_printf(T, "%s rel [#rel/#test = %ld/%ld]", phase,s,n);
     605           0 : }
     606             : 
     607             : /* Imaginary Quadratic fields */
     608             : 
     609             : static void
     610       18701 : imag_relations(struct buch_quad *B, long need, long *pc, ulong LIMC, GEN mat)
     611             : {
     612             :   pari_timer T;
     613       18701 :   long lgsub = lg(B->subFB), current = *pc, nbtest = 0, s = 0;
     614             :   long i, fpc;
     615             :   pari_sp av;
     616       18701 :   GEN col, form, ex = cgetg(lgsub, t_VECSMALL);
     617             : 
     618       18701 :   if (!current) current = 1;
     619       18701 :   if (DEBUGLEVEL>2) timer_start(&T);
     620       18701 :   av = avma;
     621             :   for(;;)
     622             :   {
     623     2813429 :     if (s >= need) break;
     624     1397364 :     set_avma(av);
     625     1397364 :     form = qfi_random(B,ex);
     626     1397364 :     form = qficomp(form, qfi_pf(B->QFR->D, B->FB[current]));
     627     1397364 :     nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
     628     1397364 :     if (!fpc)
     629             :     {
     630      287527 :       if (DEBUGLEVEL>3) err_printf(".");
     631      287527 :       if ((nbtest & 0xff) == 0 && ++current > B->KC) current = 1;
     632      287527 :       continue;
     633             :     }
     634     1109837 :     if (fpc > 1)
     635             :     {
     636      793254 :       long *fpd = largeprime(B,fpc,ex,current,0);
     637             :       ulong b1, b2, p;
     638             :       GEN form2;
     639      793254 :       if (!fpd)
     640             :       {
     641      745640 :         if (DEBUGLEVEL>3) err_printf(".");
     642      745640 :         continue;
     643             :       }
     644       47614 :       form2 = qficomp(qfi_factorback(B,fpd), qfi_pf(B->QFR->D, B->FB[fpd[-2]]));
     645       47614 :       p = fpc << 1;
     646       47614 :       b1 = umodiu(gel(form2,2), p);
     647       47614 :       b2 = umodiu(gel(form,2),  p);
     648       47614 :       if (b1 != b2 && b1+b2 != p) continue;
     649             : 
     650       47579 :       col = gel(mat,++s);
     651       47579 :       add_fact(B,col, form);
     652       47579 :       (void)factorquad(B,form2,B->KC,LIMC);
     653       47579 :       if (b1==b2)
     654             :       {
     655       23723 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
     656       23723 :         sub_fact(B, col, form2); col[fpd[-2]]++;
     657             :       }
     658             :       else
     659             :       {
     660       23856 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
     661       23856 :         add_fact(B, col, form2); col[fpd[-2]]--;
     662             :       }
     663       47579 :       if (DEBUGLEVEL>2) err_printf(" %ldP",s);
     664             :     }
     665             :     else
     666             :     {
     667      316583 :       col = gel(mat,++s);
     668      316583 :       for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
     669      316583 :       add_fact(B, col, form);
     670      316583 :       if (DEBUGLEVEL>2) err_printf(" %ld",s);
     671             :     }
     672      364162 :     col[current]--;
     673      364162 :     if (++current > B->KC) current = 1;
     674             :   }
     675       18701 :   if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
     676       18701 :   *pc = current;
     677       18701 : }
     678             : 
     679             : static int
     680           7 : imag_be_honest(struct buch_quad *B)
     681             : {
     682           7 :   long p, fpc, s = B->KC, nbtest = 0;
     683           7 :   GEN F, ex = cgetg(lg(B->subFB), t_VECSMALL);
     684           7 :   pari_sp av = avma;
     685             : 
     686         483 :   while (s<B->KC2)
     687             :   {
     688         469 :     p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
     689         469 :     F = qficomp(qfi_pf(B->QFR->D, p), qfi_random(B, ex));
     690         469 :     fpc = factorquad(B,F,s,p-1);
     691         469 :     if (fpc == 1) { nbtest=0; s++; }
     692             :     else
     693         322 :       if (++nbtest > 40) return 0;
     694         469 :     set_avma(av);
     695             :   }
     696           7 :   return 1;
     697             : }
     698             : 
     699             : /* Real Quadratic fields */
     700             : 
     701             : static void
     702        5124 : real_relations(struct buch_quad *B, long need, long *pc, long lim, ulong LIMC, GEN mat, GEN C)
     703             : {
     704             :   pari_timer T;
     705        5124 :   long lgsub = lg(B->subFB), prec = B->PRECREG, current = *pc, nbtest=0, s=0;
     706             :   long i, fpc, endcycle, rhoacc, rho;
     707             :   /* in a 2nd phase, don't include FB[current] but run along the cyle
     708             :    * ==> get more units */
     709        5124 :   int first = (current == 0);
     710             :   pari_sp av, av1;
     711        5124 :   GEN d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL);
     712             : 
     713        5124 :   if (DEBUGLEVEL>2) timer_start(&T);
     714        5124 :   if (!current) current = 1;
     715        5124 :   if (lim > need) lim = need;
     716        5124 :   av = avma;
     717             :   for(;;)
     718             :   {
     719      327348 :     if (s >= need) break;
     720      161112 :     if (first && s >= lim) {
     721        2142 :       first = 0;
     722        2142 :       if (DEBUGLEVEL>2) dbg_all(&T, "initial", s, nbtest);
     723             :     }
     724      161112 :     set_avma(av); form = qfr3_random(B, ex);
     725      161112 :     if (!first)
     726      158935 :       form = QFR3_comp(form, qfr3_pf(B->QFR, B->FB[current]), B->QFR);
     727      161112 :     av1 = avma;
     728      161112 :     form0 = form; form1 = NULL;
     729      161112 :     endcycle = rhoacc = 0;
     730      161112 :     rho = -1;
     731             : 
     732             : CYCLE:
     733     1297569 :     if (endcycle || rho > 5000)
     734             :     {
     735          21 :       if (++current > B->KC) current = 1;
     736          21 :       continue;
     737             :     }
     738     1297548 :     if (gc_needed(av,1))
     739             :     {
     740           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"real_relations");
     741           0 :       gerepileall(av1, form1? 2: 1, &form, &form1);
     742             :     }
     743     1297548 :     if (rho < 0) rho = 0; /* first time in */
     744             :     else
     745             :     {
     746     1136436 :       form = qfr3_rho(form, B->QFR); rho++;
     747     1136436 :       rhoacc++;
     748     1136436 :       if (first)
     749      444430 :         endcycle = (absequalii(gel(form,1),gel(form0,1))
     750      222215 :              && equalii(gel(form,2),gel(form0,2)));
     751             :       else
     752             :       {
     753      914221 :         if (absequalii(gel(form,1), gel(form,3))) /* a = -c */
     754             :         {
     755           0 :           if (absequalii(gel(form,1),gel(form0,1)) &&
     756           0 :                   equalii(gel(form,2),gel(form0,2))) continue;
     757           0 :           form = qfr3_rho(form, B->QFR); rho++;
     758           0 :           rhoacc++;
     759             :         }
     760             :         else
     761      914221 :           { setsigne(form[1],1); setsigne(form[3],-1); }
     762      914277 :         if (equalii(gel(form,1),gel(form0,1)) &&
     763          56 :             equalii(gel(form,2),gel(form0,2))) continue;
     764             :       }
     765             :     }
     766     1297548 :     nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
     767     1297548 :     if (!fpc)
     768             :     {
     769      385511 :       if (DEBUGLEVEL>3) err_printf(".");
     770      385511 :       goto CYCLE;
     771             :     }
     772      912037 :     if (fpc > 1)
     773             :     { /* look for Large Prime relation */
     774      764113 :       long *fpd = largeprime(B,fpc,ex,first? 0: current,rhoacc);
     775             :       ulong b1, b2, p;
     776             :       GEN form2;
     777      764113 :       if (!fpd)
     778             :       {
     779      727783 :         if (DEBUGLEVEL>3) err_printf(".");
     780      727783 :         goto CYCLE;
     781             :       }
     782       36330 :       if (!form1)
     783             :       {
     784       36330 :         form1 = qfr5_factorback(B,ex);
     785       36330 :         if (!first)
     786       36330 :           form1 = QFR5_comp(form1, qfr5_pf(B->QFR, B->FB[current], prec), B->QFR);
     787             :       }
     788       36330 :       form1 = qfr5_rho_pow(form1, rho, B->QFR);
     789       36330 :       rho = 0;
     790             : 
     791       36330 :       form2 = qfr5_factorback(B,fpd);
     792       36330 :       if (fpd[-2])
     793       23954 :         form2 = QFR5_comp(form2, qfr5_pf(B->QFR, B->FB[fpd[-2]], prec), B->QFR);
     794       36330 :       form2 = qfr5_rho_pow(form2, fpd[-3], B->QFR);
     795       36330 :       if (!absequalii(gel(form2,1),gel(form2,3)))
     796             :       {
     797       36330 :         setsigne(form2[1], 1);
     798       36330 :         setsigne(form2[3],-1);
     799             :       }
     800       36330 :       p = fpc << 1;
     801       36330 :       b1 = umodiu(gel(form2,2), p);
     802       36330 :       b2 = umodiu(gel(form1,2), p);
     803       36330 :       if (b1 != b2 && b1+b2 != p) goto CYCLE;
     804             : 
     805       36330 :       col = gel(mat,++s);
     806       36330 :       add_fact(B, col, form1);
     807       36330 :       (void)factorquad(B,form2,B->KC,LIMC);
     808       36330 :       if (b1==b2)
     809             :       {
     810       18235 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
     811       18235 :         sub_fact(B,col, form2);
     812       18235 :         if (fpd[-2]) col[fpd[-2]]++;
     813       36470 :         d = qfr5_dist(subii(gel(form1,4),gel(form2,4)),
     814       36470 :                       divrr(gel(form1,5),gel(form2,5)), prec);
     815             :       }
     816             :       else
     817             :       {
     818       18095 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
     819       18095 :         add_fact(B, col, form2);
     820       18095 :         if (fpd[-2]) col[fpd[-2]]--;
     821       36190 :         d = qfr5_dist(addii(gel(form1,4),gel(form2,4)),
     822       36190 :                       mulrr(gel(form1,5),gel(form2,5)), prec);
     823             :       }
     824       36330 :       if (DEBUGLEVEL>2) err_printf(" %ldP",s);
     825             :     }
     826             :     else
     827             :     { /* standard relation */
     828      147924 :       if (!form1)
     829             :       {
     830      124782 :         form1 = qfr5_factorback(B, ex);
     831      124782 :         if (!first)
     832      122605 :           form1 = QFR5_comp(form1, qfr5_pf(B->QFR, B->FB[current], prec), B->QFR);
     833             :       }
     834      147924 :       form1 = qfr5_rho_pow(form1, rho, B->QFR);
     835      147924 :       rho = 0;
     836             : 
     837      147924 :       col = gel(mat,++s);
     838      147924 :       for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
     839      147924 :       add_fact(B, col, form1);
     840      147924 :       d = qfr5_dist(gel(form1,4), gel(form1,5), prec);
     841      147924 :       if (DEBUGLEVEL>2) err_printf(" %ld",s);
     842             :     }
     843      184254 :     affrr(d, gel(C,s));
     844      184254 :     if (first)
     845             :     {
     846       25319 :       if (s >= lim) continue;
     847       23163 :       goto CYCLE;
     848             :     }
     849             :     else
     850             :     {
     851      158935 :       col[current]--;
     852      158935 :       if (++current > B->KC) current = 1;
     853             :     }
     854             :   }
     855        5124 :   if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
     856        5124 :   *pc = current;
     857        5124 : }
     858             : 
     859             : static int
     860           7 : real_be_honest(struct buch_quad *B)
     861             : {
     862           7 :   long p, fpc, s = B->KC, nbtest = 0;
     863           7 :   GEN F,F0, ex = cgetg(lg(B->subFB), t_VECSMALL);
     864           7 :   pari_sp av = avma;
     865             : 
     866          35 :   while (s<B->KC2)
     867             :   {
     868          21 :     p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
     869          21 :     F = QFR3_comp(qfr3_random(B, ex), qfr3_pf(B->QFR, p), B->QFR);
     870          21 :     for (F0 = F;;)
     871             :     {
     872          63 :       fpc = factorquad(B,F,s,p-1);
     873          42 :       if (fpc == 1) { nbtest=0; s++; break; }
     874          21 :       if (++nbtest > 40) return 0;
     875          21 :       F = qfr3_canon(qfr3_rho(F, B->QFR), B->QFR);
     876          21 :       if (equalii(gel(F,1),gel(F0,1))
     877           0 :        && equalii(gel(F,2),gel(F0,2))) break;
     878             :     }
     879          21 :     set_avma(av);
     880             :   }
     881           7 :   return 1;
     882             : }
     883             : 
     884             : static GEN
     885       32025 : gcdreal(GEN a,GEN b)
     886             : {
     887       32025 :   if (!signe(a)) return mpabs_shallow(b);
     888       31361 :   if (!signe(b)) return mpabs_shallow(a);
     889       31213 :   if (typ(a)==t_INT)
     890             :   {
     891           0 :     if (typ(b)==t_INT) return gcdii(a,b);
     892           0 :     a = itor(a, lg(b));
     893             :   }
     894       31213 :   else if (typ(b)==t_INT)
     895           0 :     b = itor(b, lg(a));
     896       31213 :   if (expo(a)<-5) return absr(b);
     897       11305 :   if (expo(b)<-5) return absr(a);
     898        8456 :   a = absr(a); b = absr(b);
     899       27135 :   while (expo(b) >= -5  && signe(b))
     900             :   {
     901             :     long e;
     902       10223 :     GEN r, q = gcvtoi(divrr(a,b),&e);
     903       10223 :     if (e > 0) return NULL;
     904       10223 :     r = subrr(a, mulir(q,b)); a = b; b = r;
     905             :   }
     906        8456 :   return mpabs_shallow(a);
     907             : }
     908             : 
     909             : static int
     910       19911 : get_R(struct buch_quad *B, GEN C, long sreg, GEN z, GEN *ptR)
     911             : {
     912       19911 :   GEN R = gen_1;
     913             :   double c;
     914             :   long i;
     915             : 
     916       19911 :   if (B->PRECREG)
     917             :   {
     918        2877 :     R = mpabs_shallow(gel(C,1));
     919       34902 :     for (i=2; i<=sreg; i++)
     920             :     {
     921       32025 :       R = gcdreal(gel(C,i), R);
     922       32025 :       if (!R) return fupb_PRECI;
     923             :     }
     924        2877 :     if (gexpo(R) <= -3)
     925             :     {
     926           0 :       if (DEBUGLEVEL>2) err_printf("regulator is zero.\n");
     927           0 :       return fupb_RELAT;
     928             :     }
     929        2877 :     if (DEBUGLEVEL>2) err_printf("#### Tentative regulator: %Ps\n",R);
     930             :   }
     931       19911 :   c = gtodouble(gmul(z, R));
     932       19911 :   if (c < 0.8 || c > 1.3) return fupb_RELAT;
     933       16274 :   *ptR = R; return fupb_NONE;
     934             : }
     935             : 
     936             : static int
     937       16274 : quad_be_honest(struct buch_quad *B)
     938             : {
     939             :   int r;
     940       16274 :   if (B->KC2 <= B->KC) return 1;
     941          14 :   if (DEBUGLEVEL>2)
     942           0 :     err_printf("be honest for primes from %ld to %ld\n", B->FB[B->KC+1],B->FB[B->KC2]);
     943          14 :   r = B->PRECREG? real_be_honest(B): imag_be_honest(B);
     944          14 :   if (DEBUGLEVEL>2) err_printf("\n");
     945          14 :   return r;
     946             : }
     947             : 
     948             : 
     949             : GEN
     950       16281 : Buchquad(GEN D, double cbach, double cbach2, long prec)
     951             : {
     952       16281 :   const long MAXRELSUP = 7, SFB_MAX = 3;
     953             :   pari_timer T;
     954       16281 :   pari_sp av0 = avma, av, av2;
     955       16281 :   const long RELSUP = 5;
     956             :   long i, s, current, triv, sfb_trials, nrelsup, nreldep, need, nsubFB, minSFB;
     957             :   ulong low, high, LIMC0, LIMC, LIMC2, LIMCMAX, cp;
     958       16281 :   GEN W, cyc, res, gen, dep, mat, C, extraC, B, R, invhr, h = NULL; /*-Wall*/
     959             :   double drc, sdrc, lim, LOGD, LOGD2;
     960             :   GRHcheck_t GRHcheck;
     961             :   struct qfr_data QFR;
     962             :   struct buch_quad BQ;
     963       16281 :   int FIRST = 1;
     964             : 
     965       16281 :   check_quaddisc(D, &s, /*junk*/&i, "Buchquad");
     966       16281 :   R = NULL; /* -Wall */
     967       16281 :   BQ.QFR = &QFR;
     968       16281 :   QFR.D = D;
     969       16281 :   if (s < 0)
     970             :   {
     971       14125 :     if (abscmpiu(QFR.D,4) <= 0)
     972             :     {
     973           7 :       GEN z = cgetg(5,t_VEC);
     974           7 :       gel(z,1) = gel(z,4) = gen_1; gel(z,2) = gel(z,3) = cgetg(1,t_VEC);
     975           7 :       return z;
     976             :     }
     977       14118 :     prec = BQ.PRECREG = 0;
     978             :   } else {
     979        2156 :     BQ.PRECREG = maxss(prec+EXTRAPRECWORD, nbits2prec(2*expi(QFR.D) + 128));
     980             :   }
     981       16274 :   if (DEBUGLEVEL>2) timer_start(&T);
     982       16274 :   BQ.primfact   = new_chunk(100);
     983       16274 :   BQ.exprimfact = new_chunk(100);
     984       16274 :   BQ.hashtab = (long**) new_chunk(HASHT);
     985       16274 :   for (i=0; i<HASHT; i++) BQ.hashtab[i] = NULL;
     986             : 
     987       16274 :   drc = fabs(gtodouble(QFR.D));
     988       16274 :   LOGD = log(drc);
     989       16274 :   LOGD2 = LOGD * LOGD;
     990             : 
     991       16274 :   sdrc = lim = sqrt(drc);
     992       16274 :   if (!BQ.PRECREG) lim /= sqrt(3.);
     993       16274 :   cp = (ulong)exp(sqrt(LOGD*log(LOGD)/8.0));
     994       16274 :   if (cp < 20) cp = 20;
     995       16274 :   if (cbach > 6.) {
     996           0 :     if (cbach2 < cbach) cbach2 = cbach;
     997           0 :     cbach = 6.;
     998             :   }
     999       16274 :   if (cbach < 0.)
    1000           0 :     pari_err_DOMAIN("Buchquad","Bach constant","<",gen_0,dbltor(cbach));
    1001       16274 :   av = avma;
    1002       16274 :   BQ.powsubFB = BQ.subFB = NULL;
    1003       16274 :   minSFB = (expi(D) > 15)? 3: 2;
    1004       16274 :   init_GRHcheck(&GRHcheck, 2, BQ.PRECREG? 2: 0, LOGD);
    1005       16274 :   high = low = LIMC0 = maxss((long)(cbach2*LOGD2), 1);
    1006       16274 :   LIMCMAX = (long)(6.*LOGD2);
    1007             :   /* 97/1223 below to ensure a good enough approximation of residue */
    1008       16274 :   cache_prime_quad(&GRHcheck, expi(D) < 16 ? 97: 1223, D);
    1009      141307 :   while (!quadGRHchk(D, &GRHcheck, high))
    1010             :   {
    1011      108759 :     low = high;
    1012      108759 :     high *= 2;
    1013             :   }
    1014      125103 :   while (high - low > 1)
    1015             :   {
    1016       92555 :     long test = (low+high)/2;
    1017       92555 :     if (quadGRHchk(D, &GRHcheck, test))
    1018       49196 :       high = test;
    1019             :     else
    1020       43359 :       low = test;
    1021             :   }
    1022       16274 :   if (high == LIMC0+1 && quadGRHchk(D, &GRHcheck, LIMC0))
    1023           0 :     LIMC2 = LIMC0;
    1024             :   else
    1025       16274 :     LIMC2 = high;
    1026       16274 :   if (LIMC2 > LIMCMAX) LIMC2 = LIMCMAX;
    1027       16274 :   LIMC0 = (long)(cbach*LOGD2);
    1028       16274 :   LIMC = cbach ? LIMC0 : LIMC2;
    1029       16274 :   LIMC = maxss(LIMC, nthidealquad(D, 2));
    1030             : 
    1031             : /* LIMC = Max(cbach*(log D)^2, exp(sqrt(log D loglog D) / 8)) */
    1032             : START:
    1033             :   do
    1034             :   {
    1035       16449 :     if (!FIRST) LIMC = bnf_increase_LIMC(LIMC,LIMCMAX);
    1036       16449 :     if (DEBUGLEVEL>2 && LIMC > LIMC0)
    1037           0 :       err_printf("%s*** Bach constant: %f\n", FIRST?"":"\n", LIMC/LOGD2);
    1038       16449 :     FIRST = 0; set_avma(av);
    1039       16449 :     if (BQ.subFB) gunclone(BQ.subFB);
    1040       16449 :     if (BQ.powsubFB) gunclone(BQ.powsubFB);
    1041       16449 :     clearhash(BQ.hashtab);
    1042       16449 :     if (LIMC < cp) LIMC = cp;
    1043       16449 :     if (LIMC2 < LIMC) LIMC2 = LIMC;
    1044       16449 :     if (BQ.PRECREG) qfr_data_init(QFR.D, BQ.PRECREG, &QFR);
    1045             : 
    1046       16449 :     FBquad(&BQ, LIMC2, LIMC, &GRHcheck);
    1047       16449 :     if (DEBUGLEVEL>2) timer_printf(&T, "factor base");
    1048       16449 :     BQ.subFB = subFBquad(&BQ, QFR.D, lim + 0.5, minSFB);
    1049       16449 :     if (DEBUGLEVEL>2) timer_printf(&T, "subFBquad = %Ps",
    1050             :                                  vecpermute(BQ.FB, BQ.subFB));
    1051       16449 :     nsubFB = lg(BQ.subFB) - 1;
    1052             :   }
    1053       16449 :   while (nsubFB < (expi(D) > 15 ? 3 : 2));
    1054             :   /* invhr = 2^r1 (2pi)^r2 / sqrt(D) w ~ L(chi,1) / hR */
    1055       16274 :   invhr = gmul(dbltor((BQ.PRECREG?2.:M_PI)/sdrc), compute_invresquad(&GRHcheck));
    1056       16274 :   BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
    1057       16274 :   if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
    1058       16274 :   BQ.limhash = (LIMC & HIGHMASK)? (HIGHBIT>>1): LIMC*LIMC;
    1059             : 
    1060       16274 :   need = BQ.KC + RELSUP - 2;
    1061       16274 :   current = 0;
    1062       16274 :   W = NULL;
    1063       16274 :   sfb_trials = nreldep = nrelsup = 0;
    1064       16274 :   s = nsubFB + RELSUP;
    1065       16274 :   av2 = avma;
    1066             : 
    1067             :   do
    1068             :   {
    1069       23825 :     if ((nreldep & 3) == 1 || (nrelsup & 7) == 1) {
    1070        6915 :       if (DEBUGLEVEL>2) err_printf("*** Changing sub factor base\n");
    1071        6915 :       gunclone(BQ.subFB);
    1072        6915 :       gunclone(BQ.powsubFB);
    1073        6915 :       BQ.subFB = gclone(vecslice(BQ.vperm, 1, nsubFB));
    1074        6915 :       BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
    1075        6915 :       if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
    1076        6915 :       clearhash(BQ.hashtab);
    1077             :     }
    1078       23825 :     need += 2;
    1079       23825 :     mat    = cgetg(need+1, t_MAT);
    1080       23825 :     extraC = cgetg(need+1, t_VEC);
    1081       23825 :     if (!W) { /* first time */
    1082       16274 :       C = extraC;
    1083       16274 :       triv = trivial_relations(&BQ, mat, C);
    1084       16274 :       if (DEBUGLEVEL>2) err_printf("KC = %ld, need %ld relations\n", BQ.KC, need);
    1085             :     } else {
    1086        7551 :       triv = 0;
    1087        7551 :       if (DEBUGLEVEL>2) err_printf("...need %ld more relations\n", need);
    1088             :     }
    1089       23825 :     if (BQ.PRECREG) {
    1090      189378 :       for (i = triv+1; i<=need; i++) {
    1091      184254 :         gel(mat,i) = zero_zv(BQ.KC);
    1092      184254 :         gel(extraC,i) = cgetr(BQ.PRECREG);
    1093             :       }
    1094        5124 :       real_relations(&BQ, need - triv, &current, s,LIMC,mat + triv,extraC + triv);
    1095             :     } else {
    1096      382863 :       for (i = triv+1; i<=need; i++) {
    1097      364162 :         gel(mat,i) = zero_zv(BQ.KC);
    1098      364162 :         gel(extraC,i) = gen_0;
    1099             :       }
    1100       18701 :       imag_relations(&BQ, need - triv, &current, LIMC,mat + triv);
    1101             :     }
    1102             : 
    1103       23825 :     if (!W)
    1104       16274 :       W = hnfspec_i(mat,BQ.vperm,&dep,&B,&C,nsubFB);
    1105             :     else
    1106        7551 :       W = hnfadd_i(W,BQ.vperm,&dep,&B,&C, mat,extraC);
    1107       23825 :     gerepileall(av2, 4, &W,&C,&B,&dep);
    1108       23825 :     need = BQ.KC - (lg(W)-1) - (lg(B)-1);
    1109       23825 :     if (need)
    1110             :     {
    1111        3914 :       if (++nreldep > 15 && cbach < 1) goto START;
    1112        3914 :       continue;
    1113             :     }
    1114             : 
    1115       19911 :     h = ZM_det_triangular(W);
    1116       19911 :     if (DEBUGLEVEL>2) err_printf("\n#### Tentative class number: %Ps\n", h);
    1117             : 
    1118       19911 :     switch(get_R(&BQ, C, (lg(C)-1) - (lg(B)-1) - (lg(W)-1), mulir(h,invhr), &R))
    1119             :     {
    1120             :       case fupb_PRECI:
    1121           0 :         BQ.PRECREG = precdbl(BQ.PRECREG);
    1122           0 :         FIRST = 1; goto START;
    1123             : 
    1124             :       case fupb_RELAT:
    1125        3637 :         if (++nrelsup > MAXRELSUP)
    1126             :         {
    1127          21 :           if (++sfb_trials > SFB_MAX && cbach <= 1) goto START;
    1128          21 :           if (nsubFB < minss(10,BQ.KC)) nsubFB++;
    1129             :         }
    1130        3637 :         need = minss(BQ.KC, nrelsup);
    1131             :     }
    1132             :   }
    1133       23825 :   while (need);
    1134             :   /* DONE */
    1135       16274 :   if (!quad_be_honest(&BQ)) goto START;
    1136       16274 :   if (DEBUGLEVEL>2) timer_printf(&T, "be honest");
    1137       16274 :   clearhash(BQ.hashtab);
    1138       16274 :   free_GRHcheck(&GRHcheck);
    1139             : 
    1140       16274 :   gen = get_clgp(&BQ,W,&cyc,prec);
    1141       16274 :   gunclone(BQ.subFB);
    1142       16274 :   gunclone(BQ.powsubFB);
    1143       16274 :   res = cgetg(5,t_VEC);
    1144       16274 :   gel(res,1) = h;
    1145       16274 :   gel(res,2) = cyc;
    1146       16274 :   gel(res,3) = gen;
    1147       16274 :   gel(res,4) = R; return gerepilecopy(av0,res);
    1148             : }
    1149             : 
    1150             : GEN
    1151           0 : buchimag(GEN D, GEN c, GEN c2, GEN REL)
    1152           0 : { (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),0); }
    1153             : 
    1154             : GEN
    1155           0 : buchreal(GEN D, GEN flag, GEN c, GEN c2, GEN REL, long prec) {
    1156           0 :   if (signe(flag)) pari_err_IMPL("narrow class group");
    1157           0 :   (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),prec);
    1158             : }
    1159             : 
    1160             : GEN
    1161       16281 : quadclassunit0(GEN x, long flag, GEN data, long prec)
    1162             : {
    1163             :   long lx;
    1164       16281 :   double c1 = 0.0, c2 = 0.0;
    1165             : 
    1166       16281 :   if (!data) lx=1;
    1167             :   else
    1168             :   {
    1169          28 :     lx = lg(data);
    1170          28 :     if (typ(data)!=t_VEC) pari_err_TYPE("quadclassunit", data);
    1171          28 :     if (lx > 7) pari_err_DIM("quadclassunit [tech vector]");
    1172          28 :     if (lx > 3) lx = 3;
    1173             :   }
    1174       16281 :   switch(lx)
    1175             :   {
    1176          21 :     case 3: c2 = gtodouble(gel(data,2));
    1177          28 :     case 2: c1 = gtodouble(gel(data,1));
    1178             :   }
    1179       16281 :   if (flag) pari_err_IMPL("narrow class group");
    1180       16281 :   return Buchquad(x,c1,c2,prec);
    1181             : }

Generated by: LCOV version 1.13