1 /*
2
3 classgroup.c - computations with class groups and quadratic forms
4
5 Copyright (C) 2009, 2010, 2018 Andreas Enge
6
7 This file is part of CM.
8
9 CM is free software; you can redistribute it and/or modify it under
10 the terms of the GNU General Public License as published by the
11 Free Software Foundation; either version 3 of the license, or (at your
12 option) any later version.
13
14 CM is distributed in the hope that it will be useful, but WITHOUT ANY
15 WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 for more details.
18
19 You should have received a copy of the GNU General Public License along
20 with CM; see the file COPYING. If not, write to the Free Software
21 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
22 */
23
24 #include "cm_class-impl.h"
25
26 static int_cl_t classgroup_gcdext (int_cl_t *u, int_cl_t *v, int_cl_t a,
27 int_cl_t b);
28 static int_cl_t classgroup_tonelli (int_cl_t a, int_cl_t p);
29 static int_cl_t classgroup_fundamental_discriminant_conductor (int_cl_t d,
30 uint_cl_t *cond_primes, unsigned int *cond_exp);
31
32 /*****************************************************************************/
33 /* */
34 /* functions transforming between int_cl_t and mpz_t */
35 /* */
36 /*****************************************************************************/
37
cm_classgroup_mpz_set_icl(mpz_t rop,int_cl_t op)38 void cm_classgroup_mpz_set_icl (mpz_t rop, int_cl_t op)
39 /* relies on the fact that int_cl_t has 64 bits */
40
41 {
42 int_cl_t mask = (((int_cl_t) 1) << 32) - 1;
43 int_cl_t abso = (op > 0 ? op : -op);
44
45 /* copy 32 high bits */
46 mpz_set_ui (rop, (unsigned long int) (abso >> 32));
47 /* add 32 low bits */
48 mpz_mul_2exp (rop, rop, 32);
49 mpz_add_ui (rop, rop, (unsigned long int) (abso & mask));
50
51 if (op < 0)
52 mpz_neg (rop, rop);
53 }
54
55 /*****************************************************************************/
56
cm_classgroup_mpz_get_icl(mpz_t op)57 int_cl_t cm_classgroup_mpz_get_icl (mpz_t op)
58
59 {
60 int_cl_t rop;
61 int sign = (mpz_cmp_ui (op, 0) < 0 ? -1 : 1);
62 mpz_t tmp;
63
64 mpz_init (tmp);
65
66 /* switch to absolute value */
67 if (sign < 0)
68 mpz_neg (op, op);
69
70 /* extract 32 high bits */
71 mpz_tdiv_q_2exp (tmp, op, 32);
72 rop = ((int_cl_t) (mpz_get_ui (tmp))) << 32;
73 /* extract 32 low bits */
74 mpz_tdiv_r_2exp (tmp, op, 32);
75 rop += mpz_get_ui (tmp);
76
77 if (sign < 0) {
78 rop = -rop;
79 mpz_neg (op, op);
80 }
81
82 mpz_clear (tmp);
83
84 return rop;
85 }
86
87
88 /*****************************************************************************/
89 /* */
90 /* creating and freeing variables */
91 /* */
92 /*****************************************************************************/
93
cm_classgroup_init(cm_classgroup_t * cl,int_cl_t disc,bool verbose)94 void cm_classgroup_init (cm_classgroup_t *cl, int_cl_t disc, bool verbose)
95
96 {
97 int h;
98 int length;
99 /* Nobody will ever need a class group with more than 64 components.
100 More precisely, this is the limit imposed by the class number encoded
101 in 64 bits. */
102 int_cl_t ord [64];
103 cm_form_t gen [64];
104 cm_form_t Ppow;
105 int i, j, k;
106
107 if (disc >= 0) {
108 printf ("\n*** The discriminant must be negative.\n");
109 exit (1);
110 }
111 else if (disc % 4 != 0 && (disc - 1) % 4 != 0) {
112 printf ("\n*** %"PRIicl" is not a quadratic discriminant.\n", disc);
113 exit (1);
114 }
115 else
116 cl->d = disc;
117
118 cl->h = cm_classgroup_h (cl->d);
119 cl->form = (cm_form_t *) malloc (cl->h * sizeof (cm_form_t));
120 cl->conj = (int *) malloc (cl->h * sizeof (int));
121 if (verbose)
122 printf ("Class number: h = %d\n", cl->h);
123
124 length = cm_pari_classgroup (disc, ord, gen);
125 if (verbose)
126 for (i = 0; i < length; i++)
127 printf ("%"PRIicl" [%"PRIicl", %"PRIicl"]\n",
128 ord [i], gen [i].a, gen [i].b);
129
130 /* Roll out the class group from its generators. */
131 cl->form [0].a = 1;
132 cl->form [0].b = (disc % 2 == 0 ? 0 : 1);
133 h = 1;
134 for (i = 0; i < length; i++) {
135 Ppow = gen [i];
136 for (j = 1; j < ord [i]; j++) {
137 /* Ppow contains gen [i]^j; multiply the first h forms by this. */
138 for (k = 0; k < h; k++)
139 cm_classgroup_compose (&(cl->form [h*j+k]), cl->form [k], Ppow,
140 disc);
141 cm_classgroup_compose (&Ppow, Ppow, gen [i], disc);
142 }
143 h *= ord [i];
144 }
145
146 /* Pair up inverse forms. */
147 for (i = 0; i < cl->h; i++) {
148 for (j = i;
149 j < cl->h && cl->form [j].b != -cl->form [i].b;
150 j++);
151 if (j == cl->h || cl->form [i].b == 0)
152 cl->conj [i] = i;
153 else {
154 cl->conj [i] = j;
155 cl->conj [j] = i;
156 }
157 }
158 }
159
160 /*****************************************************************************/
161
cm_classgroup_clear(cm_classgroup_t * cl)162 void cm_classgroup_clear (cm_classgroup_t *cl)
163
164 {
165 free (cl->form);
166 free (cl->conj);
167 }
168
169 /*****************************************************************************/
170 /* */
171 /* elementary number theory with int_cl_t */
172 /* */
173 /*****************************************************************************/
174
cm_classgroup_mod(int_cl_t a,uint_cl_t p)175 uint_cl_t cm_classgroup_mod (int_cl_t a, uint_cl_t p)
176 /* returns a representative of a % p in [0, p-1[ */
177
178 {
179 int_cl_t res = a % (int_cl_t) p;
180 return (res < 0 ? res + (int_cl_t) p : res);
181 }
182
183 /*****************************************************************************/
184
cm_classgroup_gcd(int_cl_t a,int_cl_t b)185 int_cl_t cm_classgroup_gcd (int_cl_t a, int_cl_t b)
186 /* returns the positive gcd of a and b */
187
188 {
189 int_cl_t r;
190
191 if (a == 0)
192 return b;
193 else if (b == 0)
194 return a;
195 else
196 {
197 if (a < 0)
198 a = -a;
199 if (b < 0)
200 b = -b;
201 r = a % b;
202 while (r > 0)
203 {
204 a = b;
205 b = r;
206 r = a % b;
207 }
208 return b;
209 }
210 }
211
212 /*****************************************************************************/
213
classgroup_gcdext(int_cl_t * u,int_cl_t * v,int_cl_t a,int_cl_t b)214 static int_cl_t classgroup_gcdext (int_cl_t *u, int_cl_t *v, int_cl_t a,
215 int_cl_t b)
216 /* returns the positive gcd d of a and b; if u and v is not NULL, */
217 /* modifies them such that u a + v b = d; it is also possible to have */
218 /* only one of u and v non NULL. */
219
220 {
221 int_cl_t r0, r1, r2, u0, u1, u2, v0, v1, v2, q;
222 int sgn_a, sgn_b;
223
224 if (a < 0) {
225 sgn_a = -1;
226 r0 = -a;
227 }
228 else {
229 sgn_a = 1;
230 r0 = a;
231 }
232 if (b < 0) {
233 sgn_b = -1;
234 r1 = -b;
235 }
236 else {
237 sgn_b = 1;
238 r1 = b;
239 }
240 /* computing u and v might be faster than permanent tests for NULL */
241 u0 = 1;
242 u1 = 0;
243 v0 = 0;
244 v1 = 1;
245
246 while (r1 != 0) {
247 q = r0 / r1;
248 r2 = r0 % r1;
249 r0 = r1;
250 r1 = r2;
251 u2 = u0 - q * u1;
252 u0 = u1;
253 u1 = u2;
254 v2 = v0 - q * v1;
255 v0 = v1;
256 v1 = v2;
257 }
258
259 if (u != NULL)
260 *u = sgn_a * u0;
261 if (v != NULL)
262 *v = sgn_b * v0;
263
264 return r0;
265 }
266
267 /*****************************************************************************/
268
cm_classgroup_kronecker(int_cl_t a,int_cl_t b)269 int cm_classgroup_kronecker (int_cl_t a, int_cl_t b)
270 /* computes the Kronecker symbol (a / b) following Algorithm 1.4.12 */
271 /* in Cohen93 (by the binary algorithm) */
272 {
273 const int tab [8] = {0, 1, 0, -1, 0, -1, 0, 1};
274 int k;
275 int_cl_t r;
276
277 /* step 1 */
278 if (b == 0) {
279 if (a == 1 || a == -1)
280 return 1;
281 else
282 return 0;
283 }
284
285 /* step 2*/
286 if (a % 2 == 0 && b % 2 == 0)
287 return 0;
288
289 while (b % 4 == 0)
290 b >>= 2;
291 if (b % 2 == 0) {
292 b >>= 1;
293 k = tab [cm_classgroup_mod (a, (uint_cl_t) 8)];
294 }
295 else
296 k = 1;
297
298 if (b < 0) {
299 b = -b;
300 if (a < 0)
301 k = -k;
302 }
303
304 /* step 3 */
305 a = cm_classgroup_mod (a, b);
306
307 /* steps 4 to 6 */
308 while (a != 0) {
309 while (a % 4 == 0)
310 a >>= 2;
311 if (a % 2 == 0) {
312 a >>= 1;
313 k *= tab [cm_classgroup_mod (b, (uint_cl_t) 8)];
314 }
315 if (b > a) {
316 r = b - a;
317 if ( cm_classgroup_mod (a, (uint_cl_t) 4) == 3
318 && cm_classgroup_mod (b, (uint_cl_t) 4) == 3)
319 k = -k;
320 b = a;
321 a = r;
322 }
323 else
324 a -= b;
325 }
326
327 if (b > 1)
328 return 0;
329 else
330 return k;
331 }
332
333 /*****************************************************************************/
334
classgroup_tonelli(int_cl_t a,int_cl_t p)335 static int_cl_t classgroup_tonelli (int_cl_t a, int_cl_t p)
336 {
337 mpz_t az, pz, rz;
338 int_cl_t r;
339
340 mpz_init (az);
341 mpz_init (pz);
342 mpz_init (rz);
343
344 cm_classgroup_mpz_set_icl (az, a);
345 cm_classgroup_mpz_set_icl (pz, p);
346 cm_nt_mpz_tonelli_z (rz, az, pz);
347 r = cm_classgroup_mpz_get_icl (rz);
348
349 mpz_clear (az);
350 mpz_clear (pz);
351 mpz_clear (rz);
352
353 return r;
354 }
355
356 /*****************************************************************************/
357 /* */
358 /* functions computing fundamental discriminants, conductors and class */
359 /* numbers */
360 /* */
361 /*****************************************************************************/
362
cm_classgroup_factor(int_cl_t d,uint_cl_t * factors,unsigned int * exponents)363 void cm_classgroup_factor (int_cl_t d, uint_cl_t *factors,
364 unsigned int *exponents)
365 /* factors the absolute value of d by trial division. The prime factors */
366 /* are stored in "factors", their multiplicities in "exponents", which */
367 /* must provide sufficiently much space. The list of prime factors is */
368 /* terminated by 0, so that 12 entries suffice for a number of 32 bits, */
369 /* and 17 entries for a number of 64 bits. */
370
371 {
372 uint_cl_t no, trial, trial2;
373 int j;
374
375 if (d < 0)
376 no = -d;
377 else
378 no = d;
379
380 j = 0;
381 trial = 0;
382 trial2 = 0;
383 while (trial2 <= no) {
384 if (trial == 0) {
385 trial = 2;
386 trial2 = 4;
387 }
388 else if (trial == 2) {
389 trial = 3;
390 trial2 = 9;
391 }
392 else {
393 trial += 2;
394 trial2 += 4 * (trial - 1);
395 }
396 if (no % trial == 0) {
397 factors [j] = trial;
398 no /= trial;
399 exponents [j] = 1;
400 while (no % trial == 0) {
401 no /= trial;
402 exponents [j]++;
403 }
404 j++;
405 }
406 }
407 if (no != 1) {
408 factors [j] = no;
409 exponents [j] = 1;
410 j++;
411 }
412 factors [j] = 0;
413 }
414
415 /*****************************************************************************/
416
classgroup_fundamental_discriminant_conductor(int_cl_t d,uint_cl_t * cond_primes,unsigned int * cond_exp)417 static int_cl_t classgroup_fundamental_discriminant_conductor (int_cl_t d,
418 uint_cl_t *cond_primes, unsigned int *cond_exp)
419 /* returns the fundamental discriminant corresponding to d, and the */
420 /* prime factorisation of the conductor via cond_primes and cond_exp. */
421
422 {
423 int i, j;
424 int_cl_t local_d;
425 int pow4;
426 uint_cl_t factors [17];
427 unsigned int exponents [17];
428 int_cl_t fundamental_d = -1;
429
430 /* handle 2 in the conductor separately */
431 local_d = d;
432 pow4 = 0;
433 while (local_d % 4 == 0) {
434 local_d /= 4;
435 if ((local_d - 1) % 4 == 0 || local_d % 4 == 0)
436 pow4++;
437 else
438 fundamental_d *= 4;
439 }
440 if (local_d % 2 == 0) {
441 local_d /= 2;
442 fundamental_d *= 2;
443 }
444
445 if (pow4 != 0) {
446 cond_primes [0] = 2;
447 cond_exp [0] = pow4;
448 j = 1;
449 }
450 else
451 j = 0;
452
453 cm_classgroup_factor (local_d, factors, exponents);
454
455 for (i = 0; factors [i] != 0; i++) {
456 if (exponents [i] >= 2) {
457 cond_primes [j] = factors [i];
458 cond_exp [j] = exponents [i] / 2;
459 j++;
460 }
461 if (exponents [i] & 1)
462 fundamental_d *= factors [i];
463 }
464 cond_primes [j] = 0;
465
466 return fundamental_d;
467 }
468
469 /*****************************************************************************/
470
cm_classgroup_fundamental_discriminant(int_cl_t d)471 int_cl_t cm_classgroup_fundamental_discriminant (int_cl_t d)
472
473 {
474 uint_cl_t cond_primes [17];
475 unsigned int cond_exp [17];
476
477 return classgroup_fundamental_discriminant_conductor (d, cond_primes,
478 cond_exp);
479 }
480
481 /*****************************************************************************/
482
cm_classgroup_h(int_cl_t d)483 int cm_classgroup_h (int_cl_t d)
484 /* computes the class number of the imaginary quadratic order with */
485 /* discriminant d using Louboutin's formula [Lou02] for the maximal part */
486 /* and the class number formula for non-maximal orders. */
487
488 {
489 int_cl_t fund;
490 uint_cl_t cond_primes [17];
491 unsigned int cond_exp [17];
492 double pi, a2, a, en, enp1, fn, deltaf, sum1, sum2;
493 /* en stands for e_n = exp (- n^2 a2), enp1 for e_{n+1}, */
494 /* deltaf for exp (-2 a2) and fn for exp (- (2n+3) a2), */
495 /* so that e_{n+2} = e_{n+1} * fn. */
496 /* sum1 contains the sum of chi (n) * e_n / n; */
497 /* sum2 contains the sum of (e_n + e_{n+1}) S_n. */
498 int m, n, S, h;
499 int chi, i;
500
501 fund = classgroup_fundamental_discriminant_conductor (d, cond_primes,
502 cond_exp);
503
504 if (fund == -3 || fund == -4)
505 h = 1;
506 else
507 {
508 pi = 2 * asin (1.0);
509 a2 = pi / (-fund);
510 a = sqrt (a2);
511 m = (int) ceil (sqrt (-fund / (2 * pi) * log (-fund / pi) + 6));
512
513 /* initialisation for n = 1 */
514 chi = 1;
515 S = 1;
516 en = exp (-a2);
517 deltaf = en * en;
518 enp1 = deltaf * deltaf;
519 fn = enp1 * en;
520 sum1 = en;
521 sum2 = en + enp1;
522
523 for (n = 2; n <= m; n++)
524 {
525 chi = cm_classgroup_kronecker (fund, (int_cl_t) n);
526 S += chi;
527 en = enp1;
528 enp1 *= fn;
529 fn *= deltaf;
530 sum1 += chi * en / (double) n;
531 sum2 += (en + enp1) * S;
532 }
533
534 h = (int) ((sum1 / a + sum2 * a) / sqrt (pi) + 0.5);
535 }
536
537 /* correct for the conductor */
538 for (i = 0; cond_primes [i] != 0; i++)
539 {
540 h *= cond_primes [i]
541 - cm_classgroup_kronecker (fund, (int_cl_t) cond_primes [i]);
542 while (cond_exp [i] > 1)
543 {
544 h *= cond_primes [i];
545 cond_exp [i]--;
546 }
547 }
548
549 /* correct for the units */
550 if (cond_primes [0] != 0) {
551 if (fund == -3)
552 h /= 3;
553 else if (fund == -4)
554 h /= 2;
555 }
556
557 return h;
558 }
559
560
561 /*****************************************************************************/
562 /* */
563 /* functions for computing in the class group */
564 /* */
565 /*****************************************************************************/
566
cm_classgroup_compute_c(int_cl_t a,int_cl_t b,int_cl_t d)567 int_cl_t cm_classgroup_compute_c (int_cl_t a, int_cl_t b, int_cl_t d)
568 /* computes c = (b^2 - d) / (4*a), potentially switching to multi- */
569 /* precision to avoid intermediate overflow */
570
571 {
572 int_cl_t c;
573
574 if ((b > 0 ? b : -b)
575 < ((int_cl_t) 1) << (4 * sizeof (int_cl_t) - 2))
576 c = (b * b - d) / a / 4;
577 else {
578 /* should happen only rarely */
579 mpz_t az, dz, cz;
580 mpz_init (az);
581 mpz_init (cz);
582 mpz_init (dz);
583 cm_classgroup_mpz_set_icl (dz, d);
584 cm_classgroup_mpz_set_icl (az, a);
585 cm_classgroup_mpz_set_icl (cz, b);
586 mpz_mul (cz, cz, cz);
587 mpz_sub (cz, cz, dz);
588 mpz_div (cz, cz, az);
589 mpz_div_2exp (cz, cz, 2);
590 c = cm_classgroup_mpz_get_icl (cz);
591 mpz_clear (az);
592 mpz_clear (cz);
593 mpz_clear (dz);
594 }
595
596 return c;
597 }
598
599 /*****************************************************************************/
600
cm_classgroup_reduce(cm_form_t * Q,int_cl_t d)601 void cm_classgroup_reduce (cm_form_t *Q, int_cl_t d)
602 /* reduces the quadratic form Q without checking if it belongs indeed to */
603 /* the discriminant d. */
604
605 {
606 int_cl_t c, a_minus_b, two_a, offset;
607 bool reduced;
608
609 reduced = false;
610 while (!reduced){
611 /* first step: obtain |b| <= a */
612 a_minus_b = Q->a - Q->b;
613 two_a = 2 * Q->a;
614 if (a_minus_b < 0) {
615 a_minus_b++;
616 /* a trick to obtain the correct rounding */
617 offset = a_minus_b / two_a;
618 /* Since a_minus_b is negative, a negative remainder is computed, */
619 /* so the quotient is effectively rounded to the nearest integer */
620 offset--;
621 /* offset is (a-b) / (2a) floored */
622 offset *= two_a;
623 Q->b += offset;
624 }
625 else if (a_minus_b >= two_a) {
626 offset = a_minus_b / two_a;
627 offset *= two_a;
628 Q->b += offset;
629 }
630
631 c = cm_classgroup_compute_c (Q->a, Q->b, d);
632
633 /* if not reduced, invert */
634 if (Q->a < c || (Q->a == c && Q->b >= 0))
635 reduced = true;
636 else {
637 Q->a = c;
638 Q->b = -Q->b;
639 }
640 }
641 }
642
643 /*****************************************************************************/
644
cm_classgroup_compose(cm_form_t * Q,cm_form_t Q1,cm_form_t Q2,int_cl_t d)645 void cm_classgroup_compose (cm_form_t *Q, cm_form_t Q1, cm_form_t Q2,
646 int_cl_t d)
647 /* computes the reduced form Q corresponding to the composition of Q1 and */
648 /* Q2. */
649
650 {
651 int_cl_t s, t, v1, v, w, a2t;
652
653 t = classgroup_gcdext (&v1, NULL, Q2.a, Q1.a);
654
655 if (t == 1) {
656 Q->a = Q1.a * Q2.a;
657 Q->b = Q2.b + Q2.a * v1 * (Q1.b - Q2.b);
658 }
659 else {
660 s = (Q1.b + Q2.b) / 2;
661 t = classgroup_gcdext (&w, &v, s, t);
662 v *= v1;
663 a2t = Q2.a / t;
664 Q->a = (Q1.a / t) * a2t;
665 Q->b = ((s - Q2.b) * v - w * (Q2.b * Q2.b - d) / (4 * Q2.a))
666 % (2 * Q->a); /* intermediate reduction to avoid overflow */
667 Q->b = (Q2.b + 2 * Q->b * a2t) % (2 * Q->a);
668 }
669 cm_classgroup_reduce (Q, d);
670 }
671
672 /*****************************************************************************/
673
cm_classgroup_prime_form(int_cl_t p,int_cl_t d)674 cm_form_t cm_classgroup_prime_form (int_cl_t p, int_cl_t d)
675 /* Assumes that p is a ramified or split prime and returns the reduction */
676 /* of the prime form above p with non-negative b. */
677
678 {
679 cm_form_t Q;
680
681 Q.a = p;
682 if (p == 2)
683 if (d % 8 == 0)
684 Q.b = 0;
685 else if ((d - 4) % 8 == 0)
686 Q.b = 2;
687 else
688 Q.b = 1;
689 else {
690 Q.b = classgroup_tonelli (d, p);
691 /* fix parity of b */
692 if ((d + Q.b) % 2 != 0)
693 Q.b += p;
694 }
695
696 cm_classgroup_reduce (&Q, d);
697
698 return Q;
699 }
700
701 /*****************************************************************************/
702 /*****************************************************************************/
703