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