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 371 : row_swap(GEN M, long i, long j)
34 : {
35 371 : long k, l = lg(M);
36 1883 : for (k = 1; k < l; k++) swap(gcoeff(M,i,k), gcoeff(M,j,k));
37 371 : }
38 :
39 : static void
40 371 : swap_step(GEN P, GEN M, long i, long j)
41 : {
42 371 : if (i == j) return;
43 371 : swap(gel(P,i), gel(P,j));
44 371 : swap(gel(M,i), gel(M,j));
45 371 : row_swap(M, i, j);
46 : }
47 :
48 : /* M <- U(i,j, [u,v; u1,v1])~ * M. In place */
49 : static void
50 658 : row_bezout(GEN M, long i, long j, GEN u, GEN v, GEN u1, GEN v1)
51 : {
52 658 : long k, l = lg(M);
53 3444 : for (k = 1; k < l; k++)
54 : {
55 2786 : GEN a = gcoeff(M,i,k), b = gcoeff(M,j,k);
56 2786 : gcoeff(M,i,k) = addii(mulii(u,a), mulii(v,b));
57 2786 : gcoeff(M,j,k) = addii(mulii(u1,a),mulii(v1,b));
58 : }
59 658 : }
60 :
61 : /* M <- M * U(i,j, [u,v; u1,v1]). In place */
62 : static void
63 1316 : col_bezout(GEN M, long i, long j, GEN u, GEN v, GEN u1, GEN v1)
64 : {
65 1316 : GEN Mi = gel(M,i), Mj = gel(M,j);
66 1316 : gel(M,i) = ZC_lincomb(u, v, Mi, Mj);
67 1316 : gel(M,j) = ZC_lincomb(u1, v1, Mi, Mj);
68 1316 : }
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 658 : bezout_apply(GEN P, GEN M, long i, long j, GEN u, GEN v, GEN u1, GEN v1)
75 : {
76 658 : col_bezout(P, i,j, u,v,u1,v1);
77 658 : row_bezout(M, i,j, u,v,u1,v1);
78 658 : col_bezout(M, i,j, u,v,u1,v1);
79 658 : }
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 658 : transvect_step(GEN P, GEN M, long i, long j, GEN q)
91 658 : { 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 742 : row_pivot(GEN m, long i)
97 : {
98 742 : long j, jx = 0, l = lg(m);
99 742 : GEN x = NULL;
100 2268 : for (j = 1; j < l; j++)
101 : {
102 2268 : GEN z = gcoeff(m,i,j);
103 2268 : if (signe(z))
104 : {
105 742 : 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 364 : ZM_symplectic_reduction(GEN M, GEN *pD)
115 : {
116 364 : long dim, n = lg(M)-1, g = n >> 1; /* n will decrease */
117 364 : GEN P = matid(n), D = zerovec(g);
118 364 : pari_sp av = avma;
119 :
120 364 : M = shallowcopy(M);
121 : /* main loop on symplectic 2-subspace */
122 1106 : for (dim = 0; dim < g; dim++)
123 : {
124 742 : long j, i = 2 * dim + 1;
125 742 : int cleared = 0;
126 : /* lines 0..2d-1 already cleared */
127 742 : 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 742 : if (j != i+1) { swap_step(P, M, j, i+1); j = i+1; }
133 742 : 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 1484 : while (!cleared)
137 : { /* clear row i */
138 : long k;
139 1526 : for (k = j + 1; k <= n; k++)
140 784 : if (signe(gcoeff(M,i,k)))
141 : {
142 294 : GEN r, q = dvmdii(gcoeff(M,i,k), gcoeff(M,i,j), &r);
143 294 : if (r == gen_0)
144 294 : 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 742 : cleared = 1;
149 : /* clear row j */
150 1526 : for (k = j + 1; k <= n; k++)
151 784 : if (signe(gcoeff(M,j,k)))
152 : {
153 364 : GEN r, q = dvmdii(gcoeff(M,j,k), gcoeff(M,i,j), &r);
154 364 : if (r == gen_0)
155 364 : 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 742 : gel(D, dim + 1) = gcoeff(M,i,j);
164 : }
165 364 : END:
166 364 : if (pD) *pD = D;
167 364 : 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 9296 : make_Aprim(GEN A, long ia, long ib)
174 : {
175 9296 : long i, j = 0, lA = lg(A);
176 9296 : GEN a = gel(A, ia), b = gel(A, ib), p = gadd(b, a), m = gsub(b, a);
177 9296 : GEN Aprim = cgetg(lA - 2, t_VEC);
178 64204 : for (i = 1; i < lA; ++i)
179 54908 : if (i != ia && i != ib)
180 36316 : gel(Aprim, ++j) = gdiv(gsub(gmulsg(2, gel(A, i)), p), m);
181 9296 : return Aprim;
182 : }
183 :
184 : static GEN
185 618534 : sqrt_affinereduction(GEN Aprim, GEN z, long prec)
186 : {
187 618534 : pari_sp av = avma;
188 618534 : GEN p = gen_1;
189 618534 : long i, l = lg(Aprim), s = 0;
190 2901528 : for (i = 1; i < l; ++i)
191 : {
192 2282994 : GEN a = gel(Aprim, i);
193 2282994 : if (signe(real_i(a)) > 0) { s++; a = gsub(a, z); } else a = gsub(z, a);
194 2282994 : p = gmul(p, gsqrt(a, prec));
195 : }
196 618534 : return gc_upto(av, gmul(p, powIs(s)));
197 : }
198 :
199 : static long
200 868 : intersection_abbd(GEN A, long ia, long ib, long id, long prec)
201 : {
202 868 : pari_sp av = avma;
203 868 : long k = lg(A)-1;
204 868 : GEN a = gel(A, ia), b = gel(A, ib), d = gel(A, id);
205 868 : GEN fbd = gmul(gpowgs(gsqrt(gsub(d, b), prec), k),
206 : sqrt_affinereduction(make_Aprim(A, ib, id), gen_m1, prec));
207 868 : GEN fab = gmul(gpowgs(gsqrt(gsub(b, a), prec), k),
208 : sqrt_affinereduction(make_Aprim(A, ia, ib), gen_1, prec));
209 868 : return gc_long(av, signe(imag_i(gdiv(fbd, fab))));
210 : }
211 :
212 : static long
213 259 : intersection_abcb(GEN A, long ia, long ib, long ic, long prec)
214 : {
215 259 : pari_sp av = avma;
216 259 : long k = lg(A)-1;
217 259 : GEN a = gel(A, ia), b = gel(A, ib), c = gel(A, ic), fab, fcb;
218 259 : fcb = gmul(gpowgs(gsqrt(gsub(b, c), prec), k),
219 : sqrt_affinereduction(make_Aprim(A, ic, ib), gen_1, prec));
220 259 : fab = gmul(gpowgs(gsqrt(gsub(b, a), prec), k),
221 : sqrt_affinereduction(make_Aprim(A, ia, ib), gen_1, prec));
222 259 : return gc_long(av, signe(imag_i(gdiv(fab, fcb))));
223 : }
224 :
225 : static long
226 182 : intersection_abad(GEN A, long ia, long ib, long id, long prec)
227 : {
228 182 : pari_sp av = avma;
229 182 : long k = lg(A)-1;
230 182 : GEN a = gel(A, ia), b = gel(A, ib), d = gel(A, id), fab, fad;
231 182 : fad = gmul(gpowgs(gsqrt(gsub(d, a), prec), k),
232 : sqrt_affinereduction(make_Aprim(A, ia, id), gen_m1, prec));
233 182 : fab = gmul(gpowgs(gsqrt(gsub(b, a), prec), k),
234 : sqrt_affinereduction(make_Aprim(A, ia, ib), gen_m1, prec));
235 182 : 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 1001 : intersection_inner(GEN A, long ia, long ib, long ic, long id, long prec)
242 : {
243 1001 : pari_sp av = avma;
244 1001 : 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 1001 : GEN bpa = gadd(b, a), bma = gsub(b, a), dpc, dmc;
247 1001 : GEN cprim = gdiv(gsub(gmulsg(2, c), bpa), bma);
248 1001 : GEN dprim = gdiv(gsub(gmulsg(2, d), bpa), bma);
249 1001 : GEN imc = imag_i(cprim), imd = imag_i(dprim);
250 : long k;
251 1001 : if (signe(imc)*signe(imd) == 1) return 0;
252 : /* on the same side */
253 : /* p the intersection */
254 476 : xp = imag_i(gmul(gconj(cprim), gsub(dprim, cprim)));
255 476 : fp = gsub(imd, imc);
256 476 : 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 2310 : intersection(GEN A, long ia, long ib, long ic, long id, long prec)
272 : {
273 2310 : if (ia == ib || ic == id) return 0; /* bad entry */
274 2310 : if (ia == ic && ib == id) return 0; /* self intersection */
275 2310 : if (ia == id && ib == ic) return 0; /* self intersection */
276 2310 : if (ia == ic) return intersection_abad(A, ia, ib, id, prec);
277 2128 : if (ib == ic) return intersection_abbd(A, ia, ib, id, prec);
278 1764 : if (ia == id) return -intersection_abbd(A, ic, id, ib, prec);
279 1260 : if (ib == id) return intersection_abcb(A, ia, ib, ic, prec);
280 1001 : return intersection_inner(A, ia, ib, ic, id, prec);
281 : }
282 :
283 : static GEN
284 364 : intersection_spanning(GEN A, GEN tree, long prec)
285 : {
286 364 : long i, j, n = lg(tree)-1;
287 364 : GEN res = cgetg(n+1, t_MAT);
288 1848 : for (i = 1; i <= n; ++i)
289 1484 : gel(res, i) = cgetg(n+1, t_VECSMALL);
290 1848 : for (i = 1; i <= n; ++i)
291 : {
292 1484 : coeff(res, i, i) = 0;
293 3794 : for (j = i+1; j <= n; ++j)
294 : {
295 2310 : long s = intersection(A, mael(tree, i, 1), mael(tree, i, 2),
296 2310 : mael(tree, j, 1), mael(tree, j, 2), prec);
297 2310 : coeff(res, i, j) = s;
298 2310 : coeff(res, j, i) = -s;
299 : }
300 : }
301 364 : return zm_to_ZM(res);
302 : }
303 :
304 : static GEN
305 1484 : int_periods_affinereduction(GEN C, GEN edge, long prec)
306 : {
307 1484 : pari_sp av = avma;
308 1484 : long g = itos(gel(C, 1)), i1 = edge[1], i2 = edge[2];
309 1484 : GEN A = gel(C, 2), a = gel(A, i1), b = gel(A, i2);
310 1484 : GEN h = gel(C, 4), int_points = gel(C, 5);
311 : GEN F, geom_factor, decprim, Aprim, res;
312 1484 : long i, j, l = lg(int_points);
313 :
314 1484 : if (gcmp(real_i(a), real_i(b)) > 0) pari_err_BUG("hyperellperiods");
315 1484 : decprim = gdiv(gadd(b, a), gsub(b, a));
316 1484 : Aprim = make_Aprim(A, i1, i2);
317 1484 : res = gpowers0(decprim, g-1, ginv(sqrt_affinereduction(Aprim, gen_0, prec)));
318 308700 : for (i = 1; i < l; i++)
319 : {
320 307216 : GEN x = gmael(int_points, i, 1), dx = gmael(int_points, i, 2);
321 307216 : GEN tp = gdiv(dx, sqrt_affinereduction(Aprim, x, prec));
322 307216 : GEN tm = gdiv(dx, sqrt_affinereduction(Aprim, gneg(x), prec));
323 307216 : GEN Tp = gpowers0(gadd(decprim,x), g-1, tp);
324 307216 : GEN Tm = gpowers0(gsub(decprim,x), g-1, tm);
325 933324 : for (j = 1; j <= g; j++)
326 626108 : gel(res,j) = gadd(gel(res,j), gadd(gel(Tp,j), gel(Tm,j)));
327 : }
328 1484 : geom_factor = gdivgs(gsub(b, a), 2);
329 1484 : F = gpowers0(geom_factor, g-1,
330 1484 : gdiv(mulcxI(h), gpowgs(gsqrt(geom_factor, prec), lg(Aprim)-1)));
331 4536 : for (j = 1; j <= g; j++) gel(res,j) = gmul(gel(res,j), gel(F,j));
332 1484 : settyp(res, t_COL); return gc_GEN(av, res);
333 : }
334 :
335 : static GEN
336 364 : periods_spanning(GEN C, long prec)
337 : {
338 364 : GEN tree = gel(C, 3);
339 364 : long k, n = lg(tree)-1;
340 364 : GEN res = cgetg(n+1, t_MAT);
341 1848 : for (k = 1; k <= n; k++)
342 1484 : gel(res, k) = gmul2n(int_periods_affinereduction(C, gel(tree,k), prec), 1);
343 364 : return res;
344 : }
345 :
346 : /* tau, lambda are t_REAL */
347 : static GEN
348 364 : phi_bound(GEN tau, GEN lambda)
349 : {
350 364 : GEN lam2 = sqrr(lambda), costau, sintau, Xtau, Ytau2;
351 364 : mpsincos(tau, &costau, &sintau);
352 364 : Ytau2 = mulrr(lam2, sintau);
353 364 : Xtau = divrr(mulrr(costau, sqrtr(subrr(Ytau2, mulrr(lam2, sqrr(sintau))))),
354 : sintau);
355 364 : return addrr(divur(2, mpcos(sqrtr(Ytau2))), invr(mpsinh(Xtau)));
356 : }
357 :
358 : /* tau and lambda are t_REAL */
359 : static GEN
360 364 : integration_parameters(GEN tau, long bit, GEN lambda, long *pnpoints)
361 : {
362 364 : GEN h = divrr(mulrr(Pi2n(1, DEFAULTPREC), tau),
363 : mplog1p(mulrr(shiftr(phi_bound(tau, lambda), 1),
364 : gexp(utoi(bit), DEFAULTPREC))));
365 364 : GEN t = mpasinh(divrr(addsr(bit, mulur(3,mplog2(DEFAULTPREC))), lambda));
366 364 : *pnpoints = itos(gceil(divrr(t, h))); return h;
367 : }
368 :
369 : static void
370 364 : integration_points_thsh(GEN h, long npoints, GEN lambda, long prec, GEN *ph, GEN *pres)
371 : {
372 364 : pari_sp av = avma;
373 : long k;
374 364 : GEN eh = gexp(h, prec), eh_inv = ginv(eh), ekh = gen_1, ekh_inv = gen_1;
375 364 : GEN res = cgetg(npoints+1, t_VEC);
376 76195 : for (k = 1; k <= npoints; ++k)
377 : {
378 : GEN sh, ch2, esh, esh_inv, chsh2_i, shsh2, thsh;
379 75831 : ekh = gmul(ekh, eh);
380 75831 : ekh_inv = gmul(ekh_inv, eh_inv);
381 75831 : sh = gdivgs(gsub(ekh, ekh_inv), 2);
382 75831 : ch2 = gadd(ekh, ekh_inv);
383 75831 : esh = gexp(gmul(lambda, sh), prec);
384 75831 : esh_inv = ginv(esh);
385 75831 : chsh2_i = ginv(gadd(esh, esh_inv));
386 75831 : shsh2 = gsub(esh, esh_inv);
387 75831 : thsh = gmul(shsh2, chsh2_i);
388 75831 : gel(res, k) = mkvec2(thsh, gmul(ch2, chsh2_i));
389 : }
390 364 : *ph = gmul(h, lambda);
391 364 : *pres = res;
392 364 : (void) gc_all(av, 2, ph, pres);
393 364 : }
394 :
395 : static GEN
396 5194 : tau_edge(GEN A, long i, long j, GEN lambda)
397 : {
398 5194 : GEN tauij = utoipos(4), Aprim = make_Aprim(A, i, j);
399 5194 : long prec = DEFAULTPREC, l = lg(Aprim), k;
400 25564 : for (k = 1; k < l; ++k)
401 : {
402 20370 : GEN xItau = gasinh(gdiv(gatanh(gel(Aprim,k), prec), lambda), prec);
403 20370 : tauij = gmin_shallow(tauij, absr(imag_i(xItau)));
404 : }
405 5194 : return tauij;
406 : }
407 :
408 : static void
409 364 : max_spanning(GEN A, long nedge, GEN lambda, GEN *ptree, GEN *ptaumin)
410 : {
411 364 : pari_sp av = avma;
412 : GEN real_A, tau_v, tau_c, per, taken, tree, taumin;
413 364 : long z, i, j, k, n = lg(A)-1, m = (n*(n-1))>>1;
414 :
415 364 : tau_v = cgetg(m+1, t_VEC);
416 364 : tau_c = cgetg(m+1, t_VEC);
417 364 : real_A = real_i(A);
418 364 : z = 1;
419 2492 : for (i = 1; i <= n; ++i)
420 7322 : for (j = i+1; j <= n; ++j)
421 : {
422 5194 : gel(tau_v, z) = tau_edge(A, i, j, lambda);
423 5194 : if (gcmp(gel(real_A, i), gel(real_A, j)) > 0)
424 1414 : gel(tau_c, z) = mkvecsmall2(j, i);
425 : else
426 3780 : gel(tau_c, z) = mkvecsmall2(i, j);
427 5194 : z++;
428 : }
429 364 : per = indexsort(tau_v);
430 364 : tau_v = vecpermute(tau_v, per);
431 364 : tau_c = vecpermute(tau_c, per);
432 364 : taken = zero_Flv(n);
433 364 : tree = cgetg(nedge+1, t_VEC);
434 364 : taken[mael(tau_c, m, 1)] = 1;
435 364 : taumin = gel(tau_v, m);
436 1848 : for (k = 1; k <= nedge; ++k)
437 : {
438 1484 : z = m;
439 4277 : while (taken[mael(tau_c, z, 1)]+taken[mael(tau_c, z, 2)] != 1) z--;
440 1484 : gel(tree, k) = gel(tau_c, z);
441 1484 : taumin = gmin_shallow(taumin, gel(tau_v, z));
442 1484 : taken[mael(tau_c, z, 1)] = taken[mael(tau_c, z, 2)] = 1;
443 : }
444 364 : *ptree = tree;
445 364 : *ptaumin = taumin;
446 364 : (void)gc_all(av, 2, ptaumin, ptree);
447 364 : }
448 :
449 : static long
450 735 : hyperellgenus(GEN H)
451 735 : { long d = degpol(H); return ((d+1)>>1)-1; }
452 : static GEN
453 364 : periodmat(GEN P, long prec)
454 : {
455 364 : pari_sp av = avma;
456 364 : GEN A = roots(P, prec), hc, lambda, tree, tau, h, coh1x_homC, IntC, ABtoC;
457 364 : long npoints, g = hyperellgenus(P);
458 :
459 364 : lambda = Pi2n(-1, DEFAULTPREC);
460 364 : max_spanning(A, 2*g, lambda, &tree, &tau);
461 364 : h = integration_parameters(tau, prec, lambda, &npoints);
462 364 : h = rtor(h, prec);
463 364 : lambda = Pi2n(-1,prec);
464 364 : hc = mkvec5(stoi(g), A, tree, gen_0, gen_0);
465 364 : integration_points_thsh(h, npoints, lambda, prec, &gel(hc, 4), &gel(hc, 5));
466 364 : coh1x_homC = periods_spanning(hc, prec);
467 364 : IntC = intersection_spanning(A, tree, prec);
468 364 : ABtoC = ZM_symplectic_reduction(IntC, NULL);
469 364 : 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 1484 : RgC_round_0(GEN x, long prec)
487 4536 : { pari_APPLY_type(t_COL, signe(gel(x,i))==0 || expo(gel(x,i))<prec ? gen_0: gel(x,i)) }
488 :
489 : static GEN
490 364 : RgM_round_0(GEN x, long prec)
491 1848 : { pari_APPLY_same(RgC_round_0(gel(x,i), prec)) }
492 :
493 : static GEN
494 364 : genus2BSDperiod(GEN C, long prec)
495 : {
496 364 : pari_sp av = avma;
497 : forsubset_t iter;
498 364 : GEN PQ, P, M, v, B = int2n(prec >> 1);
499 : long g;
500 364 : PQ = hyperellminimalmodel(C, NULL, NULL);
501 364 : P = gadd(gmulsg(4, gel(PQ, 1)), gsqr(gel(PQ, 2)));
502 364 : g = hyperellgenus(P);
503 364 : M = RgM_round_0(real_i(periodmat(P, prec)),-(prec>>1));
504 364 : forsubset_init(&iter, mkvec2s(2*g,g));
505 1106 : while ((v = forsubset_next(&iter)))
506 : {
507 1106 : GEN Om = vecpermute(M, v), Dm = det(Om);
508 1106 : if (signe(Dm) && expo(Dm) > -(prec>>1))
509 : {
510 364 : GEN r, d, Omr = bestappr(gmul(ginv(Om), M), B);
511 364 : Omr = Q_remove_denom(Omr, &d);
512 364 : r = gmul(Dm, ZM_det_triangular(ZM_hnf(Omr)));
513 364 : if (d) r = gdiv(r, gsqr(d));
514 364 : return gc_upto(av, gabs(r, prec));
515 : }
516 : }
517 0 : pari_err_BUG("genus2BSDperiod");
518 : return NULL; /* LCOV_EXCL_LINE */
519 : }
520 :
521 : GEN
522 371 : hyperellperiods(GEN C, long flag, long prec)
523 : {
524 371 : pari_sp av = avma;
525 : GEN M, H;
526 : long g;
527 371 : if (flag==2) return genus2BSDperiod(C, prec);
528 7 : H = check_hyperell(C);
529 7 : if (!H) pari_err_TYPE("hyperellperiods", C);
530 7 : if (flag<0 || flag>1) pari_err_FLAG("hyperellperiods");
531 7 : g = hyperellgenus(H); if (g < 1) pari_err_DOMAIN("hyperellperiods","genus","=",gen_0,C);
532 0 : M = periodmat(H, prec);
533 0 : if (flag==0) M = gauss(vecslice(M,1,g), vecslice(M,g+1,2*g));
534 0 : return gc_upto(av, M);
535 : }
|
