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 : #include "mpqs.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_quadclassunit
19 :
20 : /*******************************************************************/
21 : /* */
22 : /* CLASS GROUP AND REGULATOR (McCURLEY, BUCHMANN) */
23 : /* QUADRATIC FIELDS */
24 : /* */
25 : /*******************************************************************/
26 : /* For largeprime() hashtable. Note that hashed pseudoprimes are odd (unless
27 : * 2 | index), hence the low order bit is not useful. So we hash
28 : * HASHBITS bits starting at bit 1, not bit 0 */
29 : #define HASHBITS 10
30 : static const long HASHT = 1L << HASHBITS;
31 :
32 : static long
33 2565952 : hash(long q) { return (q & ((1L << (HASHBITS+1)) - 1)) >> 1; }
34 : #undef HASHBITS
35 :
36 : /* See buch2.c:
37 : * B->subFB contains split p such that \prod p > sqrt(B->Disc)
38 : * B->powsubFB contains powers of forms in B->subFB */
39 : #define RANDOM_BITS 4
40 : static const long CBUCH = (1L<<RANDOM_BITS)-1;
41 :
42 : struct buch_quad
43 : {
44 : ulong limhash;
45 : long KC, KC2, PRECREG;
46 : long *primfact, *exprimfact, **hashtab;
47 : GEN FB, numFB, prodFB;
48 : GEN powsubFB, vperm, subFB, badprim;
49 : struct qfr_data *q;
50 : };
51 :
52 : /*******************************************************************/
53 : /* */
54 : /* Routines related to binary quadratic forms (for internal use) */
55 : /* */
56 : /*******************************************************************/
57 : /* output canonical representative wrt projection Cl^+ --> Cl (a > 0) */
58 : static GEN
59 1166963 : qfr3_canon(GEN x, struct qfr_data *S)
60 : {
61 1166963 : GEN a = gel(x,1), c = gel(x,3);
62 1166963 : if (signe(a) < 0) {
63 403298 : if (absequalii(a,c)) return qfr3_rho(x, S);
64 403284 : setsigne(a, 1);
65 403284 : setsigne(c,-1);
66 : }
67 1166949 : return x;
68 : }
69 : static GEN
70 3710 : qfr3_canon_safe(GEN x, struct qfr_data *S)
71 : {
72 3710 : GEN a = gel(x,1), c = gel(x,3);
73 3710 : if (signe(a) < 0) {
74 224 : if (absequalii(a,c)) return qfr3_rho(x, S);
75 224 : gel(x,1) = negi(a);
76 224 : gel(x,3) = negi(c);
77 : }
78 3710 : return x;
79 : }
80 : static GEN
81 1952748 : qfr5_canon(GEN x, struct qfr_data *S)
82 : {
83 1952748 : GEN a = gel(x,1), c = gel(x,3);
84 1952748 : if (signe(a) < 0) {
85 723268 : if (absequalii(a,c)) return qfr5_rho(x, S);
86 723254 : setsigne(a, 1);
87 723254 : setsigne(c,-1);
88 : }
89 1952734 : return x;
90 : }
91 : static GEN
92 1735160 : QFR5_comp(GEN x,GEN y, struct qfr_data *S)
93 1735160 : { return qfr5_canon(qfr5_comp(x,y,S), S); }
94 : static GEN
95 1007545 : QFR3_comp(GEN x, GEN y, struct qfr_data *S)
96 1007545 : { return qfr3_canon(qfr3_comp(x,y,S), S); }
97 :
98 : /* compute rho^n(x) */
99 : static GEN
100 220157 : qfr5_rho_pow(GEN x, long n, struct qfr_data *S)
101 : {
102 : long i;
103 220157 : pari_sp av = avma;
104 2194360 : for (i=1; i<=n; i++)
105 : {
106 1974203 : x = qfr5_rho(x,S);
107 1974203 : if (gc_needed(av,1))
108 : {
109 0 : if(DEBUGMEM>1) pari_warn(warnmem,"qfr5_rho_pow");
110 0 : x = gc_GEN(av, x);
111 : }
112 : }
113 220157 : return gc_GEN(av, x);
114 : }
115 :
116 : static GEN
117 217588 : qfr5_pf(struct qfr_data *S, long p, long prec)
118 : {
119 217588 : GEN y = primeform_u(S->D,p);
120 217588 : return qfr5_canon(qfr5_red(qfr_to_qfr5(y,prec), S), S);
121 : }
122 :
123 : static GEN
124 159390 : qfr3_pf(struct qfr_data *S, long p)
125 : {
126 159390 : GEN y = primeform_u(S->D,p);
127 159390 : return qfr3_canon(qfr3_red(y, S), S);
128 : }
129 :
130 : #define qfi_pf primeform_u
131 :
132 : /* Warning: ex[0] not set in general */
133 : static GEN
134 4062950 : init_form(struct buch_quad *B, GEN ex,
135 : GEN (*comp)(GEN,GEN,struct qfr_data *S))
136 : {
137 4062950 : long i, l = lg(B->powsubFB);
138 4062950 : GEN F = NULL;
139 22922360 : for (i=1; i<l; i++)
140 18866302 : if (ex[i])
141 : {
142 17688905 : GEN t = gmael(B->powsubFB,i,ex[i]);
143 17688905 : F = F? comp(F, t, B->q): t;
144 : }
145 4056058 : return F;
146 : }
147 : static GEN
148 197015 : qfr5_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFR5_comp); }
149 :
150 : static GEN
151 11744917 : QFI_comp(GEN x, GEN y, struct qfr_data *S) { (void)S; return qfbcomp_i(x,y); }
152 :
153 : static GEN
154 160710 : qfi_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFI_comp); }
155 :
156 : static GEN
157 3704739 : random_form(struct buch_quad *B, GEN ex,
158 : GEN (*comp)(GEN,GEN, struct qfr_data *S))
159 : {
160 3704739 : long i, l = lg(ex);
161 3704739 : pari_sp av = avma;
162 : GEN F;
163 : for(;;)
164 : {
165 20597149 : for (i=1; i<l; i++) ex[i] = random_bits(RANDOM_BITS);
166 3705062 : if ((F = init_form(B, ex, comp))) return F;
167 512 : set_avma(av);
168 : }
169 : }
170 : static GEN
171 161567 : qfr3_random(struct buch_quad *B,GEN ex){ return random_form(B, ex, &QFR3_comp); }
172 : static GEN
173 3543167 : qfi_random(struct buch_quad *B,GEN ex) { return random_form(B, ex, &QFI_comp); }
174 :
175 : /*******************************************************************/
176 : /* */
177 : /* Common subroutines */
178 : /* */
179 : /*******************************************************************/
180 : /* Prime ideals of norm < D = 4 log^2(disc K) generate Cl(K) under GRH
181 : * (Grenie-Molteni). Failure with a factor base of primes with norm < LIMC.
182 : * Suggest new value.
183 : *
184 : * A value of c * log^2(disc K) in [0.3, 1] is OK for most fields and no
185 : * example is known where c >= 2 is necessary to generate the class group.
186 : * Nevertheless, when the bound (hence the factorbase) is small, the index
187 : * calculus algorithm may fail even if the class group is generated. We allow a
188 : * return value > D for this reason [#2653] */
189 : ulong
190 2348 : bnf_increase_LIMC(ulong LIMC, ulong D)
191 : {
192 2348 : if (LIMC <= D / 13.333) /* c <= 0.3 */
193 241 : LIMC *= 2; /* tiny c, roughly double it */
194 : else
195 2107 : LIMC += maxuu(1, D / 20); /* large c, add 0.2 to it */
196 2348 : return LIMC;
197 : }
198 :
199 : /* Is |q| <= p ? */
200 : static int
201 10416 : isless_iu(GEN q, ulong p) {
202 10416 : long l = lgefint(q);
203 10416 : return l==2 || (l == 3 && uel(q,2) <= p);
204 : }
205 :
206 : static GEN
207 5031785 : Z_isquasismooth_prod(GEN N, GEN P)
208 : {
209 5031785 : P = gcdii(P,N);
210 10036003 : while (!is_pm1(P))
211 : {
212 4999344 : N = diviiexact(N, P);
213 5007523 : P = gcdii(N, P);
214 : }
215 5026531 : return N;
216 : }
217 :
218 : static long
219 5038528 : factorquad(struct buch_quad *B, GEN f, long nFB, ulong limp)
220 : {
221 : ulong X;
222 5038528 : long i, lo = 0;
223 5038528 : GEN F, x = gel(f,1), FB = B->FB, P = B->primfact, E = B->exprimfact;
224 5038528 : if (B->badprim && !is_pm1(gcdii(x, B->badprim))) return 0;
225 5031591 : F = Z_isquasismooth_prod(x, B->prodFB);
226 5026395 : if (cmpiu(F, B->limhash) > 0) return 0;
227 4402476 : for (i=1; lgefint(x) > 3; i++)
228 : {
229 10416 : ulong p = uel(FB,i), r;
230 10416 : GEN q = absdiviu_rem(x, p, &r);
231 10416 : if (!r)
232 : {
233 1530 : long k = 0;
234 2633 : do { k++; x = q; q = absdiviu_rem(x, p, &r); } while (!r);
235 1530 : lo++; P[lo] = p; E[lo] = k;
236 : }
237 10416 : if (isless_iu(q,p)) {
238 1 : if (lgefint(x) == 3) { X = uel(x,2); goto END; }
239 34 : return 0;
240 : }
241 10415 : if (i == nFB) return 0;
242 : }
243 4392060 : X = uel(x,2);
244 4392060 : if (X == 1) { P[0] = 0; return 1; }
245 68609411 : for (;; i++)
246 68609411 : { /* single precision affair, split for efficiency */
247 72982700 : ulong p = uel(FB,i);
248 72982700 : ulong q = X / p, r = X % p; /* gcc makes a single div */
249 72982700 : if (!r)
250 : {
251 5621189 : long k = 0;
252 6789128 : do { k++; X = q; q = X / p; r = X % p; } while (!r);
253 5621189 : lo++; P[lo] = p; E[lo] = k;
254 : }
255 72982700 : if (q <= p) break;
256 68642329 : if (i == nFB) return 0;
257 : }
258 4340372 : END:
259 4340372 : if (X > B->limhash) return 0;
260 4340372 : if (X != 1 && X <= limp) { lo++; P[lo] = X; E[lo] = 1; X = 1; }
261 4340372 : P[0] = lo; return X;
262 : }
263 :
264 : /* Check for a "large prime relation" involving q; q may not be prime */
265 : static long *
266 2565918 : largeprime(struct buch_quad *B, long q, GEN ex, long np, long nrho)
267 : {
268 2565918 : const long hashv = hash(q);
269 2565956 : long *pt, i, l = lg(B->subFB);
270 :
271 2787663 : for (pt = B->hashtab[hashv]; ; pt = (long*) pt[0])
272 : {
273 2787663 : if (!pt)
274 : {
275 2366486 : pt = (long*) pari_malloc((l+3) * sizeof(long));
276 2366724 : *pt++ = nrho; /* nrho = pt[-3] */
277 2366724 : *pt++ = np; /* np = pt[-2] */
278 2366724 : *pt++ = q; /* q = pt[-1] */
279 2366724 : pt[0] = (long)B->hashtab[hashv];
280 14990634 : for (i=1; i<l; i++) pt[i]=ex[i];
281 2366724 : B->hashtab[hashv]=pt; return NULL;
282 : }
283 421177 : if (pt[-1] == q) break;
284 : }
285 238043 : for(i=1; i<l; i++)
286 234749 : if (pt[i] != ex[i]) return pt;
287 3294 : return (pt[-2]==np)? NULL: pt;
288 : }
289 :
290 : static void
291 169819 : clearhash(long **hash)
292 : {
293 : long *pt;
294 : long i;
295 174024798 : for (i=0; i<HASHT; i++) {
296 176221748 : for (pt = hash[i]; pt; ) {
297 2366769 : void *z = (void*)(pt - 3);
298 2366769 : pt = (long*) pt[0]; pari_free(z);
299 : }
300 173854979 : hash[i] = NULL;
301 : }
302 169805 : }
303 :
304 : /* last prime stored */
305 : ulong
306 0 : GRH_last_prime(GRHcheck_t *S) { return (S->primes + S->nprimes-1)->p; }
307 : /* ensure that S->primes can hold at least nb primes */
308 : void
309 402220 : GRH_ensure(GRHcheck_t *S, long nb)
310 : {
311 402220 : if (S->maxprimes <= nb)
312 : {
313 0 : do S->maxprimes *= 2; while (S->maxprimes <= nb);
314 0 : pari_realloc_ip((void**)&S->primes, S->maxprimes*sizeof(*S->primes));
315 : }
316 402220 : }
317 : /* cache data for all primes up to the LIM */
318 : static void
319 1175323 : cache_prime_quad(GRHcheck_t *S, ulong LIM, GEN D)
320 : {
321 : GRHprime_t *pr;
322 : long nb;
323 :
324 1175323 : if (S->limp >= LIM) return;
325 72563 : nb = (long)primepi_upper_bound((double)LIM); /* #{p <= LIM} <= nb */
326 72563 : GRH_ensure(S, nb+1); /* room for one extra prime */
327 72563 : for (pr = S->primes + S->nprimes;;)
328 12671724 : {
329 12744287 : ulong p = u_forprime_next(&(S->P));
330 12743856 : pr->p = p;
331 12743856 : pr->logp = log((double)p);
332 12743856 : pr->dec = (GEN)kroiu(D,p);
333 12744285 : S->nprimes++;
334 12744285 : pr++;
335 : /* store up to nextprime(LIM) included */
336 12744285 : if (p >= LIM) { S->limp = p; break; }
337 : }
338 : }
339 :
340 : static GEN
341 70438 : compute_invresquad(GRHcheck_t *S, long LIMC)
342 : {
343 70438 : pari_sp av = avma;
344 70438 : GEN invres = real_1(DEFAULTPREC);
345 70439 : double limp = log((double)LIMC) / 2;
346 : GRHprime_t *pr;
347 : long i;
348 12718968 : for (pr = S->primes, i = S->nprimes; i > 0; pr++, i--)
349 : {
350 12650967 : long s = (long)pr->dec;
351 12650967 : if (s)
352 : {
353 12535699 : ulong p = pr->p;
354 12535699 : if (s > 0 || pr->logp <= limp)
355 : /* Both p and P contribute */
356 6370611 : invres = mulur(p - s, divru(invres, p));
357 6165088 : else if (s<0)
358 : /* Only p contributes */
359 6165099 : invres = mulur(p, divru(invres, p - 1));
360 : }
361 : }
362 68001 : return gc_leaf(av, invres);
363 : }
364 :
365 : /* p | conductor of order of disc D ? */
366 : static int
367 392192 : is_bad(GEN D, ulong p)
368 : {
369 392192 : pari_sp av = avma;
370 392192 : if (p == 2)
371 : {
372 89639 : long r = mod16(D) >> 1;
373 89639 : if (r && signe(D) < 0) r = 8-r;
374 89639 : return (r < 4);
375 : }
376 302553 : return gc_bool(av, dvdii(D, sqru(p))); /* p^2 | D ? */
377 : }
378 :
379 : /* returns the n-th suitable ideal for the factorbase */
380 : static long
381 70426 : nthidealquad(GEN D, long n)
382 : {
383 70426 : pari_sp av = avma;
384 : forprime_t S;
385 : ulong p;
386 70426 : (void)u_forprime_init(&S, 2, ULONG_MAX);
387 357876 : while ((p = u_forprime_next(&S)) && n > 0)
388 287465 : if (!is_bad(D, p) && kroiu(D, p) >= 0) n--;
389 70415 : return gc_long(av, p);
390 : }
391 :
392 : static int
393 1034386 : quadGRHchk(GEN D, GRHcheck_t *S, ulong LIMC)
394 : {
395 1034386 : double logC = log((double)LIMC), SA = 0, SB = 0;
396 : long i;
397 :
398 1034386 : cache_prime_quad(S, LIMC, D);
399 1034387 : for (i = 0;; i++)
400 29533474 : {
401 30567861 : GRHprime_t *pr = S->primes+i;
402 30567861 : ulong p = pr->p;
403 : long M;
404 : double logNP, q, A, B;
405 30567861 : if (p > LIMC) break;
406 29622900 : if ((long)pr->dec < 0)
407 : {
408 14790027 : logNP = 2 * pr->logp;
409 14790027 : q = 1/(double)p;
410 : }
411 : else
412 : {
413 14832873 : logNP = pr->logp;
414 14832873 : q = 1/sqrt((double)p);
415 : }
416 29622900 : A = logNP * q; B = logNP * A; M = (long)(logC/logNP);
417 29622900 : if (M > 1)
418 : {
419 2502939 : double inv1_q = 1 / (1-q);
420 2502939 : A *= (1 - pow(q, (double) M)) * inv1_q;
421 2502939 : B *= (1 - pow(q, (double) M)*(M+1 - M*q)) * inv1_q * inv1_q;
422 : }
423 29622900 : if ((long)pr->dec>0) { SA += 2*A;SB += 2*B; } else { SA += A; SB += B; }
424 29622900 : if (p == LIMC) break;
425 : }
426 1034387 : return GRHok(S, logC, SA, SB);
427 : }
428 :
429 : /* C2 >= C1; create B->FB, B->numFB; set B->badprim. Return L(kro_D, 1) */
430 : static void
431 70517 : FBquad(struct buch_quad *B, ulong C2, ulong C1, GRHcheck_t *S)
432 : {
433 70517 : GEN D = B->q->D;
434 : long i;
435 : pari_sp av;
436 : GRHprime_t *pr;
437 :
438 70517 : cache_prime_quad(S, C2, D);
439 70517 : pr = S->primes;
440 70517 : B->numFB = cgetg(C2+1, t_VECSMALL);
441 70516 : B->FB = cgetg(C2+1, t_VECSMALL);
442 70483 : av = avma;
443 70483 : B->KC = 0; i = 0;
444 70483 : B->badprim = gen_1;
445 2813636 : for (;; pr++) /* p <= C2 */
446 2813636 : {
447 2884119 : ulong p = pr->p;
448 2884119 : if (!B->KC && p > C1) B->KC = i;
449 2884119 : if (p > C2) break;
450 2821090 : switch ((long)pr->dec)
451 : {
452 1390196 : case -1: break; /* inert */
453 104730 : case 0: /* ramified */
454 104730 : if (is_bad(D, p)) { B->badprim = muliu(B->badprim, p); break; }
455 : /* fall through */
456 : default: /* split */
457 1430865 : i++; B->numFB[p] = i; B->FB[i] = p; break;
458 : }
459 2821124 : if (p == C2)
460 : {
461 7488 : if (!B->KC) B->KC = i;
462 7488 : break;
463 : }
464 : }
465 70517 : B->KC2 = i;
466 70517 : setlg(B->FB, B->KC2+1);
467 70517 : if (B->badprim != gen_1)
468 42 : B->badprim = gc_INT(av, B->badprim);
469 : else
470 : {
471 70475 : B->badprim = NULL; set_avma(av);
472 : }
473 70517 : B->prodFB = zv_prod_Z(B->FB);
474 70515 : }
475 :
476 : /* create B->vperm, return B->subFB */
477 : static GEN
478 70514 : subFBquad(struct buch_quad *B, GEN D, double PROD, long minSFB)
479 : {
480 70514 : long i, j, lgsub = 1, ino = 1, lv = B->KC+1;
481 70514 : double prod = 1.;
482 : pari_sp av;
483 : GEN no;
484 :
485 70514 : B->vperm = cgetg(lv, t_VECSMALL);
486 70514 : av = avma;
487 70514 : no = cgetg(lv, t_VECSMALL);
488 335076 : for (j = 1; j < lv; j++)
489 : {
490 334992 : ulong p = uel(B->FB,j);
491 334992 : if (!umodiu(D, p)) no[ino++] = j; /* ramified */
492 : else
493 : {
494 255160 : B->vperm[lgsub++] = j;
495 255160 : prod *= p;
496 255160 : if (lgsub > minSFB && prod > PROD) break;
497 : }
498 : }
499 : /* lgsub >= 1 otherwise quadGRHchk is false */
500 70517 : i = lgsub;
501 150351 : for (j = 1; j < ino;i++,j++) B->vperm[i] = no[j];
502 1165958 : for ( ; i < lv; i++) B->vperm[i] = i;
503 70517 : no = gclone(vecslice(B->vperm, 1, lgsub-1));
504 70517 : set_avma(av); return no;
505 : }
506 :
507 : /* assume n >= 1, x[i][j] = B->subFB[i]^j, for j = 1..n */
508 : static GEN
509 99319 : powsubFBquad(struct buch_quad *B, long n)
510 : {
511 99319 : pari_sp av = avma;
512 99319 : long i,j, l = lg(B->subFB);
513 99319 : GEN F, y, x = cgetg(l, t_VEC), D = B->q->D;
514 :
515 99319 : if (B->PRECREG) /* real */
516 : {
517 39627 : for (i=1; i<l; i++)
518 : {
519 34510 : F = qfr5_pf(B->q, B->FB[B->subFB[i]], B->PRECREG);
520 34510 : y = cgetg(n+1, t_VEC); gel(x,i) = y;
521 34510 : gel(y,1) = F;
522 552160 : for (j=2; j<=n; j++) gel(y,j) = QFR5_comp(gel(y,j-1), F, B->q);
523 : }
524 : }
525 : else /* imaginary */
526 : {
527 427549 : for (i=1; i<l; i++)
528 : {
529 333359 : F = qfi_pf(D, B->FB[B->subFB[i]]);
530 333340 : y = cgetg(n+1, t_VEC); gel(x,i) = y;
531 333924 : gel(y,1) = F;
532 5327766 : for (j=2; j<=n; j++) gel(y,j) = qfbcomp_i(gel(y,j-1), F);
533 : }
534 : }
535 99307 : x = gclone(x); set_avma(av); return x;
536 : }
537 :
538 : static void
539 98252 : sub_fact(struct buch_quad *B, GEN col, GEN F)
540 : {
541 98252 : GEN b = gel(F,2);
542 : long i;
543 207961 : for (i=1; i<=B->primfact[0]; i++)
544 : {
545 109709 : ulong p = B->primfact[i], k = B->numFB[p];
546 109709 : long e = B->exprimfact[i];
547 109709 : if (umodiu(b, p<<1) > p) e = -e;
548 109709 : col[k] -= e;
549 : }
550 98252 : }
551 :
552 : #if 0
553 : static void
554 : dbg_fact(struct buch_quad *B)
555 : {
556 : long i;
557 : for (i=1; i<=B->primfact[0]; i++)
558 : {
559 : ulong p = B->primfact[i];
560 : long e = B->exprimfact[i];
561 : err_printf("%lu^%ld ",p,e);
562 : }
563 : }
564 :
565 : static void
566 : chk_fact(struct buch_quad *B, GEN col)
567 : {
568 : long i, l = lg(col);
569 : GEN Q = qfi_pf(B->q->D, 1);
570 : for (i=1; i< l; i++)
571 : {
572 : ulong p = B->FB[i];
573 : long k = col[i];
574 : Q = qfbcomp(qfbpowraw(qfi_pf(B->q->D, p),k),Q);
575 : }
576 : if (!gequal1(gel(Q,1))) pari_err_BUG("chk_fact");
577 : }
578 : #endif
579 :
580 : static void
581 1893320 : add_fact(struct buch_quad *B, GEN col, GEN F)
582 : {
583 1893320 : GEN b = gel(F,2);
584 : long i;
585 5932637 : for (i=1; i<=B->primfact[0]; i++)
586 : {
587 4039177 : ulong p = B->primfact[i], k = B->numFB[p];
588 4039177 : long e = B->exprimfact[i];
589 4039177 : if (umodiu(b, p<<1) > p) e = -e;
590 4039317 : col[k] += e;
591 : }
592 1893460 : }
593 :
594 : static GEN
595 70426 : get_clgp(struct buch_quad *B, GEN W, GEN *ptD)
596 : {
597 70426 : GEN res, init, u1, D = ZM_snf_group(W,NULL,&u1);
598 70425 : long i, j, l = lg(W), c = lg(D);
599 :
600 70425 : res=cgetg(c,t_VEC); init = cgetg(l,t_VEC);
601 216332 : for (i=1; i<l; i++) gel(init,i) = primeform_u(B->q->D, B->FB[B->vperm[i]]);
602 195940 : for (j=1; j<c; j++)
603 : {
604 125515 : GEN g = NULL;
605 125515 : if (signe(B->q->D) > 0)
606 : {
607 13328 : for (i=1; i<l; i++)
608 : {
609 9618 : GEN t, u = gcoeff(u1,i,j);
610 9618 : if (!signe(u)) continue;
611 4543 : t = qfr3_pow(gel(init,i), u, B->q);
612 4543 : g = g? qfr3_comp(g, t, B->q): t;
613 : }
614 3710 : g = qfr3_to_qfr(qfr3_canon_safe(qfr3_red(g, B->q), B->q), B->q->D);
615 : }
616 : else
617 : {
618 423834 : for (i=1; i<l; i++)
619 : {
620 302026 : GEN t, u = gcoeff(u1,i,j);
621 302026 : if (!signe(u)) continue;
622 203636 : t = powgi(gel(init,i), u);
623 203642 : g = g? qfbcomp_i(g, t): t;
624 : }
625 : }
626 125516 : gel(res,j) = g;
627 : }
628 70425 : *ptD = D; return res;
629 : }
630 :
631 : static long
632 70440 : trivial_relations(struct buch_quad *B, GEN mat, GEN C)
633 : {
634 70440 : long i, j = 0;
635 70440 : GEN col, D = B->q->D;
636 1500618 : for (i = 1; i <= B->KC; i++)
637 : { /* ramified prime ==> trivial relation */
638 1430176 : if (umodiu(D, B->FB[i])) continue;
639 104458 : col = zero_zv(B->KC);
640 104466 : col[i] = 2; j++;
641 104466 : gel(mat,j) = col;
642 104466 : gel(C,j) = gen_0;
643 : }
644 70442 : return j;
645 : }
646 :
647 : static void
648 0 : dbg_all(pari_timer *T, const char *phase, long s, long n)
649 : {
650 0 : err_printf("\n");
651 0 : timer_printf(T, "%s rel [#rel/#test = %ld/%ld]", phase,s,n);
652 0 : }
653 :
654 : /* Imaginary Quadratic fields */
655 :
656 : static void
657 8813 : rel_to_col(struct buch_quad *B, GEN col, GEN rel, GEN b)
658 : {
659 8813 : GEN P = gel(rel, 1), E = gel(rel, 2);
660 8813 : long i, lP = lg(P);
661 105784 : for (i=1; i<lP; i++)
662 : {
663 96971 : ulong p = uel(P, i), e = uel(E, i);
664 96971 : col[B->numFB[p]] += umodiu(b, p<<1) > p ? -e :e;
665 : }
666 8813 : }
667 :
668 : static void
669 99151 : imag_relations(struct buch_quad *B, long need, long *pc, ulong LIMC, GEN mat)
670 : {
671 : pari_timer T;
672 99151 : long lgsub = lg(B->subFB), current = *pc, nbtest = 0, s = 0;
673 : long i, fpc;
674 : pari_sp av;
675 99151 : GEN col, form, ex = cgetg(lgsub, t_VECSMALL);
676 :
677 99151 : if (!current) current = 1;
678 99151 : if (DEBUGLEVEL>2) timer_start(&T);
679 99151 : av = avma;
680 : for(;;)
681 : {
682 3641937 : if (s >= need) break;
683 3542787 : set_avma(av);
684 3542681 : form = qfi_random(B,ex);
685 3541171 : form = qfbcomp_i(form, qfi_pf(B->q->D, B->FB[current]));
686 3540993 : nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
687 3542085 : if (!fpc)
688 : {
689 291563 : if (DEBUGLEVEL>3) err_printf(".");
690 291563 : if ((nbtest & 0xff) == 0 && ++current > B->KC) current = 1;
691 291563 : continue;
692 : }
693 3250522 : if (fpc > 1)
694 : {
695 1800825 : long *fpd = largeprime(B,fpc,ex,current,0);
696 : ulong b1, b2, p;
697 : GEN form2;
698 1801271 : if (!fpd)
699 : {
700 1640583 : if (DEBUGLEVEL>3) err_printf(".");
701 1640578 : continue;
702 : }
703 160688 : form2 = qfbcomp_i(qfi_factorback(B,fpd), qfi_pf(B->q->D, B->FB[fpd[-2]]));
704 160713 : p = fpc << 1;
705 160713 : b1 = umodiu(gel(form2,2), p);
706 160713 : b2 = umodiu(gel(form,2), p);
707 160713 : if (b1 != b2 && b1+b2 != p) continue;
708 :
709 160713 : col = gel(mat,++s);
710 160713 : add_fact(B,col, form);
711 160714 : (void)factorquad(B,form2,B->KC,LIMC);
712 160719 : if (b1==b2)
713 : {
714 414762 : for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
715 80626 : sub_fact(B, col, form2); col[fpd[-2]]++;
716 : }
717 : else
718 : {
719 410975 : for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
720 80093 : add_fact(B, col, form2); col[fpd[-2]]--;
721 : }
722 160720 : if (DEBUGLEVEL>2) err_printf(" %ldP",s);
723 : }
724 : else
725 : {
726 1449697 : col = gel(mat,++s);
727 6962359 : for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
728 1449697 : add_fact(B, col, form);
729 1450075 : if (DEBUGLEVEL>2) err_printf(" %ld",s);
730 : }
731 1610645 : col[current]--;
732 1610645 : if (++current > B->KC) current = 1;
733 : }
734 99150 : if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
735 99150 : *pc = current;
736 99150 : }
737 :
738 : static void
739 336 : mpqs_relations(struct buch_quad *B, long need, long *pc, ulong LIMC, GEN mat, mpqs_handle_t *H, GEN missing_primes)
740 : {
741 : pari_timer T;
742 : long i, lV;
743 : GEN V;
744 336 : if (DEBUGLEVEL>2) timer_start(&T);
745 336 : V = mpqs_class_rels(H, need, missing_primes);
746 336 : if (!V) { imag_relations(B, need, pc, LIMC, mat); return; }
747 336 : lV = lg(V);
748 9149 : for (i = 1; i < lV && i <= need; i++)
749 : {
750 8813 : GEN formA = gel(V,i), rel = gel(formA,2), b = gel(formA,1);
751 8813 : GEN col = gel(mat,i);
752 8813 : rel_to_col(B, col, rel, b);
753 : }
754 336 : if (DEBUGLEVEL>2) timer_printf(&T, "MPQS rel [#rel = %ld]", i-1);
755 : }
756 :
757 : static int
758 7 : imag_be_honest(struct buch_quad *B)
759 : {
760 7 : long p, fpc, s = B->KC, nbtest = 0;
761 7 : GEN F, ex = cgetg(lg(B->subFB), t_VECSMALL);
762 7 : pari_sp av = avma;
763 :
764 525 : while (s<B->KC2)
765 : {
766 518 : p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
767 518 : F = qfbcomp_i(qfi_pf(B->q->D, p), qfi_random(B, ex));
768 518 : fpc = factorquad(B,F,s,p-1);
769 518 : if (fpc == 1) { nbtest=0; s++; }
770 : else
771 392 : if (++nbtest > 40) return 0;
772 518 : set_avma(av);
773 : }
774 7 : return 1;
775 : }
776 :
777 : static GEN
778 184688 : dist(GEN e, GEN d, long prec) { return mkvec2(qfr5_dist(e, d, prec), d); }
779 : /* z a pre-allocated pair [t_REAL, t_INT], D = [t,d] from dist() */
780 : static void
781 184688 : dist_set(GEN z, GEN D)
782 : {
783 184688 : affrr(gel(D,1), gel(z,1));
784 184688 : affsi(signe(gel(D,2)) < 0? 1: 0, gel(z,2));
785 184688 : }
786 :
787 : /* Real Quadratic fields */
788 :
789 : static void
790 5145 : real_relations(struct buch_quad *B, long need, long *pc, long lim, ulong LIMC, GEN mat, GEN C)
791 : {
792 : pari_timer T;
793 5145 : long lgsub = lg(B->subFB), prec = B->PRECREG, current = *pc, nbtest=0, s=0;
794 : long i, fpc, endcycle, rhoacc, rho;
795 : /* in a 2nd phase, don't include FB[current] but run along the cyle
796 : * ==> get more units */
797 5145 : int first = (current == 0);
798 : pari_sp av, av1;
799 5145 : GEN d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL);
800 :
801 5145 : if (DEBUGLEVEL>2) timer_start(&T);
802 5145 : if (!current) current = 1;
803 5145 : if (lim > need) lim = need;
804 5145 : av = avma;
805 : for(;;)
806 : {
807 166691 : if (s >= need) break;
808 161546 : if (first && s >= lim) {
809 2142 : first = 0;
810 2142 : if (DEBUGLEVEL>2) dbg_all(&T, "initial", s, nbtest);
811 : }
812 161546 : set_avma(av); form = qfr3_random(B, ex);
813 161546 : if (!first)
814 159369 : form = QFR3_comp(form, qfr3_pf(B->q, B->FB[current]), B->q);
815 161546 : av1 = avma;
816 161546 : form0 = form; form1 = NULL;
817 161546 : endcycle = rhoacc = 0;
818 161546 : rho = -1;
819 :
820 1300964 : CYCLE:
821 1300964 : if (endcycle || rho > 5000)
822 : {
823 21 : if (++current > B->KC) current = 1;
824 21 : continue;
825 : }
826 1300943 : if (gc_needed(av,1))
827 : {
828 0 : if(DEBUGMEM>1) pari_warn(warnmem,"real_relations");
829 0 : (void)gc_all(av1, form1? 2: 1, &form, &form1);
830 : }
831 1300943 : if (rho < 0) rho = 0; /* first time in */
832 : else
833 : {
834 1139397 : form = qfr3_rho(form, B->q); rho++;
835 1139397 : rhoacc++;
836 1139397 : if (first)
837 442526 : endcycle = (absequalii(gel(form,1),gel(form0,1))
838 221263 : && equalii(gel(form,2),gel(form0,2)));
839 : else
840 : {
841 918134 : if (absequalii(gel(form,1), gel(form,3))) /* a = -c */
842 : {
843 0 : if (absequalii(gel(form,1),gel(form0,1)) &&
844 0 : equalii(gel(form,2),gel(form0,2))) continue;
845 0 : form = qfr3_rho(form, B->q); rho++;
846 0 : rhoacc++;
847 : }
848 : else
849 918134 : { setsigne(form[1],1); setsigne(form[3],-1); }
850 918176 : if (equalii(gel(form,1),gel(form0,1)) &&
851 42 : equalii(gel(form,2),gel(form0,2))) continue;
852 : }
853 : }
854 1300943 : nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
855 1300943 : if (!fpc)
856 : {
857 386778 : if (DEBUGLEVEL>3) err_printf(".");
858 386778 : goto CYCLE;
859 : }
860 914165 : if (fpc > 1)
861 : { /* look for Large Prime relation */
862 764946 : long *fpd = largeprime(B,fpc,ex,first? 0: current,rhoacc);
863 : ulong b1, b2, p;
864 : GEN form2;
865 764946 : if (!fpd)
866 : {
867 729477 : if (DEBUGLEVEL>3) err_printf(".");
868 729477 : goto CYCLE;
869 : }
870 35469 : if (!form1)
871 : {
872 35469 : form1 = qfr5_factorback(B,ex);
873 35469 : if (!first)
874 35469 : form1 = QFR5_comp(form1, qfr5_pf(B->q, B->FB[current], prec), B->q);
875 : }
876 35469 : form1 = qfr5_rho_pow(form1, rho, B->q);
877 35469 : rho = 0;
878 :
879 35469 : form2 = qfr5_factorback(B,fpd);
880 35469 : if (fpd[-2])
881 23709 : form2 = QFR5_comp(form2, qfr5_pf(B->q, B->FB[fpd[-2]], prec), B->q);
882 35469 : form2 = qfr5_rho_pow(form2, fpd[-3], B->q);
883 35469 : if (!absequalii(gel(form2,1),gel(form2,3)))
884 : {
885 35469 : setsigne(form2[1], 1);
886 35469 : setsigne(form2[3],-1);
887 : }
888 35469 : p = fpc << 1;
889 35469 : b1 = umodiu(gel(form2,2), p);
890 35469 : b2 = umodiu(gel(form1,2), p);
891 35469 : if (b1 != b2 && b1+b2 != p) goto CYCLE;
892 :
893 35469 : col = gel(mat,++s);
894 35469 : add_fact(B, col, form1);
895 35469 : (void)factorquad(B,form2,B->KC,LIMC);
896 35469 : if (b1==b2)
897 : {
898 135100 : for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
899 17626 : sub_fact(B,col, form2);
900 17626 : if (fpd[-2]) col[fpd[-2]]++;
901 17626 : d = dist(subii(gel(form1,4),gel(form2,4)),
902 17626 : divrr(gel(form1,5),gel(form2,5)), prec);
903 : }
904 : else
905 : {
906 136780 : for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
907 17843 : add_fact(B, col, form2);
908 17843 : if (fpd[-2]) col[fpd[-2]]--;
909 17843 : d = dist(addii(gel(form1,4),gel(form2,4)),
910 17843 : mulrr(gel(form1,5),gel(form2,5)), prec);
911 : }
912 35469 : if (DEBUGLEVEL>2) err_printf(" %ldP",s);
913 : }
914 : else
915 : { /* standard relation */
916 149219 : if (!form1)
917 : {
918 126077 : form1 = qfr5_factorback(B, ex);
919 126077 : if (!first)
920 123900 : form1 = QFR5_comp(form1, qfr5_pf(B->q, B->FB[current], prec), B->q);
921 : }
922 149219 : form1 = qfr5_rho_pow(form1, rho, B->q);
923 149219 : rho = 0;
924 :
925 149219 : col = gel(mat,++s);
926 1147489 : for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
927 149219 : add_fact(B, col, form1);
928 149219 : d = dist(gel(form1,4), gel(form1,5), prec);
929 149219 : if (DEBUGLEVEL>2) err_printf(" %ld",s);
930 : }
931 184688 : dist_set(gel(C,s), d);
932 184688 : if (first)
933 : {
934 25319 : if (s >= lim) continue;
935 23163 : goto CYCLE;
936 : }
937 : else
938 : {
939 159369 : col[current]--;
940 159369 : if (++current > B->KC) current = 1;
941 : }
942 : }
943 5145 : if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
944 5145 : *pc = current;
945 5145 : }
946 :
947 : static int
948 7 : real_be_honest(struct buch_quad *B)
949 : {
950 7 : long p, fpc, s = B->KC, nbtest = 0;
951 7 : GEN F,F0, ex = cgetg(lg(B->subFB), t_VECSMALL);
952 7 : pari_sp av = avma;
953 :
954 28 : while (s<B->KC2)
955 : {
956 21 : p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
957 21 : F = QFR3_comp(qfr3_random(B, ex), qfr3_pf(B->q, p), B->q);
958 21 : for (F0 = F;;)
959 : {
960 49 : fpc = factorquad(B,F,s,p-1);
961 49 : if (fpc == 1) { nbtest=0; s++; break; }
962 28 : if (++nbtest > 40) return 0;
963 28 : F = qfr3_canon(qfr3_rho(F, B->q), B->q);
964 28 : if (equalii(gel(F,1),gel(F0,1))
965 0 : && equalii(gel(F,2),gel(F0,2))) break;
966 : }
967 21 : set_avma(av);
968 : }
969 7 : return 1;
970 : }
971 :
972 : static GEN
973 55146 : crabs(GEN a)
974 : {
975 55146 : return signe(real_i(a)) < 0 ? gneg(a): a;
976 : }
977 :
978 : static GEN
979 33978 : gcdreal(GEN a,GEN b)
980 : {
981 33978 : if (!signe(real_i(a))) return crabs(b);
982 33012 : if (!signe(real_i(b))) return crabs(a);
983 32897 : if (expo(real_i(a))<-5) return crabs(b);
984 12040 : if (expo(real_i(b))<-5) return crabs(a);
985 9093 : a = crabs(a); b = crabs(b);
986 20019 : while (expo(real_i(b)) >= -5 && signe(real_i(b)))
987 : {
988 : long e;
989 10926 : GEN r, q = gcvtoi(divrr(real_i(a),real_i(b)),&e);
990 10926 : if (e > 0) return NULL;
991 10926 : r = gsub(a, gmul(q,b)); a = b; b = r;
992 : }
993 9093 : return crabs(a);
994 : }
995 :
996 : static int
997 91701 : get_R(struct buch_quad *B, GEN C, long sreg, GEN z, GEN *ptR)
998 : {
999 91701 : GEN R = gen_1;
1000 : double c;
1001 : long i;
1002 91701 : if (B->PRECREG)
1003 : {
1004 2982 : R = crabs(gel(C,1));
1005 36960 : for (i=2; i<=sreg; i++)
1006 : {
1007 33978 : R = gcdreal(gel(C,i), R);
1008 33978 : if (!R) return fupb_PRECI;
1009 : }
1010 2982 : if (gexpo(real_i(R)) <= -3)
1011 : {
1012 0 : if (DEBUGLEVEL>2) err_printf("regulator is zero.\n");
1013 0 : return fupb_RELAT;
1014 : }
1015 2982 : if (DEBUGLEVEL>2) err_printf("#### Tentative regulator: %Ps\n",R);
1016 : }
1017 91701 : c = gtodouble(gmul(z, real_i(R)));
1018 91703 : if (c < 0.8 || c > 1.3) return fupb_RELAT;
1019 70424 : *ptR = R; return fupb_NONE;
1020 : }
1021 :
1022 : static int
1023 70424 : quad_be_honest(struct buch_quad *B)
1024 : {
1025 : int r;
1026 70424 : if (B->KC2 <= B->KC) return 1;
1027 14 : if (DEBUGLEVEL>2)
1028 0 : err_printf("be honest for primes from %ld to %ld\n", B->FB[B->KC+1],B->FB[B->KC2]);
1029 14 : r = B->PRECREG? real_be_honest(B): imag_be_honest(B);
1030 14 : if (DEBUGLEVEL>2) err_printf("\n");
1031 14 : return r;
1032 : }
1033 :
1034 : static GEN
1035 70599 : Buchquad_i(GEN D, double cbach, double cbach2, long prec)
1036 : {
1037 70599 : const long MAXRELSUP = 20, SFB_MAX = 3, MPQS_THRESHOLD = 60;
1038 : pari_timer T;
1039 : pari_sp av, av2;
1040 70599 : const long RELSUP = 5;
1041 : long i, s, current, triv, sfb_trials, nrelsup, nreldep, need, nsubFB, minSFB;
1042 : ulong low, high, LIMC0, LIMC, LIMC2, LIMCMAX, cp;
1043 70599 : GEN W, cyc, gen, dep, mat, C, extraC, B, R, invhr, h = NULL; /*-Wall*/
1044 : double drc, sdrc, lim, LOGD, LOGD2;
1045 : GRHcheck_t GRHcheck;
1046 : struct qfr_data q;
1047 : struct buch_quad BQ;
1048 70599 : int FIRST = 1, use_mpqs = 0;
1049 : mpqs_handle_t H;
1050 : GEN missing_primes;
1051 :
1052 70599 : check_quaddisc(D, &s, /*junk*/&i, "Buchquad");
1053 70599 : R = NULL; /* -Wall */
1054 70599 : BQ.q = &q;
1055 70599 : q.D = D;
1056 70599 : if (s < 0)
1057 : {
1058 68443 : if (abscmpiu(q.D,4) <= 0)
1059 175 : retmkvec4(gen_1, cgetg(1,t_VEC), cgetg(1,t_VEC), gen_1);
1060 68268 : prec = BQ.PRECREG = 0;
1061 68268 : if (expi(D) >= MPQS_THRESHOLD)
1062 28 : use_mpqs = 1;
1063 : } else {
1064 2156 : BQ.PRECREG = maxss(prec+EXTRAPREC64, nbits2prec(2*expi(q.D) + 128));
1065 : }
1066 70424 : if (DEBUGLEVEL>2) timer_start(&T);
1067 70424 : BQ.primfact = new_chunk(100);
1068 70424 : BQ.exprimfact = new_chunk(100);
1069 70423 : BQ.hashtab = (long**) new_chunk(HASHT);
1070 72178502 : for (i=0; i<HASHT; i++) BQ.hashtab[i] = NULL;
1071 :
1072 70426 : drc = fabs(gtodouble(q.D));
1073 70424 : LOGD = log(drc);
1074 70424 : LOGD2 = LOGD * LOGD;
1075 :
1076 70424 : sdrc = lim = sqrt(drc);
1077 70424 : if (!BQ.PRECREG) lim /= sqrt(3.);
1078 70424 : cp = (ulong)exp(sqrt(LOGD*log(LOGD)/8.0));
1079 70424 : if (cp < 20) cp = 20;
1080 70424 : if (cbach > 6.) {
1081 0 : if (cbach2 < cbach) cbach2 = cbach;
1082 0 : cbach = 6.;
1083 : }
1084 70424 : if (cbach < 0.)
1085 0 : pari_err_DOMAIN("Buchquad","Bach constant","<",gen_0,dbltor(cbach));
1086 70424 : av = avma;
1087 70424 : BQ.powsubFB = BQ.subFB = NULL;
1088 70424 : minSFB = (expi(D) > 15)? 3: 2;
1089 70424 : init_GRHcheck(&GRHcheck, 2, BQ.PRECREG? 2: 0, LOGD);
1090 70426 : high = low = LIMC0 = maxss((long)(cbach2*LOGD2), 1);
1091 70426 : LIMCMAX = (long)(4.*LOGD2);
1092 : /* 97/1223 below to ensure a good enough approximation of residue */
1093 70426 : cache_prime_quad(&GRHcheck, expi(D) < 16 ? 97: 1223, D);
1094 587589 : while (!quadGRHchk(D, &GRHcheck, high))
1095 : {
1096 517165 : low = high;
1097 517165 : high *= 2;
1098 : }
1099 517236 : while (high - low > 1)
1100 : {
1101 446811 : long test = (low+high)/2;
1102 446811 : if (quadGRHchk(D, &GRHcheck, test))
1103 234094 : high = test;
1104 : else
1105 212721 : low = test;
1106 : }
1107 70425 : if (high == LIMC0+1 && quadGRHchk(D, &GRHcheck, LIMC0))
1108 0 : LIMC2 = LIMC0;
1109 : else
1110 70425 : LIMC2 = high;
1111 70425 : if (LIMC2 > LIMCMAX) LIMC2 = LIMCMAX;
1112 70425 : LIMC0 = (long)(cbach*LOGD2);
1113 70425 : LIMC = cbach ? LIMC0 : LIMC2;
1114 70425 : LIMC = maxss(LIMC, nthidealquad(D, 2));
1115 :
1116 : /* LIMC = Max(cbach*(log D)^2, exp(sqrt(log D loglog D) / 8)) */
1117 70439 : START:
1118 70439 : missing_primes = NULL;
1119 : do
1120 : {
1121 70516 : if (!FIRST) LIMC = bnf_increase_LIMC(LIMC,LIMCMAX);
1122 70516 : if (DEBUGLEVEL>2 && LIMC > LIMC0)
1123 0 : err_printf("%s*** Bach constant: %f\n", FIRST?"":"\n", LIMC/LOGD2);
1124 70516 : FIRST = 0; set_avma(av);
1125 70516 : guncloneNULL(BQ.subFB);
1126 70516 : guncloneNULL(BQ.powsubFB);
1127 70516 : clearhash(BQ.hashtab);
1128 70517 : if (LIMC < cp) LIMC = cp;
1129 70517 : if (LIMC2 < LIMC) LIMC2 = LIMC;
1130 70517 : if (BQ.PRECREG) qfr_data_init(q.D, BQ.PRECREG, &q);
1131 :
1132 70517 : FBquad(&BQ, LIMC2, LIMC, &GRHcheck);
1133 70514 : if (DEBUGLEVEL>2) timer_printf(&T, "factor base");
1134 70514 : BQ.subFB = subFBquad(&BQ, q.D, lim + 0.5, minSFB);
1135 70516 : if (DEBUGLEVEL>2) timer_printf(&T, "subFBquad = %Ps",
1136 : vecpermute(BQ.FB, BQ.subFB));
1137 70516 : nsubFB = lg(BQ.subFB) - 1;
1138 : }
1139 70516 : while (nsubFB < (expi(D) > 15 ? 3 : 2));
1140 : /* invhr = 2^r1 (2pi)^r2 / sqrt(D) w ~ L(chi,1) / hR */
1141 70438 : invhr = gmul(dbltor((BQ.PRECREG?2.:M_PI)/sdrc),
1142 : compute_invresquad(&GRHcheck, LIMC));
1143 70438 : BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
1144 70440 : if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
1145 70440 : BQ.limhash = (LIMC & HIGHMASK)? (HIGHBIT>>1): LIMC*LIMC;
1146 :
1147 70440 : need = BQ.KC + RELSUP - 2;
1148 70440 : current = 0;
1149 70440 : W = NULL;
1150 70440 : sfb_trials = nreldep = nrelsup = 0;
1151 70440 : s = nsubFB + RELSUP;
1152 70440 : if (use_mpqs)
1153 28 : use_mpqs = mpqs_class_init(&H, D, BQ.KC);
1154 70440 : av2 = avma;
1155 :
1156 : do
1157 : {
1158 104632 : if ((nreldep & 3) == 1 || (nrelsup & 7) == 1) {
1159 28882 : if (DEBUGLEVEL>2) err_printf("*** Changing sub factor base\n");
1160 28882 : gunclone(BQ.subFB);
1161 28882 : gunclone(BQ.powsubFB);
1162 28882 : BQ.subFB = gclone(vecslice(BQ.vperm, 1, nsubFB));
1163 28882 : BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
1164 28882 : if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
1165 28882 : clearhash(BQ.hashtab);
1166 : }
1167 104632 : if (!use_mpqs) need += 2;
1168 104632 : mat = cgetg(need+1, t_MAT);
1169 104631 : extraC = cgetg(need+1, t_VEC);
1170 104631 : if (!W) { /* first time */
1171 70440 : C = extraC;
1172 70440 : triv = trivial_relations(&BQ, mat, C);
1173 70440 : if (DEBUGLEVEL>2) err_printf("KC = %ld, need %ld relations\n", BQ.KC, need);
1174 : } else {
1175 34191 : triv = 0;
1176 34191 : if (DEBUGLEVEL>2) err_printf("...need %ld more relations\n", need);
1177 : }
1178 104631 : if (BQ.PRECREG) {
1179 189833 : for (i = triv+1; i<=need; i++) {
1180 184688 : gel(mat,i) = zero_zv(BQ.KC);
1181 184688 : gel(extraC,i) = mkcomplex(cgetr(BQ.PRECREG), cgeti(3));
1182 : }
1183 5145 : real_relations(&BQ, need - triv, ¤t, s,LIMC,mat + triv,extraC + triv);
1184 : } else {
1185 1719218 : for (i = triv+1; i<=need; i++) {
1186 1619743 : gel(mat,i) = zero_zv(BQ.KC);
1187 1619732 : gel(extraC,i) = gen_0;
1188 : }
1189 99475 : if (use_mpqs)
1190 336 : mpqs_relations(&BQ, need - triv, ¤t, LIMC,mat + triv, &H, missing_primes);
1191 : else
1192 99139 : imag_relations(&BQ, need - triv, ¤t, LIMC,mat + triv);
1193 : }
1194 104631 : if (DEBUGLEVEL>2) timer_start(&T);
1195 104631 : if (!W)
1196 70439 : W = hnfspec_i(mat,BQ.vperm,&dep,&B,&C,nsubFB);
1197 : else
1198 34192 : W = hnfadd_i(W,BQ.vperm,&dep,&B,&C, mat,extraC);
1199 104632 : if (DEBUGLEVEL>2) timer_printf(&T, "hnfadd");
1200 104632 : need = BQ.KC - (lg(W)-1) - (lg(B)-1);
1201 104632 : (void)gc_all(av2, 4, &W,&C,&B,&dep);
1202 104632 : missing_primes = vecslice(BQ.vperm,1,need);
1203 104632 : if (need)
1204 : {
1205 12927 : if (++nreldep > 15 && cbach < 1) goto START;
1206 12927 : continue;
1207 : }
1208 :
1209 91705 : h = ZM_det_triangular(W);
1210 91704 : if (DEBUGLEVEL>2) err_printf("\n#### Tentative class number: %Ps\n", h);
1211 :
1212 91704 : switch(get_R(&BQ, C, (lg(C)-1) - (lg(B)-1) - (lg(W)-1), mulir(h,invhr), &R))
1213 : {
1214 0 : case fupb_PRECI:
1215 0 : BQ.PRECREG = precdbl(BQ.PRECREG);
1216 0 : FIRST = 1; goto START;
1217 :
1218 21279 : case fupb_RELAT:
1219 21279 : if (++nrelsup > MAXRELSUP)
1220 : {
1221 63 : if (++sfb_trials > SFB_MAX && cbach <= 1) goto START;
1222 49 : if (nsubFB < minss(10,BQ.KC)) nsubFB++;
1223 : }
1224 21265 : need = minss(BQ.KC, nrelsup);
1225 : }
1226 : }
1227 104616 : while (need);
1228 : /* DONE */
1229 70424 : if (!quad_be_honest(&BQ)) goto START;
1230 70424 : if (DEBUGLEVEL>2) timer_printf(&T, "be honest");
1231 70424 : clearhash(BQ.hashtab);
1232 70426 : free_GRHcheck(&GRHcheck);
1233 :
1234 70426 : gen = get_clgp(&BQ,W,&cyc);
1235 70424 : gunclone(BQ.subFB);
1236 70425 : gunclone(BQ.powsubFB);
1237 70425 : if (BQ.PRECREG)
1238 2156 : return mkvec5(h, cyc, gen, real_i(R), mpodd(imag_i(R)) ? gen_m1:gen_1);
1239 : else
1240 68269 : return mkvec4(h, cyc, gen, real_i(R));
1241 : }
1242 : GEN
1243 4424 : Buchquad(GEN D, double c, double c2, long prec)
1244 : {
1245 4424 : pari_sp av = avma;
1246 4424 : GEN z = Buchquad_i(D, c, c2, prec);
1247 4424 : return gc_GEN(av, z);
1248 : }
1249 :
1250 : GEN
1251 0 : buchimag(GEN D, GEN c, GEN c2, GEN REL)
1252 0 : { (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),0); }
1253 :
1254 : GEN
1255 0 : buchreal(GEN D, GEN flag, GEN c, GEN c2, GEN REL, long prec) {
1256 0 : if (signe(flag)) pari_err_IMPL("narrow class group");
1257 0 : (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),prec);
1258 : }
1259 :
1260 : GEN
1261 4424 : quadclassunit0(GEN x, long flag, GEN data, long prec)
1262 : {
1263 : long lx;
1264 4424 : double c1 = 0.0, c2 = 0.0;
1265 :
1266 4424 : if (!data) lx=1;
1267 : else
1268 : {
1269 28 : lx = lg(data);
1270 28 : if (typ(data)!=t_VEC) pari_err_TYPE("quadclassunit", data);
1271 28 : if (lx > 7) pari_err_DIM("quadclassunit [tech vector]");
1272 28 : if (lx > 3) lx = 3;
1273 : }
1274 4424 : switch(lx)
1275 : {
1276 21 : case 3: c2 = gtodouble(gel(data,2));
1277 28 : case 2: c1 = gtodouble(gel(data,1));
1278 : }
1279 4424 : if (flag) pari_err_IMPL("narrow class group");
1280 4424 : return Buchquad(x,c1,c2,prec);
1281 : }
1282 : GEN
1283 61164 : quadclassno(GEN D)
1284 : {
1285 61164 : pari_sp av = avma;
1286 61164 : GEN h = abgrp_get_no(Buchquad_i(D, 0, 0, 0));
1287 61169 : return icopy_avma(h, av);
1288 : }
1289 : long
1290 5008 : quadclassnos(long D)
1291 : {
1292 5008 : pari_sp av = avma;
1293 5008 : long h = itos(abgrp_get_no(Buchquad_i(stoi(D), 0, 0, 0)));
1294 5008 : return gc_long(av, h);
1295 : }
|
