1 /*
2
3 class.c - code for computing class polynomials
4
5 Copyright (C) 2009, 2010, 2011, 2012, 2015, 2016, 2017, 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 double class_get_valuation (cm_class_t c);
27 /* return some value related to heights and depending on the function */
28 static int class_get_height (cm_class_t c);
29 /* in the real case, returns the binary length of the largest */
30 /* coefficient of the minimal polynomial */
31 /* in the complex case, returns the binary length of the largest */
32 /* coefficient with respect to the decomposition over an integral basis */
33
34 static bool cm_class_compute_parameter (cm_class_t *c, bool verbose);
35 static void correct_nsystem_entry (cm_form_t *Q, int_cl_t N, int_cl_t b0,
36 int_cl_t d);
37 /* changes Q to be compatible with the N-system condition */
38 static void compute_nsystem (cm_form_t *nsystem, int *conj, cm_class_t *c,
39 cm_classgroup_t cl, bool verbose);
40 /* computes and returns an N-system of forms, together with information
41 on the embeddings in the case of real class polynomials in conj */
42 static fprec_t compute_precision (cm_class_t c, cm_classgroup_t cl,
43 bool verbose);
44
45 static void eval (cm_class_t c, cm_modclass_t mc, ctype rop, cm_form_t Q);
46 static void compute_conjugates (ctype *conjugate, cm_form_t *nsystem,
47 int *conj, cm_class_t c, cm_modclass_t mc, bool verbose);
48 static void write_conjugates (cm_class_t c, ctype *conjugates, int *conj);
49 static bool read_conjugates (cm_class_t c, ctype *conjugates, int *conj);
50
51 static bool get_quadratic (mpz_t out1, mpz_t out2, ctype in, int_cl_t d);
52
53 static void get_root_mod_P (cm_class_t c, mpz_t root, mpz_t P, bool verbose);
54 static mpz_t* cm_get_j_mod_P_from_modular (int *no, const char* modpoldir,
55 char type, int level, mpz_t root, mpz_t P);
56
57 /* the remaining functions are put separately as their code is quite long */
58 static void real_compute_minpoly (cm_class_t c, ctype *conjugate, int *conj,
59 bool print);
60 static void complex_compute_minpoly (cm_class_t c, ctype *conjugate,
61 bool print);
62 static bool doubleeta_compute_parameter (cm_class_t *c);
63 static mpz_t* weber_cm_get_j_mod_P (cm_class_t c, mpz_t root, mpz_t P,
64 int *no, bool verbose);
65 static mpz_t* simpleeta_cm_get_j_mod_P (cm_class_t c, mpz_t root, mpz_t P,
66 int *no);
67
68 /*****************************************************************************/
69 /* */
70 /* constructor and destructor */
71 /* */
72 /*****************************************************************************/
73
cm_class_init(cm_class_t * c,int_cl_t d,char inv,bool verbose)74 void cm_class_init (cm_class_t *c, int_cl_t d, char inv, bool verbose)
75
76 {
77 int i;
78
79 if (d >= 0) {
80 printf ("\n*** The discriminant must be negative.\n");
81 exit (1);
82 }
83 else if (d % 4 != 0 && (d - 1) % 4 != 0) {
84 printf ("\n*** %"PRIicl" is not a quadratic discriminant.\n", d);
85 exit (1);
86 }
87 else {
88 c->invariant = inv;
89 c->d = d;
90 if (!cm_class_compute_parameter (c, verbose))
91 exit (1);
92 if (inv == CM_INVARIANT_WEBER && c->p [1] == 1)
93 /* special case Weber with d=1 (4): compute ring class field for 4d */
94 /* If d=1 (8), this is the same as the Hilbert class field; */
95 /* if d=5 (8), it is a degree 3 relative extension, which will be */
96 /* corrected when computing a CM curve by applying a 2-isogeny/ */
97 c->d = 4 * d;
98 }
99 if (verbose)
100 printf ("\nDiscriminant %"PRIicl", invariant %c, parameter %s\n",
101 c->d, c->invariant, c->paramstr);
102
103 if (inv == CM_INVARIANT_SIMPLEETA)
104 c->field = CM_FIELD_COMPLEX;
105 else if (inv == CM_INVARIANT_MULTIETA) {
106 for (i = 0; c->p [i] != 0; i++);
107 if (i % 2 == 0)
108 c->field = CM_FIELD_REAL;
109 else
110 c->field = CM_FIELD_COMPLEX;
111 }
112 else
113 c->field = CM_FIELD_REAL;
114
115 c->h = cm_classgroup_h (c->d);
116 c->minpoly_deg = c->h;
117 c->minpoly = (mpz_t *) malloc ((c->minpoly_deg + 1) * sizeof (mpz_t));
118 for (i = 0; i <= c->minpoly_deg; i++)
119 mpz_init (c->minpoly [i]);
120 if (c->field == CM_FIELD_COMPLEX) {
121 c->minpoly_complex = (mpz_t *) malloc (c->minpoly_deg * sizeof (mpz_t));
122 for (i = 0; i < c->minpoly_deg; i++)
123 mpz_init (c->minpoly_complex [i]);
124 }
125 }
126
127 /*****************************************************************************/
128
cm_class_clear(cm_class_t * c)129 void cm_class_clear (cm_class_t *c)
130
131 {
132 int i;
133
134 for (i = 0; i <= c->minpoly_deg; i++)
135 mpz_clear (c->minpoly [i]);
136 free (c->minpoly);
137 if (c->field == CM_FIELD_COMPLEX) {
138 for (i = 0; i < c->minpoly_deg; i++)
139 mpz_clear (c->minpoly_complex [i]);
140 free (c->minpoly_complex);
141 }
142
143 ffree_cache ();
144 }
145
146 /*****************************************************************************/
147
cm_class_compute_parameter(cm_class_t * c,bool verbose)148 static bool cm_class_compute_parameter (cm_class_t *c, bool verbose)
149 /* tests whether the discriminant is suited for the chosen invariant and */
150 /* in this case computes and returns the parameter p */
151 /* otherwise, returns 0 */
152
153 {
154 int i;
155 char* pointer;
156
157 switch (c->invariant) {
158 case CM_INVARIANT_J:
159 c->p [0] = 1;
160 c->p [1] = 0;
161 c->s = 1;
162 c->e = 1;
163 break;
164 case CM_INVARIANT_GAMMA2:
165 if (c->d % 3 == 0) {
166 if (verbose) {
167 printf ("\n*** %"PRIicl" is divisible by 3, so that gamma2 ",
168 c->d);
169 printf ("cannot be used.\n");
170 }
171 return false;
172 }
173 c->p [0] = 1;
174 c->p [1] = 0;
175 c->s = 1;
176 c->e = 1;
177 break;
178 case CM_INVARIANT_GAMMA3:
179 if (c->d % 2 == 0) {
180 if (verbose) {
181 printf ("\n*** %"PRIicl" is divisible by 4, so that gamma3 ",
182 c->d);
183 printf ("cannot be used.\n");
184 }
185 return false;
186 }
187 c->p [0] = 1;
188 c->p [1] = 0;
189 c->s = 1;
190 c->e = 1;
191 break;
192 case CM_INVARIANT_WEBER:
193 if (c->d % 32 == 0) {
194 if (verbose) {
195 printf ("\n*** %"PRIicl" is divisible by 32, so that the Weber ",
196 c->d);
197 printf ("functions cannot be used.\n");
198 }
199 return false;
200 }
201
202 /* If disc is even, let m = -disc and p [0] = m % 8; if disc is */
203 /* odd, compute p [0] for 4*disc and let p [1] = 1. */
204 if (c->d % 4 == 0) {
205 c->p [0] = ((-(c->d)) / 4) % 8;
206 c->p [1] = 0;
207 }
208 else {
209 c->p [0] = (-(c->d)) % 8;
210 c->p [1] = 1;
211 c->p [2] = 0;
212 }
213 c->s = 24;
214 c->e = 1;
215 if (c->p [0] == 1 || c->p [0] == 2 || c->p [0] == 6)
216 c->e *= 2;
217 else if (c->p [0] == 4 || c->p [0] == 5)
218 c->e *= 4;
219 if (c->d % 3 == 0)
220 c->e *= 3;
221 break;
222 case CM_INVARIANT_ATKIN:
223 if (cm_classgroup_kronecker (c->d, (int_cl_t) 71) != -1)
224 /* factor 36, T_5 + T_29 + 1 */
225 c->p [0] = 71;
226 else if (cm_classgroup_kronecker (c->d, (int_cl_t) 131) != -1)
227 /* factor 33, T_61 + 1 */
228 c->p [0] = 131;
229 else if (cm_classgroup_kronecker (c->d, (int_cl_t) 59) != -1)
230 /* factor 30, T_5 + T_29 */
231 c->p [0] = 59;
232 else if (cm_classgroup_kronecker (c->d, (int_cl_t) 47) != -1)
233 /* factor 24, -T_17 */
234 c->p [0] = 47;
235 else {
236 if (verbose) {
237 printf ("\n*** 47, 59, 71 and 131 are inert for %"PRIicl, c->d);
238 printf (", so that atkin cannot be used.\n");
239 }
240 return false;
241 }
242 c->p [1] = 0;
243 c->s = 1;
244 c->e = 1;
245 break;
246 case CM_INVARIANT_MULTIETA:
247 c->p [0] = 3;
248 c->p [1] = 5;
249 c->p [2] = 7;
250 c->p [3] = 0;
251 c->s = 1;
252 c->e = 1;
253 break;
254 case CM_INVARIANT_DOUBLEETA:
255 if (!doubleeta_compute_parameter (c))
256 return false;
257 break;
258 case CM_INVARIANT_SIMPLEETA:
259 c->p [0] = 25;
260 if (cm_classgroup_kronecker (c->d, (int_cl_t) (c->p [0])) == -1) {
261 if (verbose)
262 printf ("*** Unsuited discriminant\n\n");
263 return false;
264 }
265
266 c->p [1] = 0;
267 c->s = 1;
268 c->e = 1;
269 break;
270 default: /* should not occur */
271 printf ("class_compute_parameter called for "
272 "unknown class invariant '%c'\n", c->invariant);
273 exit (1);
274 }
275
276 /* create parameter string */
277 pointer = c->paramstr;
278 i = 0;
279 do {
280 pointer += sprintf (pointer, "%i_", c->p [i]);
281 i++;
282 } while (c->p [i] != 0);
283 sprintf (pointer, "%i_%i", c->e, c->s);
284 return true;
285 }
286
287 /*****************************************************************************/
288
doubleeta_compute_parameter(cm_class_t * c)289 static bool doubleeta_compute_parameter (cm_class_t *c)
290 /* Compute p1 <= p2 prime following Cor. 3.1 of [EnSc04], that is, */
291 /* - 24 | (p1-1)(p2-1). */
292 /* - p1, p2 are not inert; */
293 /* - if p1!=p2, then p1, p2 do not divide the conductor; */
294 /* - if p1=p2=p, then either p splits or divides the conductor. */
295 /* The case p1=p2=2 of Cor. 3.1 is excluded by divisibility by 24. */
296 /* Minimise with respect to the height factor gained, which is */
297 /* 12 psi (p1*p2) / (p1-1)(p2-1); */
298 /* then p1, p2 <= the smallest split prime which is 1 (mod 24). */
299 /* Let p = 1000*p2 + p1. */
300
301 {
302 int_cl_t cond2 = c->d / cm_classgroup_fundamental_discriminant (c->d);
303 /* square of conductor */
304 const unsigned long int maxprime = 997;
305 unsigned long int primelist [169];
306 /* list of suitable primes, terminated by 0; big enough to hold all */
307 /* primes <= maxprime */
308 unsigned long int p1, p2, p1opt = 0, p2opt = 0;
309 double quality, opt;
310 bool ok;
311 int i, j;
312
313 /* determine all non-inert primes */
314 i = 0;
315 p1 = 2;
316 ok = false;
317 do {
318 int kro = cm_classgroup_kronecker (c->d, (int_cl_t) p1);
319 if (kro != -1) {
320 primelist [i] = p1;
321 i++;
322 }
323 if (kro == 1 && (p1 - 1) % 24 == 0)
324 ok = true;
325 else
326 p1 = cm_nt_next_prime (p1);
327 }
328 while (p1 <= maxprime && !ok);
329 primelist [i] = 0;
330
331 /* search for the best tuple */
332 opt = 0.0;
333 for (j = 0, p2 = primelist [j]; p2 != 0; j++, p2 = primelist [j])
334 for (i = 0, p1 = primelist [i]; i <= j; i++, p1 = primelist [i])
335 if ( ((p1 - 1)*(p2 - 1)) % 24 == 0
336 && ( (p1 != p2 && cond2 % p1 != 0 && cond2 % p2 != 0)
337 || (p1 == p2 && ((-c->d) % p1 != 0 || cond2 % p1 == 0)))) {
338 quality = (p1 == p2 ? p1 : p1 + 1) * (p2 + 1) / (double) (p1 - 1)
339 / (double) (p2 - 1);
340 if (quality > opt) {
341 p1opt = p1;
342 p2opt = p2;
343 opt = quality;
344 ok = true;
345 }
346 }
347 c->p [0] = p1opt;
348 c->p [1] = p2opt;
349 c->p [2] = 0;
350 c->s = 1;
351 c->e = 1;
352 return ok;
353 }
354
355 /*****************************************************************************/
356 /* */
357 /* functions handling files */
358 /* */
359 /*****************************************************************************/
360
cm_class_write(cm_class_t c)361 void cm_class_write (cm_class_t c)
362 /* writes the class polynomial to the file */
363 /* CM_CLASS_DATADIR + "/cp_" + d + "_" + invariant + "_" + paramstr */
364 /* + ".dat" */
365
366 {
367 char filename [400];
368 FILE *f;
369 int i;
370
371 sprintf (filename, "%s/cp_%"PRIicl"_%c_%s.dat", CM_CLASS_DATADIR, -c.d,
372 c.invariant, c.paramstr);
373
374 if (!cm_file_open_write (&f, filename))
375 exit (1);
376
377 fprintf (f, "%"PRIicl"\n", -c.d);
378 fprintf (f, "%c\n", c.invariant);
379 fprintf (f, "%s\n", c.paramstr);
380 fprintf (f, "%i\n", c.minpoly_deg);
381 for (i = c.minpoly_deg; i >= 0; i--) {
382 mpz_out_str (f, 10, c.minpoly [i]);
383 if (c.field == CM_FIELD_COMPLEX) {
384 fprintf (f, " ");
385 mpz_out_str (f, 10, c.minpoly_complex [i]);
386 }
387 fprintf (f, "\n");
388 }
389
390 cm_file_close (f);
391 }
392
393 /*****************************************************************************/
394
cm_class_read(cm_class_t c)395 bool cm_class_read (cm_class_t c)
396 /* reads the class polynomial from a file written by cm_class_write */
397 /* If an error occurs, the return value is false. */
398
399 {
400 char filename [400];
401 FILE* f;
402 int i;
403 char inv;
404 char pars [255];
405 int_cl_t disc;
406
407 sprintf (filename, "%s/cp_%"PRIicl"_%c_%s.dat", CM_CLASS_DATADIR, -c.d,
408 c.invariant, c.paramstr);
409
410 if (!cm_file_open_read (&f, filename))
411 return false;
412
413 if (!fscanf (f, "%"SCNicl"\n", &disc))
414 return false;
415 if (-disc != c.d) {
416 printf ("*** Inconsistency between file '%s' ", filename);
417 printf ("and internal data:\n");
418 printf ("*** discriminant %"PRIicl" instead of %"PRIicl"\n", -disc, c.d);
419 exit (1);
420 }
421 if (!fscanf (f, "%c", &inv))
422 return false;
423 if (inv != c.invariant) {
424 printf ("*** Inconsistency between file '%s' ", filename);
425 printf ("and internal data:\n");
426 printf ("*** invariant '%c' instead of '%c'\n", inv, c.invariant);
427 exit (1);
428 }
429 if (!fscanf (f, "%254s", pars))
430 return false;
431 if (strcmp (pars, c.paramstr)) {
432 printf ("*** Inconsistency between file '%s' ", filename);
433 printf ("and internal data:\n");
434 printf ("*** parameter %s instead of %s\n", pars, c.paramstr);
435 exit (1);
436 }
437 if (!fscanf (f, "%i", &i))
438 return false;
439 if (i != c.minpoly_deg) {
440 printf ("*** Inconsistency between file '%s' ", filename);
441 printf ("and internal data:\n");
442 printf ("*** degree %i instead of %i\n", i, c.minpoly_deg);
443 exit (1);
444 }
445
446 for (i = c.minpoly_deg; i >= 0; i--) {
447 mpz_inp_str (c.minpoly [i], f, 10);
448 if (c.field == CM_FIELD_COMPLEX)
449 mpz_inp_str (c.minpoly_complex [i], f, 10);
450 }
451
452 cm_file_close (f);
453
454 return true;
455 }
456
457 /*****************************************************************************/
458
write_conjugates(cm_class_t c,ctype * conjugate,int * conj)459 static void write_conjugates (cm_class_t c, ctype *conjugate, int *conj)
460 /* writes the conjugates to the file
461 CM_CLASS_TMPDIR + "/tmp_" + d + "_" + invariant + "_" + paramstr + "_"
462 + prec + "_conjugates.dat" */
463
464 {
465 char filename [400];
466 FILE *f;
467 int i;
468
469 sprintf (filename, "%s/tmp_%"PRIicl"_%c_%s_%i_conjugates.dat",
470 CM_CLASS_TMPDIR, -c.d, c.invariant, c.paramstr,
471 (int) cget_prec (conjugate [0]));
472
473 if (!cm_file_open_write (&f, filename))
474 exit (1);
475
476 for (i = 0; i < c.h; i++)
477 if (conj [i] >= i) {
478 cout_str (f, 16, 0, conjugate [i]);
479 fprintf (f, "\n");
480 }
481
482 cm_file_close (f);
483 }
484
485 /*****************************************************************************/
486
read_conjugates(cm_class_t c,ctype * conjugate,int * conj)487 static bool read_conjugates (cm_class_t c, ctype *conjugate, int *conj)
488 /* reads the conjugates from a file written by write_conjugates
489 If the file could not be openend, the return value is false. */
490 {
491 char filename [400];
492 FILE *f;
493 int i;
494
495 sprintf (filename, "%s/tmp_%"PRIicl"_%c_%s_%i_conjugates.dat",
496 CM_CLASS_TMPDIR, -c.d, c.invariant, c.paramstr,
497 (int) cget_prec (conjugate [0]));
498
499 if (!cm_file_open_read (&f, filename))
500 return false;
501
502 for (i = 0; i < c.h; i++)
503 if (conj [i] >= i)
504 cinp_str (conjugate [i], f, NULL, 16);
505
506 cm_file_close (f);
507
508 return true;
509 }
510
511 /*****************************************************************************/
512 /* */
513 /* valuation at infinity and height */
514 /* */
515 /*****************************************************************************/
516
class_get_valuation(cm_class_t c)517 static double class_get_valuation (cm_class_t c)
518 /* returns the (negative) order of the modular function at infinity */
519
520 {
521 double result;
522
523 switch (c.invariant) {
524 case CM_INVARIANT_J:
525 result = 1;
526 break;
527 case CM_INVARIANT_GAMMA2:
528 result = 1.0 / 3;
529 break;
530 case CM_INVARIANT_GAMMA3:
531 result = 0.5;
532 break;
533 case CM_INVARIANT_ATKIN:
534 if (c.p [0] == 47)
535 result = 1.0 / 24;
536 else if (c.p [0] == 59)
537 result = 1.0 / 30;
538 else if (c.p [0] == 71)
539 result = 1.0 / 36;
540 else /* 131 */
541 result = 1.0 / 33;
542 break;
543 case CM_INVARIANT_WEBER:
544 result = 1.0 / 72;
545 break;
546 case CM_INVARIANT_DOUBLEETA:
547 case CM_INVARIANT_MULTIETA:
548 {
549 int num = 1, den = 1, i;
550 for (i = 0; c.p [i] != 0; i++) {
551 num *= c.p [i] - 1;
552 den *= c.p [i] + 1;
553 }
554 if (i == 2)
555 result = num / (double) (12 * den);
556 else if (i == 3)
557 result = num / (double) (6 * den);
558 else /* i == 4 */
559 result = num / (double) (3 * den);
560 }
561 break;
562 case CM_INVARIANT_SIMPLEETA:
563 result = (c.p [0] - 1) / (double) (24 * (c.p [0] + 1));
564 break;
565 default: /* should not occur */
566 printf ("class_get_valuation called for unknown class ");
567 printf ("invariant\n");
568 exit (1);
569 }
570
571 result *= c.e;
572
573 return result;
574 }
575
576 /*****************************************************************************/
577
class_get_height(cm_class_t c)578 static int class_get_height (cm_class_t c)
579 /* in the real case, returns the binary length of the largest */
580 /* coefficient of the minimal polynomial */
581 /* in the complex case, returns the binary length of the largest */
582 /* coefficient with respect to the decomposition over an integral basis */
583
584 {
585 int i, height, cand;
586
587 height = -1;
588 for (i = 0; i < c.minpoly_deg; i++) {
589 cand = mpz_sizeinbase (c.minpoly [i], 2);
590 if (cand > height)
591 height = cand;
592 }
593
594 return height;
595 }
596
597 /*****************************************************************************/
598 /* */
599 /* computing the class polynomial */
600 /* */
601 /*****************************************************************************/
602
correct_nsystem_entry(cm_form_t * Q,int_cl_t N,int_cl_t b0,int_cl_t d)603 static void correct_nsystem_entry (cm_form_t *Q, int_cl_t N, int_cl_t b0,
604 int_cl_t d)
605 /* Changes the form Q by a unimodular transformation so that Q.a is
606 coprime to N and Q.b is congruent to b0 modulo 2*N. */
607
608 {
609 int_cl_t c, tmp;
610
611 /* First achieve gcd (Q->a, N) = 1, which is likely to hold already. */
612 c = (Q->b * Q->b - d) / (4 * Q->a) ;
613 if (cm_classgroup_gcd (Q->a, N) != 1) {
614 /* Translation by k yields C' = A k^2 + B k + C; we wish to reach */
615 /* gcd (C', N) = 1, so for each prime p dividing N, this excludes */
616 /* at most two values modulo p. For p = 2, A and B odd and C even, */
617 /* there is no solution; in this case, we first apply S to exchange */
618 /* A and C. */
619 if (N % 2 == 0 && Q->a % 2 != 0 && Q->b % 2 != 0 && c % 2 == 0) {
620 tmp = Q->a;
621 Q->a = c;
622 c = tmp;
623 Q->b = -Q->b;
624 }
625 while (cm_classgroup_gcd (c, N) != 1) {
626 /* Translate by 1 */
627 c += Q->a + Q->b;
628 Q->b += 2 * Q->a;
629 }
630 /* Apply S */
631 tmp = Q->a;
632 Q->a = c;
633 c = tmp;
634 Q->b = -Q->b;
635 }
636 /* Translate so that Q->b = b0 mod (2 N). */
637 while ((Q->b - b0) % (2*N) != 0) {
638 c += Q->a + Q->b;
639 Q->b += 2 * Q->a;
640 }
641 }
642
643 /*****************************************************************************/
644
compute_nsystem(cm_form_t * nsystem,int * conj,cm_class_t * c,cm_classgroup_t cl,bool verbose)645 static void compute_nsystem (cm_form_t *nsystem, int *conj, cm_class_t *c,
646 cm_classgroup_t cl, bool verbose)
647 /* Compute an N-system, or to be more precise, some part of an N-system
648 that yields all different conjugates up to complex conjugation;
649 this information is passed in conj as follows:
650 If the class polynomial is real, then conj [i] == j if form j in
651 If the class polynomial is complex, then conj [i] = i. */
652
653 {
654 int_cl_t b0, N;
655 cm_form_t neutral, inverse;
656 int h1, h2;
657 int i, j;
658
659 /* Compute the targeted b0 for the N-system and (in the real case) the
660 neutral form. */
661 if (c->invariant == CM_INVARIANT_SIMPLEETA) {
662 bool ok = false;
663
664 if (c->p[0] != 4) {
665 b0 = c->d % 2;
666 while (!ok) {
667 b0 += 2;
668 if ((b0*b0 - c->d) % (4*c->p[0]) != 0)
669 ok = false;
670 else if (c->p[0] != 2 && (c->s/c->e) % 2 == 0 && (b0 - 1) % 4 != 0)
671 ok = false;
672 else if (c->p[0] != 2 && (c->s/c->e) % 3 == 0 && b0 % 3 != 0)
673 ok = false;
674 else
675 ok = true;
676 }
677 }
678 else {
679 b0 = (int_cl_t) -7;
680 }
681 if (verbose)
682 printf ("N %i\ns %i\ne %i\nb0 %"PRIicl"\n", c->p[0], c->s, c->e, b0);
683 N = c->p[0] * c->s / c->e;
684 }
685
686 else {
687 neutral.a = 1;
688 if (c->d % 2 == 0)
689 neutral.b = 0;
690 else
691 neutral.b = 1;
692
693 switch (c->invariant) {
694 case CM_INVARIANT_J:
695 b0 = c->d % 2;
696 /* An even N makes c even if 2 is split so that during the
697 evaluation of eta(z/2) for f1 all forms can be taken from
698 the precomputed ones. */
699 N = 2;
700 break;
701 case CM_INVARIANT_GAMMA2:
702 b0 = 3 * (c->d % 2);
703 /* Use an even N as for j. */
704 N = 6;
705 break;
706 case CM_INVARIANT_GAMMA3:
707 b0 = 1;
708 N = 2;
709 break;
710 case CM_INVARIANT_ATKIN:
711 N = c->p [0];
712 if (c->d % 2 == 0)
713 b0 = 0;
714 else
715 b0 = 1;
716 while ((b0*b0 - c->d) % N != 0)
717 b0 += 2;
718 neutral.a = N;
719 neutral.b = -b0;
720 cm_classgroup_reduce (&neutral, c->d);
721 break;
722 case CM_INVARIANT_WEBER:
723 neutral.a = 1;
724 neutral.b = 0;
725 b0 = 0;
726 N = 48;
727 break;
728 case CM_INVARIANT_DOUBLEETA:
729 case CM_INVARIANT_MULTIETA:
730 {
731 int_cl_t C;
732 N = 1;
733 for (i = 0; c->p [i] != 0; i++)
734 N *= c->p [i];
735 if (c->d % 2 == 0)
736 b0 = 2;
737 else
738 b0 = 1;
739 while (true) {
740 C = (b0*b0 - c->d) / 4;
741 if (C % N == 0 && cm_nt_gcd (C / N, N) == 1)
742 break;
743 b0 += 2;
744 }
745 neutral.a = N;
746 neutral.b = -b0;
747 cm_classgroup_reduce (&neutral, c->d);
748 /* The neutral form corresponds to the product of the primes, */
749 /* but the n-system needs to take s/e into account */
750 N *= c->s / c->e;
751 break;
752 }
753 default: /* should not occur */
754 printf ("compute_nsystem called for unknown class invariant\n");
755 exit (1);
756 }
757 }
758
759 for (i = 0; i < cl.h; i++) {
760 nsystem [i] = cl.form [i];
761 conj [i] = -1;
762 }
763
764 for (i = 0; i < cl.h; i++)
765 /* Pair forms yielding complex conjugate roots. */
766 if (c->field == CM_FIELD_REAL) {
767 if (conj [i] == -1) {
768 /* form did not yet occur in a pair; look for its inverse
769 with respect to neutral_class */
770 nsystem [i].b = -nsystem [i].b;
771 cm_classgroup_compose (&inverse, neutral, nsystem [i], c->d);
772 nsystem [i].b = -nsystem [i].b;
773 j = 0;
774 /* So far, nsystem still contains the reduced forms, so we may look
775 for the inverse form. */
776 while (nsystem [j].a != inverse.a || nsystem [j].b != inverse.b)
777 j++;
778 conj [i] = j;
779 conj [j] = i;
780 }
781 }
782 else /* c->field == CM_FIELD_COMPLEX */
783 conj [i] = i;
784
785 /* Now modify the entries of nsystem. */
786 for (i = 0; i < cl.h; i++)
787 if (conj [i] >= i)
788 correct_nsystem_entry (&(nsystem [i]), N, b0, c->d);
789
790 /* Compute h1 and h2, only for printing. */
791 if (verbose) {
792 h1 = 0;
793 h2 = 0;
794 for (i = 0; i < cl.h; i++)
795 if (conj [i] == i)
796 h1++;
797 else if (conj [i] > i)
798 h2++;
799 printf ("h = %i, h1 = %i, h2 = %i\n", c->h, h1, h2);
800 }
801 }
802
803 /*****************************************************************************/
804
compute_precision(cm_class_t c,cm_classgroup_t cl,bool verbose)805 static fprec_t compute_precision (cm_class_t c, cm_classgroup_t cl,
806 bool verbose) {
807 /* returns an approximation of the required precision */
808
809 const double C = 2114.567;
810 const double pisqrtd = 3.14159265358979323846 * sqrt ((double) (-cl.d));
811 const double cf = class_get_valuation (c);
812 double x, binom = 1.0, simpleprec = 0, prec = 0, M;
813 int_cl_t amax;
814 int i, m;
815 fprec_t precision;
816
817 /* heuristic formula: log (height) = pi * sqrt (|d|) * \sum 1/A */
818 for (i = 0; i < cl.h; i++)
819 simpleprec += 1.0 / cl.form [i].a;
820 simpleprec = ceil (simpleprec * pisqrtd / log (2.0) * cf);
821
822 /* formula of Lemma 8 of [Sutherland11] */
823 amax = 0;
824 for (i = 0; i < cl.h; i++) {
825 x = pisqrtd / cl.form [i].a;
826 if (x < 42)
827 M = log (exp (x) + C);
828 else /* prevent overflow in exponential without changing the result */
829 M = x;
830 prec += M;
831 if (cl.form [i].a > amax)
832 amax = cl.form [i].a;
833 }
834 M = exp (pisqrtd / amax) + C;
835 m = (int) ((cl.h + 1) / (M + 1));
836 for (i = 1; i <= m; i++)
837 binom *= (double) (cl.h - 1 + i) / i / M;
838 prec = ceil ((prec + log (binom)) / log (2.0) * cf);
839
840 if (verbose) {
841 printf ("Heuristic precision bound: %ld\n", (long int) simpleprec);
842 printf ("Less heuristic precision bound: %ld\n", (long int) prec);
843 }
844
845 if (c.invariant == CM_INVARIANT_GAMMA3) {
846 /* increase height estimate by bit size of sqrt (|D|)^h in constant */
847 /* coefficient */
848 prec += (int) (log ((double) (-cl.d)) / log (2.0) / 2.0 * cl.h);
849 if (verbose)
850 printf ("Corrected bound for gamma3: %ld\n", (long int) prec);
851 }
852
853 /* add a security margin */
854 precision = (fprec_t) (prec + 256);
855
856 if (verbose)
857 printf ("Precision: %ld\n", (long int) precision);
858
859 return precision;
860 }
861 /*****************************************************************************/
862
eval(cm_class_t c,cm_modclass_t mc,ctype rop,cm_form_t Q)863 static void eval (cm_class_t c, cm_modclass_t mc, ctype rop, cm_form_t Q)
864
865 {
866 switch (c.invariant) {
867 case CM_INVARIANT_J:
868 cm_modclass_j_eval_quad (mc, rop, Q.a, Q.b);
869 break;
870 case CM_INVARIANT_GAMMA2:
871 cm_modclass_gamma2_eval_quad (mc, rop, Q.a, Q.b);
872 break;
873 case CM_INVARIANT_GAMMA3:
874 cm_modclass_gamma3_eval_quad (mc, rop, Q.a, Q.b);
875 break;
876 case CM_INVARIANT_ATKIN:
877 cm_modclass_atkinhecke_level_eval_quad (mc, rop, Q.a, Q.b, c.p [0]);
878 break;
879 case CM_INVARIANT_SIMPLEETA:
880 case CM_INVARIANT_DOUBLEETA:
881 case CM_INVARIANT_MULTIETA:
882 cm_modclass_multieta_eval_quad (mc, rop, Q.a, Q.b, c.p, c.e);
883 break;
884 case CM_INVARIANT_WEBER:
885 if (c.p [0] == 1) {
886 cm_modclass_f_eval_quad (mc, rop, Q.a, Q.b, 2);
887 cmul_fr (rop, rop, mc.sqrt2_over2);
888 }
889 else if (c.p [0] == 3)
890 cm_modclass_f_eval_quad (mc, rop, Q.a, Q.b, 1);
891 else if (c.p [0] == 5) {
892 cm_modclass_f_eval_quad (mc, rop, Q.a, Q.b, 2);
893 csqr (rop, rop);
894 cdiv_ui (rop, rop, 2ul);
895 }
896 else if (c.p [0] == 7) {
897 cm_modclass_f_eval_quad (mc, rop, Q.a, Q.b, 1);
898 cmul_fr (rop, rop, mc.sqrt2_over2);
899 }
900 else if (c.p [0] == 2 || c.p [0] == 6) {
901 cm_modclass_f1_eval_quad (mc, rop, Q.a, Q.b, 1);
902 csqr (rop, rop);
903 cmul_fr (rop, rop, mc.sqrt2_over2);
904 }
905 else {
906 /* c.p [0] == 4 */
907 cm_modclass_f1_eval_quad (mc, rop, Q.a, Q.b, 1);
908 cpow_ui (rop, rop, 4ul);
909 cmul_fr (rop, rop, mc.sqrt2_over4);
910 }
911
912 if (c.d % 3 == 0)
913 cpow_ui (rop, rop, 3ul);
914
915 if (c.p [0] != 3 && c.p [0] != 5)
916 if (cm_classgroup_kronecker ((int_cl_t) 2, Q.a) == -1)
917 cneg (rop, rop);
918
919 break;
920 default: /* should not occur */
921 printf ("class_eval called for unknown class invariant\n");
922 exit (1);
923 }
924 }
925
926 /*****************************************************************************/
927
compute_conjugates(ctype * conjugate,cm_form_t * nsystem,int * conj,cm_class_t c,cm_modclass_t mc,bool verbose)928 static void compute_conjugates (ctype *conjugate, cm_form_t *nsystem,
929 int *conj, cm_class_t c, cm_modclass_t mc, bool verbose)
930 /* computes the conjugates of the singular value over Q */
931
932 {
933 int i;
934
935 for (i = 0; i < c.h; i++) {
936 if (conj [i] >= i)
937 eval (c, mc, conjugate [i], nsystem [i]);
938 if (verbose && i % 200 == 0) {
939 printf (".");
940 fflush (stdout);
941 }
942 }
943 if (verbose)
944 printf ("\n");
945 }
946
947 /*****************************************************************************/
948
cm_class_compute_minpoly(cm_class_t c,bool checkpoints,bool disk,bool print,bool verbose)949 void cm_class_compute_minpoly (cm_class_t c, bool checkpoints, bool disk,
950 bool print, bool verbose)
951 /* checkpoints indicates whether intermediate results are to be kept in */
952 /* and potentially read from files. */
953 /* disk indicates whether the result should be written to disk. */
954 /* print indicates whether the result should be printed on screen. */
955 {
956 cm_classgroup_t cl, cl2;
957 cm_form_t *nsystem;
958 int *conj;
959 fprec_t prec;
960 cm_modclass_t mc;
961 ctype *conjugate;
962 int i;
963 cm_timer clock_global, clock_local;
964
965 cm_timer_start (clock_global);
966
967 cm_timer_start (clock_local);
968 cm_classgroup_init (&cl, c.d, verbose);
969 cm_timer_stop (clock_local);
970 if (verbose)
971 printf ("--- Time for class group: %.1f\n", cm_timer_get (clock_local));
972
973 if ( (c.invariant == CM_INVARIANT_WEBER
974 && ( (c.p [0] / 10) % 10 == 1
975 || c.p [0] % 10 == 3 || c.p [0] % 10 == 4 || c.p [0] % 10 == 7))
976 || ( (c.invariant == CM_INVARIANT_J || c.invariant == CM_INVARIANT_GAMMA2)
977 && c.d % 4 == 0
978 && ((c.d / 4) % 4 == 0 || ((c.d / 4) - 1) % 4 == 0))) {
979 /* also compute class group for order of conductor less by a factor */
980 /* of 2 */
981 cm_timer_start (clock_local);
982 cm_classgroup_init (&cl2, c.d / 4, verbose);
983 cm_timer_stop (clock_local);
984 if (verbose)
985 printf ("--- Time for class group2: %.1f\n", cm_timer_get (clock_local));
986 }
987 else
988 cl2.d = 0;
989
990 nsystem = (cm_form_t *) malloc (c.h * sizeof (cm_form_t));
991 conj = (int *) malloc (c.h * sizeof (int));
992 cm_timer_start (clock_local);
993 compute_nsystem (nsystem, conj, &c, cl, verbose);
994 cm_timer_stop (clock_local);
995 if (verbose)
996 printf ("--- Time for N-system: %.1f\n", cm_timer_get (clock_local));
997 if (c.invariant == CM_INVARIANT_WEBER && c.p [0] % 100 == 17)
998 prec = compute_precision (c, cl2, verbose);
999 else
1000 prec = compute_precision (c, cl, verbose);
1001
1002 conjugate = (ctype *) malloc (c.h * sizeof (ctype));
1003 for (i = 0; i < c.h; i++)
1004 if (conj [i] >= i)
1005 cinit (conjugate [i], prec);
1006 cm_timer_start (clock_local);
1007 if (!checkpoints || !read_conjugates (c, conjugate, conj)) {
1008 cm_modclass_init (&mc, cl, cl2, prec, checkpoints, verbose);
1009 compute_conjugates (conjugate, nsystem, conj, c, mc, verbose);
1010 if (checkpoints)
1011 write_conjugates (c, conjugate, conj);
1012 cm_modclass_clear (&mc);
1013 }
1014 cm_timer_stop (clock_local);
1015 if (verbose)
1016 printf ("--- Time for conjugates: %.1f\n", cm_timer_get (clock_local));
1017
1018 cm_timer_start (clock_local);
1019 if (c.field == CM_FIELD_REAL)
1020 real_compute_minpoly (c, conjugate, conj, print);
1021 else
1022 complex_compute_minpoly (c, conjugate, print);
1023 cm_timer_stop (clock_local);
1024 if (verbose)
1025 printf ("--- Time for polynomial reconstruction: %.1f\n",
1026 cm_timer_get (clock_local));
1027
1028 for (i = 0; i < c.h; i++)
1029 if (conj [i] >= i)
1030 cclear (conjugate [i]);
1031 free (conjugate);
1032 free (nsystem);
1033 free (conj);
1034 cm_classgroup_clear (&cl);
1035
1036 cm_timer_stop (clock_global);
1037 if (verbose) {
1038 printf ("--- Total time for minimal polynomial: %.1f\n",
1039 cm_timer_get (clock_global));
1040 printf ("Height of minimal polynomial: %d\n", class_get_height (c));
1041 }
1042 if (disk)
1043 cm_class_write (c);
1044 }
1045
1046 /*****************************************************************************/
1047
real_compute_minpoly(cm_class_t c,ctype * conjugate,int * conj,bool print)1048 static void real_compute_minpoly (cm_class_t c, ctype *conjugate, int *conj,
1049 bool print)
1050 /* computes the minimal polynomial of the function over Q */
1051
1052 {
1053 mpfrx_t mpol;
1054 fprec_t prec;
1055 int i;
1056
1057 for (i = 0; conj [i] < i; i++);
1058 prec = cget_prec (conjugate [i]);
1059 mpfrx_init (mpol, c.minpoly_deg + 1, prec);
1060 mpfrcx_reconstruct_from_roots (mpol, conjugate, conj, c.minpoly_deg);
1061
1062 /* the minimal polynomial is now in mpol, rounding to integral polynomial */
1063 for (i = 0; i < c.minpoly_deg; i++)
1064 if (!cm_nt_fget_z (c.minpoly [i], mpol->coeff [i])) {
1065 printf ("*** Accuracy not sufficient for coefficient of X^%d = ", i);
1066 fout_str (stdout, 10, 0ul, mpol->coeff [i]);
1067 printf ("\n");
1068 exit (1);
1069 }
1070 mpz_set_ui (c.minpoly [c.minpoly_deg], 1ul);
1071
1072 if (print) {
1073 printf ("x^%i", c.minpoly_deg);
1074 for (i = c.minpoly_deg - 1; i >= 0; i--) {
1075 int cmp = mpz_cmp_ui (c.minpoly [i], 0);
1076 if (cmp != 0) {
1077 if (cmp > 0) {
1078 printf (" +");
1079 mpz_out_str (stdout, 10, c.minpoly [i]);
1080 }
1081 else
1082 mpz_out_str (stdout, 10, c.minpoly [i]);
1083 if (i >= 2)
1084 printf ("*x^%i", i);
1085 else if (i == 1)
1086 printf ("*x");
1087 }
1088 }
1089 printf ("\n");
1090 }
1091
1092 mpfrx_clear (mpol);
1093 }
1094
1095 /*****************************************************************************/
1096
get_quadratic(mpz_t out1,mpz_t out2,ctype in,int_cl_t d)1097 static bool get_quadratic (mpz_t out1, mpz_t out2, ctype in, int_cl_t d)
1098 /* tries to round the complex number to an integer in the quadratic order */
1099 /* of discriminant d by decomposing it over the integral basis [1, omega] */
1100 /* with omega = sqrt (d/4) or omega = (1 + sqrt (d)) / 2 */
1101 /* The return value reflects the success of the operation. */
1102 /* out1 and out2 are changed. */
1103
1104 {
1105 ftype omega_i, tmp;
1106 bool div4, ok;
1107
1108 finit (omega_i, cget_prec (in));
1109 finit (tmp, cget_prec (in));
1110
1111 div4 = (cm_classgroup_mod (d, (uint_cl_t) 4) == 0);
1112 fsqrt_ui (omega_i, -d);
1113 fdiv_2ui (omega_i, omega_i, 1ul);
1114
1115 fdiv (tmp, in->im, omega_i);
1116 ok = cm_nt_fget_z (out2, tmp);
1117
1118 if (ok) {
1119 if (div4)
1120 fset (tmp, in->re);
1121 else {
1122 fset_z (tmp, out2);
1123 fdiv_2ui (tmp, tmp, 1ul);
1124 fsub (tmp, in->re, tmp);
1125 }
1126 ok = cm_nt_fget_z (out1, tmp);
1127 }
1128
1129 fclear (omega_i);
1130 fclear (tmp);
1131
1132 return ok;
1133 }
1134
1135 /*****************************************************************************/
1136
complex_compute_minpoly(cm_class_t c,ctype * conjugate,bool print)1137 static void complex_compute_minpoly (cm_class_t c, ctype *conjugate,
1138 bool print)
1139 /* computes the minimal polynomial of the function over Q (sqrt D) */
1140
1141 {
1142 int_cl_t fund;
1143 int i;
1144 fprec_t prec;
1145 mpcx_t mpol;
1146
1147 prec = cget_prec (conjugate [0]);
1148 mpcx_init (mpol, c.minpoly_deg + 1, prec);
1149 mpcx_reconstruct_from_roots (mpol, conjugate, c.minpoly_deg);
1150
1151 /* the minimal polynomial is now in mpol; round to integral polynomial */
1152 fund = cm_classgroup_fundamental_discriminant (c.d);
1153 for (i = 0; i < c.h; i++) {
1154 if (!get_quadratic (c.minpoly [i], c.minpoly_complex [i],
1155 mpol->coeff[i], fund)) {
1156 printf ("*** accuracy not sufficient for coefficient of X^%d = ", i);
1157 cout_str (stdout, 10, 0ul, mpol->coeff [i]);
1158 printf ("\n");
1159 exit (1);
1160 }
1161 }
1162
1163 mpcx_clear (mpol);
1164
1165 if (print) {
1166 printf ("Minimal polynomial:\nx^%i", c.minpoly_deg);
1167 for (i = c.minpoly_deg - 1; i >= 0; i--) {
1168 printf (" + (");
1169 mpz_out_str (stdout, 10, c.minpoly [i]);
1170 if (mpz_cmp_ui (c.minpoly_complex [i], 0) >= 0)
1171 printf ("+");
1172 mpz_out_str (stdout, 10, c.minpoly_complex [i]);
1173 printf ("*omega) * x^%i", i);
1174 }
1175 printf ("\n");
1176 }
1177 }
1178
1179 /*****************************************************************************/
1180 /* */
1181 /* computing j-invariants */
1182 /* */
1183 /*****************************************************************************/
1184
get_root_mod_P(cm_class_t c,mpz_t root,mpz_t P,bool verbose)1185 static void get_root_mod_P (cm_class_t c, mpz_t root, mpz_t P, bool verbose)
1186 /* returns a root of the minimal polynomial modulo P */
1187 /* root is changed */
1188
1189 {
1190 /* The local variables are needed only for the complex case. */
1191 int_cl_t fund;
1192 cm_timer clock;
1193 mpz_t *minpoly_P;
1194 int i;
1195 mpz_t omega, tmp;
1196 /* the second element of the integral basis of the maximal order */
1197 /* modulo P */
1198
1199 if (c.field == CM_FIELD_REAL)
1200 cm_pari_oneroot (root, c.minpoly, c.minpoly_deg, P, verbose);
1201 else {
1202 mpz_init (omega);
1203 mpz_init (tmp);
1204
1205 minpoly_P = (mpz_t*) malloc ((c.minpoly_deg + 1) * sizeof (mpz_t));
1206 for (i = 0; i <= c.minpoly_deg; i++)
1207 mpz_init (minpoly_P [i]);
1208
1209 /* compute the second element of the integral basis modulo P */
1210 cm_timer_start (clock);
1211 fund = cm_classgroup_fundamental_discriminant (c.d);
1212 cm_nt_mpz_tonelli (omega, fund, P);
1213 cm_timer_stop (clock);
1214 if (verbose)
1215 printf ("--- Time for square root: %.1f\n", cm_timer_get (clock));
1216 if (fund % 4 != 0)
1217 mpz_add_ui (omega, omega, 1ul);
1218 mpz_set_ui (tmp, 2ul);
1219 mpz_invert (tmp, tmp, P);
1220 mpz_mul (omega, omega, tmp);
1221 mpz_mod (omega, omega, P);
1222
1223 /* create the minimal polynomial modulo P */
1224 for (i = 0; i < c.minpoly_deg; i++) {
1225 mpz_mod (minpoly_P [i], c.minpoly_complex [i], P);
1226 mpz_mul (tmp, minpoly_P [i], omega);
1227 mpz_mod (minpoly_P [i], tmp, P);
1228 mpz_add (tmp, minpoly_P [i], c.minpoly [i]);
1229 mpz_mod (minpoly_P [i], tmp, P);
1230 }
1231 mpz_set_ui (minpoly_P [c.minpoly_deg], 1ul);
1232
1233 cm_pari_oneroot (root, minpoly_P, c.minpoly_deg, P, verbose);
1234
1235 mpz_clear (omega);
1236 mpz_clear (tmp);
1237 for (i = 0; i <= c.minpoly_deg; i++)
1238 mpz_clear (minpoly_P [i]);
1239 free (minpoly_P);
1240 }
1241
1242 if (verbose) {
1243 printf ("Root: ");
1244 mpz_out_str (stdout, 10, root);
1245 printf ("\n");
1246 }
1247 }
1248
1249 /*****************************************************************************/
1250
cm_class_get_j_mod_P(int_cl_t d,char inv,mpz_t P,int * no,const char * modpoldir,bool readwrite,bool verbose)1251 mpz_t* cm_class_get_j_mod_P (int_cl_t d, char inv, mpz_t P, int *no,
1252 const char* modpoldir, bool readwrite, bool verbose)
1253 /* If readwrite is true, tries to read the polynomial from a file, and */
1254 /* computes it and writes it to the file if it is not yet present. */
1255
1256 {
1257 cm_class_t c;
1258 mpz_t *j;
1259 mpz_t root, d_mpz, tmp, tmp2;
1260 cm_timer clock;
1261
1262 cm_class_init (&c, d, inv, verbose);
1263 if (!readwrite || !cm_class_read (c))
1264 cm_class_compute_minpoly (c, false, readwrite, false, verbose);
1265
1266 cm_timer_start (clock);
1267 mpz_init (root);
1268 if (inv != CM_INVARIANT_WEBER || c.p [0] != 3 || c.p [1] != 1)
1269 /* avoid special case of Weber polynomial factoring over extension */
1270 /* of degree 3; handled in weber_cm_get_j_mod_P */
1271 get_root_mod_P (c, root, P, verbose);
1272 switch (inv)
1273 {
1274 case CM_INVARIANT_J:
1275 j = (mpz_t*) malloc (sizeof (mpz_t));
1276 mpz_init_set (j [0], root);
1277 *no = 1;
1278 break;
1279 case CM_INVARIANT_GAMMA2:
1280 j = (mpz_t*) malloc (sizeof (mpz_t));
1281 mpz_init_set (j [0], root);
1282 mpz_powm_ui (j [0], j [0], 3, P);
1283 *no = 1;
1284 break;
1285 case CM_INVARIANT_GAMMA3:
1286 j = (mpz_t*) malloc (sizeof (mpz_t));
1287 mpz_init (j [0]);
1288
1289 mpz_init_set_si (d_mpz, d);
1290 mpz_init (tmp);
1291 mpz_init (tmp2);
1292
1293 mpz_powm_ui (tmp, root, 2, P);
1294 mpz_invert (tmp2, d_mpz, P);
1295 mpz_mul (root, tmp, tmp2);
1296 mpz_add_ui (root, root, 1728);
1297 mpz_mod (j [0], root, P);
1298
1299 *no = 1;
1300
1301 mpz_clear (d_mpz);
1302 mpz_clear (tmp);
1303 mpz_clear (tmp2);
1304 break;
1305 case CM_INVARIANT_WEBER:
1306 j = weber_cm_get_j_mod_P (c, root, P, no, verbose);
1307 break;
1308 case CM_INVARIANT_DOUBLEETA:
1309 #if 0
1310 case CM_INVARIANT_RAMIFIED:
1311 #endif
1312 if (c.s != c.e)
1313 mpz_powm_ui (root, root, (unsigned long int) (c.s / c.e), P);
1314 j = cm_get_j_mod_P_from_modular (no, modpoldir, CM_MODPOL_DOUBLEETA,
1315 c.p [0] * c.p [1], root, P);
1316 break;
1317 case CM_INVARIANT_MULTIETA:
1318 {
1319 int N = 1, i;
1320 for (i = 0; c.p [i] != 0; i++)
1321 N *= c.p [i];
1322 if (c.s != c.e)
1323 mpz_powm_ui (root, root, (unsigned long int) (c.s / c.e), P);
1324 j = cm_get_j_mod_P_from_modular (no, modpoldir,
1325 CM_MODPOL_MULTIETA, N, root, P);
1326 break;
1327 }
1328 case CM_INVARIANT_ATKIN:
1329 j = cm_get_j_mod_P_from_modular (no, modpoldir, CM_MODPOL_ATKIN,
1330 c.p [0], root, P);
1331 break;
1332 case CM_INVARIANT_SIMPLEETA:
1333 j = simpleeta_cm_get_j_mod_P (c, root, P, no);
1334 break;
1335 default: /* should not occur */
1336 printf ("class_cm_get_j_mod_P called for unknown class ");
1337 printf ("invariant\n");
1338 exit (1);
1339 }
1340 mpz_clear (root);
1341 cm_timer_stop (clock);
1342 if (verbose) {
1343 int i;
1344
1345 printf ("%i candidate", *no);
1346 if (*no > 1)
1347 printf ("s");
1348 printf (" for j:");
1349 for (i = 0; i < *no; i++) {
1350 printf ("\n ");
1351 mpz_out_str (stdout, 10, j [i]);
1352 }
1353 printf ("\n");
1354 printf ("--- Time for j: %.1f\n", cm_timer_get (clock));
1355 }
1356
1357 cm_class_clear (&c);
1358
1359 return j;
1360 }
1361
1362 /*****************************************************************************/
1363
cm_get_j_mod_P_from_modular(int * no,const char * modpoldir,char type,int level,mpz_t root,mpz_t P)1364 static mpz_t* cm_get_j_mod_P_from_modular (int *no, const char* modpoldir,
1365 char type, int level, mpz_t root, mpz_t P)
1366 /* computes the possible j-values as roots of the modular polynomial */
1367 /* specified by type and level */
1368
1369 {
1370 mpz_t *j;
1371 mpz_t* poly_j;
1372 /* a polynomial one of whose roots is j mod P */
1373 int n, i;
1374
1375 poly_j = cm_modpol_read_specialised_mod (&n, level, type, P, root,
1376 modpoldir);
1377 j = cm_pari_find_roots (poly_j, n, P, no);
1378 for (i = 0; i <= n; i++)
1379 mpz_clear (poly_j [i]);
1380 free (poly_j);
1381
1382 return j;
1383 }
1384
1385 /*****************************************************************************/
1386
weber_cm_get_j_mod_P(cm_class_t c,mpz_t root,mpz_t P,int * no,bool verbose)1387 static mpz_t* weber_cm_get_j_mod_P (cm_class_t c, mpz_t root, mpz_t P, int *no,
1388 bool verbose)
1389
1390 {
1391 mpz_t* j = (mpz_t*) malloc (sizeof (mpz_t));
1392 mpz_t f24, tmp;
1393 mpfp_t one;
1394
1395 mpz_init (j [0]);
1396 mpz_init (f24);
1397 mpz_init (tmp);
1398
1399 if (c.p [0] != 3 || c.p [1] != 1) {
1400 if (c.d % 3 == 0)
1401 mpz_powm_ui (f24, root, 2ul, P);
1402 else
1403 mpz_powm_ui (f24, root, 6ul, P);
1404
1405 if (c.p [0] != 7 || c.p [1] != 1) {
1406 if (c.p [0] == 1) {
1407 mpz_mul_2exp (tmp, f24, 3ul);
1408 mpz_mod (tmp, tmp, P);
1409 mpz_powm_ui (f24, tmp, 2ul, P);
1410 mpz_sub_ui (j [0], f24, 16ul);
1411 mpz_powm_ui (j [0], j [0], 3ul, P);
1412 mpz_invert (tmp, f24, P);
1413 mpz_mul (j [0], j [0], tmp);
1414 mpz_mod (j [0], j [0], P);
1415 }
1416 else if (c.p [0] == 3) {
1417 mpz_powm_ui (tmp, f24, 4ul, P);
1418 mpz_set (f24, tmp);
1419 mpz_sub_ui (j [0], f24, 16ul);
1420 mpz_powm_ui (j [0], j [0], 3ul, P);
1421 mpz_invert (tmp, f24, P);
1422 mpz_mul (j [0], j [0], tmp);
1423 mpz_mod (j [0], j [0], P);
1424 }
1425 else if (c.p [0] == 5) {
1426 mpz_mul_2exp (tmp, f24, 6ul);
1427 mpz_set (f24, tmp);
1428 mpz_sub_ui (j [0], f24, 16ul);
1429 mpz_powm_ui (j [0], j [0], 3ul, P);
1430 mpz_invert (tmp, f24, P);
1431 mpz_mul (j [0], j [0], tmp);
1432 mpz_mod (j [0], j [0], P);
1433 }
1434 else if (c.p [0] == 7) {
1435 mpz_mul_2exp (tmp, f24, 3ul);
1436 mpz_mod (tmp, tmp, P);
1437 mpz_powm_ui (f24, tmp, 4ul, P);
1438 mpz_sub_ui (j [0], f24, 16ul);
1439 mpz_powm_ui (j [0], j [0], 3ul, P);
1440 mpz_invert (tmp, f24, P);
1441 mpz_mul (j [0], j [0], tmp);
1442 mpz_mod (j [0], j [0], P);
1443 }
1444 else if (c.p [0] == 2 || c.p [0] == 6) {
1445 mpz_mul_2exp (tmp, f24, 3ul);
1446 mpz_mod (tmp, tmp, P);
1447 mpz_powm_ui (f24, tmp, 2ul, P);
1448 mpz_add_ui (j [0], f24, 16ul);
1449 mpz_powm_ui (j [0], j [0], 3ul, P);
1450 mpz_invert (tmp, f24, P);
1451 mpz_mul (j [0], j [0], tmp);
1452 mpz_mod (j [0], j [0], P);
1453 }
1454 else {
1455 /* c.p [0] == 4 */
1456 mpz_mul_2exp (tmp, f24, 9ul);
1457 mpz_set (f24, tmp);
1458 mpz_add_ui (j [0], f24, 16ul);
1459 mpz_powm_ui (j [0], j [0], 3ul, P);
1460 mpz_invert (tmp, f24, P);
1461 mpz_mul (j [0], j [0], tmp);
1462 mpz_mod (j [0], j [0], P);
1463 }
1464 }
1465 else {
1466 /* d/4 = 1 (mod 8), where d/4 is the real complex multiplication */
1467 /* discriminant */
1468 mpz_invert (tmp, f24, P);
1469 mpz_powm_ui (f24, tmp, 4ul, P);
1470 mpz_sub_ui (j [0], f24, 16ul);
1471 mpz_powm_ui (j [0], j [0], 3ul, P);
1472 mpz_invert (tmp, f24, P);
1473 mpz_mul (j [0], j [0], tmp);
1474 mpz_mod (j [0], j [0], P);
1475 }
1476 }
1477 else {
1478 /* d/4 = 5 (mod 8), we need to factor over an extension of degree 3 */
1479 mpz_t factor [4];
1480 int i;
1481 mpfpx_t root_poly, f24_poly, tmp_poly, tmp2_poly;
1482 const unsigned long int x2 [] = {0, 0, 1};
1483 const unsigned long int x6 [] = {0, 0, 0, 0, 0, 0, 1};
1484
1485 for (i = 0; i <= 3; i++)
1486 mpz_init (factor [i]);
1487 mpfpx_type_init (P, MPFPX_KARATSUBA);
1488 mpfpx_init (root_poly);
1489 mpfpx_init (f24_poly);
1490 mpfpx_init (tmp_poly);
1491 mpfpx_init (tmp2_poly);
1492 mpfp_init (one);
1493 mpfp_set_ui (one, 1ul);
1494
1495 cm_pari_onefactor (factor, c.minpoly, c.minpoly_deg, 3, P, verbose);
1496 mpfpx_set_z_array (root_poly, factor, 4);
1497 if (verbose) {
1498 printf ("Root: ");
1499 mpfpx_out (root_poly);
1500 }
1501
1502 if (c.d % 3 == 0)
1503 mpfpx_set_ui_array (f24_poly, x2, 3);
1504 else {
1505 mpfpx_set_ui_array (f24_poly, x6, 7);
1506 mpfpx_div_r (f24_poly, f24_poly, root_poly, one);
1507 }
1508
1509 mpfpx_invert (f24_poly, f24_poly, root_poly, one);
1510 mpfpx_mul_ui (f24_poly, f24_poly, 8ul);
1511 mpfpx_sqr (f24_poly, f24_poly);
1512 mpfpx_div_r (f24_poly, f24_poly, root_poly, one);
1513 mpfpx_sqr (f24_poly, f24_poly);
1514 mpfpx_div_r (f24_poly, f24_poly, root_poly, one);
1515 mpfpx_invert (tmp_poly, f24_poly, root_poly, one);
1516 mpfpx_sub_ui (f24_poly, f24_poly, 16ul);
1517 mpfpx_sqr (tmp2_poly, f24_poly);
1518 mpfpx_div_r (tmp2_poly, tmp2_poly, root_poly, one);
1519 mpfpx_mul (f24_poly, tmp2_poly, f24_poly);
1520 mpfpx_div_r (f24_poly, f24_poly, root_poly, one);
1521 mpfpx_mul (f24_poly, f24_poly, tmp_poly);
1522 mpfpx_div_r (f24_poly, f24_poly, root_poly, one);
1523 mpfpx_coeff_get_z (j [0], f24_poly, 0);
1524
1525 for (i = 0; i <= 3; i++)
1526 mpz_clear (factor [i]);
1527 mpfpx_clear (root_poly);
1528 mpfpx_clear (f24_poly);
1529 mpfpx_clear (tmp_poly);
1530 mpfpx_clear (tmp2_poly);
1531 mpfp_clear (one);
1532 }
1533
1534 *no = 1;
1535
1536 mpz_clear (f24);
1537 mpz_clear (tmp);
1538
1539 return j;
1540 }
1541
1542 /*****************************************************************************/
1543
simpleeta_cm_get_j_mod_P(cm_class_t c,mpz_t root,mpz_t P,int * no)1544 static mpz_t* simpleeta_cm_get_j_mod_P (cm_class_t c, mpz_t root, mpz_t P,
1545 int *no)
1546
1547 {
1548 mpz_t* j = (mpz_t*) malloc (sizeof (mpz_t));
1549 mpz_t f3, tmp;
1550
1551 mpz_init (j [0]);
1552 mpz_init (f3);
1553 mpz_init (tmp);
1554
1555 mpz_powm_ui (root, root, (unsigned long int) (c.s / c.e), P);
1556
1557 if (c.p [0] == 3) {
1558 mpz_add_ui (f3, root, 3ul);
1559 mpz_powm_ui (f3, f3, 3ul, P);
1560 mpz_invert (tmp, root, P);
1561 mpz_mul_ui (tmp, tmp, 27ul);
1562 mpz_add_ui (tmp, tmp, 1ul);
1563 mpz_mul (j [0], f3, tmp);
1564 mpz_mod (j [0], j [0], P);
1565 }
1566 else if (c.p [0] == 5) {
1567 mpz_add_ui (tmp, root, 10ul);
1568 mpz_mul (f3, root, tmp);
1569 mpz_mod (f3, f3, P);
1570 mpz_add_ui (f3, f3, 5ul);
1571 mpz_powm_ui (f3, f3, 3ul, P);
1572 mpz_invert (tmp, root, P);
1573 mpz_mul (j [0], f3, tmp);
1574 mpz_mod (j [0], j [0], P);
1575 }
1576 else if (c.p [0] == 7) {
1577 mpz_add_ui (tmp, root, 5ul);
1578 mpz_mul (f3, root, tmp);
1579 mpz_mod (f3, f3, P);
1580 mpz_add_ui (f3, f3, 1ul);
1581 mpz_powm_ui (f3, f3, 3ul, P);
1582 mpz_add_ui (j [0], root, 13ul);
1583 mpz_mul (j [0], j [0], root);
1584 mpz_mod (j [0], j [0], P);
1585 mpz_add_ui (j [0], j [0], 49ul);
1586 mpz_mul (j [0], j [0], f3);
1587 mpz_mod (j [0], j [0], P);
1588 mpz_invert (tmp, root, P);
1589 mpz_mul (j [0], j [0], tmp);
1590 mpz_mod (j [0], j [0], P);
1591 }
1592 else if (c.p [0] == 13) {
1593 mpz_add_ui (f3, root, 7ul);
1594 mpz_mul (f3, f3, root);
1595 mpz_mod (f3, f3, P);
1596 mpz_add_ui (f3, f3, 20ul);
1597 mpz_mul (f3, f3, root);
1598 mpz_mod (f3, f3, P);
1599 mpz_add_ui (f3, f3, 19ul);
1600 mpz_mul (f3, f3, root);
1601 mpz_mod (f3, f3, P);
1602 mpz_add_ui (f3, f3, 1ul);
1603 mpz_powm_ui (f3, f3, 3ul, P);
1604 mpz_add_ui (j [0], root, 5ul);
1605 mpz_mul (j [0], j [0], root);
1606 mpz_mod (j [0], j [0], P);
1607 mpz_add_ui (j [0], j [0], 13ul);
1608 mpz_mul (j [0], j [0], f3);
1609 mpz_mod (j [0], j [0], P);
1610 mpz_invert (tmp, root, P);
1611 mpz_mul (j [0], j [0], tmp);
1612 mpz_mod (j [0], j [0], P);
1613 }
1614 else if (c.p [0] == 4) {
1615 mpz_add_ui (f3, root, 16ul);
1616 mpz_mul (f3, f3, root);
1617 mpz_mod (f3, f3, P);
1618 mpz_invert (tmp, f3, P);
1619 mpz_add_ui (f3, f3, 16ul);
1620 mpz_powm_ui (f3, f3, 3ul, P);
1621 mpz_mul (j [0], tmp, f3);
1622 mpz_mod (j [0], j [0], P);
1623 }
1624 else if (c.p [0] == 9) {
1625 mpz_add_ui (tmp, root, 9ul);
1626 mpz_mul (tmp, tmp, root);
1627 mpz_mod (tmp, tmp, P);
1628 mpz_add_ui (tmp, tmp, 27ul);
1629 mpz_mul (f3, f3, root);
1630 mpz_mod (f3, f3, P);
1631 mpz_invert (j [0], tmp, P);
1632 mpz_add_ui (f3, tmp, 3ul);
1633 mpz_add_ui (tmp, root, 3ul);
1634 mpz_mul (f3, f3, tmp);
1635 mpz_mod (f3, f3, P);
1636 mpz_powm_ui (f3, f3, 3ul, P);
1637 mpz_mul (j [0], j [0], f3);
1638 mpz_mod (j [0], j [0], P);
1639 }
1640 else if (c.p [0] == 25) {
1641 mpz_add_ui (f3, root, 10ul);
1642 mpz_mul (f3, f3, root);
1643 mpz_mod (f3, f3, P);
1644 mpz_add_ui (f3, f3, 55ul);
1645 mpz_mul (f3, f3, root);
1646 mpz_mod (f3, f3, P);
1647 mpz_add_ui (f3, f3, 200ul);
1648 mpz_mul (f3, f3, root);
1649 mpz_mod (f3, f3, P);
1650 mpz_add_ui (f3, f3, 525ul);
1651 mpz_mul (f3, f3, root);
1652 mpz_mod (f3, f3, P);
1653 mpz_add_ui (f3, f3, 1010ul);
1654 mpz_mul (f3, f3, root);
1655 mpz_mod (f3, f3, P);
1656 mpz_add_ui (f3, f3, 1425ul);
1657 mpz_mul (f3, f3, root);
1658 mpz_mod (f3, f3, P);
1659 mpz_add_ui (f3, f3, 1400ul);
1660 mpz_mul (f3, f3, root);
1661 mpz_mod (f3, f3, P);
1662 mpz_add_ui (f3, f3, 875ul);
1663 mpz_mul (f3, f3, root);
1664 mpz_mod (f3, f3, P);
1665 mpz_add_ui (f3, f3, 250ul);
1666 mpz_mul (f3, f3, root);
1667 mpz_mod (f3, f3, P);
1668 mpz_add_ui (f3, f3, 5ul);
1669 mpz_powm_ui (f3, f3, 3ul, P);
1670 mpz_add_ui (tmp, root, 5ul);
1671 mpz_mul (tmp, tmp, root);
1672 mpz_mod (tmp, tmp, P);
1673 mpz_add_ui (tmp, tmp, 15ul);
1674 mpz_mul (tmp, tmp, root);
1675 mpz_mod (tmp, tmp, P);
1676 mpz_add_ui (tmp, tmp, 25ul);
1677 mpz_mul (tmp, tmp, root);
1678 mpz_mod (tmp, tmp, P);
1679 mpz_add_ui (tmp, tmp, 25ul);
1680 mpz_mul (tmp, tmp, root);
1681 mpz_mod (tmp, tmp, P);
1682 mpz_invert (j [0], tmp, P);
1683 mpz_mul (j [0], j [0], f3);
1684 mpz_mod (j [0], j [0], P);
1685 }
1686
1687 *no = 1;
1688
1689 mpz_clear (f3);
1690 mpz_clear (tmp);
1691
1692 return j;
1693 }
1694
1695 /*****************************************************************************/
1696 /*****************************************************************************/
1697