Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
10 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 : #include "pari.h"
14 : #include "paripriv.h"
15 :
16 : /*********************************************************************/
17 : /** **/
18 : /** PERIODS OF HYPERELLIPTIC CURVES **/
19 : /** contributed by Pascal Molin (2019) **/
20 : /** **/
21 : /*********************************************************************/
22 :
23 : /*********************************************************************/
24 : /* */
25 : /* Symplectic pairing and basis */
26 : /* */
27 : /*********************************************************************/
28 :
29 : /* compute symplectic homology basis */
30 :
31 : /* exchange rows and columns i,j, in place */
32 : static void
33 357 : row_swap(GEN m, long i1, long i2)
34 : {
35 357 : long k, l = lg(m);
36 1785 : for (k = 1; k < l; k++)
37 1428 : swap(gcoeff(m,i1,k),gcoeff(m,i2,k));
38 357 : }
39 :
40 : static void
41 714 : col_swap(GEN m, long j1, long j2)
42 : {
43 714 : swap(gel(m,j1),gel(m,j2));
44 714 : }
45 :
46 : static void
47 357 : swap_step(GEN p, GEN m, long i1, long i2)
48 : {
49 357 : if (i1==i2)
50 0 : return;
51 357 : col_swap(p, i1, i2);
52 357 : row_swap(m, i1, i2);
53 357 : col_swap(m, i1, i2);
54 : }
55 :
56 : /* ci1 <- u ci1 + v ci2, ci2 <- u1 ci1 + v1 ci2, in place */
57 : static void
58 581 : row_bezout(GEN m, long i1, long i2, GEN u, GEN v, GEN u1, GEN v1)
59 : {
60 581 : long j, l = lg(m);
61 2905 : for (j = 1; j < l; j++)
62 : {
63 2324 : GEN mi1j = gcoeff(m,i1,j), mi2j = gcoeff(m,i2,j);
64 2324 : gcoeff(m,i1,j) = gadd(gmul(u,mi1j),gmul(v,mi2j));
65 2324 : gcoeff(m,i2,j) = gadd(gmul(u1,mi1j),gmul(v1,mi2j));
66 : }
67 581 : }
68 :
69 : static void
70 1162 : col_bezout(GEN m, long j1, long j2, GEN u, GEN v, GEN u1, GEN v1)
71 : {
72 1162 : GEN mj1 = gel(m,j1), mj2 = gel(m,j2);
73 1162 : gel(m,j1) = gadd(gmul(u,mj1),gmul(v,mj2));
74 1162 : gel(m,j2) = gadd(gmul(u1,mj1),gmul(v1,mj2));
75 1162 : }
76 :
77 : /* (i,k) <- (u i + v * k, u1 i + v1 k)*/
78 : static void
79 0 : bezout_step(GEN p, GEN m, long i, long k, GEN a, GEN b)
80 : {
81 : GEN d,u,v,u1,v1;
82 0 : d = gbezout(a,b,&u,&v);
83 0 : u1 = gneg(gdiv(b,d));
84 0 : v1 = gdiv(a,d);
85 0 : col_bezout(p,i,k,u,v,u1,v1);
86 0 : row_bezout(m,i,k,u,v,u1,v1);
87 0 : col_bezout(m,i,k,u,v,u1,v1);
88 0 : }
89 :
90 : /* i <- i + q * k */
91 : static void
92 581 : transvect_step(GEN p, GEN m, long i, long k, GEN q)
93 : {
94 581 : GEN u = gen_1, v = q, u1 = gen_0, v1 = gen_1;
95 581 : col_bezout(p,i,k,u,v,u1,v1);
96 581 : row_bezout(m,i,k,u,v,u1,v1);
97 581 : col_bezout(m,i,k,u,v,u1,v1);
98 581 : }
99 :
100 : /* choose +/-1 or smallest element */
101 : static long
102 700 : pivot_line(GEN m, long i, long len)
103 : {
104 700 : long j, jx = 0;
105 700 : GEN x = NULL;
106 2100 : for (j = 1; j <= len; j++)
107 : {
108 2100 : GEN z = gcoeff(m,i,j);
109 2100 : if (signe(z)==0)
110 1400 : continue;
111 700 : if (is_pm1(z))
112 700 : return j;
113 0 : if (x == NULL || abscmpii(x, z) > 0)
114 0 : x = z, jx = j;
115 : }
116 0 : return jx;
117 : }
118 :
119 : /* returns p s.t. p~*m*p = J_g(d), d vector of diagonal coefficients */
120 : static GEN
121 350 : symplectic_reduction(GEN m, long g, GEN * d)
122 : {
123 : long dim, i, j, k, len;
124 : GEN p, diag;
125 350 : pari_sp av = avma;
126 :
127 350 : len = lg(m)-1;
128 350 : if (2 * g > len)
129 0 : pari_err_DOMAIN("symplectic_reduction","dimension",">", stoi(g), stoi(len));
130 :
131 350 : p = matid(len);
132 350 : diag = zerovec(g);
133 350 : m = shallowcopy(m);
134 : /* main loop on symplectic 2-subspace */
135 1050 : for (dim = 0; dim < g; dim++)
136 : {
137 700 : int cleared = 0;
138 700 : i = 2 * dim + 1;
139 : /* lines 0..2d-1 already cleared */
140 700 : while ((j = pivot_line(m, i, len)) == 0)
141 : {
142 : /* no intersection -> move ci to end */
143 0 : swap_step(p, m, i, len);
144 0 : len--;
145 0 : if (len == 2*dim)
146 : {
147 0 : if (d) *d = diag;
148 0 : return gc_all(av, d ? 2 : 1, &p, d);
149 : }
150 : }
151 :
152 : /* move j to i + 1 */
153 700 : if (j != i+1)
154 0 : swap_step(p, m, j, i + 1);
155 700 : j = i + 1;
156 :
157 : /* make sure gcoeff(m,i,j) > 0 */
158 700 : if (signe(gcoeff(m,i,j)) < 0)
159 357 : swap_step(p, m, i, j);
160 :
161 1400 : while(!cleared)
162 : {
163 : /* clear row i */
164 1400 : for (k = j + 1; k <= len; k++)
165 700 : if (signe(gcoeff(m,i,k)))
166 : {
167 273 : if (gequal0(remii(gcoeff(m,i,k), gcoeff(m,i,j))))
168 : {
169 273 : GEN q = gneg(gdiv(gcoeff(m,i,k), gcoeff(m,i,j)));
170 273 : transvect_step(p, m, k, j, q);
171 : }
172 : else
173 0 : bezout_step(p, m, j, k, gcoeff(m,i,j), gcoeff(m,i,k));
174 : }
175 700 : cleared = 1;
176 : /* clear row j */
177 1400 : for (k = j + 1; k <= len; k++)
178 700 : if (signe(gcoeff(m,j,k)))
179 : {
180 308 : if (gequal0(remii(gcoeff(m,j,k), gcoeff(m,i,j))))
181 : {
182 308 : GEN q = gdiv(gcoeff(m,j,k), gcoeff(m,i,j));
183 308 : transvect_step(p, m, k, i, q);
184 : }
185 : else
186 : {
187 0 : bezout_step(p, m, i, k, gcoeff(m,j,i), gcoeff(m,j,k));
188 : /* warning: row i now contains some ck.cl... */
189 0 : cleared = 0;
190 : }
191 : }
192 : }
193 700 : gel(diag, dim + 1) = gcoeff(m,i,j);
194 : }
195 350 : if (d) *d = diag;
196 350 : return gc_all(av, d ? 2: 1, &p, d);
197 : }
198 :
199 : /* below is GP code from Pascal converted to C by Bill. */
200 :
201 : static GEN
202 3794 : make_Aprim(GEN A, long ia, long ib)
203 : {
204 3794 : long i, j = 0, lA = lg(A);
205 : GEN Aprim;
206 3794 : if (ib)
207 : {
208 3794 : GEN a = gel(A, ia), b = gel(A, ib);
209 3794 : Aprim = cgetg(lA-2, t_VEC);
210 25788 : for (i = 1; i < lA; ++i)
211 21994 : if (i != ia && i != ib)
212 14406 : gel(Aprim, ++j) = gdiv(gsub(gsub(gmulsg(2, gel(A, i)), a), b), gsub(b, a));
213 : }
214 : else
215 : {
216 0 : GEN a = gel(A, ia);
217 0 : Aprim = cgetg(lA-1, t_VEC);
218 0 : for (i = 1; i < lA; ++i)
219 0 : if (i != ia)
220 0 : gel(Aprim, ++j) = gdiv(gsub(a, gmulsg(2, gel(A, i))), a);
221 : }
222 3794 : return Aprim;
223 : }
224 :
225 : static GEN
226 586642 : sqrt_affinereduction(GEN Aprim, GEN z, long prec)
227 : {
228 586642 : pari_sp av = avma;
229 586642 : GEN p = gen_1;
230 586642 : long i, l = lg(Aprim), s = 0;
231 2724176 : for (i = 1; i < l; ++i)
232 : {
233 2137534 : if (signe(real_i(gel(Aprim, i))) > 0)
234 : {
235 917756 : s++;
236 917756 : p = gmul(p, gsqrt(gsub(gel(Aprim, i), z), prec));
237 : }
238 : else
239 1219778 : p = gmul(p, gsqrt(gsub(z, gel(Aprim, i)), prec));
240 : }
241 586642 : return gc_upto(av, gmul(p, powIs(s)));
242 : }
243 :
244 : static long
245 805 : intersection_abbd(GEN A, long ia, long ib, long ic, long id, long prec)
246 : {
247 805 : pari_sp av = avma;
248 805 : long k = lg(A)-1;
249 805 : GEN a = gel(A, ia), b = gel(A, ib), d = gel(A, id);
250 805 : GEN fbd = gmul(gpowgs(gsqrt(gsub(d, b), prec), k),
251 : sqrt_affinereduction(make_Aprim(A, ib, id), gen_m1, prec));
252 805 : GEN fab = gmul(gpowgs(gsqrt(gsub(b, a), prec), k),
253 : sqrt_affinereduction(make_Aprim(A, ia, ib), gen_1, prec));
254 805 : return gc_long(av, gcmpgs(imag_i(gdiv(fbd, fab)), 0) < 0 ? -1: 1);
255 : }
256 :
257 : static long
258 217 : intersection_abcb(GEN A, long ia, long ib, long ic, long id, long prec)
259 : {
260 217 : pari_sp av = avma;
261 217 : long k = lg(A)-1;
262 217 : GEN a = gel(A, ia), b = gel(A, ib), c = gel(A, ic);
263 : GEN Aprim, fab, fcb;
264 217 : Aprim = make_Aprim(A, ic, ib);
265 217 : fcb = gmul(gpowgs(gsqrt(gsub(b, c), prec), k), sqrt_affinereduction(Aprim, stoi(1), prec));
266 217 : Aprim = make_Aprim(A, ia, ib);
267 217 : fab = gmul(gpowgs(gsqrt(gsub(b, a), prec), k), sqrt_affinereduction(Aprim, stoi(1), prec));
268 217 : return gc_long(av, gcmpgs(imag_i(gdiv(fab, fcb)), 0) < 0 ? -1: 1);
269 : }
270 :
271 : static long
272 175 : intersection_abad(GEN A, long ia, long ib, long ic, long id, long prec)
273 : {
274 175 : pari_sp av = avma;
275 175 : long k = lg(A)-1;
276 175 : GEN a = gel(A, ia), b = gel(A, ib), d = gel(A, id);
277 : GEN Aprim, fab, fad;
278 175 : long i, j = 0;
279 175 : Aprim = make_Aprim(A, ia, id);
280 175 : fad = gmul(gpowgs(gsqrt(gsub(d, a), prec), k), sqrt_affinereduction(Aprim, gen_m1, prec));
281 175 : Aprim = make_Aprim(A, ia, ib);
282 1190 : for (i = 1; i <= k; ++i)
283 : {
284 1015 : if (!((i == ia) || (i == ib)))
285 665 : gel(Aprim, ++j) = gdiv(gsub(gsub(gmulsg(2, gel(A, i)), a), b), gsub(b, a));
286 : }
287 175 : fab = gmul(gpowgs(gsqrt(gsub(b, a), prec), k), sqrt_affinereduction(Aprim, gen_m1, prec));
288 175 : return gc_long(av, gcmpgs(imag_i(gdiv(fab, fad)), 0) < 0 ? -1 : 1);
289 : }
290 :
291 : /* inner intersection I[ab].I[cd] */
292 : /* assume different end points */
293 : static long
294 903 : intersection_inner(GEN A, long ia, long ib, long ic, long id, long prec)
295 : {
296 903 : pari_sp av = avma;
297 903 : GEN a = gel(A, ia), b = gel(A, ib), c = gel(A, ic), d = gel(A, id);
298 : GEN xp, fp, xpab, pb, xpcd, fpab, fpcd, p;
299 903 : long i, lA = lg(A);
300 903 : GEN cprim = gdiv(gsub(gsub(gmulsg(2, c), a), b), gsub(b, a));
301 903 : GEN dprim = gdiv(gsub(gsub(gmulsg(2, d), a), b), gsub(b, a));
302 903 : GEN imc = imag_i(cprim), imd = imag_i(dprim);
303 903 : if (signe(imc)*signe(imd) == 1)
304 434 : return 0;
305 : /* on the same side */
306 : /* p the intersection */
307 469 : xp = imag_i(gmul(gconj(cprim), gsub(dprim, cprim)));
308 469 : fp = gsub(imd, imc);
309 469 : if (gcmp(gabs(xp, prec), gabs(fp, prec)) >= 0)
310 469 : return 0;
311 : /* discard if xp not in ]-1,1[ */
312 0 : xpab = gdiv(xp, fp);
313 0 : pb = gadd(gmul(gdivgs(gsub(b, a), 2), xpab), gdivgs(gadd(b, a), 2));
314 0 : xpcd = gdiv(gsub(gsub(gmulsg(2, pb), c), d), gsub(d, c));
315 0 : p = gen_1;
316 0 : for (i = 1; i < lA; ++i)
317 0 : if (i != ia && i != ib)
318 0 : p = gmul(p, gsqrt(gsub(xpab, gdiv(gsub(gsub(gmulsg(2, gel(A, i)), a), b), gsub(b, a))), prec));
319 : /* should be in ]-1,1[ */
320 0 : fpab = gmul(gpowgs(gsqrt(gsub(b, a), prec), lA-1), p);
321 0 : p = gen_1;
322 0 : for (i = 1; i < lA; ++i)
323 0 : if (i != ic && i != id)
324 0 : p = gmul(p, gsqrt(gsub(xpcd, gdiv(gsub(gsub(gmulsg(2, gel(A, i)), c), d), gsub(d, c))), prec));
325 0 : fpcd = gmul(gpowgs(gsqrt(gsub(d, c), prec), lA-1), p);
326 0 : return gc_long(av, 2*signe(fp)*signe(real_i(gdiv(fpab, fpcd))));
327 : }
328 :
329 : static long
330 2100 : intersection(GEN A, long ia, long ib, long ic, long id, long prec)
331 : {
332 2100 : if (ia == ib || ic == id)
333 0 : return 0; /* bad entry */
334 2100 : if ((ia == ic && ib == id) || (ia == id && ib == ic))
335 0 : return 0; /* self intersection */
336 2100 : if (ia == ic)
337 175 : return intersection_abad(A, ia, ib, ic, id, prec);
338 1925 : if (ib == ic)
339 364 : return intersection_abbd(A, ia, ib, ic, id, prec);
340 1561 : if (ia == id)
341 441 : return -intersection_abbd(A, ic, id, ia, ib, prec);
342 1120 : if (ib == id)
343 217 : return intersection_abcb(A, ia, ib, ic, id, prec);
344 903 : return intersection_inner(A, ia, ib, ic, id, prec);
345 : }
346 :
347 : static GEN
348 350 : intersection_spanning(GEN A, GEN tree, long prec)
349 : {
350 350 : long i, j, n = lg(tree)-1;
351 350 : GEN res = cgetg(n+1, t_MAT);
352 1750 : for (i = 1; i <= n; ++i)
353 1400 : gel(res, i) = cgetg(n+1, t_VECSMALL);
354 1750 : for (i = 1; i <= n; ++i)
355 : {
356 1400 : coeff(res, i, i) = 0;
357 3500 : for (j = i+1; j <= n; ++j)
358 : {
359 2100 : long s = intersection(A, mael(tree, i, 1), mael(tree, i, 2),
360 2100 : mael(tree, j, 1), mael(tree, j, 2), prec);
361 2100 : coeff(res, i, j) = s;
362 2100 : coeff(res, j, i) = -s;
363 : }
364 : }
365 350 : return zm_to_ZM(res);
366 : }
367 :
368 : static GEN
369 1400 : int_periods_affinereduction(GEN C, long i1, long i2, long prec) /* vec */
370 : {
371 1400 : pari_sp av = avma;
372 1400 : long g = itos(gel(C, 1));
373 1400 : GEN A = gel(C, 2), a = gel(A, i1), b = gel(A, i2);
374 1400 : GEN h = gel(C, 4), int_points = gel(C, 5);
375 : GEN geom_factor, decprim, Aprim, res, fact;
376 1400 : long i, j, l = lg(int_points);
377 1400 : if (gcmp(real_i(a), real_i(b)) > 0)
378 0 : pari_err_BUG("hyperellperiods");
379 1400 : geom_factor = gdivgs(gsub(b, a), 2);
380 1400 : decprim = gdiv(gadd(b, a), gsub(b, a));
381 1400 : Aprim = make_Aprim(A, i1, i2);
382 1400 : res = cgetg(g+1, t_COL);
383 1400 : gel(res, 1) = ginv(sqrt_affinereduction(Aprim, gen_0, prec));
384 2800 : for (j = 1; j < g; ++j)
385 1400 : gel(res, j + 1) = gmul(gel(res, j), decprim);
386 292824 : for (i = 1; i < l; ++i)
387 : {
388 291424 : GEN x = gmael(int_points, i, 1), dx = gmael(int_points, i, 2);
389 291424 : GEN tp = gdiv(dx, sqrt_affinereduction(Aprim, x, prec));
390 291424 : GEN tm = gdiv(dx, sqrt_affinereduction(Aprim, gneg(x), prec));
391 291424 : gel(res, 1) = gadd(gel(res, 1), gadd(tp, tm));
392 582848 : for (j = 2; j <= g; ++j)
393 : {
394 291424 : tp = gmul(tp, gadd(decprim, x));
395 291424 : tm = gmul(tm, gsub(decprim, x));
396 291424 : gel(res, j) = gadd(gel(res, j), gadd(tp, tm));
397 : }
398 : }
399 1400 : fact = gdiv(mulcxI(h), gpowgs(gsqrt(geom_factor, prec), lg(Aprim)-1));
400 1400 : gel(res, 1) = gmul(gel(res, 1), fact);
401 2800 : for (j = 2; j <= g; ++j)
402 : {
403 1400 : fact = gmul(fact, geom_factor);
404 1400 : gel(res, j) = gmul(gel(res, j), fact);
405 : }
406 1400 : return gc_GEN(av, res);
407 : }
408 :
409 : static GEN
410 350 : periods_spanning(GEN C, long prec)
411 : {
412 350 : GEN tree = gel(C, 3);
413 350 : long k, n = lg(tree)-1;
414 350 : GEN res = cgetg(n+1, t_MAT);
415 1750 : for (k = 1; k <= n; ++k)
416 : {
417 1400 : long i = mael(tree, k, 1), j = mael(tree, k, 2);
418 1400 : gel(res, k) = gmul2n(int_periods_affinereduction(C, i, j, prec), 1);
419 : }
420 350 : return res;
421 : }
422 :
423 : static GEN
424 350 : t_opt(GEN D, GEN M1, GEN lambda, long prec)
425 : {
426 350 : return gasinh(gdiv(gadd(D, glog(gmulsg(4, M1), prec)), lambda), prec);
427 : }
428 :
429 : static GEN
430 350 : phi_bound(GEN tau, GEN lambda, long prec)
431 : {
432 350 : GEN lam2 = gsqr(lambda);
433 350 : GEN Ytau = gsqrt(gmul(lam2, gsin(tau, prec)), prec);
434 350 : GEN Xtau = gdiv(gsqrt(gsub(gsqr(Ytau), gmul(lam2, gsqr(gsin(tau, prec)))), prec), gtan(tau, prec));
435 350 : return gadd(gdiv(gen_2, gcos(Ytau, prec)), ginv(gsinh(Xtau, prec)));
436 : }
437 :
438 : static GEN
439 350 : h_opt(GEN D, GEN tau, GEN M2, GEN lambda, long prec)
440 : {
441 350 : return gdiv(gmul(mulsr(2, mppi(prec)), tau),
442 : glog(gaddsg(1, gmul(gmul(gmulsg(2, M2), phi_bound(tau, lambda, prec)), gexp(D, prec))), prec));
443 : }
444 :
445 : static GEN
446 350 : integration_parameters(GEN A, GEN tree, GEN tau, GEN D, GEN lambda, long prec, long *pnpoints)
447 : {
448 350 : GEN h = h_opt(D, tau, gen_1, lambda, prec);
449 350 : *pnpoints = itos(gceil(gdiv(t_opt(D, gen_2, lambda, prec), h)));
450 350 : return h;
451 : }
452 :
453 : static void
454 350 : integration_points_thsh(GEN h, long npoints, GEN lambda, long prec, GEN *ph, GEN *pres)
455 : {
456 350 : pari_sp av = avma;
457 : long k;
458 350 : GEN eh = gexp(h, prec), eh_inv = ginv(eh), ekh = gen_1, ekh_inv = gen_1;
459 350 : GEN res = cgetg(npoints+1, t_VEC);
460 73206 : for (k = 1; k <= npoints; ++k)
461 : {
462 : GEN sh, ch2, esh, esh_inv, chsh2_i, shsh2, thsh;
463 72856 : ekh = gmul(ekh, eh);
464 72856 : ekh_inv = gmul(ekh_inv, eh_inv);
465 72856 : sh = gdivgs(gsub(ekh, ekh_inv), 2);
466 72856 : ch2 = gadd(ekh, ekh_inv);
467 72856 : esh = gexp(gmul(lambda, sh), prec);
468 72856 : esh_inv = ginv(esh);
469 72856 : chsh2_i = ginv(gadd(esh, esh_inv));
470 72856 : shsh2 = gsub(esh, esh_inv);
471 72856 : thsh = gmul(shsh2, chsh2_i);
472 72856 : gel(res, k) = mkvec2(thsh, gmul(ch2, chsh2_i));
473 : }
474 350 : *ph = gmul(h, lambda);
475 350 : *pres = res;
476 350 : (void) gc_all(av, 2, ph, pres);
477 350 : }
478 :
479 : static GEN
480 18900 : tau_3(GEN ai, GEN aj, GEN ak, GEN lambda, long prec)
481 : {
482 18900 : GEN zk = gdiv(gsub(gsub(gmulsg(2, ak), ai), aj), gsub(aj, ai));
483 18900 : GEN xItau = gasinh(gdiv(gatanh(zk, prec), lambda), prec);
484 18900 : return gabs(imag_i(xItau), prec);
485 : }
486 :
487 : static GEN
488 4900 : tau_edge(GEN A, long i, long j, GEN lambda, long prec)
489 : {
490 4900 : GEN tauij = utoipos(4);
491 4900 : long l = lg(A), k;
492 33600 : for (k = 1; k < l; ++k)
493 : {
494 28700 : if (k == i || k == j)
495 9800 : continue;
496 18900 : tauij = gmin(tauij, tau_3(gel(A, i), gel(A, j), gel(A, k), lambda, prec));
497 : }
498 4900 : return tauij;
499 : }
500 :
501 : static void
502 350 : max_spanning(GEN A, long nedges, GEN lambda, long prec, GEN *ptree, GEN *ptaumin)
503 : {
504 350 : pari_sp av = avma;
505 : GEN tau_v, tau_c, per;
506 : long z;
507 : GEN taken, tree, taumin;
508 : long i, j, k;
509 350 : long n = lg(A)-1, m = (n*(n-1))>>1;
510 350 : tau_v = cgetg(m+1, t_VEC);
511 350 : tau_c = cgetg(m+1, t_VEC);
512 350 : z = 1;
513 2380 : for (i = 1; i <= n; ++i)
514 6930 : for (j = i+1; j <= n; ++j)
515 : {
516 4900 : gel(tau_v, z) = tau_edge(A, i, j, lambda, prec);
517 4900 : if (gcmp(real_i(gel(A, i)), real_i(gel(A, j))) > 0)
518 1260 : gel(tau_c, z) = mkvecsmall2(j, i);
519 : else
520 3640 : gel(tau_c, z) = mkvecsmall2(i, j);
521 4900 : z++;
522 : }
523 350 : per = indexsort(tau_v);
524 350 : tau_v = vecpermute(tau_v, per);
525 350 : tau_c = vecpermute(tau_c, per);
526 350 : taken = zero_Flv(n);
527 350 : tree = cgetg(nedges+1, t_VEC);
528 350 : taken[mael(tau_c, m, 1)] = 1;
529 350 : taumin = gel(tau_v, m);
530 1750 : for (k = 1; k <= nedges; ++k)
531 : {
532 1400 : z = m;
533 3955 : while(taken[mael(tau_c, z, 1)]+taken[mael(tau_c, z, 2)] != 1)
534 2555 : z--;
535 1400 : gel(tree, k) = gel(tau_c, z);
536 1400 : taumin = gmin_shallow(taumin, gel(tau_v, z));
537 1400 : taken[mael(tau_c, z, 1)] = taken[mael(tau_c, z, 2)] = 1;
538 : }
539 350 : *ptree = tree;
540 350 : *ptaumin = taumin;
541 350 : (void)gc_all(av, 2, ptaumin, ptree);
542 350 : }
543 :
544 : static GEN
545 350 : periodmat(GEN P, long prec)
546 : {
547 350 : pari_sp av = avma;
548 350 : GEN l = leading_coeff(P), A = roots(P, prec);
549 : GEN hc, lambda;
550 : GEN tree, tau, h, coh1x_homC, IntC, ABtoC;
551 350 : long g = (lg(A)-2)/2;
552 : long npoints;
553 350 : lambda = Pi2n(-1, DEFAULTPREC);
554 350 : max_spanning(A, 2*g, lambda, DEFAULTPREC, &tree, &tau);
555 350 : h = integration_parameters(A, tree, tau, stoi(prec), lambda,
556 : DEFAULTPREC, &npoints);
557 350 : h = rtor(h, prec);
558 350 : lambda = Pi2n(-1,prec);
559 350 : hc = mkvec5(stoi(g), A, tree, gen_0, gen_0);
560 350 : integration_points_thsh(h, npoints, lambda, prec, &gel(hc, 4), &gel(hc, 5));
561 350 : coh1x_homC = periods_spanning(hc, prec);
562 350 : IntC = intersection_spanning(A, tree, prec);
563 350 : ABtoC = symplectic_reduction(IntC, g, NULL);
564 350 : return gc_upto(av, gdiv(gmul(coh1x_homC, ABtoC), gmul2n(gsqrt(l, prec),-1)));
565 : }
566 : static long
567 357 : hyperellgenus(GEN H)
568 357 : { long d = degpol(H); return ((d+1)>>1)-1; }
569 :
570 : /* return 4P + Q^2 */
571 : static GEN
572 7 : check_hyperell(GEN PQ)
573 : {
574 : GEN H;
575 7 : if (is_vec_t(typ(PQ)) && lg(PQ)==3)
576 0 : H = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
577 : else
578 7 : H = gmul2n(PQ, 2);
579 7 : return typ(H) == t_POL? H: NULL;
580 : }
581 :
582 : static GEN
583 350 : genus2BSDperiod(GEN C, long prec)
584 : {
585 350 : pari_sp av = avma;
586 : forsubset_t iter;
587 350 : GEN PQ, P, M, Om = NULL, Omr, d, r, v;
588 : long g;
589 350 : PQ = hyperellminimalmodel(C, NULL, NULL);
590 350 : P = gadd(gmulsg(4, gel(PQ, 1)), gsqr(gel(PQ, 2)));
591 350 : g = hyperellgenus(P);
592 350 : M = real_i(periodmat(P, prec));
593 350 : forsubset_init(&iter, mkvec2s(2*g,g));
594 938 : while ((v = forsubset_next(&iter)))
595 : {
596 938 : Om = extract0(M, identity_perm(g), v);
597 938 : if (expo(gabs(det(Om), prec))>-(prec>>1))
598 350 : break;
599 : }
600 350 : Omr = bestappr(gmul(ginv(Om), M), int2n(prec>>1));
601 350 : d = denominator(Omr, NULL);
602 350 : r = gabs(gdiv(gmul(det(Om), det(mathnf0(gmul(Omr, d), 0))), gsqr(d)), prec);
603 350 : return gc_upto(av, r);
604 : }
605 :
606 : GEN
607 357 : hyperellperiods(GEN C, long flag, long prec)
608 : {
609 357 : pari_sp av = avma;
610 : GEN M, H;
611 : long g;
612 357 : if (flag==2) return genus2BSDperiod(C, prec);
613 7 : H = check_hyperell(C);
614 7 : if (!H) pari_err_TYPE("hyperellperiods", C);
615 7 : if (flag<0 || flag>1) pari_err_FLAG("hyperellperiods");
616 7 : g = hyperellgenus(H); if (g < 1) pari_err_DOMAIN("hyperellperiods","genus","=",gen_0,C);
617 0 : M = periodmat(H, prec);
618 0 : if (flag==0) M = gauss(vecslice(M,1,g), vecslice(M,g+1,2*g));
619 0 : return gc_upto(av, M);
620 : }
|
