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 i,j, in place */
32 : static void
33 357 : row_swap(GEN M, long i, long j)
34 : {
35 357 : long k, l = lg(M);
36 1785 : for (k = 1; k < l; k++) swap(gcoeff(M,i,k), gcoeff(M,j,k));
37 357 : }
38 :
39 : static void
40 357 : swap_step(GEN P, GEN M, long i, long j)
41 : {
42 357 : if (i == j) return;
43 357 : swap(gel(P,i), gel(P,j));
44 357 : swap(gel(M,i), gel(M,j));
45 357 : row_swap(M, i, j);
46 : }
47 :
48 : /* M <- U(i,j, [u,v; u1,v1])~ * M. In place */
49 : static void
50 581 : row_bezout(GEN M, long i, long j, GEN u, GEN v, GEN u1, GEN v1)
51 : {
52 581 : long k, l = lg(M);
53 2905 : for (k = 1; k < l; k++)
54 : {
55 2324 : GEN a = gcoeff(M,i,k), b = gcoeff(M,j,k);
56 2324 : gcoeff(M,i,k) = addii(mulii(u,a), mulii(v,b));
57 2324 : gcoeff(M,j,k) = addii(mulii(u1,a),mulii(v1,b));
58 : }
59 581 : }
60 :
61 : /* M <- M * U(i,j, [u,v; u1,v1]). In place */
62 : static void
63 1162 : col_bezout(GEN M, long i, long j, GEN u, GEN v, GEN u1, GEN v1)
64 : {
65 1162 : GEN Mi = gel(M,i), Mj = gel(M,j);
66 1162 : gel(M,i) = ZC_lincomb(u, v, Mi, Mj);
67 1162 : gel(M,j) = ZC_lincomb(u1, v1, Mi, Mj);
68 1162 : }
69 :
70 : /* P <- P*U(i,j, [u,v;u1,v1]); where U[k,k] = 1 for k != i,j,
71 : * U[i,i] = u, U[i,j] = v, U[j,i] = u1, U[j,j] = v1
72 : * M <- U~ * M * U */
73 : static void
74 581 : bezout_apply(GEN P, GEN M, long i, long j, GEN u, GEN v, GEN u1, GEN v1)
75 : {
76 581 : col_bezout(P, i,j, u,v,u1,v1);
77 581 : row_bezout(M, i,j, u,v,u1,v1);
78 581 : col_bezout(M, i,j, u,v,u1,v1);
79 581 : }
80 : /* (i,j) <- (u i + v j, u1 i + v1 j)*/
81 : static void
82 0 : bezout_step(GEN P, GEN M, long i, long j, GEN a, GEN b)
83 : {
84 0 : GEN u, v, d = bezout(a,b,&u,&v);
85 0 : if (!is_pm1(d)) { a = diviiexact(a, d); b = diviiexact(b, d); }
86 0 : bezout_apply(P,M, i,j, u,v,negi(b),a);
87 0 : }
88 : /* i <- i + q * j */
89 : static void
90 581 : transvect_step(GEN P, GEN M, long i, long j, GEN q)
91 581 : { bezout_apply(P,M, i,j, gen_1,q,gen_0,gen_1); }
92 :
93 : /* index of non-zero element in m[i,]; return 0 if none exist else index
94 : * of smallest element in absolute value */
95 : static long
96 700 : row_pivot(GEN m, long i)
97 : {
98 700 : long j, jx = 0, l = lg(m);
99 700 : GEN x = NULL;
100 2100 : for (j = 1; j < l; j++)
101 : {
102 2100 : GEN z = gcoeff(m,i,j);
103 2100 : if (signe(z))
104 : {
105 700 : if (is_pm1(z)) return j; /* minimal */
106 0 : if (!x || abscmpii(x, z) > 0) { x = z; jx = j; }
107 : }
108 : }
109 0 : return jx;
110 : }
111 :
112 : /* M symplectic in M_2g(Z). Returns P such that P~*M*P = J_g(D), D a ZV */
113 : static GEN
114 350 : ZM_symplectic_reduction(GEN M, GEN *pD)
115 : {
116 350 : long dim, n = lg(M)-1, g = n >> 1; /* n will decrease */
117 350 : GEN P = matid(n), D = zerovec(g);
118 350 : pari_sp av = avma;
119 :
120 350 : M = shallowcopy(M);
121 : /* main loop on symplectic 2-subspace */
122 1050 : for (dim = 0; dim < g; dim++)
123 : {
124 700 : long j, i = 2 * dim + 1;
125 700 : int cleared = 0;
126 : /* lines 0..2d-1 already cleared */
127 700 : while ((j = row_pivot(M, i)) == 0)
128 : { /* no intersection: move M[,i] to end and decrease n */
129 0 : swap_step(P, M, i, n);
130 0 : if (--n == 2*dim) goto END;
131 : }
132 700 : if (j != i+1) { swap_step(P, M, j, i+1); j = i+1; }
133 700 : if (signe(gcoeff(M,i,j)) < 0) swap_step(P, M, i, j);
134 : /* now j = i+1 and M[i,j] > 0 */
135 :
136 1400 : while (!cleared)
137 : { /* clear row i */
138 : long k;
139 1400 : for (k = j + 1; k <= n; k++)
140 700 : if (signe(gcoeff(M,i,k)))
141 : {
142 273 : GEN r, q = dvmdii(gcoeff(M,i,k), gcoeff(M,i,j), &r);
143 273 : if (r == gen_0)
144 273 : transvect_step(P, M, k, j, negi(q));
145 : else
146 0 : bezout_step(P, M, j, k, gcoeff(M,i,j), gcoeff(M,i,k));
147 : }
148 700 : cleared = 1;
149 : /* clear row j */
150 1400 : for (k = j + 1; k <= n; k++)
151 700 : if (signe(gcoeff(M,j,k)))
152 : {
153 308 : GEN r, q = dvmdii(gcoeff(M,j,k), gcoeff(M,i,j), &r);
154 308 : if (r == gen_0)
155 308 : transvect_step(P, M, k, i, q);
156 : else
157 : {
158 0 : bezout_step(P, M, i, k, gcoeff(M,j,i), gcoeff(M,j,k));
159 0 : cleared = 0; /* M[i,] now contains some ck.cl ! */
160 : }
161 : }
162 : }
163 700 : gel(D, dim + 1) = gcoeff(M,i,j);
164 : }
165 350 : END:
166 350 : if (pD) *pD = D;
167 350 : return gc_all(av, pD ? 2: 1, &P, pD);
168 : }
169 :
170 : /* below is GP code from Pascal converted to C by Bill. */
171 :
172 : static GEN
173 8694 : make_Aprim(GEN A, long ia, long ib)
174 : {
175 8694 : long i, j = 0, lA = lg(A);
176 8694 : GEN a = gel(A, ia), b = gel(A, ib), p = gadd(b, a), m = gsub(b, a);
177 8694 : GEN Aprim = cgetg(lA - 2, t_VEC);
178 59388 : for (i = 1; i < lA; ++i)
179 50694 : if (i != ia && i != ib)
180 33306 : gel(Aprim, ++j) = gdiv(gsub(gmulsg(2, gel(A, i)), p), m);
181 8694 : return Aprim;
182 : }
183 :
184 : static GEN
185 586642 : sqrt_affinereduction(GEN Aprim, GEN z, long prec)
186 : {
187 586642 : pari_sp av = avma;
188 586642 : GEN p = gen_1;
189 586642 : long i, l = lg(Aprim), s = 0;
190 2724176 : for (i = 1; i < l; ++i)
191 : {
192 2137534 : GEN a = gel(Aprim, i);
193 2137534 : if (signe(real_i(a)) > 0) { s++; a = gsub(a, z); } else a = gsub(z, a);
194 2137534 : p = gmul(p, gsqrt(a, prec));
195 : }
196 586642 : return gc_upto(av, gmul(p, powIs(s)));
197 : }
198 :
199 : static long
200 805 : intersection_abbd(GEN A, long ia, long ib, long id, long prec)
201 : {
202 805 : pari_sp av = avma;
203 805 : long k = lg(A)-1;
204 805 : GEN a = gel(A, ia), b = gel(A, ib), d = gel(A, id);
205 805 : GEN fbd = gmul(gpowgs(gsqrt(gsub(d, b), prec), k),
206 : sqrt_affinereduction(make_Aprim(A, ib, id), gen_m1, prec));
207 805 : GEN fab = gmul(gpowgs(gsqrt(gsub(b, a), prec), k),
208 : sqrt_affinereduction(make_Aprim(A, ia, ib), gen_1, prec));
209 805 : return gc_long(av, signe(imag_i(gdiv(fbd, fab))));
210 : }
211 :
212 : static long
213 217 : intersection_abcb(GEN A, long ia, long ib, long ic, long prec)
214 : {
215 217 : pari_sp av = avma;
216 217 : long k = lg(A)-1;
217 217 : GEN a = gel(A, ia), b = gel(A, ib), c = gel(A, ic), fab, fcb;
218 217 : fcb = gmul(gpowgs(gsqrt(gsub(b, c), prec), k),
219 : sqrt_affinereduction(make_Aprim(A, ic, ib), gen_1, prec));
220 217 : fab = gmul(gpowgs(gsqrt(gsub(b, a), prec), k),
221 : sqrt_affinereduction(make_Aprim(A, ia, ib), gen_1, prec));
222 217 : return gc_long(av, signe(imag_i(gdiv(fab, fcb))));
223 : }
224 :
225 : static long
226 175 : intersection_abad(GEN A, long ia, long ib, long id, long prec)
227 : {
228 175 : pari_sp av = avma;
229 175 : long k = lg(A)-1;
230 175 : GEN a = gel(A, ia), b = gel(A, ib), d = gel(A, id), fab, fad;
231 175 : fad = gmul(gpowgs(gsqrt(gsub(d, a), prec), k),
232 : sqrt_affinereduction(make_Aprim(A, ia, id), gen_m1, prec));
233 175 : fab = gmul(gpowgs(gsqrt(gsub(b, a), prec), k),
234 : sqrt_affinereduction(make_Aprim(A, ia, ib), gen_m1, prec));
235 175 : return gc_long(av, signe(imag_i(gdiv(fab, fad))));
236 : }
237 :
238 : /* inner intersection I[ab].I[cd] */
239 : /* assume different end points */
240 : static long
241 903 : intersection_inner(GEN A, long ia, long ib, long ic, long id, long prec)
242 : {
243 903 : pari_sp av = avma;
244 903 : GEN a = gel(A, ia), b = gel(A, ib), c = gel(A, ic), d = gel(A, id);
245 : GEN xp, fp, xpab, pb, xpcd, fpab, fpcd;
246 903 : GEN bpa = gadd(b, a), bma = gsub(b, a), dpc, dmc;
247 903 : GEN cprim = gdiv(gsub(gmulsg(2, c), bpa), bma);
248 903 : GEN dprim = gdiv(gsub(gmulsg(2, d), bpa), bma);
249 903 : GEN imc = imag_i(cprim), imd = imag_i(dprim);
250 : long k;
251 903 : if (signe(imc)*signe(imd) == 1) return 0;
252 : /* on the same side */
253 : /* p the intersection */
254 469 : xp = imag_i(gmul(gconj(cprim), gsub(dprim, cprim)));
255 469 : fp = gsub(imd, imc);
256 469 : if (gcmp(gabs(xp, prec), gabs(fp, prec)) >= 0) return 0;
257 : /* discard if xp not in ]-1,1[ */
258 0 : xpab = gdiv(xp, fp); dpc = gadd(d, c); dmc = gsub(d, c);
259 0 : pb = gadd(gmul(gdivgs(bma, 2), xpab), gdivgs(bpa, 2));
260 0 : xpcd = gdiv(gsub(gmulsg(2, pb), dpc), dmc);
261 : /* should be in ]-1,1[ */
262 0 : k = lg(A)-1;
263 0 : fpab = gmul(gpowgs(gsqrt(bma, prec), k),
264 : sqrt_affinereduction(make_Aprim(A, ia, ib), xpab, prec));
265 0 : fpcd = gmul(gpowgs(gsqrt(dmc, prec), k),
266 : sqrt_affinereduction(make_Aprim(A, ic, id), xpcd, prec));
267 0 : return gc_long(av, 2*signe(fp)*signe(real_i(gdiv(fpab, fpcd))));
268 : }
269 :
270 : static long
271 2100 : intersection(GEN A, long ia, long ib, long ic, long id, long prec)
272 : {
273 2100 : if (ia == ib || ic == id) return 0; /* bad entry */
274 2100 : if (ia == ic && ib == id) return 0; /* self intersection */
275 2100 : if (ia == id && ib == ic) return 0; /* self intersection */
276 2100 : if (ia == ic) return intersection_abad(A, ia, ib, id, prec);
277 1925 : if (ib == ic) return intersection_abbd(A, ia, ib, id, prec);
278 1561 : if (ia == id) return -intersection_abbd(A, ic, id, ib, prec);
279 1120 : if (ib == id) return intersection_abcb(A, ia, ib, ic, prec);
280 903 : return intersection_inner(A, ia, ib, ic, id, prec);
281 : }
282 :
283 : static GEN
284 350 : intersection_spanning(GEN A, GEN tree, long prec)
285 : {
286 350 : long i, j, n = lg(tree)-1;
287 350 : GEN res = cgetg(n+1, t_MAT);
288 1750 : for (i = 1; i <= n; ++i)
289 1400 : gel(res, i) = cgetg(n+1, t_VECSMALL);
290 1750 : for (i = 1; i <= n; ++i)
291 : {
292 1400 : coeff(res, i, i) = 0;
293 3500 : for (j = i+1; j <= n; ++j)
294 : {
295 2100 : long s = intersection(A, mael(tree, i, 1), mael(tree, i, 2),
296 2100 : mael(tree, j, 1), mael(tree, j, 2), prec);
297 2100 : coeff(res, i, j) = s;
298 2100 : coeff(res, j, i) = -s;
299 : }
300 : }
301 350 : return zm_to_ZM(res);
302 : }
303 :
304 : static GEN
305 1400 : int_periods_affinereduction(GEN C, GEN edge, long prec)
306 : {
307 1400 : pari_sp av = avma;
308 1400 : long g = itos(gel(C, 1)), i1 = edge[1], i2 = edge[2];
309 1400 : GEN A = gel(C, 2), a = gel(A, i1), b = gel(A, i2);
310 1400 : GEN h = gel(C, 4), int_points = gel(C, 5);
311 : GEN F, geom_factor, decprim, Aprim, res;
312 1400 : long i, j, l = lg(int_points);
313 :
314 1400 : if (gcmp(real_i(a), real_i(b)) > 0) pari_err_BUG("hyperellperiods");
315 1400 : decprim = gdiv(gadd(b, a), gsub(b, a));
316 1400 : Aprim = make_Aprim(A, i1, i2);
317 1400 : res = gpowers0(decprim, g-1, ginv(sqrt_affinereduction(Aprim, gen_0, prec)));
318 292824 : for (i = 1; i < l; i++)
319 : {
320 291424 : GEN x = gmael(int_points, i, 1), dx = gmael(int_points, i, 2);
321 291424 : GEN tp = gdiv(dx, sqrt_affinereduction(Aprim, x, prec));
322 291424 : GEN tm = gdiv(dx, sqrt_affinereduction(Aprim, gneg(x), prec));
323 291424 : GEN Tp = gpowers0(gadd(decprim,x), g-1, tp);
324 291424 : GEN Tm = gpowers0(gsub(decprim,x), g-1, tm);
325 874272 : for (j = 1; j <= g; j++)
326 582848 : gel(res,j) = gadd(gel(res,j), gadd(gel(Tp,j), gel(Tm,j)));
327 : }
328 1400 : geom_factor = gdivgs(gsub(b, a), 2);
329 1400 : F = gpowers0(geom_factor, g-1,
330 1400 : gdiv(mulcxI(h), gpowgs(gsqrt(geom_factor, prec), lg(Aprim)-1)));
331 4200 : for (j = 1; j <= g; j++) gel(res,j) = gmul(gel(res,j), gel(F,j));
332 1400 : settyp(res, t_COL); return gc_GEN(av, res);
333 : }
334 :
335 : static GEN
336 350 : periods_spanning(GEN C, long prec)
337 : {
338 350 : GEN tree = gel(C, 3);
339 350 : long k, n = lg(tree)-1;
340 350 : GEN res = cgetg(n+1, t_MAT);
341 1750 : for (k = 1; k <= n; k++)
342 1400 : gel(res, k) = gmul2n(int_periods_affinereduction(C, gel(tree,k), prec), 1);
343 350 : return res;
344 : }
345 :
346 : static GEN
347 350 : t_opt(GEN D, GEN M1, GEN lambda, long prec)
348 : {
349 350 : return gasinh(gdiv(gadd(D, glog(gmulsg(4, M1), prec)), lambda), prec);
350 : }
351 :
352 : static GEN
353 350 : phi_bound(GEN tau, GEN lambda, long prec)
354 : {
355 350 : GEN lam2 = gsqr(lambda);
356 350 : GEN Ytau = gsqrt(gmul(lam2, gsin(tau, prec)), prec);
357 350 : GEN Xtau = gdiv(gsqrt(gsub(gsqr(Ytau), gmul(lam2, gsqr(gsin(tau, prec)))), prec), gtan(tau, prec));
358 350 : return gadd(gdiv(gen_2, gcos(Ytau, prec)), ginv(gsinh(Xtau, prec)));
359 : }
360 :
361 : static GEN
362 350 : h_opt(GEN D, GEN tau, GEN M2, GEN lambda, long prec)
363 : {
364 350 : return gdiv(gmul(mulsr(2, mppi(prec)), tau),
365 : glog(gaddsg(1, gmul(gmul(gmulsg(2, M2), phi_bound(tau, lambda, prec)), gexp(D, prec))), prec));
366 : }
367 :
368 : static GEN
369 350 : integration_parameters(GEN tau, GEN D, GEN lambda, long prec, long *pnpoints)
370 : {
371 350 : GEN h = h_opt(D, tau, gen_1, lambda, prec);
372 350 : *pnpoints = itos(gceil(gdiv(t_opt(D, gen_2, lambda, prec), h)));
373 350 : return h;
374 : }
375 :
376 : static void
377 350 : integration_points_thsh(GEN h, long npoints, GEN lambda, long prec, GEN *ph, GEN *pres)
378 : {
379 350 : pari_sp av = avma;
380 : long k;
381 350 : GEN eh = gexp(h, prec), eh_inv = ginv(eh), ekh = gen_1, ekh_inv = gen_1;
382 350 : GEN res = cgetg(npoints+1, t_VEC);
383 73206 : for (k = 1; k <= npoints; ++k)
384 : {
385 : GEN sh, ch2, esh, esh_inv, chsh2_i, shsh2, thsh;
386 72856 : ekh = gmul(ekh, eh);
387 72856 : ekh_inv = gmul(ekh_inv, eh_inv);
388 72856 : sh = gdivgs(gsub(ekh, ekh_inv), 2);
389 72856 : ch2 = gadd(ekh, ekh_inv);
390 72856 : esh = gexp(gmul(lambda, sh), prec);
391 72856 : esh_inv = ginv(esh);
392 72856 : chsh2_i = ginv(gadd(esh, esh_inv));
393 72856 : shsh2 = gsub(esh, esh_inv);
394 72856 : thsh = gmul(shsh2, chsh2_i);
395 72856 : gel(res, k) = mkvec2(thsh, gmul(ch2, chsh2_i));
396 : }
397 350 : *ph = gmul(h, lambda);
398 350 : *pres = res;
399 350 : (void) gc_all(av, 2, ph, pres);
400 350 : }
401 :
402 : static GEN
403 4900 : tau_edge(GEN A, long i, long j, GEN lambda, long prec)
404 : {
405 4900 : GEN tauij = utoipos(4), Aprim = make_Aprim(A, i, j);
406 4900 : long l = lg(Aprim), k;
407 23800 : for (k = 1; k < l; ++k)
408 : {
409 18900 : GEN xItau = gasinh(gdiv(gatanh(gel(Aprim,k), prec), lambda), prec);
410 18900 : tauij = gmin_shallow(tauij, absr(imag_i(xItau)));
411 : }
412 4900 : return tauij;
413 : }
414 :
415 : static void
416 350 : max_spanning(GEN A, long nedge, GEN lambda, long prec, GEN *ptree, GEN *ptaumin)
417 : {
418 350 : pari_sp av = avma;
419 : GEN real_A, tau_v, tau_c, per, taken, tree, taumin;
420 350 : long z, i, j, k, n = lg(A)-1, m = (n*(n-1))>>1;
421 :
422 350 : tau_v = cgetg(m+1, t_VEC);
423 350 : tau_c = cgetg(m+1, t_VEC);
424 350 : real_A = real_i(A);
425 350 : z = 1;
426 2380 : for (i = 1; i <= n; ++i)
427 6930 : for (j = i+1; j <= n; ++j)
428 : {
429 4900 : gel(tau_v, z) = tau_edge(A, i, j, lambda, prec);
430 4900 : if (gcmp(gel(real_A, i), gel(real_A, j)) > 0)
431 1260 : gel(tau_c, z) = mkvecsmall2(j, i);
432 : else
433 3640 : gel(tau_c, z) = mkvecsmall2(i, j);
434 4900 : z++;
435 : }
436 350 : per = indexsort(tau_v);
437 350 : tau_v = vecpermute(tau_v, per);
438 350 : tau_c = vecpermute(tau_c, per);
439 350 : taken = zero_Flv(n);
440 350 : tree = cgetg(nedge+1, t_VEC);
441 350 : taken[mael(tau_c, m, 1)] = 1;
442 350 : taumin = gel(tau_v, m);
443 1750 : for (k = 1; k <= nedge; ++k)
444 : {
445 1400 : z = m;
446 3955 : while (taken[mael(tau_c, z, 1)]+taken[mael(tau_c, z, 2)] != 1) z--;
447 1400 : gel(tree, k) = gel(tau_c, z);
448 1400 : taumin = gmin_shallow(taumin, gel(tau_v, z));
449 1400 : taken[mael(tau_c, z, 1)] = taken[mael(tau_c, z, 2)] = 1;
450 : }
451 350 : *ptree = tree;
452 350 : *ptaumin = taumin;
453 350 : (void)gc_all(av, 2, ptaumin, ptree);
454 350 : }
455 :
456 : static long
457 707 : hyperellgenus(GEN H)
458 707 : { long d = degpol(H); return ((d+1)>>1)-1; }
459 : static GEN
460 350 : periodmat(GEN P, long prec)
461 : {
462 350 : pari_sp av = avma;
463 350 : GEN A = roots(P, prec), hc, lambda, tree, tau, h, coh1x_homC, IntC, ABtoC;
464 350 : long npoints, g = hyperellgenus(P);
465 :
466 350 : lambda = Pi2n(-1, DEFAULTPREC);
467 350 : max_spanning(A, 2*g, lambda, DEFAULTPREC, &tree, &tau);
468 350 : h = integration_parameters(tau, stoi(prec), lambda, DEFAULTPREC, &npoints);
469 350 : h = rtor(h, prec);
470 350 : lambda = Pi2n(-1,prec);
471 350 : hc = mkvec5(stoi(g), A, tree, gen_0, gen_0);
472 350 : integration_points_thsh(h, npoints, lambda, prec, &gel(hc, 4), &gel(hc, 5));
473 350 : coh1x_homC = periods_spanning(hc, prec);
474 350 : IntC = intersection_spanning(A, tree, prec);
475 350 : ABtoC = ZM_symplectic_reduction(IntC, NULL);
476 350 : return gc_upto(av, gdiv(gmul(coh1x_homC, ABtoC),
477 : gmul2n(gsqrt(leading_coeff(P), prec),-1)));
478 : }
479 :
480 : /* return 4P + Q^2 */
481 : static GEN
482 7 : check_hyperell(GEN PQ)
483 : {
484 : GEN H;
485 7 : if (is_vec_t(typ(PQ)) && lg(PQ)==3)
486 0 : H = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
487 : else
488 7 : H = gmul2n(PQ, 2);
489 7 : return typ(H) == t_POL? H: NULL;
490 : }
491 :
492 : static GEN
493 350 : genus2BSDperiod(GEN C, long prec)
494 : {
495 350 : pari_sp av = avma;
496 : forsubset_t iter;
497 350 : GEN PQ, P, M, v, B = int2n(prec >> 1);
498 : long g;
499 350 : PQ = hyperellminimalmodel(C, NULL, NULL);
500 350 : P = gadd(gmulsg(4, gel(PQ, 1)), gsqr(gel(PQ, 2)));
501 350 : g = hyperellgenus(P);
502 350 : M = real_i(periodmat(P, prec));
503 350 : forsubset_init(&iter, mkvec2s(2*g,g));
504 938 : while ((v = forsubset_next(&iter)))
505 : {
506 938 : GEN Om = vecpermute(M, v), Dm = det(Om);
507 938 : if (expo(Dm) > -(prec>>1))
508 : {
509 350 : GEN r, d, Omr = bestappr(gmul(ginv(Om), M), B);
510 350 : Omr = Q_remove_denom(Omr, &d);
511 350 : r = gmul(Dm, ZM_det_triangular(ZM_hnf(Omr)));
512 350 : if (d) r = gdiv(r, gsqr(d));
513 350 : return gc_upto(av, gabs(r, prec));
514 : }
515 : }
516 0 : pari_err_BUG("genus2BSDperiod");
517 : return NULL; /* LCOV_EXCL_LINE */
518 : }
519 :
520 : GEN
521 357 : hyperellperiods(GEN C, long flag, long prec)
522 : {
523 357 : pari_sp av = avma;
524 : GEN M, H;
525 : long g;
526 357 : if (flag==2) return genus2BSDperiod(C, prec);
527 7 : H = check_hyperell(C);
528 7 : if (!H) pari_err_TYPE("hyperellperiods", C);
529 7 : if (flag<0 || flag>1) pari_err_FLAG("hyperellperiods");
530 7 : g = hyperellgenus(H); if (g < 1) pari_err_DOMAIN("hyperellperiods","genus","=",gen_0,C);
531 0 : M = periodmat(H, prec);
532 0 : if (flag==0) M = gauss(vecslice(M,1,g), vecslice(M,g+1,2*g));
533 0 : return gc_upto(av, M);
534 : }
|
