1 /* jinvariants.c -- computation of invariants from theta constants
2 *
3 * Copyright (C) 2006, 2011, 2012 INRIA
4 *
5 * This file is part of CMH.
6 *
7 * CMH is free software; you can redistribute it and/or modify it under
8 * the terms of the GNU General Public License as published by the
9 * Free Software Foundation; either version 3 of the License, or (at your
10 * option) any later version.
11 *
12 * CMH is distributed in the hope that it will be useful, but WITHOUT ANY
13 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for
15 * more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see http://www.gnu.org/licenses/ .
19 */
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <assert.h>
24
25 #include <mpc.h>
26 #include <time.h>
27
28 #include "params.h"
29 #include "macros.h"
30 #include "jinvariants.h"
31
product_of_six(mpc_t v,mpc_t * t,int i0,int i1,int i2,int i3,int i4,int i5)32 static void product_of_six(mpc_t v, mpc_t * t, int i0, int i1, int i2, int i3, int i4, int i5)
33 {
34 cset_ui(v, 1);
35 cmul(v, v, t[i0]);
36 cmul(v, v, t[i1]);
37 cmul(v, v, t[i2]);
38 cmul(v, v, t[i3]);
39 cmul(v, v, t[i4]);
40 cmul(v, v, t[i5]);
41 }
42
sum_of_four(mpc_t v,mpc_t * t,int i0,int i1,int i2,int i3)43 static void sum_of_four(mpc_t v, mpc_t * t, int i0, int i1, int i2, int i3)
44 {
45 cset_ui(v, 0);
46 cadd(v, v, t[i0]);
47 cadd(v, v, t[i1]);
48 cadd(v, v, t[i2]);
49 cadd(v, v, t[i3]);
50 }
51
52 /* {{{ compute h12 and h16 by enumerating all Göpel quadruples */
h12h16_from_th4th8(mpc_t * h,mpc_t * th4,mpc_t * th8)53 static void h12h16_from_th4th8(mpc_t * h, mpc_t * th4, mpc_t * th8)
54 {
55 mpc_ptr h12 = h[0];
56 mpc_ptr h16 = h[1];
57
58 mpc_t v[15];
59 int prec = cprec(h12);
60 mpc_t aux;
61
62 cinit(aux, prec);
63
64 for (int i=0; i<15; i++)
65 cinit(v[i], prec);
66
67 cset_ui(h12, 0);
68 cset_ui(h16, 0);
69
70 product_of_six(v[0], th4, 0, 1, 4, 5, 8, 9);
71 sum_of_four(aux, th8, 2, 3, 6, 7);
72 cadd(h12, h12, v[0]);
73 cmul(aux, aux, v[0]);
74 cadd(h16, h16, aux);
75
76 product_of_six(v[1], th4, 0, 1, 2, 5, 7, 8);
77 sum_of_four(aux, th8, 3, 4, 6, 9);
78 cadd(h12, h12, v[1]);
79 cmul(aux, aux, v[1]);
80 cadd(h16, h16, aux);
81
82 product_of_six(v[2], th4, 4, 5, 6, 7, 8, 9);
83 sum_of_four(aux, th8, 0, 1, 2, 3);
84 cadd(h12, h12, v[2]);
85 cmul(aux, aux, v[2]);
86 cadd(h16, h16, aux);
87
88 product_of_six(v[3], th4, 1, 2, 4, 5, 6, 7);
89 sum_of_four(aux, th8, 0, 3, 8, 9);
90 cadd(h12, h12, v[3]);
91 cmul(aux, aux, v[3]);
92 cadd(h16, h16, aux);
93
94 product_of_six(v[4], th4, 0, 1, 2, 4, 6, 9);
95 sum_of_four(aux, th8, 3, 5, 7, 8);
96 cadd(h12, h12, v[4]);
97 cmul(aux, aux, v[4]);
98 cadd(h16, h16, aux);
99
100 product_of_six(v[5], th4, 0, 2, 6, 7, 8, 9);
101 sum_of_four(aux, th8, 1, 3, 4, 5);
102 cadd(h12, h12, v[5]);
103 cmul(aux, aux, v[5]);
104 cadd(h16, h16, aux);
105
106 product_of_six(v[6], th4, 0, 1, 3, 4, 7, 9);
107 sum_of_four(aux, th8, 2, 5, 6, 8);
108 cadd(h12, h12, v[6]);
109 cmul(aux, aux, v[6]);
110 cadd(h16, h16, aux);
111
112 product_of_six(v[7], th4, 0, 2, 3, 4, 7, 8);
113 sum_of_four(aux, th8, 1, 5, 6, 9);
114 cadd(h12, h12, v[7]);
115 cmul(aux, aux, v[7]);
116 cadd(h16, h16, aux);
117
118 product_of_six(v[8], th4, 2, 3, 4, 5, 8, 9);
119 sum_of_four(aux, th8, 0, 1, 6, 7);
120 cadd(h12, h12, v[8]);
121 cmul(aux, aux, v[8]);
122 cadd(h16, h16, aux);
123
124 product_of_six(v[9], th4, 1, 2, 3, 5, 7, 9);
125 sum_of_four(aux, th8, 0, 4, 6, 8);
126 cadd(h12, h12, v[9]);
127 cmul(aux, aux, v[9]);
128 cadd(h16, h16, aux);
129
130 product_of_six(v[10], th4, 0, 1, 3, 5, 6, 8);
131 sum_of_four(aux, th8, 2, 4, 7, 9);
132 cadd(h12, h12, v[10]);
133 cmul(aux, aux, v[10]);
134 cadd(h16, h16, aux);
135
136 product_of_six(v[11], th4, 0, 3, 4, 5, 6, 7);
137 sum_of_four(aux, th8, 1, 2, 8, 9);
138 cadd(h12, h12, v[11]);
139 cmul(aux, aux, v[11]);
140 cadd(h16, h16, aux);
141
142 product_of_six(v[12], th4, 1, 3, 6, 7, 8, 9);
143 sum_of_four(aux, th8, 0, 2, 4, 5);
144 cadd(h12, h12, v[12]);
145 cmul(aux, aux, v[12]);
146 cadd(h16, h16, aux);
147
148 product_of_six(v[13], th4, 1, 2, 3, 4, 6, 8);
149 sum_of_four(aux, th8, 0, 5, 7, 9);
150 cadd(h12, h12, v[13]);
151 cmul(aux, aux, v[13]);
152 cadd(h16, h16, aux);
153
154 product_of_six(v[14], th4, 0, 2, 3, 5, 6, 9);
155 sum_of_four(aux, th8, 1, 4, 7, 8);
156 cadd(h12, h12, v[14]);
157 cmul(aux, aux, v[14]);
158 cadd(h16, h16, aux);
159
160 for (int i=0; i<15; i++) cclear(v[i]);
161 cclear(aux);
162 }
163 /* }}} */
164
h4h10h12h16_from_th2(mpc_t * h,mpc_t * th2)165 void h4h10h12h16_from_th2(mpc_t * h, mpc_t * th2)/*{{{*/
166 {
167 mpc_ptr h4 = h[0];
168 mpc_ptr h10 = h[1];
169 // mpc_ptr h12 = h[2];
170 // mpc_ptr h16 = h[3];
171
172 mpc_t th4[10], th8[10];
173 int prec = cprec(h4);
174
175 for (int i=0; i<10; i++)
176 {
177 cinit(th4[i], prec);
178 cinit(th8[i], prec);
179 }
180 for (int i=0; i<10; i++) {
181 csqr(th4[i], th2[i]);
182 csqr(th8[i], th4[i]);
183 }
184
185 // sum of ^8
186 cset_ui(h4, 0);
187 for (int i=0; i<10; i++)
188 cadd(h4, h4, th8[i]);
189
190 // product of ^2
191 cset_ui(h10, 1);
192 for (int i=0; i<10; i++)
193 cmul(h10, h10, th2[i]);
194
195 h12h16_from_th4th8(h + 2, th4, th8);
196
197 for (int i=0; i<10; i++)
198 {
199 cclear(th4[i]);
200 cclear(th8[i]);
201 }
202 }/*}}}*/
203
I2I4I6I10_from_h4h10h12h16(mpc_t * I,mpc_t * h)204 void I2I4I6I10_from_h4h10h12h16(mpc_t * I, mpc_t * h)/*{{{*/
205 {
206 mpc_ptr h4 = h[0];
207 mpc_ptr h10 = h[1];
208 mpc_ptr h12 = h[2];
209 mpc_ptr h16 = h[3];
210 mpc_ptr I2 = I[0];
211 mpc_ptr I4 = I[1];
212 mpc_ptr I6 = I[2];
213 mpc_ptr I10 = I[3];
214
215 int prec = cprec(I2);
216
217 mpc_t aux;
218 cinit(aux, prec);
219
220 cinv(aux, h10);
221 cmul(I2, h12, aux); // I2 == h12/h10
222 cset(I4, h4); // I4 == h4
223 cmul(I6, h16, aux); // I16 = h16 / h10
224 cset(I10, h10); // I10 = h10
225
226 cclear(aux);
227 }/*}}}*/
228
gorenlauter_i1i2i3_from_I2I4I6I10(mpc_t * absolute_invariants,mpc_t * I)229 void gorenlauter_i1i2i3_from_I2I4I6I10(mpc_t * absolute_invariants, mpc_t * I)
230 {
231 mpc_ptr i1 = absolute_invariants[0];
232 mpc_ptr i2 = absolute_invariants[1];
233 mpc_ptr i3 = absolute_invariants[2];
234 mpc_ptr I2 = I[0];
235 mpc_ptr I4 = I[1];
236 mpc_ptr I6 = I[2];
237 mpc_ptr I10 = I[3];
238
239 int prec = cprec(I2); // not i1 !!!
240
241 mpc_t aux;
242 mpc_t I10_inv;
243 mpc_t I2_2;
244 mpc_t I2_3;
245
246 cinit(aux, prec);
247 cinit(I10_inv, prec);
248 cinit(I2_2, prec);
249 cinit(I2_3, prec);
250
251 cinv(I10_inv, I10);
252 csqr(I2_2, I2);
253 cmul(I2_3, I2_2, I2);
254
255 cmul(aux, I2_3, I2_2);
256 cmul(i1, aux, I10_inv); // jinv[0] == I2^5/I10 == i1 in streng notes
257
258 cmul(aux, I2_3, I4);
259 cmul(i2, aux, I10_inv); // jinv[1] == I2^3*I4/I10 == i2 in streng notes
260
261 cmul(aux, I6, I2_2);
262 cmul(i3, aux, I10_inv); // jinv[2] == I2^2*I6/I10 == i3 in streng notes
263
264 cclear(aux);
265 cclear(I2_3);
266 cclear(I2_2);
267 cclear(I10_inv);
268 }
269
kohel_j1j2j3_from_I2I4I6I10(mpc_t * absolute_invariants,mpc_t * I)270 void kohel_j1j2j3_from_I2I4I6I10(mpc_t * absolute_invariants, mpc_t * I)
271 {
272 mpc_ptr j1 = absolute_invariants[0];
273 mpc_ptr j2 = absolute_invariants[1];
274 mpc_ptr j3 = absolute_invariants[2];
275 mpc_ptr I2 = I[0];
276 mpc_ptr I4 = I[1];
277 mpc_ptr I6 = I[2];
278 mpc_ptr I10 = I[3];
279
280 int prec = cprec(I2); // not i1 !!!
281
282 mpc_t aux;
283 mpc_t I10_inv;
284 mpc_t I2_2;
285 mpc_t I2_3;
286
287 cinit(aux, prec);
288 cinit(I10_inv, prec);
289 cinit(I2_2, prec);
290 cinit(I2_3, prec);
291
292 cinv(I10_inv, I10);
293 csqr(I2_2, I2);
294 cmul(I2_3, I2_2, I2);
295
296 cmul(aux, I4, I10_inv);
297 cmul(j2, aux, I2_3);
298
299 cmul(j1, aux, I6);
300
301 cmul(aux, I6, I10_inv);
302 cmul(j3, aux, I2_2);
303
304 cclear(aux);
305 cclear(I2_3);
306 cclear(I2_2);
307 cclear(I10_inv);
308 }
309
streng_i1i2i3_from_I2I4I6I10(mpc_t * absolute_invariants,mpc_t * I)310 void streng_i1i2i3_from_I2I4I6I10(mpc_t * absolute_invariants, mpc_t * I)
311 {
312 mpc_ptr i1 = absolute_invariants[0];
313 mpc_ptr i2 = absolute_invariants[1];
314 mpc_ptr i3 = absolute_invariants[2];
315 mpc_ptr I2 = I[0];
316 mpc_ptr I4 = I[1];
317 mpc_ptr I6 = I[2];
318 mpc_ptr I10 = I[3];
319
320 // old-Streng -2^10I4I6'/I10, 2^18I4^5/I10^2,2^28/3I2I4^2/I10
321 // Streng I4I6'/I10, I2I4^2/I10, I4^5/I10^2
322 //
323 int prec = cprec(I2); // not i1 !!!
324
325 mpc_t aux, aux2, I4_div_I10;
326
327 cinit(aux, prec);
328 cinit(aux2, prec);
329 cinit(I4_div_I10, prec);
330
331 cinv(I4_div_I10, I10);
332 cmul(I4_div_I10, I4_div_I10, I4);
333
334 cmul_ui(aux, I6, 3);
335 cmul(aux2, I2, I4);
336 csub(aux, aux2, aux);
337 cdiv_2ui(aux, aux, 1);
338 // aux = I6' = (I2I4-3I6)/2
339
340 // i1 = I4I6'/I10 = I6' * I4/I10
341 cmul(i1, aux, I4_div_I10);
342
343 // i2 = I2I4^2/I10 = I2I4 * I4/I10
344 cmul(i2, aux2, I4_div_I10);
345
346 // i3 = I4^5/I10^2 = I4^2 * I4 * (I4/I10)^2
347 csqr(aux, I4);
348 cmul(aux, aux, I4);
349 csqr(aux2, I4_div_I10);
350 cmul(i3, aux, aux2);
351
352 cclear(aux);
353 cclear(aux2);
354 cclear(I4_div_I10);
355 }
356
invariants_from_th2(mpc_t * jinv,mpc_t * th2,int invariant_set)357 void invariants_from_th2(mpc_t * jinv, mpc_t * th2, int invariant_set)
358 {
359 mpc_t h[4];
360 mpc_t I[4];
361
362 int prec = cprec(jinv[0])+3*SEC_MARGIN;
363
364 for(int i = 0 ; i < 4 ; i++) cinit(h[i], prec);
365 for(int i = 0 ; i < 4 ; i++) cinit(I[i], prec);
366
367 // since most invariants are computed from I, factor it out. We could
368 // also make it appear in the switch statement, of course.
369 h4h10h12h16_from_th2(h, th2);
370 I2I4I6I10_from_h4h10h12h16(I, h);
371
372 switch(invariant_set) {
373 case 0: // Goren-Lauter / van Wamelen (I2^5/I10,I2^3I4/I10,I2^2I6/I10)
374 gorenlauter_i1i2i3_from_I2I4I6I10(jinv, I);
375 break;
376 case 1: // Kohel I4I6/I10,I2^3I4/I10,I2^2I6/I10
377 kohel_j1j2j3_from_I2I4I6I10(jinv, I);
378 break;
379 case 2:
380 // old-Streng -2^10I4I6'/I10, 2^18I4^5/I10^2,2^28/3I2I4^2/I10
381 // Streng I4I6'/I10, I2I4^2/I10, I4^5/I10^2
382 // ( I6' = (I2I4-3I6)/2)
383 streng_i1i2i3_from_I2I4I6I10(jinv, I);
384 break;
385 default:
386 fprintf(stderr, "No such invariant set defined !!!\n");
387 abort();
388 }
389
390 for(int i = 0 ; i < 4 ; i++) cclear(h[i]);
391 for(int i = 0 ; i < 4 ; i++) cclear(I[i]);
392
393 return;
394 }
395
396