1 /*
2
3 modpol.c - code for handling modular polynomials
4
5 Copyright (C) 2009, 2012 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 unsigned long int read_gz_ui (gzFile f);
27 static void read_gz_mpz (mpz_t rop, gzFile f);
28
29 /*****************************************************************************/
30
read_gz_ui(gzFile f)31 static unsigned long int read_gz_ui (gzFile f)
32 /* reads and returns an unsigned long integer in base 10 from the gzipped */
33 /* file f */
34 {
35 char c;
36 unsigned long int res = 0ul;
37
38 do {
39 c = gzgetc (f);
40 if (c == EOF) {
41 printf ("*** EOF in 'modpol_read_gz_ui'\n");
42 exit (1);
43 }
44 } while (isspace (c));
45 do {
46 if (c < '0' || c > '9') {
47 printf ("*** Non-digit in 'modpol_read_gz_ui'\n");
48 exit (1);
49 }
50 res = 10 * res + (c - '0');
51 c = gzgetc (f);
52 } while (!(isspace (c) || c == EOF));
53
54 return res;
55 }
56
57 /*****************************************************************************/
58
read_gz_mpz(mpz_t rop,gzFile f)59 static void read_gz_mpz (mpz_t rop, gzFile f)
60 /* reads a coefficient of type mpz in base 10 from the gzipped file f */
61 {
62 char c;
63 int sign = 1;
64
65 mpz_set_ui (rop, 0);
66
67 do {
68 c = gzgetc (f);
69 if (c == EOF) {
70 printf ("*** EOF in 'modpol_read_gz_mpz'\n");
71 exit (1);
72 }
73 } while (isspace (c));
74 if (c == '-') {
75 sign = -1;
76 c = gzgetc (f);
77 if (c == EOF) {
78 printf ("*** EOF in 'modpol_read_gz_mpz'\n");
79 exit (1);
80 }
81 }
82 do {
83 if (c < '0' || c > '9') {
84 printf ("*** Non-digit in 'modpol_read_gz_mpz'\n");
85 exit (1);
86 }
87 mpz_mul_ui (rop, rop, 10ul);
88 mpz_add_ui (rop, rop, (unsigned long int) (c - '0'));
89 c = gzgetc (f);
90 } while (!(isspace (c) || c == EOF));
91
92 if (sign == -1)
93 mpz_neg (rop, rop);
94 }
95
96 /*****************************************************************************/
97
cm_modpol_read_specialised_mod(int * n,int level,char type,mpz_t p,mpz_t x,const char * datadir)98 mpz_t* cm_modpol_read_specialised_mod (int* n, int level, char type, mpz_t p,
99 mpz_t x, const char* datadir)
100 /* returns the modular polynomial of the given type and level modulo p */
101 /* and evaluated in the first argument x; n is replaced by its degree */
102
103 {
104 char filename [255];
105 gzFile f;
106 int lev, N, i, k;
107 char c;
108 mpz_t* xpow;
109 mpz_t tmp;
110 mpz_t* res;
111
112 sprintf (filename, "%s/%cf/%cf_%.4i.dat.gz", datadir,
113 type, type, level);
114 cm_file_gzopen_read (&f, filename);
115
116 lev = read_gz_ui (f);
117 if (lev != level) {
118 printf ("*** Trying to read modular polynomial of level %i ", level);
119 printf ("in a file for the level %i!\n", lev);
120 exit (1);
121 }
122 c = gzgetc (f);
123 if (c != type) {
124 printf ("*** Trying to read modular polynomial of type %c ", type);
125 printf ("in a file for the type %c!\n", c);
126 exit (1);
127 }
128
129 N = read_gz_ui (f);
130 xpow = (mpz_t *) malloc ((N + 1) * sizeof (mpz_t));
131 mpz_init_set_ui (xpow [0], 1ul);
132 for (i = 1; i <= N; i++) {
133 mpz_init (xpow [i]);
134 mpz_mul (xpow [i], xpow [i-1], x);
135 mpz_mod (xpow [i], xpow [i], p);
136 }
137 mpz_init (tmp);
138
139 *n = read_gz_ui (f);
140 res = (mpz_t*) malloc ((*n+1) * sizeof (mpz_t));
141 for (i = 0; i <= *n; i++)
142 mpz_init_set_ui (res [i], 0ul);
143 do {
144 i = read_gz_ui (f);
145 k = read_gz_ui (f);
146 read_gz_mpz (tmp, f);
147 mpz_mul (tmp, tmp, xpow [i]);
148 mpz_add (res [k], res [k], tmp);
149 mpz_mod (res [k], res [k], p);
150 } while (k != 0 || i != 0);
151 /* we assume that the last entry in the file is the constant one */
152
153 for (i = 0; i <= N; i++)
154 mpz_clear (xpow [i]);
155 free (xpow);
156 mpz_clear (tmp);
157
158 cm_file_gzclose (f);
159
160 return res;
161 }
162
163 /*****************************************************************************/
164
cm_modpol_print_pari(int level,char type,const char * datadir)165 void cm_modpol_print_pari (int level, char type, const char* datadir)
166 /* prints the modular polynomial of the given type and level in the */
167 /* pari seadata format */
168
169 {
170 char filename [255];
171 gzFile f;
172 int lev, i_old, i, k;
173 char c;
174 mpz_t tmp;
175
176 if (type != 'a') {
177 printf ("*** Trying to read modular polynomial of type %c ", type);
178 printf ("instead of 'a'!\n");
179 exit (1);
180 }
181 sprintf (filename, "%s/%cf/%cf_%.4i.dat.gz", datadir,
182 type, type, level);
183 cm_file_gzopen_read (&f, filename);
184
185 lev = read_gz_ui (f);
186 if (lev != level) {
187 printf ("*** Trying to read modular polynomial of level %i ", level);
188 printf ("in a file for the level %i!\n", lev);
189 exit (1);
190 }
191 c = gzgetc (f);
192 if (c != type) {
193 printf ("*** Trying to read modular polynomial of type '%c' ", type);
194 printf ("in a file for the type %c!\n", c);
195 exit (1);
196 }
197
198 /* skip N and n */
199 read_gz_ui (f);
200 read_gz_ui (f);
201 mpz_init (tmp);
202
203 printf ("[%i, \"A\", [", level);
204 i_old = level + 2;
205 do {
206 i = read_gz_ui (f);
207 k = read_gz_ui (f);
208 read_gz_mpz (tmp, f);
209 if (i != i_old && k != 0)
210 printf ("[");
211 mpz_out_str (stdout, 10, tmp);
212 if (k != 0)
213 printf (", ");
214 else {
215 if (i == i_old)
216 printf ("]");
217 if (i != 0)
218 printf (", ");
219 }
220 i_old = i;
221 } while (k != 0 || i != 0);
222 /* we assume that the last entry in the file is the constant one */
223 printf ("]]");
224
225 cm_file_gzclose (f);
226 }
227
228 /*****************************************************************************/
229
cm_modpol_print_magma(int level,char type,const char * datadir)230 void cm_modpol_print_magma (int level, char type, const char* datadir)
231 /* prints the modular polynomial of the given type and level in magma */
232 /* format */
233
234 {
235 char filename [255];
236 gzFile f;
237 int lev, i, k;
238 char c;
239 mpz_t tmp;
240
241 sprintf (filename, "%s/%cf/%cf_%.4i.dat.gz", datadir,
242 type, type, level);
243 cm_file_gzopen_read (&f, filename);
244
245 lev = read_gz_ui (f);
246 if (lev != level) {
247 printf ("*** Trying to read modular polynomial of level %i ", level);
248 printf ("in a file for the level %i!\n", lev);
249 exit (1);
250 }
251 c = gzgetc (f);
252 if (c != type) {
253 printf ("*** Trying to read modular polynomial of type %c ", type);
254 printf ("in a file for the type %c!\n", c);
255 exit (1);
256 }
257
258 /* skip N and n */
259 read_gz_ui (f);
260 read_gz_ui (f);
261 mpz_init (tmp);
262
263 printf ("phi :=");
264 do {
265 i = read_gz_ui (f);
266 k = read_gz_ui (f);
267 read_gz_mpz (tmp, f);
268 printf (" +(");
269 mpz_out_str (stdout, 10, tmp);
270 printf (")*F^%i*J^%i", i, k);
271 } while (k != 0 || i != 0);
272 /* we assume that the last entry in the file is the constant one */
273 printf (";\n");
274
275 cm_file_gzclose (f);
276 }
277
278 /*****************************************************************************/
279