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 594874 : sqrt_affinereduction(GEN Aprim, GEN z, long prec)
186 : {
187 594874 : pari_sp av = avma;
188 594874 : GEN p = gen_1;
189 594874 : long i, l = lg(Aprim), s = 0;
190 2759568 : for (i = 1; i < l; ++i)
191 : {
192 2164694 : GEN a = gel(Aprim, i);
193 2164694 : if (signe(real_i(a)) > 0) { s++; a = gsub(a, z); } else a = gsub(z, a);
194 2164694 : p = gmul(p, gsqrt(a, prec));
195 : }
196 594874 : 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 296940 : for (i = 1; i < l; i++)
319 : {
320 295540 : GEN x = gmael(int_points, i, 1), dx = gmael(int_points, i, 2);
321 295540 : GEN tp = gdiv(dx, sqrt_affinereduction(Aprim, x, prec));
322 295540 : GEN tm = gdiv(dx, sqrt_affinereduction(Aprim, gneg(x), prec));
323 295540 : GEN Tp = gpowers0(gadd(decprim,x), g-1, tp);
324 295540 : GEN Tm = gpowers0(gsub(decprim,x), g-1, tm);
325 886620 : for (j = 1; j <= g; j++)
326 591080 : 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 : /* tau, lambda are t_REAL */
347 : static GEN
348 350 : phi_bound(GEN tau, GEN lambda)
349 : {
350 350 : GEN lam2 = sqrr(lambda), costau, sintau, Xtau, Ytau2;
351 350 : mpsincos(tau, &costau, &sintau);
352 350 : Ytau2 = mulrr(lam2, sintau);
353 350 : Xtau = divrr(mulrr(costau, sqrtr(subrr(Ytau2, mulrr(lam2, sqrr(sintau))))),
354 : sintau);
355 350 : return addrr(divur(2, mpcos(sqrtr(Ytau2))), invr(mpsinh(Xtau)));
356 : }
357 :
358 : /* tau and lambda are t_REAL */
359 : static GEN
360 350 : integration_parameters(GEN tau, long bit, GEN lambda, long *pnpoints)
361 : {
362 350 : GEN h = divrr(mulrr(Pi2n(1, DEFAULTPREC), tau),
363 : mplog1p(mulrr(shiftr(phi_bound(tau, lambda), 1),
364 : gexp(utoi(bit), DEFAULTPREC))));
365 350 : GEN t = mpasinh(divrr(addsr(bit, mulur(3,mplog2(DEFAULTPREC))), lambda));
366 350 : *pnpoints = itos(gceil(divrr(t, h))); return h;
367 : }
368 :
369 : static void
370 350 : integration_points_thsh(GEN h, long npoints, GEN lambda, long prec, GEN *ph, GEN *pres)
371 : {
372 350 : pari_sp av = avma;
373 : long k;
374 350 : GEN eh = gexp(h, prec), eh_inv = ginv(eh), ekh = gen_1, ekh_inv = gen_1;
375 350 : GEN res = cgetg(npoints+1, t_VEC);
376 74235 : for (k = 1; k <= npoints; ++k)
377 : {
378 : GEN sh, ch2, esh, esh_inv, chsh2_i, shsh2, thsh;
379 73885 : ekh = gmul(ekh, eh);
380 73885 : ekh_inv = gmul(ekh_inv, eh_inv);
381 73885 : sh = gdivgs(gsub(ekh, ekh_inv), 2);
382 73885 : ch2 = gadd(ekh, ekh_inv);
383 73885 : esh = gexp(gmul(lambda, sh), prec);
384 73885 : esh_inv = ginv(esh);
385 73885 : chsh2_i = ginv(gadd(esh, esh_inv));
386 73885 : shsh2 = gsub(esh, esh_inv);
387 73885 : thsh = gmul(shsh2, chsh2_i);
388 73885 : gel(res, k) = mkvec2(thsh, gmul(ch2, chsh2_i));
389 : }
390 350 : *ph = gmul(h, lambda);
391 350 : *pres = res;
392 350 : (void) gc_all(av, 2, ph, pres);
393 350 : }
394 :
395 : static GEN
396 4900 : tau_edge(GEN A, long i, long j, GEN lambda)
397 : {
398 4900 : GEN tauij = utoipos(4), Aprim = make_Aprim(A, i, j);
399 4900 : long prec = DEFAULTPREC, l = lg(Aprim), k;
400 23800 : for (k = 1; k < l; ++k)
401 : {
402 18900 : GEN xItau = gasinh(gdiv(gatanh(gel(Aprim,k), prec), lambda), prec);
403 18900 : tauij = gmin_shallow(tauij, absr(imag_i(xItau)));
404 : }
405 4900 : return tauij;
406 : }
407 :
408 : static void
409 350 : max_spanning(GEN A, long nedge, GEN lambda, GEN *ptree, GEN *ptaumin)
410 : {
411 350 : pari_sp av = avma;
412 : GEN real_A, tau_v, tau_c, per, taken, tree, taumin;
413 350 : long z, i, j, k, n = lg(A)-1, m = (n*(n-1))>>1;
414 :
415 350 : tau_v = cgetg(m+1, t_VEC);
416 350 : tau_c = cgetg(m+1, t_VEC);
417 350 : real_A = real_i(A);
418 350 : z = 1;
419 2380 : for (i = 1; i <= n; ++i)
420 6930 : for (j = i+1; j <= n; ++j)
421 : {
422 4900 : gel(tau_v, z) = tau_edge(A, i, j, lambda);
423 4900 : if (gcmp(gel(real_A, i), gel(real_A, j)) > 0)
424 1260 : gel(tau_c, z) = mkvecsmall2(j, i);
425 : else
426 3640 : gel(tau_c, z) = mkvecsmall2(i, j);
427 4900 : z++;
428 : }
429 350 : per = indexsort(tau_v);
430 350 : tau_v = vecpermute(tau_v, per);
431 350 : tau_c = vecpermute(tau_c, per);
432 350 : taken = zero_Flv(n);
433 350 : tree = cgetg(nedge+1, t_VEC);
434 350 : taken[mael(tau_c, m, 1)] = 1;
435 350 : taumin = gel(tau_v, m);
436 1750 : for (k = 1; k <= nedge; ++k)
437 : {
438 1400 : z = m;
439 3955 : while (taken[mael(tau_c, z, 1)]+taken[mael(tau_c, z, 2)] != 1) z--;
440 1400 : gel(tree, k) = gel(tau_c, z);
441 1400 : taumin = gmin_shallow(taumin, gel(tau_v, z));
442 1400 : taken[mael(tau_c, z, 1)] = taken[mael(tau_c, z, 2)] = 1;
443 : }
444 350 : *ptree = tree;
445 350 : *ptaumin = taumin;
446 350 : (void)gc_all(av, 2, ptaumin, ptree);
447 350 : }
448 :
449 : static long
450 707 : hyperellgenus(GEN H)
451 707 : { long d = degpol(H); return ((d+1)>>1)-1; }
452 : static GEN
453 350 : periodmat(GEN P, long prec)
454 : {
455 350 : pari_sp av = avma;
456 350 : GEN A = roots(P, prec), hc, lambda, tree, tau, h, coh1x_homC, IntC, ABtoC;
457 350 : long npoints, g = hyperellgenus(P);
458 :
459 350 : lambda = Pi2n(-1, DEFAULTPREC);
460 350 : max_spanning(A, 2*g, lambda, &tree, &tau);
461 350 : h = integration_parameters(tau, prec, lambda, &npoints);
462 350 : h = rtor(h, prec);
463 350 : lambda = Pi2n(-1,prec);
464 350 : hc = mkvec5(stoi(g), A, tree, gen_0, gen_0);
465 350 : integration_points_thsh(h, npoints, lambda, prec, &gel(hc, 4), &gel(hc, 5));
466 350 : coh1x_homC = periods_spanning(hc, prec);
467 350 : IntC = intersection_spanning(A, tree, prec);
468 350 : ABtoC = ZM_symplectic_reduction(IntC, NULL);
469 350 : return gc_upto(av, gdiv(gmul(coh1x_homC, ABtoC),
470 : gmul2n(gsqrt(leading_coeff(P), prec),-1)));
471 : }
472 :
473 : /* return 4P + Q^2 */
474 : static GEN
475 7 : check_hyperell(GEN PQ)
476 : {
477 : GEN H;
478 7 : if (is_vec_t(typ(PQ)) && lg(PQ)==3)
479 0 : H = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
480 : else
481 7 : H = gmul2n(PQ, 2);
482 7 : return typ(H) == t_POL? H: NULL;
483 : }
484 :
485 : static GEN
486 350 : genus2BSDperiod(GEN C, long prec)
487 : {
488 350 : pari_sp av = avma;
489 : forsubset_t iter;
490 350 : GEN PQ, P, M, v, B = int2n(prec >> 1);
491 : long g;
492 350 : PQ = hyperellminimalmodel(C, NULL, NULL);
493 350 : P = gadd(gmulsg(4, gel(PQ, 1)), gsqr(gel(PQ, 2)));
494 350 : g = hyperellgenus(P);
495 350 : M = real_i(periodmat(P, prec));
496 350 : forsubset_init(&iter, mkvec2s(2*g,g));
497 938 : while ((v = forsubset_next(&iter)))
498 : {
499 938 : GEN Om = vecpermute(M, v), Dm = det(Om);
500 938 : if (expo(Dm) > -(prec>>1))
501 : {
502 350 : GEN r, d, Omr = bestappr(gmul(ginv(Om), M), B);
503 350 : Omr = Q_remove_denom(Omr, &d);
504 350 : r = gmul(Dm, ZM_det_triangular(ZM_hnf(Omr)));
505 350 : if (d) r = gdiv(r, gsqr(d));
506 350 : return gc_upto(av, gabs(r, prec));
507 : }
508 : }
509 0 : pari_err_BUG("genus2BSDperiod");
510 : return NULL; /* LCOV_EXCL_LINE */
511 : }
512 :
513 : GEN
514 357 : hyperellperiods(GEN C, long flag, long prec)
515 : {
516 357 : pari_sp av = avma;
517 : GEN M, H;
518 : long g;
519 357 : if (flag==2) return genus2BSDperiod(C, prec);
520 7 : H = check_hyperell(C);
521 7 : if (!H) pari_err_TYPE("hyperellperiods", C);
522 7 : if (flag<0 || flag>1) pari_err_FLAG("hyperellperiods");
523 7 : g = hyperellgenus(H); if (g < 1) pari_err_DOMAIN("hyperellperiods","genus","=",gen_0,C);
524 0 : M = periodmat(H, prec);
525 0 : if (flag==0) M = gauss(vecslice(M,1,g), vecslice(M,g+1,2*g));
526 0 : return gc_upto(av, M);
527 : }
|
