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 2557036 : 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;
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 4035311 : init_form(struct buch_quad *B, GEN ex,
134 : GEN (*comp)(GEN,GEN,struct qfr_data *S))
135 : {
136 4035311 : long i, l = lg(B->powsubFB);
137 4035311 : GEN F = NULL;
138 22818424 : for (i=1; i<l; i++)
139 18791510 : if (ex[i])
140 : {
141 17619698 : GEN t = gmael(B->powsubFB,i,ex[i]);
142 17619698 : F = F? comp(F, t, B->q): t;
143 : }
144 4026914 : 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 11703431 : QFI_comp(GEN x, GEN y, struct qfr_data *S) { (void)S; return qfbcomp_i(x,y); }
151 :
152 : static GEN
153 158628 : qfi_factorback(struct buch_quad *B, GEN ex) { return init_form(B, ex, &QFI_comp); }
154 :
155 : static GEN
156 3678114 : random_form(struct buch_quad *B, GEN ex,
157 : GEN (*comp)(GEN,GEN, struct qfr_data *S))
158 : {
159 3678114 : long i, l = lg(ex);
160 3678114 : pari_sp av = avma;
161 : GEN F;
162 : for(;;)
163 : {
164 20499675 : for (i=1; i<l; i++) ex[i] = random_bits(RANDOM_BITS);
165 3679264 : if ((F = init_form(B, ex, comp))) return F;
166 401 : 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 3517068 : qfi_random(struct buch_quad *B,GEN ex) { return random_form(B, ex, &QFI_comp); }
173 :
174 : /*******************************************************************/
175 : /* */
176 : /* Common subroutines */
177 : /* */
178 : /*******************************************************************/
179 : long
180 2469 : bnf_increase_LIMC(long LIMC, long LIMCMAX)
181 : {
182 2469 : if (LIMC >= LIMCMAX) pari_err_BUG("Buchmann's algorithm");
183 2469 : if (LIMC <= LIMCMAX/40) /* cbach <= 0.3 */
184 208 : LIMC *= 2;
185 2261 : else if (LIMCMAX < 60) /* \Delta_K <= 9 */
186 0 : LIMC++;
187 : else
188 2261 : LIMC += LIMCMAX / 60; /* cbach += 0.2 */
189 2469 : if (LIMC > LIMCMAX) LIMC = LIMCMAX;
190 2469 : return LIMC;
191 : }
192 :
193 : /* Is |q| <= p ? */
194 : static int
195 73056 : isless_iu(GEN q, ulong p) {
196 73056 : long l = lgefint(q);
197 73056 : return l==2 || (l == 3 && uel(q,2) <= p);
198 : }
199 :
200 : static long
201 5007886 : factorquad(struct buch_quad *B, GEN f, long nFB, ulong limp)
202 : {
203 : ulong X;
204 5007886 : long i, lo = 0;
205 5007886 : GEN x = gel(f,1), FB = B->FB, P = B->primfact, E = B->exprimfact;
206 :
207 5080550 : for (i=1; lgefint(x) > 3; i++)
208 : {
209 73056 : ulong p = uel(FB,i), r;
210 73056 : GEN q = absdiviu_rem(x, p, &r);
211 73056 : if (!r)
212 : {
213 3822 : long k = 0;
214 6154 : do { k++; x = q; q = absdiviu_rem(x, p, &r); } while (!r);
215 3822 : lo++; P[lo] = p; E[lo] = k;
216 : }
217 73056 : if (isless_iu(q,p)) {
218 1 : if (lgefint(x) == 3) { X = uel(x,2); goto END; }
219 391 : return 0;
220 : }
221 73055 : if (i == nFB) return 0;
222 : }
223 5007494 : X = uel(x,2);
224 5007494 : if (X == 1) { P[0] = 0; return 1; }
225 123368751 : for (;; i++)
226 123368751 : { /* single precision affair, split for efficiency */
227 128361299 : ulong p = uel(FB,i);
228 128361299 : ulong q = X / p, r = X % p; /* gcc makes a single div */
229 128361299 : if (!r)
230 : {
231 5987314 : long k = 0;
232 7213536 : do { k++; X = q; q = X / p; r = X % p; } while (!r);
233 5987314 : lo++; P[lo] = p; E[lo] = k;
234 : }
235 128361299 : if (q <= p) break;
236 124038157 : if (i == nFB) return 0;
237 : }
238 4323143 : END:
239 4323143 : if (X > B->limhash) return 0;
240 4323115 : if (X != 1 && X <= limp)
241 : {
242 1434055 : if (B->badprim && ugcdui(X, B->badprim) > 1) return 0;
243 1433740 : lo++; P[lo] = X; E[lo] = 1;
244 1433740 : X = 1;
245 : }
246 4322800 : P[0] = lo; return X;
247 : }
248 :
249 : /* Check for a "large prime relation" involving q; q may not be prime */
250 : static long *
251 2557071 : largeprime(struct buch_quad *B, long q, GEN ex, long np, long nrho)
252 : {
253 2557071 : const long hashv = hash(q);
254 2557009 : long *pt, i, l = lg(B->subFB);
255 :
256 2557009 : for (pt = B->hashtab[hashv]; ; pt = (long*) pt[0])
257 : {
258 2779374 : if (!pt)
259 : {
260 2358208 : pt = (long*) pari_malloc((l+3) * sizeof(long));
261 2358387 : *pt++ = nrho; /* nrho = pt[-3] */
262 2358387 : *pt++ = np; /* np = pt[-2] */
263 2358387 : *pt++ = q; /* q = pt[-1] */
264 2358387 : pt[0] = (long)B->hashtab[hashv];
265 14948468 : for (i=1; i<l; i++) pt[i]=ex[i];
266 2358387 : B->hashtab[hashv]=pt; return NULL;
267 : }
268 421166 : if (pt[-1] == q) break;
269 : }
270 242237 : for(i=1; i<l; i++)
271 238381 : if (pt[i] != ex[i]) return pt;
272 3856 : return (pt[-2]==np)? NULL: pt;
273 : }
274 :
275 : static void
276 167780 : clearhash(long **hash)
277 : {
278 : long *pt;
279 : long i;
280 171932585 : for (i=0; i<HASHT; i++) {
281 174123318 : for (pt = hash[i]; pt; ) {
282 2358513 : void *z = (void*)(pt - 3);
283 2358513 : pt = (long*) pt[0]; pari_free(z);
284 : }
285 171764805 : hash[i] = NULL;
286 : }
287 168865 : }
288 :
289 : /* last prime stored */
290 : ulong
291 0 : GRH_last_prime(GRHcheck_t *S) { return (S->primes + S->nprimes-1)->p; }
292 : /* ensure that S->primes can hold at least nb primes */
293 : void
294 398203 : GRH_ensure(GRHcheck_t *S, long nb)
295 : {
296 398203 : if (S->maxprimes <= nb)
297 : {
298 0 : do S->maxprimes *= 2; while (S->maxprimes <= nb);
299 0 : pari_realloc_ip((void**)&S->primes, S->maxprimes*sizeof(*S->primes));
300 : }
301 398203 : }
302 : /* cache data for all primes up to the LIM */
303 : static void
304 1166315 : cache_prime_quad(GRHcheck_t *S, ulong LIM, GEN D)
305 : {
306 : GRHprime_t *pr;
307 : long nb;
308 :
309 1166315 : if (S->limp >= LIM) return;
310 71490 : nb = (long)primepi_upper_bound((double)LIM); /* #{p <= LIM} <= nb */
311 71490 : GRH_ensure(S, nb+1); /* room for one extra prime */
312 71490 : for (pr = S->primes + S->nprimes;;)
313 12627974 : {
314 12699464 : ulong p = u_forprime_next(&(S->P));
315 12698893 : pr->p = p;
316 12698893 : pr->logp = log((double)p);
317 12698893 : pr->dec = (GEN)kroiu(D,p);
318 12699466 : S->nprimes++;
319 12699466 : pr++;
320 : /* store up to nextprime(LIM) included */
321 12699466 : if (p >= LIM) { S->limp = p; break; }
322 : }
323 : }
324 :
325 : static GEN
326 69411 : compute_invresquad(GRHcheck_t *S, long LIMC)
327 : {
328 69411 : pari_sp av = avma;
329 69411 : GEN invres = real_1(DEFAULTPREC);
330 69411 : double limp = log((double)LIMC) / 2;
331 : GRHprime_t *pr;
332 : long i;
333 12677616 : for (pr = S->primes, i = S->nprimes; i > 0; pr++, i--)
334 : {
335 12609055 : long s = (long)pr->dec;
336 12609055 : if (s)
337 : {
338 12480850 : ulong p = pr->p;
339 12480850 : if (s > 0 || pr->logp <= limp)
340 : /* Both p and P contribute */
341 6344289 : invres = mulur(p - s, divru(invres, p));
342 6136561 : else if (s<0)
343 : /* Only p contributes */
344 6136649 : invres = mulur(p, divru(invres, p - 1));
345 : }
346 : }
347 68561 : return gerepileuptoleaf(av, invres);
348 : }
349 :
350 : /* p | conductor of order of disc D ? */
351 : static int
352 388365 : is_bad(GEN D, ulong p)
353 : {
354 388365 : pari_sp av = avma;
355 388365 : if (p == 2)
356 : {
357 88514 : long r = mod16(D) >> 1;
358 88513 : if (r && signe(D) < 0) r = 8-r;
359 88513 : return (r < 4);
360 : }
361 299851 : return gc_bool(av, dvdii(D, sqru(p))); /* p^2 | D ? */
362 : }
363 :
364 : /* returns the n-th suitable ideal for the factorbase */
365 : static long
366 69410 : nthidealquad(GEN D, long n)
367 : {
368 69410 : pari_sp av = avma;
369 : forprime_t S;
370 : ulong p;
371 69410 : (void)u_forprime_init(&S, 2, ULONG_MAX);
372 353782 : while ((p = u_forprime_next(&S)) && n > 0)
373 284378 : if (!is_bad(D, p) && kroiu(D, p) >= 0) n--;
374 69399 : return gc_long(av, p);
375 : }
376 :
377 : static int
378 1027258 : quadGRHchk(GEN D, GRHcheck_t *S, ulong LIMC)
379 : {
380 1027258 : double logC = log((double)LIMC), SA = 0, SB = 0;
381 : long i;
382 :
383 1027258 : cache_prime_quad(S, LIMC, D);
384 1027250 : for (i = 0;; i++)
385 29318447 : {
386 30345697 : GRHprime_t *pr = S->primes+i;
387 30345697 : ulong p = pr->p;
388 : long M;
389 : double logNP, q, A, B;
390 30345697 : if (p > LIMC) break;
391 29405958 : if ((long)pr->dec < 0)
392 : {
393 14687915 : logNP = 2 * pr->logp;
394 14687915 : q = 1/(double)p;
395 : }
396 : else
397 : {
398 14718043 : logNP = pr->logp;
399 14718043 : q = 1/sqrt((double)p);
400 : }
401 29405958 : A = logNP * q; B = logNP * A; M = (long)(logC/logNP);
402 29405958 : if (M > 1)
403 : {
404 2493095 : double inv1_q = 1 / (1-q);
405 2493095 : A *= (1 - pow(q, (double) M)) * inv1_q;
406 2493095 : B *= (1 - pow(q, (double) M)*(M+1 - M*q)) * inv1_q * inv1_q;
407 : }
408 29405958 : if ((long)pr->dec>0) { SA += 2*A;SB += 2*B; } else { SA += A; SB += B; }
409 29405958 : if (p == LIMC) break;
410 : }
411 1027250 : return GRHok(S, logC, SA, SB);
412 : }
413 :
414 : /* C2 >= C1; create B->FB, B->numFB; set B->badprim. Return L(kro_D, 1) */
415 : static void
416 69656 : FBquad(struct buch_quad *B, ulong C2, ulong C1, GRHcheck_t *S)
417 : {
418 69656 : GEN D = B->q->D;
419 : long i;
420 : pari_sp av;
421 : GRHprime_t *pr;
422 :
423 69656 : cache_prime_quad(S, C2, D);
424 69656 : pr = S->primes;
425 69656 : B->numFB = cgetg(C2+1, t_VECSMALL);
426 69656 : B->FB = cgetg(C2+1, t_VECSMALL);
427 69656 : av = avma;
428 69656 : B->KC = 0; i = 0;
429 69656 : B->badprim = gen_1;
430 2793076 : for (;; pr++) /* p <= C2 */
431 2793076 : {
432 2862732 : ulong p = pr->p;
433 2862732 : if (!B->KC && p > C1) B->KC = i;
434 2862732 : if (p > C2) break;
435 2800607 : switch ((long)pr->dec)
436 : {
437 1380633 : case -1: break; /* inert */
438 103989 : case 0: /* ramified */
439 103989 : if (is_bad(D, p)) { B->badprim = muliu(B->badprim, p); break; }
440 : /* fall through */
441 : default: /* split */
442 1419946 : i++; B->numFB[p] = i; B->FB[i] = p; break;
443 : }
444 2800607 : if (p == C2)
445 : {
446 7531 : if (!B->KC) B->KC = i;
447 7531 : break;
448 : }
449 : }
450 69656 : B->KC2 = i;
451 69656 : setlg(B->FB, B->KC2+1);
452 69656 : if (B->badprim != gen_1)
453 21 : B->badprim = gerepileuptoint(av, B->badprim);
454 : else
455 : {
456 69635 : B->badprim = NULL; set_avma(av);
457 : }
458 69656 : }
459 :
460 : /* create B->vperm, return B->subFB */
461 : static GEN
462 69656 : subFBquad(struct buch_quad *B, GEN D, double PROD, long minSFB)
463 : {
464 69656 : long i, j, lgsub = 1, ino = 1, lv = B->KC+1;
465 69656 : double prod = 1.;
466 : pari_sp av;
467 : GEN no;
468 :
469 69656 : B->vperm = cgetg(lv, t_VECSMALL);
470 69656 : av = avma;
471 69656 : no = cgetg(lv, t_VECSMALL);
472 332436 : for (j = 1; j < lv; j++)
473 : {
474 332184 : ulong p = uel(B->FB,j);
475 332184 : if (!umodiu(D, p)) no[ino++] = j; /* ramified */
476 : else
477 : {
478 252878 : B->vperm[lgsub++] = j;
479 252878 : prod *= p;
480 252878 : if (lgsub > minSFB && prod > PROD) break;
481 : }
482 : }
483 : /* lgsub >= 1 otherwise quadGRHchk is false */
484 69656 : i = lgsub;
485 148962 : for (j = 1; j < ino;i++,j++) B->vperm[i] = no[j];
486 1156929 : for ( ; i < lv; i++) B->vperm[i] = i;
487 69656 : no = gclone(vecslice(B->vperm, 1, lgsub-1));
488 69656 : set_avma(av); return no;
489 : }
490 :
491 : /* assume n >= 1, x[i][j] = B->subFB[i]^j, for j = 1..n */
492 : static GEN
493 98127 : powsubFBquad(struct buch_quad *B, long n)
494 : {
495 98127 : pari_sp av = avma;
496 98127 : long i,j, l = lg(B->subFB);
497 98127 : GEN F, y, x = cgetg(l, t_VEC), D = B->q->D;
498 :
499 98128 : if (B->PRECREG) /* real */
500 : {
501 39088 : for (i=1; i<l; i++)
502 : {
503 34055 : F = qfr5_pf(B->q, B->FB[B->subFB[i]], B->PRECREG);
504 34055 : y = cgetg(n+1, t_VEC); gel(x,i) = y;
505 34055 : gel(y,1) = F;
506 544880 : for (j=2; j<=n; j++) gel(y,j) = QFR5_comp(gel(y,j-1), F, B->q);
507 : }
508 : }
509 : else /* imaginary */
510 : {
511 423941 : for (i=1; i<l; i++)
512 : {
513 330867 : F = qfi_pf(D, B->FB[B->subFB[i]]);
514 330860 : y = cgetg(n+1, t_VEC); gel(x,i) = y;
515 331861 : gel(y,1) = F;
516 5285216 : for (j=2; j<=n; j++) gel(y,j) = qfbcomp_i(gel(y,j-1), F);
517 : }
518 : }
519 98107 : x = gclone(x); set_avma(av); return x;
520 : }
521 :
522 : static void
523 97650 : sub_fact(struct buch_quad *B, GEN col, GEN F)
524 : {
525 97650 : GEN b = gel(F,2);
526 : long i;
527 207793 : for (i=1; i<=B->primfact[0]; i++)
528 : {
529 110143 : ulong p = B->primfact[i], k = B->numFB[p];
530 110143 : long e = B->exprimfact[i];
531 110143 : if (umodiu(b, p<<1) > p) e = -e;
532 110143 : col[k] -= e;
533 : }
534 97650 : }
535 : static void
536 1877196 : add_fact(struct buch_quad *B, GEN col, GEN F)
537 : {
538 1877196 : GEN b = gel(F,2);
539 : long i;
540 5903992 : for (i=1; i<=B->primfact[0]; i++)
541 : {
542 4026718 : ulong p = B->primfact[i], k = B->numFB[p];
543 4026718 : long e = B->exprimfact[i];
544 4026718 : if (umodiu(b, p<<1) > p) e = -e;
545 4026796 : col[k] += e;
546 : }
547 1877274 : }
548 :
549 : static GEN
550 69411 : get_clgp(struct buch_quad *B, GEN W, GEN *ptD)
551 : {
552 69411 : GEN res, init, u1, D = ZM_snf_group(W,NULL,&u1);
553 69410 : long i, j, l = lg(W), c = lg(D);
554 :
555 69410 : res=cgetg(c,t_VEC); init = cgetg(l,t_VEC);
556 214513 : for (i=1; i<l; i++) gel(init,i) = primeform_u(B->q->D, B->FB[B->vperm[i]]);
557 194167 : for (j=1; j<c; j++)
558 : {
559 124756 : GEN g = NULL;
560 124756 : if (signe(B->q->D) > 0)
561 : {
562 13293 : for (i=1; i<l; i++)
563 : {
564 9583 : GEN t, u = gcoeff(u1,i,j);
565 9583 : if (!signe(u)) continue;
566 4501 : t = qfr3_pow(gel(init,i), u, B->q);
567 4501 : g = g? qfr3_comp(g, t, B->q): t;
568 : }
569 3710 : g = qfr3_to_qfr(qfr3_canon_safe(qfr3_red(g, B->q), B->q), B->q->D);
570 : }
571 : else
572 : {
573 421524 : for (i=1; i<l; i++)
574 : {
575 300477 : GEN t, u = gcoeff(u1,i,j);
576 300477 : if (!signe(u)) continue;
577 202686 : t = powgi(gel(init,i), u);
578 202689 : g = g? qfbcomp_i(g, t): t;
579 : }
580 : }
581 124757 : gel(res,j) = g;
582 : }
583 69411 : *ptD = D; return res;
584 : }
585 :
586 : static long
587 69411 : trivial_relations(struct buch_quad *B, GEN mat, GEN C)
588 : {
589 69411 : long i, j = 0;
590 69411 : GEN col, D = B->q->D;
591 1487989 : for (i = 1; i <= B->KC; i++)
592 : { /* ramified prime ==> trivial relation */
593 1418570 : if (umodiu(D, B->FB[i])) continue;
594 103325 : col = zero_zv(B->KC);
595 103331 : col[i] = 2; j++;
596 103331 : gel(mat,j) = col;
597 103331 : gel(C,j) = gen_0;
598 : }
599 69419 : return j;
600 : }
601 :
602 : static void
603 0 : dbg_all(pari_timer *T, const char *phase, long s, long n)
604 : {
605 0 : err_printf("\n");
606 0 : timer_printf(T, "%s rel [#rel/#test = %ld/%ld]", phase,s,n);
607 0 : }
608 :
609 : /* Imaginary Quadratic fields */
610 :
611 : static void
612 97761 : imag_relations(struct buch_quad *B, long need, long *pc, ulong LIMC, GEN mat)
613 : {
614 : pari_timer T;
615 97761 : long lgsub = lg(B->subFB), current = *pc, nbtest = 0, s = 0;
616 : long i, fpc;
617 : pari_sp av;
618 97761 : GEN col, form, ex = cgetg(lgsub, t_VECSMALL);
619 :
620 97761 : if (!current) current = 1;
621 97761 : if (DEBUGLEVEL>2) timer_start(&T);
622 97761 : av = avma;
623 : for(;;)
624 : {
625 3615012 : if (s >= need) break;
626 3517251 : set_avma(av);
627 3516719 : form = qfi_random(B,ex);
628 3514995 : form = qfbcomp_i(form, qfi_pf(B->q->D, B->FB[current]));
629 3514970 : nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
630 3517256 : if (!fpc)
631 : {
632 287099 : if (DEBUGLEVEL>3) err_printf(".");
633 287099 : if ((nbtest & 0xff) == 0 && ++current > B->KC) current = 1;
634 287099 : continue;
635 : }
636 3230157 : if (fpc > 1)
637 : {
638 1793010 : long *fpd = largeprime(B,fpc,ex,current,0);
639 : ulong b1, b2, p;
640 : GEN form2;
641 1793135 : if (!fpd)
642 : {
643 1634536 : if (DEBUGLEVEL>3) err_printf(".");
644 1634529 : continue;
645 : }
646 158599 : form2 = qfbcomp_i(qfi_factorback(B,fpd), qfi_pf(B->q->D, B->FB[fpd[-2]]));
647 158625 : p = fpc << 1;
648 158625 : b1 = umodiu(gel(form2,2), p);
649 158626 : b2 = umodiu(gel(form,2), p);
650 158632 : if (b1 != b2 && b1+b2 != p) continue;
651 :
652 158597 : col = gel(mat,++s);
653 158597 : add_fact(B,col, form);
654 158596 : (void)factorquad(B,form2,B->KC,LIMC);
655 158597 : if (b1==b2)
656 : {
657 408950 : for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
658 79415 : sub_fact(B, col, form2); col[fpd[-2]]++;
659 : }
660 : else
661 : {
662 408008 : for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
663 79182 : add_fact(B, col, form2); col[fpd[-2]]--;
664 : }
665 158597 : if (DEBUGLEVEL>2) err_printf(" %ldP",s);
666 : }
667 : else
668 : {
669 1437147 : col = gel(mat,++s);
670 6919341 : for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
671 1437147 : add_fact(B, col, form);
672 1437260 : if (DEBUGLEVEL>2) err_printf(" %ld",s);
673 : }
674 1595588 : col[current]--;
675 1595588 : if (++current > B->KC) current = 1;
676 : }
677 97761 : if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
678 97761 : *pc = current;
679 97761 : }
680 :
681 : static int
682 7 : imag_be_honest(struct buch_quad *B)
683 : {
684 7 : long p, fpc, s = B->KC, nbtest = 0;
685 7 : GEN F, ex = cgetg(lg(B->subFB), t_VECSMALL);
686 7 : pari_sp av = avma;
687 :
688 476 : while (s<B->KC2)
689 : {
690 469 : p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
691 469 : F = qfbcomp_i(qfi_pf(B->q->D, p), qfi_random(B, ex));
692 469 : fpc = factorquad(B,F,s,p-1);
693 469 : if (fpc == 1) { nbtest=0; s++; }
694 : else
695 322 : if (++nbtest > 40) return 0;
696 469 : set_avma(av);
697 : }
698 7 : return 1;
699 : }
700 :
701 : /* Real Quadratic fields */
702 :
703 : static void
704 5124 : real_relations(struct buch_quad *B, long need, long *pc, long lim, ulong LIMC, GEN mat, GEN C)
705 : {
706 : pari_timer T;
707 5124 : long lgsub = lg(B->subFB), prec = B->PRECREG, current = *pc, nbtest=0, s=0;
708 : long i, fpc, endcycle, rhoacc, rho;
709 : /* in a 2nd phase, don't include FB[current] but run along the cyle
710 : * ==> get more units */
711 5124 : int first = (current == 0);
712 : pari_sp av, av1;
713 5124 : GEN d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL);
714 :
715 5124 : if (DEBUGLEVEL>2) timer_start(&T);
716 5124 : if (!current) current = 1;
717 5124 : if (lim > need) lim = need;
718 5124 : av = avma;
719 : for(;;)
720 : {
721 166236 : if (s >= need) break;
722 161112 : if (first && s >= lim) {
723 2142 : first = 0;
724 2142 : if (DEBUGLEVEL>2) dbg_all(&T, "initial", s, nbtest);
725 : }
726 161112 : set_avma(av); form = qfr3_random(B, ex);
727 161112 : if (!first)
728 158935 : form = QFR3_comp(form, qfr3_pf(B->q, B->FB[current]), B->q);
729 161112 : av1 = avma;
730 161112 : form0 = form; form1 = NULL;
731 161112 : endcycle = rhoacc = 0;
732 161112 : rho = -1;
733 :
734 1297569 : CYCLE:
735 1297569 : if (endcycle || rho > 5000)
736 : {
737 21 : if (++current > B->KC) current = 1;
738 21 : continue;
739 : }
740 1297548 : if (gc_needed(av,1))
741 : {
742 0 : if(DEBUGMEM>1) pari_warn(warnmem,"real_relations");
743 0 : gerepileall(av1, form1? 2: 1, &form, &form1);
744 : }
745 1297548 : if (rho < 0) rho = 0; /* first time in */
746 : else
747 : {
748 1136436 : form = qfr3_rho(form, B->q); rho++;
749 1136436 : rhoacc++;
750 1136436 : if (first)
751 444430 : endcycle = (absequalii(gel(form,1),gel(form0,1))
752 222215 : && equalii(gel(form,2),gel(form0,2)));
753 : else
754 : {
755 914221 : if (absequalii(gel(form,1), gel(form,3))) /* a = -c */
756 : {
757 0 : if (absequalii(gel(form,1),gel(form0,1)) &&
758 0 : equalii(gel(form,2),gel(form0,2))) continue;
759 0 : form = qfr3_rho(form, B->q); rho++;
760 0 : rhoacc++;
761 : }
762 : else
763 914221 : { setsigne(form[1],1); setsigne(form[3],-1); }
764 914277 : if (equalii(gel(form,1),gel(form0,1)) &&
765 56 : equalii(gel(form,2),gel(form0,2))) continue;
766 : }
767 : }
768 1297548 : nbtest++; fpc = factorquad(B,form,B->KC,LIMC);
769 1297548 : if (!fpc)
770 : {
771 385511 : if (DEBUGLEVEL>3) err_printf(".");
772 385511 : goto CYCLE;
773 : }
774 912037 : if (fpc > 1)
775 : { /* look for Large Prime relation */
776 764113 : long *fpd = largeprime(B,fpc,ex,first? 0: current,rhoacc);
777 : ulong b1, b2, p;
778 : GEN form2;
779 764113 : if (!fpd)
780 : {
781 727783 : if (DEBUGLEVEL>3) err_printf(".");
782 727783 : goto CYCLE;
783 : }
784 36330 : if (!form1)
785 : {
786 36330 : form1 = qfr5_factorback(B,ex);
787 36330 : if (!first)
788 36330 : form1 = QFR5_comp(form1, qfr5_pf(B->q, B->FB[current], prec), B->q);
789 : }
790 36330 : form1 = qfr5_rho_pow(form1, rho, B->q);
791 36330 : rho = 0;
792 :
793 36330 : form2 = qfr5_factorback(B,fpd);
794 36330 : if (fpd[-2])
795 23954 : form2 = QFR5_comp(form2, qfr5_pf(B->q, B->FB[fpd[-2]], prec), B->q);
796 36330 : form2 = qfr5_rho_pow(form2, fpd[-3], B->q);
797 36330 : if (!absequalii(gel(form2,1),gel(form2,3)))
798 : {
799 36330 : setsigne(form2[1], 1);
800 36330 : setsigne(form2[3],-1);
801 : }
802 36330 : p = fpc << 1;
803 36330 : b1 = umodiu(gel(form2,2), p);
804 36330 : b2 = umodiu(gel(form1,2), p);
805 36330 : if (b1 != b2 && b1+b2 != p) goto CYCLE;
806 :
807 36330 : col = gel(mat,++s);
808 36330 : add_fact(B, col, form1);
809 36330 : (void)factorquad(B,form2,B->KC,LIMC);
810 36330 : if (b1==b2)
811 : {
812 139377 : for (i=1; i<lgsub; i++) col[B->subFB[i]] += fpd[i]-ex[i];
813 18235 : sub_fact(B,col, form2);
814 18235 : if (fpd[-2]) col[fpd[-2]]++;
815 18235 : d = qfr5_dist(subii(gel(form1,4),gel(form2,4)),
816 18235 : divrr(gel(form1,5),gel(form2,5)), prec);
817 : }
818 : else
819 : {
820 138670 : for (i=1; i<lgsub; i++) col[B->subFB[i]] += -fpd[i]-ex[i];
821 18095 : add_fact(B, col, form2);
822 18095 : if (fpd[-2]) col[fpd[-2]]--;
823 18095 : d = qfr5_dist(addii(gel(form1,4),gel(form2,4)),
824 18095 : mulrr(gel(form1,5),gel(form2,5)), prec);
825 : }
826 36330 : if (DEBUGLEVEL>2) err_printf(" %ldP",s);
827 : }
828 : else
829 : { /* standard relation */
830 147924 : if (!form1)
831 : {
832 124782 : form1 = qfr5_factorback(B, ex);
833 124782 : if (!first)
834 122605 : form1 = QFR5_comp(form1, qfr5_pf(B->q, B->FB[current], prec), B->q);
835 : }
836 147924 : form1 = qfr5_rho_pow(form1, rho, B->q);
837 147924 : rho = 0;
838 :
839 147924 : col = gel(mat,++s);
840 1138494 : for (i=1; i<lgsub; i++) col[B->subFB[i]] = -ex[i];
841 147924 : add_fact(B, col, form1);
842 147924 : d = qfr5_dist(gel(form1,4), gel(form1,5), prec);
843 147924 : if (DEBUGLEVEL>2) err_printf(" %ld",s);
844 : }
845 184254 : affrr(d, gel(C,s));
846 184254 : if (first)
847 : {
848 25319 : if (s >= lim) continue;
849 23163 : goto CYCLE;
850 : }
851 : else
852 : {
853 158935 : col[current]--;
854 158935 : if (++current > B->KC) current = 1;
855 : }
856 : }
857 5124 : if (DEBUGLEVEL>2) dbg_all(&T, "random", s, nbtest);
858 5124 : *pc = current;
859 5124 : }
860 :
861 : static int
862 7 : real_be_honest(struct buch_quad *B)
863 : {
864 7 : long p, fpc, s = B->KC, nbtest = 0;
865 7 : GEN F,F0, ex = cgetg(lg(B->subFB), t_VECSMALL);
866 7 : pari_sp av = avma;
867 :
868 28 : while (s<B->KC2)
869 : {
870 21 : p = B->FB[s+1]; if (DEBUGLEVEL>2) err_printf(" %ld",p);
871 21 : F = QFR3_comp(qfr3_random(B, ex), qfr3_pf(B->q, p), B->q);
872 21 : for (F0 = F;;)
873 : {
874 42 : fpc = factorquad(B,F,s,p-1);
875 42 : if (fpc == 1) { nbtest=0; s++; break; }
876 21 : if (++nbtest > 40) return 0;
877 21 : F = qfr3_canon(qfr3_rho(F, B->q), B->q);
878 21 : if (equalii(gel(F,1),gel(F0,1))
879 0 : && equalii(gel(F,2),gel(F0,2))) break;
880 : }
881 21 : set_avma(av);
882 : }
883 7 : return 1;
884 : }
885 :
886 : static GEN
887 32025 : gcdreal(GEN a,GEN b)
888 : {
889 32025 : if (!signe(a)) return mpabs_shallow(b);
890 31361 : if (!signe(b)) return mpabs_shallow(a);
891 31213 : if (typ(a)==t_INT)
892 : {
893 0 : if (typ(b)==t_INT) return gcdii(a,b);
894 0 : a = itor(a, lg(b));
895 : }
896 31213 : else if (typ(b)==t_INT)
897 0 : b = itor(b, lg(a));
898 31213 : if (expo(a)<-5) return absr(b);
899 11305 : if (expo(b)<-5) return absr(a);
900 8456 : a = absr(a); b = absr(b);
901 18679 : while (expo(b) >= -5 && signe(b))
902 : {
903 : long e;
904 10223 : GEN r, q = gcvtoi(divrr(a,b),&e);
905 10223 : if (e > 0) return NULL;
906 10223 : r = subrr(a, mulir(q,b)); a = b; b = r;
907 : }
908 8456 : return mpabs_shallow(a);
909 : }
910 :
911 : static int
912 90005 : get_R(struct buch_quad *B, GEN C, long sreg, GEN z, GEN *ptR)
913 : {
914 90005 : GEN R = gen_1;
915 : double c;
916 : long i;
917 :
918 90005 : if (B->PRECREG)
919 : {
920 2877 : R = mpabs_shallow(gel(C,1));
921 34902 : for (i=2; i<=sreg; i++)
922 : {
923 32025 : R = gcdreal(gel(C,i), R);
924 32025 : if (!R) return fupb_PRECI;
925 : }
926 2877 : if (gexpo(R) <= -3)
927 : {
928 0 : if (DEBUGLEVEL>2) err_printf("regulator is zero.\n");
929 0 : return fupb_RELAT;
930 : }
931 2877 : if (DEBUGLEVEL>2) err_printf("#### Tentative regulator: %Ps\n",R);
932 : }
933 90005 : c = gtodouble(gmul(z, R));
934 90005 : if (c < 0.8 || c > 1.3) return fupb_RELAT;
935 69409 : *ptR = R; return fupb_NONE;
936 : }
937 :
938 : static int
939 69409 : quad_be_honest(struct buch_quad *B)
940 : {
941 : int r;
942 69409 : if (B->KC2 <= B->KC) return 1;
943 14 : if (DEBUGLEVEL>2)
944 0 : err_printf("be honest for primes from %ld to %ld\n", B->FB[B->KC+1],B->FB[B->KC2]);
945 14 : r = B->PRECREG? real_be_honest(B): imag_be_honest(B);
946 14 : if (DEBUGLEVEL>2) err_printf("\n");
947 14 : return r;
948 : }
949 :
950 : static GEN
951 69523 : Buchquad_i(GEN D, double cbach, double cbach2, long prec)
952 : {
953 69523 : const long MAXRELSUP = 7, SFB_MAX = 3;
954 : pari_timer T;
955 : pari_sp av, av2;
956 69523 : const long RELSUP = 5;
957 : long i, s, current, triv, sfb_trials, nrelsup, nreldep, need, nsubFB, minSFB;
958 : ulong low, high, LIMC0, LIMC, LIMC2, LIMCMAX, cp;
959 69523 : GEN W, cyc, gen, dep, mat, C, extraC, B, R, invhr, h = NULL; /*-Wall*/
960 : double drc, sdrc, lim, LOGD, LOGD2;
961 : GRHcheck_t GRHcheck;
962 : struct qfr_data q;
963 : struct buch_quad BQ;
964 69523 : int FIRST = 1;
965 :
966 69523 : check_quaddisc(D, &s, /*junk*/&i, "Buchquad");
967 69523 : R = NULL; /* -Wall */
968 69523 : BQ.q = &q;
969 69523 : q.D = D;
970 69523 : if (s < 0)
971 : {
972 67367 : if (abscmpiu(q.D,4) <= 0)
973 112 : retmkvec4(gen_1, cgetg(1,t_VEC), cgetg(1,t_VEC), gen_1);
974 67255 : prec = BQ.PRECREG = 0;
975 : } else {
976 2156 : BQ.PRECREG = maxss(prec+EXTRAPRECWORD, nbits2prec(2*expi(q.D) + 128));
977 : }
978 69411 : if (DEBUGLEVEL>2) timer_start(&T);
979 69411 : BQ.primfact = new_chunk(100);
980 69411 : BQ.exprimfact = new_chunk(100);
981 69411 : BQ.hashtab = (long**) new_chunk(HASHT);
982 71141298 : for (i=0; i<HASHT; i++) BQ.hashtab[i] = NULL;
983 :
984 69411 : drc = fabs(gtodouble(q.D));
985 69410 : LOGD = log(drc);
986 69410 : LOGD2 = LOGD * LOGD;
987 :
988 69410 : sdrc = lim = sqrt(drc);
989 69410 : if (!BQ.PRECREG) lim /= sqrt(3.);
990 69410 : cp = (ulong)exp(sqrt(LOGD*log(LOGD)/8.0));
991 69410 : if (cp < 20) cp = 20;
992 69410 : if (cbach > 6.) {
993 0 : if (cbach2 < cbach) cbach2 = cbach;
994 0 : cbach = 6.;
995 : }
996 69410 : if (cbach < 0.)
997 0 : pari_err_DOMAIN("Buchquad","Bach constant","<",gen_0,dbltor(cbach));
998 69410 : av = avma;
999 69410 : BQ.powsubFB = BQ.subFB = NULL;
1000 69410 : minSFB = (expi(D) > 15)? 3: 2;
1001 69410 : init_GRHcheck(&GRHcheck, 2, BQ.PRECREG? 2: 0, LOGD);
1002 69410 : high = low = LIMC0 = maxss((long)(cbach2*LOGD2), 1);
1003 69410 : LIMCMAX = (long)(6.*LOGD2);
1004 : /* 97/1223 below to ensure a good enough approximation of residue */
1005 69410 : cache_prime_quad(&GRHcheck, expi(D) < 16 ? 97: 1223, D);
1006 583014 : while (!quadGRHchk(D, &GRHcheck, high))
1007 : {
1008 513603 : low = high;
1009 513603 : high *= 2;
1010 : }
1011 513673 : while (high - low > 1)
1012 : {
1013 444273 : long test = (low+high)/2;
1014 444273 : if (quadGRHchk(D, &GRHcheck, test))
1015 232609 : high = test;
1016 : else
1017 211661 : low = test;
1018 : }
1019 69400 : if (high == LIMC0+1 && quadGRHchk(D, &GRHcheck, LIMC0))
1020 0 : LIMC2 = LIMC0;
1021 : else
1022 69400 : LIMC2 = high;
1023 69400 : if (LIMC2 > LIMCMAX) LIMC2 = LIMCMAX;
1024 69400 : LIMC0 = (long)(cbach*LOGD2);
1025 69400 : LIMC = cbach ? LIMC0 : LIMC2;
1026 69400 : LIMC = maxss(LIMC, nthidealquad(D, 2));
1027 :
1028 : /* LIMC = Max(cbach*(log D)^2, exp(sqrt(log D loglog D) / 8)) */
1029 245 : START:
1030 : do
1031 : {
1032 69652 : if (!FIRST) LIMC = bnf_increase_LIMC(LIMC,LIMCMAX);
1033 69652 : if (DEBUGLEVEL>2 && LIMC > LIMC0)
1034 0 : err_printf("%s*** Bach constant: %f\n", FIRST?"":"\n", LIMC/LOGD2);
1035 69652 : FIRST = 0; set_avma(av);
1036 69652 : guncloneNULL(BQ.subFB);
1037 69652 : guncloneNULL(BQ.powsubFB);
1038 69652 : clearhash(BQ.hashtab);
1039 69656 : if (LIMC < cp) LIMC = cp;
1040 69656 : if (LIMC2 < LIMC) LIMC2 = LIMC;
1041 69656 : if (BQ.PRECREG) qfr_data_init(q.D, BQ.PRECREG, &q);
1042 :
1043 69656 : FBquad(&BQ, LIMC2, LIMC, &GRHcheck);
1044 69656 : if (DEBUGLEVEL>2) timer_printf(&T, "factor base");
1045 69656 : BQ.subFB = subFBquad(&BQ, q.D, lim + 0.5, minSFB);
1046 69656 : if (DEBUGLEVEL>2) timer_printf(&T, "subFBquad = %Ps",
1047 : vecpermute(BQ.FB, BQ.subFB));
1048 69656 : nsubFB = lg(BQ.subFB) - 1;
1049 : }
1050 69656 : while (nsubFB < (expi(D) > 15 ? 3 : 2));
1051 : /* invhr = 2^r1 (2pi)^r2 / sqrt(D) w ~ L(chi,1) / hR */
1052 69411 : invhr = gmul(dbltor((BQ.PRECREG?2.:M_PI)/sdrc),
1053 : compute_invresquad(&GRHcheck, LIMC));
1054 69406 : BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
1055 69411 : if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
1056 69411 : BQ.limhash = (LIMC & HIGHMASK)? (HIGHBIT>>1): LIMC*LIMC;
1057 :
1058 69411 : need = BQ.KC + RELSUP - 2;
1059 69411 : current = 0;
1060 69411 : W = NULL;
1061 69411 : sfb_trials = nreldep = nrelsup = 0;
1062 69411 : s = nsubFB + RELSUP;
1063 69411 : av2 = avma;
1064 :
1065 : do
1066 : {
1067 102885 : if ((nreldep & 3) == 1 || (nrelsup & 7) == 1) {
1068 28720 : if (DEBUGLEVEL>2) err_printf("*** Changing sub factor base\n");
1069 28720 : gunclone(BQ.subFB);
1070 28720 : gunclone(BQ.powsubFB);
1071 28720 : BQ.subFB = gclone(vecslice(BQ.vperm, 1, nsubFB));
1072 28720 : BQ.powsubFB = powsubFBquad(&BQ,CBUCH+1);
1073 28720 : if (DEBUGLEVEL>2) timer_printf(&T, "powsubFBquad");
1074 28720 : clearhash(BQ.hashtab);
1075 : }
1076 102885 : need += 2;
1077 102885 : mat = cgetg(need+1, t_MAT);
1078 102885 : extraC = cgetg(need+1, t_VEC);
1079 102885 : if (!W) { /* first time */
1080 69411 : C = extraC;
1081 69411 : triv = trivial_relations(&BQ, mat, C);
1082 69411 : if (DEBUGLEVEL>2) err_printf("KC = %ld, need %ld relations\n", BQ.KC, need);
1083 : } else {
1084 33474 : triv = 0;
1085 33474 : if (DEBUGLEVEL>2) err_printf("...need %ld more relations\n", need);
1086 : }
1087 102885 : if (BQ.PRECREG) {
1088 189378 : for (i = triv+1; i<=need; i++) {
1089 184254 : gel(mat,i) = zero_zv(BQ.KC);
1090 184254 : gel(extraC,i) = cgetr(BQ.PRECREG);
1091 : }
1092 5124 : real_relations(&BQ, need - triv, ¤t, s,LIMC,mat + triv,extraC + triv);
1093 : } else {
1094 1693783 : for (i = triv+1; i<=need; i++) {
1095 1596014 : gel(mat,i) = zero_zv(BQ.KC);
1096 1596022 : gel(extraC,i) = gen_0;
1097 : }
1098 97769 : imag_relations(&BQ, need - triv, ¤t, LIMC,mat + triv);
1099 : }
1100 :
1101 102885 : if (!W)
1102 69411 : W = hnfspec_i(mat,BQ.vperm,&dep,&B,&C,nsubFB);
1103 : else
1104 33474 : W = hnfadd_i(W,BQ.vperm,&dep,&B,&C, mat,extraC);
1105 102885 : gerepileall(av2, 4, &W,&C,&B,&dep);
1106 102885 : need = BQ.KC - (lg(W)-1) - (lg(B)-1);
1107 102885 : if (need)
1108 : {
1109 12878 : if (++nreldep > 15 && cbach < 1) goto START;
1110 12878 : continue;
1111 : }
1112 :
1113 90007 : h = ZM_det_triangular(W);
1114 90006 : if (DEBUGLEVEL>2) err_printf("\n#### Tentative class number: %Ps\n", h);
1115 :
1116 90006 : switch(get_R(&BQ, C, (lg(C)-1) - (lg(B)-1) - (lg(W)-1), mulir(h,invhr), &R))
1117 : {
1118 0 : case fupb_PRECI:
1119 0 : BQ.PRECREG = precdbl(BQ.PRECREG);
1120 0 : FIRST = 1; goto START;
1121 :
1122 20596 : case fupb_RELAT:
1123 20596 : if (++nrelsup > MAXRELSUP)
1124 : {
1125 21 : if (++sfb_trials > SFB_MAX && cbach <= 1) goto START;
1126 21 : if (nsubFB < minss(10,BQ.KC)) nsubFB++;
1127 : }
1128 20596 : need = minss(BQ.KC, nrelsup);
1129 : }
1130 : }
1131 102883 : while (need);
1132 : /* DONE */
1133 69409 : if (!quad_be_honest(&BQ)) goto START;
1134 69409 : if (DEBUGLEVEL>2) timer_printf(&T, "be honest");
1135 69409 : clearhash(BQ.hashtab);
1136 69411 : free_GRHcheck(&GRHcheck);
1137 :
1138 69411 : gen = get_clgp(&BQ,W,&cyc);
1139 69411 : gunclone(BQ.subFB);
1140 69411 : gunclone(BQ.powsubFB);
1141 69411 : return mkvec4(h, cyc, gen, R);
1142 : }
1143 : GEN
1144 4396 : Buchquad(GEN D, double c, double c2, long prec)
1145 : {
1146 4396 : pari_sp av = avma;
1147 4396 : GEN z = Buchquad_i(D, c, c2, prec);
1148 4396 : return gerepilecopy(av, z);
1149 : }
1150 :
1151 : GEN
1152 0 : buchimag(GEN D, GEN c, GEN c2, GEN REL)
1153 0 : { (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),0); }
1154 :
1155 : GEN
1156 0 : buchreal(GEN D, GEN flag, GEN c, GEN c2, GEN REL, long prec) {
1157 0 : if (signe(flag)) pari_err_IMPL("narrow class group");
1158 0 : (void)REL; return Buchquad(D,gtodouble(c),gtodouble(c2),prec);
1159 : }
1160 :
1161 : GEN
1162 4396 : quadclassunit0(GEN x, long flag, GEN data, long prec)
1163 : {
1164 : long lx;
1165 4396 : double c1 = 0.0, c2 = 0.0;
1166 :
1167 4396 : if (!data) lx=1;
1168 : else
1169 : {
1170 28 : lx = lg(data);
1171 28 : if (typ(data)!=t_VEC) pari_err_TYPE("quadclassunit", data);
1172 28 : if (lx > 7) pari_err_DIM("quadclassunit [tech vector]");
1173 28 : if (lx > 3) lx = 3;
1174 : }
1175 4396 : switch(lx)
1176 : {
1177 21 : case 3: c2 = gtodouble(gel(data,2));
1178 28 : case 2: c1 = gtodouble(gel(data,1));
1179 : }
1180 4396 : if (flag) pari_err_IMPL("narrow class group");
1181 4396 : return Buchquad(x,c1,c2,prec);
1182 : }
1183 : GEN
1184 61169 : quadclassno(GEN D)
1185 : {
1186 61169 : pari_sp av = avma;
1187 61169 : GEN h = abgrp_get_no(Buchquad_i(D, 0, 0, 0));
1188 61168 : return icopy_avma(h, av);
1189 : }
1190 : long
1191 3958 : quadclassnos(long D)
1192 : {
1193 3958 : pari_sp av = avma;
1194 3958 : long h = itos(abgrp_get_no(Buchquad_i(stoi(D), 0, 0, 0)));
1195 3958 : return gc_long(av, h);
1196 : }
|
