1 /*
2
3 modular.c - code for evaluating modular functions in floating point arguments
4
5 Copyright (C) 2009, 2010, 2015, 2016 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_common-impl.h"
25
26 static void modular_fundamental_matrix (csrcptr z,cm_matrix_t *M);
27 static void modular_fundamental_domain_matrix (ctype z, cm_matrix_t *M);
28 static void atkin_eval (cm_modular_t m, ctype rop, ctype op,
29 unsigned long int l);
30
31 /*****************************************************************************/
32 /* */
33 /* internal functions */
34 /* */
35 /*****************************************************************************/
36
modular_fundamental_matrix(csrcptr z,cm_matrix_t * M)37 static void modular_fundamental_matrix (csrcptr z,cm_matrix_t *M)
38 /* computes the transformation matrix M such that Mz is in the */
39 /* fundamental domain */
40
41 {
42 bool ok = false;
43 long int tmp_int;
44 ftype tmp_fr;
45 ctype local_z;
46
47 finit (tmp_fr, (mp_exp_t) 100);
48 cinit (local_z, (mp_exp_t) 100);
49
50 M->a = 1;
51 M->b = 0;
52 M->c = 0;
53 M->d = 1;
54
55 /* determine the matrix from a low precision approximation of z */
56 cset (local_z, z);
57 while (!ok) {
58 ok = true;
59 /* obtain -0.5 <= real part <= 0.5 */
60 fround (tmp_fr, crealref (local_z));
61 tmp_int = -fget_si (tmp_fr);
62 cadd_si (local_z, local_z, tmp_int);
63 /* multiply M from the left by T^tmp_int */
64 M->a += M->c * tmp_int;
65 M->b += M->d * tmp_int;
66
67 cnorm (tmp_fr, local_z);
68 if (fcmp_d (tmp_fr, 0.999) < 0) {
69 /* apply S */
70 cneg (local_z, local_z);
71 cui_div (local_z, 1ul, local_z);
72 tmp_int = M->a;
73 M->a = -M->c;
74 M->c = tmp_int;
75 tmp_int = M->b;
76 M->b = -M->d;
77 M->d = tmp_int;
78 ok = false;
79 }
80 }
81
82 /* normalise the matrix */
83 if (M->c < 0 || (M->c == 0 && M->d < 0)) {
84 M->a = -M->a;
85 M->b = -M->b;
86 M->c = -M->c;
87 M->d = -M->d;
88 }
89
90 fclear (tmp_fr);
91 cclear (local_z);
92 }
93
94 /*****************************************************************************/
95
modular_fundamental_domain_matrix(ctype z,cm_matrix_t * M)96 static void modular_fundamental_domain_matrix (ctype z, cm_matrix_t *M)
97 /* transforms z into the fundamental domain and returns the inverse */
98 /* transformation matrix M */
99
100 {
101 long int tmp_int;
102 ctype tmp_c1;
103
104 cinit (tmp_c1, cget_prec (z));
105
106 modular_fundamental_matrix (z, M);
107
108 /* apply the matrix to z */
109 cmul_si (tmp_c1, z, M->a);
110 cadd_si (tmp_c1, tmp_c1, M->b);
111 cmul_si (z, z, M->c);
112 cadd_si (z, z, M->d);
113 cdiv (z, tmp_c1, z);
114
115 /* invert and normalize the matrix */
116 tmp_int = M->a;
117 M->a = M->d;
118 M->d = tmp_int;
119 M->b = -M->b;
120 M->c = -M->c;
121 if (M->c < 0 || (M->c == 0 && M->d < 0)) {
122 M->a = -M->a;
123 M->b = -M->b;
124 M->c = -M->c;
125 M->d = -M->d;
126 }
127
128 cclear (tmp_c1);
129 }
130
131 /*****************************************************************************/
132
cm_modular_fundamental_domain(ctype z)133 void cm_modular_fundamental_domain (ctype z)
134 /* transforms z into the fundamental domain */
135
136 {
137 ctype tmp_c1;
138 cm_matrix_t M;
139
140 cinit (tmp_c1, cget_prec (z));
141
142 modular_fundamental_matrix (z, &M);
143
144 /* apply the matrix to z */
145 cmul_si (tmp_c1, z, M.a);
146 cadd_si (tmp_c1, tmp_c1, M.b);
147 cmul_si (z, z, M.c);
148 cadd_si (z, z, M.d);
149 cdiv (z, tmp_c1, z);
150
151 cclear (tmp_c1);
152 }
153
154 /*****************************************************************************/
155 /* */
156 /* creating and freeing variables */
157 /* */
158 /*****************************************************************************/
159
cm_modular_init(cm_modular_t * m,fprec_t prec)160 void cm_modular_init (cm_modular_t *m, fprec_t prec)
161
162 {
163 int i;
164
165 m->prec = prec;
166
167 finit (m->pi, prec);
168 cinit (m->twopii, prec);
169 cinit (m->log_zeta24, prec);
170 cinit (m->zeta48inv, prec);
171
172 fconst_pi (m->pi);
173 cset_ui_ui (m->twopii, 0ul, 0ul);
174 fmul_2ui (cimagref (m->twopii), m->pi, 1ul);
175 cset_ui_ui (m->log_zeta24, 0ul, 0ul);
176 fdiv_ui (cimagref (m->log_zeta24), m->pi, 12ul);
177 cdiv_ui (m->zeta48inv, m->log_zeta24, 2ul);
178 cneg (m->zeta48inv, m->zeta48inv);
179 cexp (m->zeta48inv, m->zeta48inv);
180
181 cinit (m->zeta24 [0], prec);
182 cset_ui_ui (m->zeta24 [0], 1ul, 0ul);
183 cinit (m->zeta24 [1], prec);
184 cexp (m->zeta24 [1], m->log_zeta24);
185 for (i = 2; i < 24; i++)
186 {
187 cinit (m->zeta24 [i], prec);
188 cmul (m->zeta24 [i], m->zeta24 [i-1], m->zeta24 [1]);
189 }
190
191 finit (m->sqrt2, prec);
192 fsqrt_ui (m->sqrt2, 2ul);
193
194 cm_qdev_init (&(m->eta), prec);
195 }
196
197 /*****************************************************************************/
198
cm_modular_clear(cm_modular_t * m)199 void cm_modular_clear (cm_modular_t *m)
200
201 {
202 int i;
203
204 cclear (m->log_zeta24);
205 cclear (m->twopii);
206 fclear (m->pi);
207 cclear (m->zeta48inv);
208 for (i = 0; i < 24; i++)
209 cclear (m->zeta24 [i]);
210 fclear (m->sqrt2);
211 cm_qdev_clear (&(m->eta));
212 }
213
214 /*****************************************************************************/
215 /* */
216 /* external functions */
217 /* */
218 /*****************************************************************************/
219
cm_modular_eta_transform(long int * e,ctype czplusd,ctype z,cm_matrix_t M)220 int cm_modular_eta_transform (long int *e, ctype czplusd, ctype z,
221 cm_matrix_t M)
222 /* Compute the transformation between eta (Mz) and eta (z), that is,
223 the exponent e and the value c z + d such that
224 eta (Mz) = zeta_24^e * sqrt (c z + d) * eta (z).
225 A previous version computed the complex value eta (Mz) / eta (z);
226 but for quotients of powers of eta functions, it may be more
227 efficient to handle all transformation values together.
228 The return value is 0 if M is the identity matrix and 1 otherwise
229 so that the transformation quotient is different from 1. */
230
231 {
232 long int c1, lambda;
233 int a, b, c, d;
234
235 /* Check the case that M is the identity matrix, which happens quite
236 often (particularly for j). */
237 if (M.a == 1 && M.b == 0 && M.c == 0 && M.d == 1) {
238 *e = 0;
239 cset_ui (czplusd, 1);
240 return (0);
241 }
242 else {
243 /* Reduce the coefficients for the computation of zeta_exp to avoid
244 an overflow.
245 M.a is reduced modulo 48 since the formula contains (M.a*M.a-1)/2. */
246 a = M.a % 48;
247 b = M.b % 24;
248 c = M.c % 24;
249 d = M.d % 24;
250
251 /* Compute the exponent. */
252 if (M.c == 0) {
253 c1 = 1;
254 lambda = 1;
255 }
256 else {
257 c1 = M.c;
258 lambda = 0;
259 while (c1 % 2 == 0) {
260 c1 /= 2;
261 lambda++;
262 }
263 }
264 *e = a * b + c * (d * (1 - a*a) - a) + 3 * (c1 % 8) * (a - 1)
265 + (3 * lambda * (a*a - 1)) / 2;
266 if (cm_nt_kronecker (M.a, c1) == -1)
267 *e += 12;
268 *e %= 24;
269 if (*e < 0)
270 *e += 24;
271
272 /* Compute M.c * z + M.d. */
273 cmul_si (czplusd, z, M.c);
274 cadd_si (czplusd, czplusd, M.d);
275
276 return (1);
277 }
278 }
279
280 /*****************************************************************************/
281
cm_modular_eta_series(cm_modular_t m,ctype rop,ctype q_24)282 void cm_modular_eta_series (cm_modular_t m, ctype rop, ctype q_24)
283 /* evaluates the power series for eta with q_24 standing for the 24th */
284 /* root of q; uses an automatically optimised addition chain. */
285 /* All computations are carried out with the precision of rop, that must */
286 /* have the same precision for its real and its imaginary part. */
287
288 {
289 ctype factor;
290
291 cinit (factor, cget_prec (rop));
292
293 if (fzero_p (cimagref (q_24))) {
294 /* avoid cpow_ui, since it calls cpow to determine the sign of 0 */
295 fpow_ui (crealref (factor), crealref (q_24), 24ul);
296 fset_ui (cimagref (factor), 0ul);
297 }
298 else
299 cpow_ui (factor, q_24, 24ul);
300 cm_qdev_eval (factor, m.eta, factor);
301 cmul (rop, q_24, factor);
302
303 cclear (factor);
304 }
305
306 /*****************************************************************************/
307
cm_modular_eta_series_fr(cm_modular_t m,ftype rop,ftype q_24)308 void cm_modular_eta_series_fr (cm_modular_t m, ftype rop, ftype q_24)
309 /* evaluates the power series for eta with q_24 standing for the 24th */
310 /* root of q; uses an automatically optimised addition chain. */
311 /* All computations are carried out with the precision of rop. */
312
313 {
314 ftype factor;
315
316 finit (factor, fget_prec (rop));
317
318 fpow_ui (factor, q_24, 24ul);
319 cm_qdev_eval_fr (factor, m.eta, factor);
320 fmul (rop, q_24, factor);
321
322 fclear (factor);
323 }
324
325 /*****************************************************************************/
326
cm_modular_eta_eval(cm_modular_t m,ctype rop,ctype op)327 void cm_modular_eta_eval (cm_modular_t m, ctype rop, ctype op)
328 /* evaluates the eta function at op */
329
330 {
331 cm_matrix_t M;
332 long int e;
333 int transformed;
334
335 ctype q24, op_local;
336 /* Trick: The transformation into the fundamental domain is carried out */
337 /* with the precision of op, the series evaluation with the precision of */
338 /* rop. In this way, the first computation, which has a tendency to lose */
339 /* digits, can be carried out at a higher precision. */
340 cinit (q24, cget_prec (rop));
341 cinit (op_local, cget_prec (op));
342
343 cset (op_local, op);
344 modular_fundamental_domain_matrix (op_local, &M);
345 transformed = cm_modular_eta_transform (&e, rop, op_local, M);
346 /* Now rop contains cz+d. */
347 if (transformed) {
348 csqrt (rop, rop);
349 cmul (rop, rop, m.zeta24 [e]);
350 }
351 /* workaround to efficiently handle almost real arguments; here, cexp */
352 /* cannot be improved, since the almost zero imaginary part does have an */
353 /* influence on the imaginary part of the result. */
354 if (!fzero_p (crealref (op_local)) &&
355 fget_exp (crealref (op_local))
356 < - 0.8 * ((double) fget_prec (crealref (op_local))))
357 {
358 cm_modular_eta_eval_fr (m, crealref (q24), cimagref (op_local));
359 cmul_fr (rop, rop, crealref (q24));
360 }
361 else
362 {
363 cmul (q24, m.log_zeta24, op_local);
364 cexp (q24, q24);
365 cm_modular_eta_series (m, q24, q24);
366 cmul (rop, rop, q24);
367 }
368
369 cclear (q24);
370 cclear (op_local);
371 }
372
373 /*****************************************************************************/
374
cm_modular_eta_eval_fr(cm_modular_t m,ftype rop,ftype op)375 void cm_modular_eta_eval_fr (cm_modular_t m, ftype rop, ftype op)
376 /* evaluates the eta function at the purely imaginary argument op*I, */
377 /* of course without transforming into the fundamental domain */
378
379 {
380 ftype q24;
381
382 finit (q24, fget_prec (rop));
383
384 fmul (q24, op, cimagref (m.log_zeta24));
385 fneg (q24, q24);
386 fexp (q24, q24);
387 cm_modular_eta_series_fr (m, rop, q24);
388
389 fclear (q24);
390 }
391
392 /*****************************************************************************/
393
atkin_eval(cm_modular_t m,ctype rop,ctype op,unsigned long int l)394 static void atkin_eval (cm_modular_t m, ctype rop, ctype op,
395 unsigned long int l)
396 /* evaluates eta (z)*eta (lz) */
397
398 {
399 ctype tmp, rop_local, op_local;
400
401 cinit (tmp, cget_prec (rop));
402 cinit (rop_local, cget_prec (rop));
403 cinit (op_local, cget_prec (op));
404
405 cmul_ui (op_local, op, l);
406 cm_modular_eta_eval (m, rop_local, op_local);
407 cm_modular_eta_eval (m, tmp, op);
408 cmul (rop, rop_local, tmp);
409
410 cclear (tmp);
411 cclear (rop_local);
412 cclear (op_local);
413 }
414
415 /*****************************************************************************/
416
cm_modular_atkinhecke_eval(cm_modular_t m,ctype rop,ctype op,unsigned long int l,unsigned long int r)417 void cm_modular_atkinhecke_eval (cm_modular_t m, ctype rop, ctype op,
418 unsigned long int l, unsigned long int r)
419 /* evaluates the quotient of the r-th Hecke operator (for r prime), */
420 /* applied to eta (z)*eta (lz), and the function itself, in the argument */
421 /* -1/z (to obtain a function for Gamma^0 (l) instead of Gamma_0 (l)) */
422 /* Expresses the numerator as a in transformed arguments. */
423
424 {
425 ctype Mz, tmp, op_local;
426 ctype rop_local;
427 unsigned long int i;
428
429 cinit (op_local, cget_prec (op));
430 cinit (Mz, cget_prec (op));
431 cinit (tmp, cget_prec (rop));
432 cinit (rop_local, cget_prec (rop));
433
434 cui_div (op_local, 1ul, op);
435 cneg (op_local, op_local);
436 cset_ui_ui (rop_local, 0ul, 0ul);
437 for (i = 0; i < r; i++) {
438 cadd_ui (Mz, op_local, 24*i);
439 cdiv_ui (Mz, Mz, r);
440 atkin_eval (m, tmp, Mz, l);
441 cadd (rop_local, rop_local, tmp);
442 }
443 cdiv_ui (rop_local, rop_local, r);
444 cmul_ui (Mz, op_local, r);
445 atkin_eval (m, tmp, Mz, l);
446 cadd (rop_local, rop_local, tmp);
447
448 atkin_eval (m, tmp, op_local, l);
449 cdiv (rop, rop_local, tmp);
450
451 cclear (op_local);
452 cclear (Mz);
453 cclear (tmp);
454 cclear (rop_local);
455 }
456
457 /*****************************************************************************/
458
cm_modular_atkinhecke_level_eval(cm_modular_t m,ctype rop,ctype op,unsigned long int l)459 void cm_modular_atkinhecke_level_eval (cm_modular_t m, ctype rop, ctype op,
460 unsigned long int l)
461 /* evaluates Atkin's optimised function for Gamma^0^* (l) */
462
463 {
464 if (l == 47) {
465 cm_modular_atkinhecke_eval (m, rop, op, 47, 17);
466 cneg (rop, rop);
467 }
468 else if (l == 59) {
469 ctype z, tmp;
470 cinit (z, cget_prec (op));
471 cinit (tmp, cget_prec (rop));
472 cset (z, op);
473 cm_modular_atkinhecke_eval (m, rop, z, 59, 5);
474 cm_modular_atkinhecke_eval (m, tmp, z, 59, 29);
475 cadd (rop, rop, tmp);
476 cadd_ui (rop, rop, 1ul);
477 cclear (z);
478 cclear (tmp);
479 }
480 else if (l == 71) {
481 ctype z, tmp;
482 cinit (z, cget_prec (op));
483 cinit (tmp, cget_prec (rop));
484 cset (z, op);
485 cm_modular_atkinhecke_eval (m, rop, z, 71, 5);
486 cm_modular_atkinhecke_eval (m, tmp, z, 71, 29);
487 cadd (rop, rop, tmp);
488 cadd_ui (rop, rop, 1ul);
489 cclear (z);
490 cclear (tmp);
491 }
492 else if (l == 131) {
493 cm_modular_atkinhecke_eval (m, rop, op, 131, 61);
494 cadd_ui (rop, rop, 1ul);
495 }
496 else {
497 printf ("*** Called cm_modular_atkinhecke_level_eval with level %li, "
498 "for which the optimal Atkin invariant is not implemented.\n", l);
499 exit (1);
500 }
501
502 }
503
504 /*****************************************************************************/
505