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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 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.8.0 lcov report (development 16624-25b9976) Lines: 640 671 95.4 %
Date: 2014-06-24 Functions: 46 47 97.9 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 428 535 80.0 %

           Branch data     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                 :    1069528 : 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                 :     832486 : qfr3_canon(GEN x, struct qfr_data *S)
      56                 :            : {
      57                 :     832486 :   GEN a = gel(x,1), c = gel(x,3);
      58         [ +  + ]:     832486 :   if (signe(a) < 0) {
      59         [ +  + ]:     358694 :     if (absi_equal(a,c)) return qfr3_rho(x, S);
      60                 :     358690 :     setsigne(a, 1);
      61                 :     358690 :     setsigne(c,-1);
      62                 :            :   }
      63                 :     832486 :   return x;
      64                 :            : }
      65                 :            : static GEN
      66                 :       2655 : qfr3_canon_safe(GEN x, struct qfr_data *S)
      67                 :            : {
      68                 :       2655 :   GEN a = gel(x,1), c = gel(x,3);
      69         [ +  + ]:       2655 :   if (signe(a) < 0) {
      70         [ -  + ]:        136 :     if (absi_equal(a,c)) return qfr3_rho(x, S);
      71                 :        136 :     gel(x,1) = negi(a);
      72                 :        136 :     gel(x,3) = negi(c);
      73                 :            :   }
      74                 :       2655 :   return x;
      75                 :            : }
      76                 :            : static GEN
      77                 :    1392056 : qfr5_canon(GEN x, struct qfr_data *S)
      78                 :            : {
      79                 :    1392056 :   GEN a = gel(x,1), c = gel(x,3);
      80         [ +  + ]:    1392056 :   if (signe(a) < 0) {
      81         [ +  + ]:     620797 :     if (absi_equal(a,c)) return qfr5_rho(x, S);
      82                 :     620793 :     setsigne(a, 1);
      83                 :     620793 :     setsigne(c,-1);
      84                 :            :   }
      85                 :    1392056 :   return x;
      86                 :            : }
      87                 :            : static GEN
      88                 :    1237031 : QFR5_comp(GEN x,GEN y, struct qfr_data *S)
      89                 :    1237031 : { return qfr5_canon(qfr5_comp(x,y,S), S); }
      90                 :            : static GEN
      91                 :     718742 : QFR3_comp(GEN x, GEN y, struct qfr_data *S)
      92                 :     718742 : { return qfr3_canon(qfr3_comp(x,y,S), S); }
      93                 :            : 
      94                 :            : /* compute rho^n(x) */
      95                 :            : static GEN
      96                 :     157604 : qfr5_rho_pow(GEN x, long n, struct qfr_data *S)
      97                 :            : {
      98                 :            :   long i;
      99                 :     157604 :   pari_sp av = avma, lim = stack_lim(av, 1);
     100         [ +  + ]:    1583681 :   for (i=1; i<=n; i++)
     101                 :            :   {
     102                 :    1426077 :     x = qfr5_rho(x,S);
     103         [ -  + ]:    1426077 :     if (low_stack(lim, stack_lim(av,1)))
     104                 :            :     {
     105         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"qfr5_rho_pow");
     106                 :          0 :       x = gerepilecopy(av, x);
     107                 :            :     }
     108                 :            :   }
     109                 :     157604 :   return gerepilecopy(av, x);
     110                 :            : }
     111                 :            : 
     112                 :            : static GEN
     113                 :     155025 : qfr5_pf(struct qfr_data *S, long p, long prec)
     114                 :            : {
     115                 :     155025 :   GEN y = primeform_u(S->D,p);
     116                 :     155025 :   return qfr5_canon(qfr5_red(qfr_to_qfr5(y,prec), S), S);
     117                 :            : }
     118                 :            : 
     119                 :            : static GEN
     120                 :     113708 : qfr3_pf(struct qfr_data *S, long p)
     121                 :            : {
     122                 :     113708 :   GEN y = primeform_u(S->D,p);
     123                 :     113708 :   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                 :    1128869 : init_form(struct buch_quad *B, GEN ex,
     131                 :            :           GEN (*comp)(GEN,GEN,struct qfr_data *S))
     132                 :            : {
     133                 :    1128869 :   long i, l = lg(B->powsubFB);
     134                 :    1128869 :   GEN F = NULL;
     135         [ +  + ]:    8427527 :   for (i=1; i<l; i++)
     136         [ +  + ]:    7298658 :     if (ex[i])
     137                 :            :     {
     138                 :    6842677 :       GEN t = gmael(B->powsubFB,i,ex[i]);
     139         [ +  + ]:    6842677 :       F = F? comp(F, t, B->QFR): t;
     140                 :            :     }
     141                 :    1128869 :   return F;
     142                 :            : }
     143                 :            : static GEN
     144                 :     141030 : qfr5_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFR5_comp); }
     145                 :            : 
     146                 :            : static GEN
     147                 :    4368723 : QFI_comp(GEN x, GEN y, struct qfr_data *S) { (void)S; return qficomp(x,y); }
     148                 :            : 
     149                 :            : static GEN
     150                 :      27222 : qfi_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFI_comp); }
     151                 :            : 
     152                 :            : static GEN
     153                 :     960612 : random_form(struct buch_quad *B, GEN ex,
     154                 :            :             GEN (*comp)(GEN,GEN, struct qfr_data *S))
     155                 :            : {
     156                 :     960612 :   long i, l = lg(ex);
     157                 :     960612 :   pari_sp av = avma;
     158                 :            :   GEN F;
     159                 :            :   for(;;)
     160                 :            :   {
     161         [ +  + ]:    7144627 :     for (i=1; i<l; i++) ex[i] = random_bits(RANDOM_BITS);
     162         [ +  + ]:     960617 :     if ((F = init_form(B, ex, comp))) return F;
     163                 :          5 :     avma = av;
     164                 :          5 :   }
     165                 :            : }
     166                 :            : static GEN
     167                 :     115269 : qfr3_random(struct buch_quad *B,GEN ex){ return random_form(B, ex, &QFR3_comp); }
     168                 :            : static GEN
     169                 :     845343 : 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                 :        130 : check_LIMC(long LIMC, long LIMCMAX)
     178                 :            : {
     179         [ -  + ]:        130 :   if (LIMC >= LIMCMAX) pari_err_BUG("Buchmann's algorithm");
     180         [ +  + ]:        130 :   if (LIMC <= LIMCMAX/40) /* cbach <= 0.3 */
     181                 :          5 :     LIMC *= 2;
     182         [ -  + ]:        125 :   else if (LIMCMAX < 60) /* \Delta_K <= 9 */
     183                 :          0 :     LIMC++;
     184                 :            :   else
     185                 :        125 :     LIMC += LIMCMAX / 60; /* cbach += 0.2 */
     186         [ -  + ]:        130 :   if (LIMC > LIMCMAX) LIMC = LIMCMAX;
     187                 :        130 :   return LIMC;
     188                 :            : }
     189                 :            : 
     190                 :            : /* Is |q| <= p ? */
     191                 :            : static int
     192                 :      67575 : isless_iu(GEN q, ulong p) {
     193                 :      67575 :   long l = lgefint(q);
     194 [ +  - ][ +  + ]:      67575 :   return l==2 || (l == 3 && uel(q,2) <= p);
                 [ +  + ]
     195                 :            : }
     196                 :            : 
     197                 :            : static long
     198                 :    1826354 : factorquad(struct buch_quad *B, GEN f, long nFB, ulong limp)
     199                 :            : {
     200                 :            :   ulong X;
     201                 :    1826354 :   long i, lo = 0;
     202                 :    1826354 :   GEN x = gel(f,1), FB = B->FB, P = B->primfact, E = B->exprimfact;
     203                 :            : 
     204         [ +  + ]:    1893927 :   for (i=1; lgefint(x) > 3; i++)
     205                 :            :   {
     206                 :      67575 :     ulong p = uel(FB,i), r;
     207                 :      67575 :     GEN q = diviu_rem(x, p, &r);
     208         [ +  + ]:      67575 :     if (!r)
     209                 :            :     {
     210                 :       3617 :       long k = 0;
     211         [ +  + ]:       5759 :       do { k++; x = q; q = diviu_rem(x, p, &r); } while (!r);
     212                 :       3617 :       lo++; P[lo] = p; E[lo] = k;
     213                 :            :     }
     214         [ +  + ]:      67575 :     if (isless_iu(q,p)) {
     215         [ +  - ]:          2 :       if (lgefint(x) == 3) { X = uel(x,2); goto END; }
     216                 :          0 :       return 0;
     217                 :            :     }
     218         [ +  + ]:      67573 :     if (i == nFB) return 0;
     219                 :            :   }
     220                 :    1825993 :   X = uel(x,2);
     221         [ +  + ]:    1825993 :   if (X == 1) { P[0] = 0; return 1; }
     222                 :   82995834 :   for (;; i++)
     223                 :            :   { /* single precision affair, split for efficiency */
     224                 :   84821357 :     ulong p = uel(FB,i);
     225                 :   84821357 :     ulong q = X / p, r = X % p; /* gcc makes a single div */
     226         [ +  + ]:   84821357 :     if (!r)
     227                 :            :     {
     228                 :    3011037 :       long k = 0;
     229         [ +  + ]:    3646077 :       do { k++; X = q; q = X / p; r = X % p; } while (!r);
     230                 :    3011037 :       lo++; P[lo] = p; E[lo] = k;
     231                 :            :     }
     232         [ +  + ]:   84821357 :     if (q <= p) break;
     233         [ +  + ]:   83476508 :     if (i == nFB) return 0;
     234                 :   82995834 :   }
     235                 :            : END:
     236         [ +  + ]:    1344851 :   if (X > B->limhash) return 0;
     237 [ +  + ][ +  + ]:    1344833 :   if (X != 1 && X <= limp)
     238                 :            :   {
     239 [ +  + ][ +  + ]:     216542 :     if (B->badprim && ugcd(X, umodiu(B->badprim,X)) > 1) return 0;
     240                 :     216285 :     lo++; P[lo] = X; E[lo] = 1;
     241                 :     216285 :     X = 1;
     242                 :            :   }
     243                 :    1826354 :   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                 :    1069528 : largeprime(struct buch_quad *B, long q, GEN ex, long np, long nrho)
     249                 :            : {
     250                 :    1069528 :   const long hashv = hash(q);
     251                 :    1069528 :   long *pt, i, l = lg(B->subFB);
     252                 :            : 
     253                 :    1069528 :   for (pt = B->hashtab[hashv]; ; pt = (long*) pt[0])
     254                 :            :   {
     255         [ +  + ]:    1227232 :     if (!pt)
     256                 :            :     {
     257                 :    1013952 :       pt = (long*) pari_malloc((l+3) * sizeof(long));
     258                 :    1013952 :       *pt++ = nrho; /* nrho = pt[-3] */
     259                 :    1013952 :       *pt++ = np;   /* np   = pt[-2] */
     260                 :    1013952 :       *pt++ = q;    /* q    = pt[-1] */
     261                 :    1013952 :       pt[0] = (long)B->hashtab[hashv];
     262         [ +  + ]:    7638544 :       for (i=1; i<l; i++) pt[i]=ex[i];
     263                 :    1013952 :       B->hashtab[hashv]=pt; return NULL;
     264                 :            :     }
     265         [ +  + ]:     213280 :     if (pt[-1] == q) break;
     266                 :     157704 :   }
     267         [ +  + ]:      79645 :   for(i=1; i<l; i++)
     268         [ +  + ]:      77082 :     if (pt[i] != ex[i]) return pt;
     269         [ -  + ]:    1069528 :   return (pt[-2]==np)? NULL: pt;
     270                 :            : }
     271                 :            : 
     272                 :            : static void
     273                 :       9890 : clearhash(long **hash)
     274                 :            : {
     275                 :            :   long *pt;
     276                 :            :   long i;
     277         [ +  + ]:   10137250 :   for (i=0; i<HASHT; i++) {
     278         [ +  + ]:   11141312 :     for (pt = hash[i]; pt; ) {
     279                 :    1013952 :       void *z = (void*)(pt - 3);
     280                 :    1013952 :       pt = (long*) pt[0]; pari_free(z);
     281                 :            :     }
     282                 :   10127360 :     hash[i] = NULL;
     283                 :            :   }
     284                 :       9890 : }
     285                 :            : 
     286                 :            : /* last prime stored */
     287                 :            : ulong
     288                 :       4490 : 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                 :       5460 : GRH_ensure(GRHcheck_t *S, long nb)
     292                 :            : {
     293         [ -  + ]:       5460 :   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                 :       5460 : }
     300                 :            : /* cache data for all primes up to the LIM */
     301                 :            : static void
     302                 :      72040 : cache_prime_quad(GRHcheck_t *S, ulong LIM, GEN D)
     303                 :            : {
     304                 :            :   GRHprime_t *pr;
     305                 :            :   long nb;
     306                 :            : 
     307         [ +  + ]:      76440 :   if (S->limp >= LIM) return;
     308                 :       4400 :   nb = (long)primepi_upper_bound((double)LIM); /* #{p <= LIM} <= nb */
     309                 :       4400 :   GRH_ensure(S, nb+1); /* room for one extra prime */
     310                 :       4400 :   for (pr = S->primes + S->nprimes;;)
     311                 :            :   {
     312                 :     742135 :     ulong p = u_forprime_next(&(S->P));
     313                 :     742135 :     pr->p = p;
     314                 :     742135 :     pr->logp = log((double)p);
     315                 :     742135 :     pr->dec = (GEN)kroiu(D,p);
     316                 :     742135 :     S->nprimes++;
     317                 :     742135 :     pr++;
     318                 :            :     /* store up to nextprime(LIM) included */
     319         [ +  + ]:     742135 :     if (p >= LIM) { S->limp = p; break; }
     320                 :     737735 :   }
     321                 :            : }
     322                 :            : 
     323                 :            : static GEN
     324                 :       3365 : compute_invresquad(GRHcheck_t *S)
     325                 :            : {
     326                 :       3365 :   pari_sp av = avma;
     327                 :       3365 :   GEN invres = real_1(DEFAULTPREC);
     328                 :       3365 :   GRHprime_t *pr = S->primes;
     329                 :       3365 :   long i = S->nprimes, LIMC = GRH_last_prime(S)+diffptr[i]-1; /* nextprime(p+1)-1*/
     330                 :       3365 :   double limp = log(LIMC) / 2;
     331         [ +  + ]:     745500 :   for (; i > 0; pr++, i--)
     332                 :            :   {
     333                 :     742135 :     long s = (long)pr->dec;
     334         [ +  + ]:     742135 :     if (s)
     335                 :            :     {
     336                 :     735915 :       ulong p = pr->p;
     337 [ +  + ][ +  + ]:     735915 :       if (s>0 || pr->logp <= limp)
     338                 :            :         /* Both p and P contribute */
     339                 :     384815 :         invres = mulur(p - s, divru(invres, p));
     340         [ +  - ]:     351100 :       else if (s<0)
     341                 :            :         /* Only p contributes */
     342                 :     351100 :         invres = mulur(p, divru(invres, p - 1));
     343                 :            :     }
     344                 :            :   }
     345                 :       3365 :   return gerepileuptoleaf(av, invres);
     346                 :            : }
     347                 :            : 
     348                 :            : /* p | conductor of order of disc D ? */
     349                 :            : static int
     350                 :      18135 : is_bad(GEN D, ulong p)
     351                 :            : {
     352                 :      18135 :   pari_sp av = avma;
     353                 :            :   int r;
     354         [ +  + ]:      18135 :   if (p == 2)
     355                 :            :   {
     356                 :       4470 :     r = mod16(D) >> 1;
     357 [ +  + ][ +  + ]:       4470 :     if (r && signe(D) < 0) r = 8-r;
     358                 :       4470 :     return (r < 4);
     359                 :            :   }
     360                 :      13665 :   r = (remii(D, sqru(p)) == gen_0); /* p^2 | D ? */
     361                 :      18135 :   avma = av; return r;
     362                 :            : }
     363                 :            : 
     364                 :            : /* returns the n-th suitable ideal for the factorbase */
     365                 :            : static long
     366                 :       3365 : nthidealquad(GEN D, long n)
     367                 :            : {
     368                 :       3365 :   pari_sp av = avma;
     369                 :            :   forprime_t S;
     370                 :            :   ulong p;
     371                 :       3365 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
     372 [ +  - ][ +  + ]:      15460 :   while ((p = u_forprime_next(&S)) && n > 0)
     373 [ +  + ][ +  + ]:      12095 :     if (!is_bad(D, p) && kroiu(D, p) >= 0) n--;
     374                 :       3365 :   avma = av; return p;
     375                 :            : }
     376                 :            : 
     377                 :            : static int
     378                 :      65280 : quadGRHchk(GEN D, GRHcheck_t *S, ulong LIMC)
     379                 :            : {
     380                 :      65280 :   double logC = log((double)LIMC), SA = 0, SB = 0;
     381                 :            :   long i;
     382                 :            : 
     383                 :      65280 :   cache_prime_quad(S, LIMC, D);
     384                 :      65280 :   for (i = 0;; i++)
     385                 :            :   {
     386                 :    6246760 :     GRHprime_t *pr = S->primes+i;
     387                 :    6246760 :     ulong p = pr->p;
     388                 :            :     long M;
     389                 :            :     double logNP, q, A, B;
     390         [ +  + ]:    6246760 :     if (p > LIMC) break;
     391         [ +  + ]:    6181480 :     if ((long)pr->dec < 0)
     392                 :            :     {
     393                 :    3063145 :       logNP = 2 * pr->logp;
     394                 :    3063145 :       q = 1/(double)p;
     395                 :            :     }
     396                 :            :     else
     397                 :            :     {
     398                 :    3118335 :       logNP = pr->logp;
     399                 :    3118335 :       q = 1/sqrt((double)p);
     400                 :            :     }
     401                 :    6181480 :     A = logNP * q; B = logNP * A; M = (long)(logC/logNP);
     402         [ +  + ]:    6181480 :     if (M > 1)
     403                 :            :     {
     404                 :     285610 :       double inv1_q = 1 / (1-q);
     405                 :     285610 :       A *= (1 - pow(q, M)) * inv1_q;
     406                 :     285610 :       B *= (1 - pow(q, M)*(M+1 - M*q)) * inv1_q * inv1_q;
     407                 :            :     }
     408         [ +  + ]:    6181480 :     if ((long)pr->dec>0) { SA += 2*A;SB += 2*B; } else { SA += A; SB += B; }
     409                 :    6181480 :   }
     410                 :      65280 :   return GRHok(S, logC, SA, SB);
     411                 :            : }
     412                 :            : 
     413                 :            : /* create B->FB, B->numFB; set B->badprim. Return L(kro_D, 1) */
     414                 :            : static void
     415                 :       3395 : FBquad(struct buch_quad *B, ulong C2, ulong C1, GRHcheck_t *S)
     416                 :            : {
     417                 :       3395 :   GEN D = B->QFR->D;
     418                 :            :   long i;
     419                 :            :   pari_sp av;
     420                 :            :   GRHprime_t *pr;
     421                 :            : 
     422                 :       3395 :   cache_prime_quad(S, C2, D);
     423                 :       3395 :   pr = S->primes;
     424                 :       3395 :   B->numFB = cgetg(C2+1, t_VECSMALL);
     425                 :       3395 :   B->FB    = cgetg(C2+1, t_VECSMALL);
     426                 :       3395 :   av = avma;
     427                 :       3395 :   B->KC = 0; i = 0;
     428                 :       3395 :   B->badprim = gen_1;
     429                 :     494055 :   for (;; pr++) /* p <= C2 */
     430                 :            :   {
     431                 :     497450 :     ulong p = pr->p;
     432 [ +  + ][ +  + ]:     497450 :     if (!B->KC && p > C1) B->KC = i;
     433         [ +  + ]:     497450 :     if (p > C2) break;
     434      [ +  +  + ]:     494055 :     switch ((long)pr->dec)
     435                 :            :     {
     436                 :     244625 :       case -1: break; /* inert */
     437                 :            :       case  0: /* ramified */
     438         [ +  + ]:       6040 :         if (is_bad(D, p)) { B->badprim = muliu(B->badprim, p); break; }
     439                 :            :         /* fall through */
     440                 :            :       default:  /* split */
     441                 :     249405 :         i++; B->numFB[p] = i; B->FB[i] = p; break;
     442                 :            :     }
     443                 :     494055 :   }
     444                 :       3395 :   B->KC2 = i;
     445                 :       3395 :   setlg(B->FB, B->KC2+1);
     446         [ +  + ]:       3395 :   if (B->badprim != gen_1)
     447                 :         20 :     B->badprim = gerepileuptoint(av, B->badprim);
     448                 :            :   else
     449                 :            :   {
     450                 :       3375 :     B->badprim = NULL; avma = av;
     451                 :            :   }
     452                 :       3395 : }
     453                 :            : 
     454                 :            : /* create B->vperm, return B->subFB */
     455                 :            : static GEN
     456                 :       3395 : subFBquad(struct buch_quad *B, GEN D, double PROD, long minSFB)
     457                 :            : {
     458                 :       3395 :   long i, j, lgsub = 1, ino = 1, lv = B->KC+1;
     459                 :       3395 :   double prod = 1.;
     460                 :            :   pari_sp av;
     461                 :            :   GEN no;
     462                 :            : 
     463                 :       3395 :   B->vperm = cgetg(lv, t_VECSMALL);
     464                 :       3395 :   av = avma;
     465                 :       3395 :   no    = cgetg(lv, t_VECSMALL);
     466         [ +  + ]:      25710 :   for (j = 1; j < lv; j++)
     467                 :            :   {
     468                 :      25680 :     ulong p = uel(B->FB,j);
     469         [ +  + ]:      25680 :     if (!umodiu(D, p)) no[ino++] = j; /* ramified */
     470                 :            :     else
     471                 :            :     {
     472                 :      21340 :       B->vperm[lgsub++] = j;
     473                 :      21340 :       prod *= p;
     474 [ +  + ][ +  + ]:      21340 :       if (lgsub > minSFB && prod > PROD) break;
     475                 :            :     }
     476                 :            :   }
     477                 :            :   /* lgsub >= 1 otherwise quadGRHchk is false */
     478                 :       3395 :   i = lgsub;
     479         [ +  + ]:       7735 :   for (j = 1; j < ino;i++,j++) B->vperm[i] = no[j];
     480         [ +  + ]:     226730 :   for (     ; i < lv; i++)     B->vperm[i] = i;
     481                 :       3395 :   no = gclone(vecslice(B->vperm, 1, lgsub-1));
     482                 :       3395 :   avma = av; return no;
     483                 :            : }
     484                 :            : 
     485                 :            : /* assume n >= 1, x[i][j] = B->subFB[i]^j, for j = 1..n */
     486                 :            : static GEN
     487                 :       6495 : powsubFBquad(struct buch_quad *B, long n)
     488                 :            : {
     489                 :       6495 :   pari_sp av = avma;
     490                 :       6495 :   long i,j, l = lg(B->subFB);
     491                 :       6495 :   GEN F, y, x = cgetg(l, t_VEC), D = B->QFR->D;
     492                 :            : 
     493         [ +  + ]:       6495 :   if (B->PRECREG) /* real */
     494                 :            :   {
     495         [ +  + ]:      28040 :     for (i=1; i<l; i++)
     496                 :            :     {
     497                 :      24425 :       F = qfr5_pf(B->QFR, B->FB[B->subFB[i]], B->PRECREG);
     498                 :      24425 :       y = cgetg(n+1, t_VEC); gel(x,i) = y;
     499                 :      24425 :       gel(y,1) = F;
     500         [ +  + ]:     390800 :       for (j=2; j<=n; j++) gel(y,j) = QFR5_comp(gel(y,j-1), F, B->QFR);
     501                 :            :     }
     502                 :            :   }
     503                 :            :   else /* imaginary */
     504                 :            :   {
     505         [ +  + ]:      20378 :     for (i=1; i<l; i++)
     506                 :            :     {
     507                 :      17498 :       F = qfi_pf(D, B->FB[B->subFB[i]]);
     508                 :      17498 :       y = cgetg(n+1, t_VEC); gel(x,i) = y;
     509                 :      17498 :       gel(y,1) = F;
     510         [ +  + ]:     279968 :       for (j=2; j<=n; j++) gel(y,j) = qficomp(gel(y,j-1), F);
     511                 :            :     }
     512                 :            :   }
     513                 :       6495 :   x = gclone(x); avma = av; return x;
     514                 :            : }
     515                 :            : 
     516                 :            : static void
     517                 :      26543 : sub_fact(struct buch_quad *B, GEN col, GEN F)
     518                 :            : {
     519                 :      26543 :   GEN b = gel(F,2);
     520                 :            :   long i;
     521         [ +  + ]:      86927 :   for (i=1; i<=B->primfact[0]; i++)
     522                 :            :   {
     523                 :      60384 :     ulong p = B->primfact[i], k = B->numFB[p];
     524                 :      60384 :     long e = B->exprimfact[i];
     525         [ +  + ]:      60384 :     if (umodiu(b, p<<1) > p) e = -e;
     526                 :      60384 :     col[k] -= e;
     527                 :            :   }
     528                 :      26543 : }
     529                 :            : static void
     530                 :     301525 : add_fact(struct buch_quad *B, GEN col, GEN F)
     531                 :            : {
     532                 :     301525 :   GEN b = gel(F,2);
     533                 :            :   long i;
     534         [ +  + ]:    1331880 :   for (i=1; i<=B->primfact[0]; i++)
     535                 :            :   {
     536                 :    1030355 :     ulong p = B->primfact[i], k = B->numFB[p];
     537                 :    1030355 :     long e = B->exprimfact[i];
     538         [ +  + ]:    1030355 :     if (umodiu(b, p<<1) > p) e = -e;
     539                 :    1030355 :     col[k] += e;
     540                 :            :   }
     541                 :     301525 : }
     542                 :            : 
     543                 :            : static GEN
     544                 :       3365 : get_clgp(struct buch_quad *B, GEN W, GEN *ptD, long prec)
     545                 :            : {
     546         [ +  + ]:       3365 :   GEN res, init, u1, D = ZM_snf_group(W,NULL,&u1), Z = prec? real_0(prec): NULL;
     547                 :       3365 :   long i, j, l = lg(W), c = lg(D);
     548                 :            : 
     549                 :       3365 :   res=cgetg(c,t_VEC); init = cgetg(l,t_VEC);
     550         [ +  + ]:      11099 :   for (i=1; i<l; i++) gel(init,i) = primeform_u(B->QFR->D, B->FB[B->vperm[i]]);
     551         [ +  + ]:      10335 :   for (j=1; j<c; j++)
     552                 :            :   {
     553                 :       6970 :     GEN g = NULL;
     554         [ +  + ]:       6970 :     if (prec)
     555                 :            :     {
     556         [ +  + ]:       9508 :       for (i=1; i<l; i++)
     557                 :            :       {
     558                 :       6853 :         GEN t, u = gcoeff(u1,i,j);
     559         [ +  + ]:       6853 :         if (!signe(u)) continue;
     560                 :       3140 :         t = qfr3_pow(gel(init,i), u, B->QFR);
     561         [ +  + ]:       3140 :         g = g? qfr3_comp(g, t, B->QFR): t;
     562                 :            :       }
     563                 :       2655 :       g = qfr3_to_qfr(qfr3_canon_safe(qfr3_red(g, B->QFR), B->QFR), Z);
     564                 :            :     }
     565                 :            :     else
     566                 :            :     {
     567         [ +  + ]:      18590 :       for (i=1; i<l; i++)
     568                 :            :       {
     569                 :      14275 :         GEN t, u = gcoeff(u1,i,j);
     570         [ +  + ]:      14275 :         if (!signe(u)) continue;
     571                 :       8289 :         t = powgi(gel(init,i), u);
     572         [ +  + ]:       8289 :         g = g? qficomp(g, t): t;
     573                 :            :       }
     574                 :            :     }
     575                 :       6970 :     gel(res,j) = g;
     576                 :            :   }
     577                 :       3365 :   *ptD = D; return res;
     578                 :            : }
     579                 :            : 
     580                 :            : static long
     581                 :       3365 : trivial_relations(struct buch_quad *B, GEN mat, GEN C)
     582                 :            : {
     583                 :       3365 :   long i, j = 0;
     584                 :       3365 :   GEN col, D = B->QFR->D;
     585         [ +  + ]:     252270 :   for (i = 1; i <= B->KC; i++)
     586                 :            :   { /* ramified prime ==> trivial relation */
     587         [ +  + ]:     248905 :     if (umodiu(D, B->FB[i])) continue;
     588                 :       5920 :     col = zero_zv(B->KC);
     589                 :       5920 :     col[i] = 2; j++;
     590                 :       5920 :     gel(mat,j) = col;
     591                 :       5920 :     gel(C,j) = gen_0;
     592                 :            :   }
     593                 :       3365 :   return j;
     594                 :            : }
     595                 :            : 
     596                 :            : static void
     597                 :          0 : dbg_all(pari_timer *T, const char *phase, long s, long n)
     598                 :            : {
     599                 :          0 :   err_printf("\n");
     600                 :          0 :   timer_printf(T, "%s rel [#rel/#test = %ld/%ld]", phase,s,n);
     601                 :          0 : }
     602                 :            : 
     603                 :            : /* Imaginary Quadratic fields */
     604                 :            : 
     605                 :            : static void
     606                 :       2958 : imag_relations(struct buch_quad *B, long need, long *pc, ulong LIMC, GEN mat)
     607                 :            : {
     608                 :            :   pari_timer T;
     609                 :       2958 :   long lgsub = lg(B->subFB), current = *pc, nbtest = 0, s = 0;
     610                 :            :   long i, fpc;
     611                 :            :   pari_sp av;
     612                 :       2958 :   GEN col, form, ex = cgetg(lgsub, t_VECSMALL);
     613                 :            : 
     614         [ +  + ]:       2958 :   if (!current) current = 1;
     615         [ -  + ]:       2958 :   if (DEBUGLEVEL>2) timer_start(&T);
     616                 :       2958 :   av = avma;
     617                 :            :   for(;;)
     618                 :            :   {
     619         [ +  + ]:     847932 :     if (s >= need) break;
     620                 :     844974 :     avma = av;
     621                 :     844974 :     form = qfi_random(B,ex);
     622                 :     844974 :     form = qficomp(form, qfi_pf(B->QFR->D, B->FB[current]));
     623                 :     844974 :     nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
     624         [ +  + ]:     844974 :     if (!fpc)
     625                 :            :     {
     626         [ -  + ]:     205568 :       if (DEBUGLEVEL>3) err_printf(".");
     627 [ +  + ][ +  + ]:     205568 :       if ((nbtest & 0xff) == 0 && ++current > B->KC) current = 1;
     628                 :     205568 :       continue;
     629                 :            :     }
     630         [ +  + ]:     639406 :     if (fpc > 1)
     631                 :            :     {
     632                 :     523330 :       long *fpd = largeprime(B,fpc,ex,current,0);
     633                 :            :       ulong b1, b2, p;
     634                 :            :       GEN form2;
     635         [ +  + ]:     523330 :       if (!fpd)
     636                 :            :       {
     637         [ -  + ]:     496108 :         if (DEBUGLEVEL>3) err_printf(".");
     638                 :     496108 :         continue;
     639                 :            :       }
     640                 :      27222 :       form2 = qficomp(qfi_factorback(B,fpd), qfi_pf(B->QFR->D, B->FB[fpd[-2]]));
     641                 :      27222 :       p = fpc << 1;
     642                 :      27222 :       b1 = umodiu(gel(form2,2), p);
     643                 :      27222 :       b2 = umodiu(gel(form,2),  p);
     644 [ +  + ][ +  + ]:      27222 :       if (b1 != b2 && b1+b2 != p) continue;
     645                 :            : 
     646                 :      27194 :       col = gel(mat,++s);
     647                 :      27194 :       add_fact(B,col, form);
     648                 :      27194 :       (void)factorquad(B,form2,B->KC,LIMC);
     649         [ +  + ]:      27194 :       if (b1==b2)
     650                 :            :       {
     651         [ +  + ]:     100930 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
     652                 :      13585 :         sub_fact(B, col, form2); col[fpd[-2]]++;
     653                 :            :       }
     654                 :            :       else
     655                 :            :       {
     656         [ +  + ]:     101268 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
     657                 :      13609 :         add_fact(B, col, form2); col[fpd[-2]]--;
     658                 :            :       }
     659         [ -  + ]:      27194 :       if (DEBUGLEVEL>2) err_printf(" %ldP",s);
     660                 :            :     }
     661                 :            :     else
     662                 :            :     {
     663                 :     116076 :       col = gel(mat,++s);
     664         [ +  + ]:     860088 :       for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
     665                 :     116076 :       add_fact(B, col, form);
     666         [ -  + ]:     116076 :       if (DEBUGLEVEL>2) err_printf(" %ld",s);
     667                 :            :     }
     668                 :     143270 :     col[current]--;
     669         [ +  + ]:     143270 :     if (++current > B->KC) current = 1;
     670                 :     844974 :   }
     671         [ -  + ]:       2958 :   if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
     672                 :       2958 :   *pc = current;
     673                 :       2958 : }
     674                 :            : 
     675                 :            : static int
     676                 :          5 : imag_be_honest(struct buch_quad *B)
     677                 :            : {
     678                 :          5 :   long p, fpc, s = B->KC, nbtest = 0;
     679                 :          5 :   GEN F, ex = cgetg(lg(B->subFB), t_VECSMALL);
     680                 :          5 :   pari_sp av = avma;
     681                 :            : 
     682         [ +  + ]:        374 :   while (s<B->KC2)
     683                 :            :   {
     684         [ -  + ]:        369 :     p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
     685                 :        369 :     F = qficomp(qfi_pf(B->QFR->D, p), qfi_random(B, ex));
     686                 :        369 :     fpc = factorquad(B,F,s,p-1);
     687         [ +  + ]:        369 :     if (fpc == 1) { nbtest=0; s++; }
     688                 :            :     else
     689         [ -  + ]:        264 :       if (++nbtest > 40) return 0;
     690                 :        369 :     avma = av;
     691                 :            :   }
     692                 :          5 :   return 1;
     693                 :            : }
     694                 :            : 
     695                 :            : /* Real Quadratic fields */
     696                 :            : 
     697                 :            : static void
     698                 :       3677 : real_relations(struct buch_quad *B, long need, long *pc, long lim, ulong LIMC, GEN mat, GEN C)
     699                 :            : {
     700                 :            :   pari_timer T;
     701                 :       3677 :   long lgsub = lg(B->subFB), prec = B->PRECREG, current = *pc, nbtest=0, s=0;
     702                 :            :   long i, fpc, endcycle, rhoacc, rho;
     703                 :            :   /* in a 2nd phase, don't include FB[current] but run along the cyle
     704                 :            :    * ==> get more units */
     705                 :       3677 :   int first = (current == 0);
     706                 :            :   pari_sp av, av1, limstack;
     707                 :       3677 :   GEN d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL);
     708                 :            : 
     709         [ -  + ]:       3677 :   if (DEBUGLEVEL>2) timer_start(&T);
     710         [ +  + ]:       3677 :   if (!current) current = 1;
     711         [ +  + ]:       3677 :   if (lim > need) lim = need;
     712                 :       3677 :   av = avma; limstack = stack_lim(av,1);
     713                 :            :   for(;;)
     714                 :            :   {
     715         [ +  + ]:     118916 :     if (s >= need) break;
     716 [ +  + ][ +  + ]:     115239 :     if (first && s >= lim) {
     717                 :       1535 :       first = 0;
     718         [ -  + ]:       1535 :       if (DEBUGLEVEL>2) dbg_all(&T, "initial", s, nbtest);
     719                 :            :     }
     720                 :     115239 :     avma = av; form = qfr3_random(B, ex);
     721         [ +  + ]:     115239 :     if (!first)
     722                 :     113678 :       form = QFR3_comp(form, qfr3_pf(B->QFR, B->FB[current]), B->QFR);
     723                 :     115239 :     av1 = avma;
     724                 :     115239 :     form0 = form; form1 = NULL;
     725                 :     115239 :     endcycle = rhoacc = 0;
     726                 :     115239 :     rho = -1;
     727                 :            : 
     728                 :            : CYCLE:
     729 [ +  + ][ -  + ]:     927976 :     if (endcycle || rho > 5000)
     730                 :            :     {
     731         [ -  + ]:         16 :       if (++current > B->KC) current = 1;
     732                 :         16 :       continue;
     733                 :            :     }
     734         [ -  + ]:     927960 :     if (low_stack(limstack, stack_lim(av,1)))
     735                 :            :     {
     736         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"real_relations");
     737         [ #  # ]:          0 :       gerepileall(av1, form1? 2: 1, &form, &form1);
     738                 :            :     }
     739         [ +  + ]:     927960 :     if (rho < 0) rho = 0; /* first time in */
     740                 :            :     else
     741                 :            :     {
     742                 :     812721 :       form = qfr3_rho(form, B->QFR); rho++;
     743                 :     812721 :       rhoacc++;
     744         [ +  + ]:     812721 :       if (first)
     745         [ +  + ]:     286014 :         endcycle = (absi_equal(gel(form,1),gel(form0,1))
     746         [ +  + ]:     127117 :              && equalii(gel(form,2),gel(form0,2)));
           [ +  -  +  - ]
     747                 :            :       else
     748                 :            :       {
     749         [ -  + ]:     653824 :         if (absi_equal(gel(form,1), gel(form,3))) /* a = -c */
     750                 :            :         {
     751   [ #  #  #  # ]:          0 :           if (absi_equal(gel(form,1),gel(form0,1)) &&
     752                 :          0 :                   equalii(gel(form,2),gel(form0,2))) continue;
     753                 :          0 :           form = qfr3_rho(form, B->QFR); rho++;
     754                 :          0 :           rhoacc++;
     755                 :            :         }
     756                 :            :         else
     757                 :     653824 :           { setsigne(form[1],1); setsigne(form[3],-1); }
     758   [ +  +  -  + ]:     653860 :         if (equalii(gel(form,1),gel(form0,1)) &&
     759                 :         36 :             equalii(gel(form,2),gel(form0,2))) continue;
     760                 :            :       }
     761                 :            :     }
     762                 :     927960 :     nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
     763         [ +  + ]:     927960 :     if (!fpc)
     764                 :            :     {
     765         [ -  + ]:     275740 :       if (DEBUGLEVEL>3) err_printf(".");
     766                 :     275740 :       goto CYCLE;
     767                 :            :     }
     768         [ +  + ]:     652220 :     if (fpc > 1)
     769                 :            :     { /* look for Large Prime relation */
     770         [ +  + ]:     546198 :       long *fpd = largeprime(B,fpc,ex,first? 0: current,rhoacc);
     771                 :            :       ulong b1, b2, p;
     772                 :            :       GEN form2;
     773         [ +  + ]:     546198 :       if (!fpd)
     774                 :            :       {
     775         [ -  + ]:     520407 :         if (DEBUGLEVEL>3) err_printf(".");
     776                 :     520407 :         goto CYCLE;
     777                 :            :       }
     778         [ +  - ]:      25791 :       if (!form1)
     779                 :            :       {
     780                 :      25791 :         form1 = qfr5_factorback(B,ex);
     781         [ +  - ]:      25791 :         if (!first)
     782                 :      25791 :           form1 = QFR5_comp(form1, qfr5_pf(B->QFR, B->FB[current], prec), B->QFR);
     783                 :            :       }
     784                 :      25791 :       form1 = qfr5_rho_pow(form1, rho, B->QFR);
     785                 :      25791 :       rho = 0;
     786                 :            : 
     787                 :      25791 :       form2 = qfr5_factorback(B,fpd);
     788         [ +  + ]:      25791 :       if (fpd[-2])
     789                 :      16922 :         form2 = QFR5_comp(form2, qfr5_pf(B->QFR, B->FB[fpd[-2]], prec), B->QFR);
     790                 :      25791 :       form2 = qfr5_rho_pow(form2, fpd[-3], B->QFR);
     791         [ +  - ]:      25791 :       if (!absi_equal(gel(form2,1),gel(form2,3)))
     792                 :            :       {
     793                 :      25791 :         setsigne(form2[1], 1);
     794                 :      25791 :         setsigne(form2[3],-1);
     795                 :            :       }
     796                 :      25791 :       p = fpc << 1;
     797                 :      25791 :       b1 = umodiu(gel(form2,2), p);
     798                 :      25791 :       b2 = umodiu(gel(form1,2), p);
     799 [ +  + ][ -  + ]:      25791 :       if (b1 != b2 && b1+b2 != p) goto CYCLE;
     800                 :            : 
     801                 :      25791 :       col = gel(mat,++s);
     802                 :      25791 :       add_fact(B, col, form1);
     803                 :      25791 :       (void)factorquad(B,form2,B->KC,LIMC);
     804         [ +  + ]:      25791 :       if (b1==b2)
     805                 :            :       {
     806         [ +  + ]:      99019 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
     807                 :      12958 :         sub_fact(B,col, form2);
     808         [ +  + ]:      12958 :         if (fpd[-2]) col[fpd[-2]]++;
     809                 :      12958 :         d = qfr5_dist(subii(gel(form1,4),gel(form2,4)),
     810                 :      25916 :                       divrr(gel(form1,5),gel(form2,5)), prec);
     811                 :            :       }
     812                 :            :       else
     813                 :            :       {
     814         [ +  + ]:      98324 :         for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
     815                 :      12833 :         add_fact(B, col, form2);
     816         [ +  + ]:      12833 :         if (fpd[-2]) col[fpd[-2]]--;
     817                 :      12833 :         d = qfr5_dist(addii(gel(form1,4),gel(form2,4)),
     818                 :      25666 :                       mulrr(gel(form1,5),gel(form2,5)), prec);
     819                 :            :       }
     820         [ -  + ]:      25791 :       if (DEBUGLEVEL>2) err_printf(" %ldP",s);
     821                 :            :     }
     822                 :            :     else
     823                 :            :     { /* standard relation */
     824         [ +  + ]:     106022 :       if (!form1)
     825                 :            :       {
     826                 :      89448 :         form1 = qfr5_factorback(B, ex);
     827         [ +  + ]:      89448 :         if (!first)
     828                 :      87887 :           form1 = QFR5_comp(form1, qfr5_pf(B->QFR, B->FB[current], prec), B->QFR);
     829                 :            :       }
     830                 :     106022 :       form1 = qfr5_rho_pow(form1, rho, B->QFR);
     831                 :     106022 :       rho = 0;
     832                 :            : 
     833                 :     106022 :       col = gel(mat,++s);
     834         [ +  + ]:     815584 :       for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
     835                 :     106022 :       add_fact(B, col, form1);
     836                 :     106022 :       d = qfr5_dist(gel(form1,4), gel(form1,5), prec);
     837         [ -  + ]:     106022 :       if (DEBUGLEVEL>2) err_printf(" %ld",s);
     838                 :            :     }
     839                 :     131813 :     affrr(d, gel(C,s));
     840         [ +  + ]:     131813 :     if (first)
     841                 :            :     {
     842         [ +  + ]:      18135 :       if (s >= lim) continue;
     843                 :      16590 :       goto CYCLE;
     844                 :            :     }
     845                 :            :     else
     846                 :            :     {
     847                 :     113678 :       col[current]--;
     848         [ +  + ]:     113678 :       if (++current > B->KC) current = 1;
     849                 :            :     }
     850                 :     115239 :   }
     851         [ -  + ]:       3677 :   if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
     852                 :       3677 :   *pc = current;
     853                 :       3677 : }
     854                 :            : 
     855                 :            : static int
     856                 :         10 : real_be_honest(struct buch_quad *B)
     857                 :            : {
     858                 :         10 :   long p, fpc, s = B->KC, nbtest = 0;
     859                 :         10 :   GEN F,F0, ex = cgetg(lg(B->subFB), t_VECSMALL);
     860                 :         10 :   pari_sp av = avma;
     861                 :            : 
     862         [ +  + ]:         40 :   while (s<B->KC2)
     863                 :            :   {
     864         [ -  + ]:         30 :     p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
     865                 :         30 :     F = QFR3_comp(qfr3_random(B, ex), qfr3_pf(B->QFR, p), B->QFR);
     866                 :         30 :     for (F0 = F;;)
     867                 :            :     {
     868                 :         66 :       fpc = factorquad(B,F,s,p-1);
     869         [ +  + ]:         66 :       if (fpc == 1) { nbtest=0; s++; break; }
     870         [ -  + ]:         36 :       if (++nbtest > 40) return 0;
     871                 :         36 :       F = qfr3_canon(qfr3_rho(F, B->QFR), B->QFR);
     872         [ -  + ]:         36 :       if (equalii(gel(F,1),gel(F0,1))
     873         [ #  # ]:          0 :        && equalii(gel(F,2),gel(F0,2))) break;
     874                 :         36 :     }
     875                 :         30 :     avma = av;
     876                 :            :   }
     877                 :         10 :   return 1;
     878                 :            : }
     879                 :            : 
     880                 :            : static GEN
     881                 :      23120 : gcdreal(GEN a,GEN b)
     882                 :            : {
     883         [ +  + ]:      23120 :   if (!signe(a)) return mpabs(b);
     884         [ +  + ]:      22591 :   if (!signe(b)) return mpabs(a);
     885         [ -  + ]:      22488 :   if (typ(a)==t_INT)
     886                 :            :   {
     887         [ #  # ]:          0 :     if (typ(b)==t_INT) return gcdii(a,b);
     888                 :          0 :     a = itor(a, lg(b));
     889                 :            :   }
     890         [ -  + ]:      22488 :   else if (typ(b)==t_INT)
     891                 :          0 :     b = itor(b, lg(a));
     892         [ +  + ]:      22488 :   if (expo(a)<-5) return absr(b);
     893         [ +  + ]:       8192 :   if (expo(b)<-5) return absr(a);
     894                 :       6140 :   a = absr(a); b = absr(b);
     895 [ +  + ][ +  - ]:      13589 :   while (expo(b) >= -5  && signe(b))
     896                 :            :   {
     897                 :            :     long e;
     898                 :       7449 :     GEN r, q = gcvtoi(divrr(a,b),&e);
     899         [ -  + ]:       7449 :     if (e > 0) return NULL;
     900                 :       7449 :     r = subrr(a, mulir(q,b)); a = b; b = r;
     901                 :            :   }
     902                 :      23120 :   return absr(a);
     903                 :            : }
     904                 :            : 
     905                 :            : static int
     906                 :       4511 : get_R(struct buch_quad *B, GEN C, long sreg, GEN z, GEN *ptR)
     907                 :            : {
     908                 :       4511 :   GEN R = gen_1;
     909                 :            :   double c;
     910                 :            :   long i;
     911                 :            : 
     912         [ +  + ]:       4511 :   if (B->PRECREG)
     913                 :            :   {
     914                 :       2076 :     R = mpabs(gel(C,1));
     915         [ +  + ]:      25196 :     for (i=2; i<=sreg; i++)
     916                 :            :     {
     917                 :      23120 :       R = gcdreal(gel(C,i), R);
     918         [ -  + ]:      23120 :       if (!R) return fupb_PRECI;
     919                 :            :     }
     920         [ -  + ]:       2076 :     if (gexpo(R) <= -3)
     921                 :            :     {
     922         [ #  # ]:          0 :       if (DEBUGLEVEL>2) err_printf("regulator is zero.\n");
     923                 :          0 :       return fupb_RELAT;
     924                 :            :     }
     925         [ -  + ]:       2076 :     if (DEBUGLEVEL>2) err_printf("#### Tentative regulator: %Ps\n",R);
     926                 :            :   }
     927                 :       4511 :   c = gtodouble(gmul(z, R));
     928 [ +  - ][ +  + ]:       4511 :   if (c < 0.8 || c > 1.3) return fupb_RELAT;
     929                 :       4511 :   *ptR = R; return fupb_NONE;
     930                 :            : }
     931                 :            : 
     932                 :            : static int
     933                 :       3365 : quad_be_honest(struct buch_quad *B)
     934                 :            : {
     935                 :            :   int r;
     936         [ +  + ]:       3365 :   if (B->KC2 <= B->KC) return 1;
     937         [ -  + ]:         15 :   if (DEBUGLEVEL>2)
     938                 :          0 :     err_printf("be honest for primes from %ld to %ld\n", B->FB[B->KC+1],B->FB[B->KC2]);
     939         [ +  + ]:         15 :   r = B->PRECREG? real_be_honest(B): imag_be_honest(B);
     940         [ -  + ]:         15 :   if (DEBUGLEVEL>2) err_printf("\n");
     941                 :       3365 :   return r;
     942                 :            : }
     943                 :            : 
     944                 :            : 
     945                 :            : GEN
     946                 :       3370 : Buchquad(GEN D, double cbach, double cbach2, long prec)
     947                 :            : {
     948                 :       3370 :   const long MAXRELSUP = 7, SFB_MAX = 3;
     949                 :            :   pari_timer T;
     950                 :       3370 :   pari_sp av0 = avma, av, av2;
     951                 :       3370 :   const long RELSUP = 5;
     952                 :            :   long i, s, current, triv, sfb_trials, nrelsup, nreldep, need, nsubFB, minSFB;
     953                 :            :   ulong low, high, LIMC0, LIMC, LIMC2, LIMCMAX, cp;
     954                 :       3370 :   GEN W, cyc, res, gen, dep, mat, C, extraC, B, R, invhr, h = NULL; /*-Wall*/
     955                 :            :   double drc, sdrc, lim, LOGD, LOGD2;
     956                 :            :   GRHcheck_t GRHcheck;
     957                 :            :   struct qfr_data QFR;
     958                 :            :   struct buch_quad BQ;
     959                 :       3370 :   int FIRST = 1;
     960                 :            : 
     961                 :       3370 :   check_quaddisc(D, &s, /*junk*/&i, "Buchquad");
     962                 :       3370 :   R = NULL; /* -Wall */
     963                 :       3370 :   BQ.QFR = &QFR;
     964                 :       3370 :   QFR.D = D;
     965         [ +  + ]:       3370 :   if (s < 0)
     966                 :            :   {
     967         [ +  + ]:       1825 :     if (cmpiu(QFR.D,4) <= 0)
     968                 :            :     {
     969                 :          5 :       GEN z = cgetg(5,t_VEC);
     970                 :          5 :       gel(z,1) = gel(z,4) = gen_1; gel(z,2) = gel(z,3) = cgetg(1,t_VEC);
     971                 :          5 :       return z;
     972                 :            :     }
     973                 :       1820 :     BQ.PRECREG = 0;
     974                 :            :   } else {
     975                 :       1545 :     BQ.PRECREG = maxss(prec+EXTRAPRECWORD, nbits2prec(2*expi(QFR.D) + 128));
     976                 :            :   }
     977         [ -  + ]:       3365 :   if (DEBUGLEVEL>2) timer_start(&T);
     978                 :       3365 :   BQ.primfact   = new_chunk(100);
     979                 :       3365 :   BQ.exprimfact = new_chunk(100);
     980                 :       3365 :   BQ.hashtab = (long**) new_chunk(HASHT);
     981         [ +  + ]:    3449125 :   for (i=0; i<HASHT; i++) BQ.hashtab[i] = NULL;
     982                 :            : 
     983                 :       3365 :   drc = fabs(gtodouble(QFR.D));
     984                 :       3365 :   LOGD = log(drc);
     985                 :       3365 :   LOGD2 = LOGD * LOGD;
     986                 :            : 
     987                 :       3365 :   sdrc = lim = sqrt(drc);
     988         [ +  + ]:       3365 :   if (!BQ.PRECREG) lim /= sqrt(3.);
     989                 :       3365 :   cp = (ulong)exp(sqrt(LOGD*log(LOGD)/8.0));
     990         [ +  + ]:       3365 :   if (cp < 20) cp = 20;
     991         [ -  + ]:       3365 :   if (cbach > 6.) {
     992         [ #  # ]:          0 :     if (cbach2 < cbach) cbach2 = cbach;
     993                 :          0 :     cbach = 6.;
     994                 :            :   }
     995         [ -  + ]:       3365 :   if (cbach < 0.)
     996                 :          0 :     pari_err_DOMAIN("Buchquad","Bach constant","<",gen_0,dbltor(cbach));
     997                 :       3365 :   av = avma;
     998                 :       3365 :   BQ.powsubFB = BQ.subFB = NULL;
     999         [ +  + ]:       3365 :   minSFB = (expi(D) > 15)? 3: 2;
    1000         [ +  + ]:       3365 :   init_GRHcheck(&GRHcheck, 2, BQ.PRECREG? 2: 0, LOGD);
    1001                 :       3365 :   high = low = LIMC0 = maxss((long)(cbach2*LOGD2), 1);
    1002                 :       3365 :   LIMCMAX = (long)(6.*LOGD2);
    1003                 :            :   /* 97/1223 below to ensure a good enough approximation of residue */
    1004         [ +  + ]:       3365 :   cache_prime_quad(&GRHcheck, expi(D) < 16 ? 97: 1223, D);
    1005         [ +  + ]:      35960 :   while (!quadGRHchk(D, &GRHcheck, high))
    1006                 :            :   {
    1007                 :      32595 :     low = high;
    1008                 :      32595 :     high *= 2;
    1009                 :            :   }
    1010         [ +  + ]:      32685 :   while (high - low > 1)
    1011                 :            :   {
    1012                 :      29320 :     long test = (low+high)/2;
    1013         [ +  + ]:      29320 :     if (quadGRHchk(D, &GRHcheck, test))
    1014                 :      15120 :       high = test;
    1015                 :            :     else
    1016                 :      14200 :       low = test;
    1017                 :            :   }
    1018 [ -  + ][ #  # ]:       3365 :   if (high == LIMC0+1 && quadGRHchk(D, &GRHcheck, LIMC0))
    1019                 :          0 :     LIMC2 = LIMC0;
    1020                 :            :   else
    1021                 :       3365 :     LIMC2 = high;
    1022         [ -  + ]:       3365 :   if (LIMC2 > LIMCMAX) LIMC2 = LIMCMAX;
    1023                 :       3365 :   LIMC0 = (long)(cbach*LOGD2);
    1024         [ +  + ]:       3365 :   LIMC = cbach ? LIMC0 : LIMC2;
    1025                 :       3365 :   LIMC = maxss(LIMC, nthidealquad(D, 2));
    1026                 :            : 
    1027                 :            : /* LIMC = Max(cbach*(log D)^2, exp(sqrt(log D loglog D) / 8)) */
    1028                 :            : START:
    1029                 :            :   do
    1030                 :            :   {
    1031         [ +  + ]:       3395 :     if (!FIRST) LIMC = check_LIMC(LIMC,LIMCMAX);
    1032 [ -  + ][ #  # ]:       3395 :     if (DEBUGLEVEL>2 && LIMC > LIMC0)
    1033         [ #  # ]:          0 :       err_printf("%s*** Bach constant: %f\n", FIRST?"":"\n", LIMC/LOGD2);
    1034                 :       3395 :     FIRST = 0; avma = av;
    1035         [ +  + ]:       3395 :     if (BQ.subFB) gunclone(BQ.subFB);
    1036         [ -  + ]:       3395 :     if (BQ.powsubFB) gunclone(BQ.powsubFB);
    1037                 :       3395 :     clearhash(BQ.hashtab);
    1038         [ +  + ]:       3395 :     if (LIMC < cp) LIMC = cp;
    1039         [ +  + ]:       3395 :     if (LIMC2 < LIMC) LIMC2 = LIMC;
    1040         [ +  + ]:       3395 :     if (BQ.PRECREG) qfr_data_init(QFR.D, BQ.PRECREG, &QFR);
    1041                 :            : 
    1042                 :       3395 :     FBquad(&BQ, LIMC2, LIMC, &GRHcheck);
    1043         [ -  + ]:       3395 :     if (DEBUGLEVEL>2) timer_printf(&T, "factor base");
    1044                 :       3395 :     BQ.subFB = subFBquad(&BQ, QFR.D, lim + 0.5, minSFB);
    1045         [ -  + ]:       3395 :     if (DEBUGLEVEL>2) timer_printf(&T, "subFBquad = %Ps",
    1046                 :            :                                  vecpermute(BQ.FB, BQ.subFB));
    1047                 :       3395 :     nsubFB = lg(BQ.subFB) - 1;
    1048                 :            :   }
    1049 [ +  + ][ +  + ]:       3395 :   while (nsubFB < (expi(D) > 15 ? 3 : 2));
    1050                 :            :   /* invhr = 2^r1 (2pi)^r2 / sqrt(D) w ~ L(chi,1) / hR */
    1051         [ +  + ]:       3365 :   invhr = gmul(dbltor((BQ.PRECREG?2.:PI)/sdrc), compute_invresquad(&GRHcheck));
    1052                 :       3365 :   BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
    1053         [ -  + ]:       3365 :   if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
    1054         [ +  - ]:       3365 :   BQ.limhash = (LIMC & HIGHMASK)? (HIGHBIT>>1): LIMC*LIMC;
    1055                 :            : 
    1056                 :       3365 :   need = BQ.KC + RELSUP - 2;
    1057                 :       3365 :   current = 0;
    1058                 :       3365 :   W = NULL;
    1059                 :       3365 :   sfb_trials = nreldep = nrelsup = 0;
    1060                 :       3365 :   s = nsubFB + RELSUP;
    1061                 :       3365 :   av2 = avma;
    1062                 :            : 
    1063                 :            :   do
    1064                 :            :   {
    1065 [ +  + ][ +  + ]:       6635 :     if ((nreldep & 3) == 1 || (nrelsup & 7) == 1) {
    1066         [ -  + ]:       3130 :       if (DEBUGLEVEL>2) err_printf("*** Changing sub factor base\n");
    1067                 :       3130 :       gunclone(BQ.subFB);
    1068                 :       3130 :       gunclone(BQ.powsubFB);
    1069                 :       3130 :       BQ.subFB = gclone(vecslice(BQ.vperm, 1, nsubFB));
    1070                 :       3130 :       BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
    1071         [ -  + ]:       3130 :       if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
    1072                 :       3130 :       clearhash(BQ.hashtab);
    1073                 :            :     }
    1074                 :       6635 :     need += 2;
    1075                 :       6635 :     mat    = cgetg(need+1, t_MAT);
    1076                 :       6635 :     extraC = cgetg(need+1, t_VEC);
    1077         [ +  + ]:       6635 :     if (!W) { /* first time */
    1078                 :       3365 :       C = extraC;
    1079                 :       3365 :       triv = trivial_relations(&BQ, mat, C);
    1080         [ -  + ]:       3365 :       if (DEBUGLEVEL>2) err_printf("KC = %ld, need %ld relations\n", BQ.KC, need);
    1081                 :            :     } else {
    1082                 :       3270 :       triv = 0;
    1083         [ -  + ]:       3270 :       if (DEBUGLEVEL>2) err_printf("...need %ld more relations\n", need);
    1084                 :            :     }
    1085         [ +  + ]:       6635 :     if (BQ.PRECREG) {
    1086         [ +  + ]:     135490 :       for (i = triv+1; i<=need; i++) {
    1087                 :     131813 :         gel(mat,i) = zero_zv(BQ.KC);
    1088                 :     131813 :         gel(extraC,i) = cgetr(BQ.PRECREG);
    1089                 :            :       }
    1090                 :       3677 :       real_relations(&BQ, need - triv, &current, s,LIMC,mat + triv,extraC + triv);
    1091                 :            :     } else {
    1092         [ +  + ]:     146228 :       for (i = triv+1; i<=need; i++) {
    1093                 :     143270 :         gel(mat,i) = zero_zv(BQ.KC);
    1094                 :     143270 :         gel(extraC,i) = gen_0;
    1095                 :            :       }
    1096                 :       2958 :       imag_relations(&BQ, need - triv, &current, LIMC,mat + triv);
    1097                 :            :     }
    1098                 :            : 
    1099         [ +  + ]:       6635 :     if (!W)
    1100                 :       3365 :       W = hnfspec_i(mat,BQ.vperm,&dep,&B,&C,nsubFB);
    1101                 :            :     else
    1102                 :       3270 :       W = hnfadd_i(W,BQ.vperm,&dep,&B,&C, mat,extraC);
    1103                 :       6635 :     gerepileall(av2, 4, &W,&C,&B,&dep);
    1104                 :       6635 :     need = BQ.KC - (lg(W)-1) - (lg(B)-1);
    1105         [ +  + ]:       6635 :     if (need)
    1106                 :            :     {
    1107 [ -  + ][ #  # ]:       2124 :       if (++nreldep > 15 && cbach < 1) goto START;
    1108                 :       2124 :       continue;
    1109                 :            :     }
    1110                 :            : 
    1111                 :       4511 :     h = ZM_det_triangular(W);
    1112         [ -  + ]:       4511 :     if (DEBUGLEVEL>2) err_printf("\n#### Tentative class number: %Ps\n", h);
    1113                 :            : 
    1114      [ -  +  + ]:       4511 :     switch(get_R(&BQ, C, (lg(C)-1) - (lg(B)-1) - (lg(W)-1), mulir(h,invhr), &R))
    1115                 :            :     {
    1116                 :            :       case fupb_PRECI:
    1117                 :          0 :         BQ.PRECREG = precdbl(BQ.PRECREG);
    1118                 :          0 :         FIRST = 1; goto START;
    1119                 :            : 
    1120                 :            :       case fupb_RELAT:
    1121         [ +  + ]:       1146 :         if (++nrelsup > MAXRELSUP)
    1122                 :            :         {
    1123 [ -  + ][ #  # ]:         15 :           if (++sfb_trials > SFB_MAX && cbach <= 1) goto START;
    1124         [ +  - ]:         15 :           if (nsubFB < minss(10,BQ.KC)) nsubFB++;
    1125                 :            :         }
    1126                 :       1146 :         need = minss(BQ.KC, nrelsup);
    1127                 :            :     }
    1128                 :            :   }
    1129         [ +  + ]:       6635 :   while (need);
    1130                 :            :   /* DONE */
    1131         [ -  + ]:       3365 :   if (!quad_be_honest(&BQ)) goto START;
    1132         [ -  + ]:       3365 :   if (DEBUGLEVEL>2) timer_printf(&T, "be honest");
    1133                 :       3365 :   clearhash(BQ.hashtab);
    1134                 :       3365 :   free_GRHcheck(&GRHcheck);
    1135                 :            : 
    1136                 :       3365 :   gen = get_clgp(&BQ,W,&cyc,BQ.PRECREG);
    1137                 :       3365 :   gunclone(BQ.subFB);
    1138                 :       3365 :   gunclone(BQ.powsubFB);
    1139                 :       3365 :   res = cgetg(5,t_VEC);
    1140                 :       3365 :   gel(res,1) = h;
    1141                 :       3365 :   gel(res,2) = cyc;
    1142                 :       3365 :   gel(res,3) = gen;
    1143                 :       3370 :   gel(res,4) = R; return gerepilecopy(av0,res);
    1144                 :            : }
    1145                 :            : 
    1146                 :            : GEN
    1147                 :          5 : buchimag(GEN D, GEN c, GEN c2, GEN REL)
    1148                 :          5 : { (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),0); }
    1149                 :            : 
    1150                 :            : GEN
    1151                 :          5 : buchreal(GEN D, GEN flag, GEN c, GEN c2, GEN REL, long prec) {
    1152         [ -  + ]:          5 :   if (signe(flag)) pari_err_IMPL("narrow class group");
    1153                 :          5 :   (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),prec);
    1154                 :            : }
    1155                 :            : 
    1156                 :            : GEN
    1157                 :       3360 : quadclassunit0(GEN x, long flag, GEN data, long prec)
    1158                 :            : {
    1159                 :            :   long lx;
    1160                 :       3360 :   double c1 = 0.0, c2 = 0.0;
    1161                 :            : 
    1162         [ +  + ]:       3360 :   if (!data) lx=1;
    1163                 :            :   else
    1164                 :            :   {
    1165                 :         15 :     lx = lg(data);
    1166         [ -  + ]:         15 :     if (typ(data)!=t_VEC) pari_err_TYPE("quadclassunit", data);
    1167         [ -  + ]:         15 :     if (lx > 7) pari_err_DIM("quadclassunit [tech vector]");
    1168         [ -  + ]:         15 :     if (lx > 3) lx = 3;
    1169                 :            :   }
    1170      [ +  +  + ]:       3360 :   switch(lx)
    1171                 :            :   {
    1172                 :         10 :     case 3: c2 = gtodouble(gel(data,2));
    1173                 :         15 :     case 2: c1 = gtodouble(gel(data,1));
    1174                 :            :   }
    1175         [ -  + ]:       3360 :   if (flag) pari_err_IMPL("narrow class group");
    1176                 :       3360 :   return Buchquad(x,c1,c2,prec);
    1177                 :            : }

Generated by: LCOV version 1.9