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; either version 2 of the License, or (at your option) any later
8 version. It is distributed in the hope that it will be useful, but WITHOUT
9 ANY WARRANTY WHATSOEVER.
10 
11 Check the License for details. You should have received a copy of it, along
12 with the package; see the file 'COPYING'. If not, write to the Free Software
13 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
14 #include "pari.h"
15 #include "paripriv.h"
16 
17 /********************************************************************/
18 /**                                                                **/
19 /**                       IGUSA INVARIANTS                         **/
20 /**                       (GP2C-generated)                         **/
21 /**                                                                **/
22 /********************************************************************/
23 /*
24 j2(a0,a1,a2,a3,a4,a5,a6) = (-120*a0*a6+20*a1*a5-8*a2*a4+3*a3^2) / 4;
25 */
26 static GEN
igusaj2(GEN a0,GEN a1,GEN a2,GEN a3,GEN a4,GEN a5,GEN a6)27 igusaj2(GEN a0, GEN a1, GEN a2, GEN a3, GEN a4, GEN a5, GEN a6)
28 {
29   pari_sp av = avma;
30   return gerepileupto(av, gmul2n(gadd(gsub(gadd(gmul(gmulsg(-120, a0), a6), gmul(gmulsg(20, a1), a5)), gmul(gmulsg(8, a2), a4)), gmulsg(3, gsqr(a3))), -2));
31 }
32 
33 /*
34 j4(a0,a1,a2,a3,a4,a5,a6) = (240*(a0*a3*a4*a5+a1*a2*a3*a6)-400*(a0*a2*a5^2+a1^2*a4*a6)-64*(a0*a4^3+a2^3*a6)+16*(a1*a3*a4^2+a2^2*a3*a5)-672*a0*a3^2*a6+240*a1^2*a5^2-112*a1*a2*a4*a5-8*a1*a3^2*a5+16*a2^2*a4^2-16*a2*a3^2*a4+3*a3^4+2640*a0^2*a6^2-880*a0*a1*a5*a6+1312*a0*a2*a4*a6) / 2^7
35 */
36 static GEN
igusaj4(GEN a0,GEN a1,GEN a2,GEN a3,GEN a4,GEN a5,GEN a6)37 igusaj4(GEN a0, GEN a1, GEN a2, GEN a3, GEN a4, GEN a5, GEN a6)
38 {
39   pari_sp av = avma;
40   return gerepileupto(av,
41 gmul2n(gadd(gsub(gadd(gadd(gsub(gadd(gsub(gsub(gadd(gsub(gadd(gsub(gsub(gmulsg(240,
42 gadd(gmul(gmul(gmul(a0, a3), a4), a5), gmul(gmul(gmul(a1, a2), a3), a6))),
43 gmulsg(400, gadd(gmul(gmul(a0, a2), gsqr(a5)), gmul(gmul(gsqr(a1), a4), a6)))),
44 gmulsg(64, gadd(gmul(a0, gpowgs(a4, 3)), gmul(gpowgs(a2, 3), a6)))), gmulsg(16,
45 gadd(gmul(gmul(a1, a3), gsqr(a4)), gmul(gmul(gsqr(a2), a3), a5)))),
46 gmul(gmul(gmulsg(672, a0), gsqr(a3)), a6)), gmul(gmulsg(240, gsqr(a1)),
47 gsqr(a5))), gmul(gmul(gmul(gmulsg(112, a1), a2), a4), a5)), gmul(gmul(gmulsg(8,
48 a1), gsqr(a3)), a5)), gmul(gmulsg(16, gsqr(a2)), gsqr(a4))),
49 gmul(gmul(gmulsg(16, a2), gsqr(a3)), a4)), gmulsg(3, gpowgs(a3, 4))),
50 gmul(gmulsg(2640, gsqr(a0)), gsqr(a6))), gmul(gmul(gmul(gmulsg(880, a0), a1),
51 a5), a6)), gmul(gmul(gmul(gmulsg(1312, a0), a2), a4), a6)), -7));
52 }
53 
54 /*
55 j6(a0,a1,a2,a3,a4,a5,a6) = (1600*(a0^2*a4^2*a5^2+a1^2*a2^2*a6^2)+1600*(a0*a1*a2*a5^3+a1^3*a4*a5*a6)+640*(a0*a1*a3*a4*a5^2+a1^2*a2*a3*a5*a6)-4000*(a0^2*a3*a5^3+a1^3*a3*a6^2)-384*(a0*a1*a4^3*a5+a1*a2^3*a5*a6)-640*(a0*a2^2*a4*a5^2+a1^2*a2*a4^2*a6)+80*(a0*a2*a3^2*a5^2+a1^2*a3^2*a4*a6)+192*(a0*a2*a3*a4^2*a5+a1*a2^2*a3*a4*a6)-48*(a0*a3^3*a4*a5+a1*a2*a3^3*a6)-224*(a1^2*a3*a4^2*a5+a1*a2^2*a3*a5^2)+64*(a1^2*a4^4+a2^4*a5^2)-64*(a1*a2*a3*a4^3+a2^3*a3*a4*a5)+16*(a1*a3^3*a4^2+a2^2*a3^3*a5)-4096*(a0^2*a4^3*a6+a0*a2^3*a6^2)+6400*(a0^2*a2*a5^2*a6+a0*a1^2*a4*a6^2)+10560*(a0^2*a3*a4*a5*a6+a0*a1*a2*a3*a6^2)+2624*(a0*a1*a3*a4^2*a6+a0*a2^2*a3*a5*a6)-4432*a0*a1*a3^2*a5*a6-8*a2*a3^4*a4+a3^6-320*a1^3*a5^3+64*a1^2*a2*a4*a5^2+176*a1^2*a3^2*a5^2+128*a1*a2^2*a4^2*a5+112*a1*a2*a3^2*a4*a5-28*a1*a3^4*a5+16*a2^2*a3^2*a4^2+5120*a0^3*a6^3-2544*a0^2*a3^2*a6^2+312*a0*a3^4*a6-14336*a0^2*a2*a4*a6^2+1024*a0*a2^2*a4^2*a6-2560*a0^2*a1*a5*a6^2-2240*a0*a1^2*a5^2*a6-6528*a0*a1*a2*a4*a5*a6-1568*a0*a2*a3^2*a4*a6) / 2^10
56 */
57 static GEN
igusaj6(GEN a0,GEN a1,GEN a2,GEN a3,GEN a4,GEN a5,GEN a6)58 igusaj6(GEN a0, GEN a1, GEN a2, GEN a3, GEN a4, GEN a5, GEN a6)
59 {
60   pari_sp av = avma;
61   return gerepileupto(av,
62 gmul2n(gsub(gsub(gsub(gsub(gadd(gsub(gadd(gsub(gadd(gadd(gsub(gadd(gadd(gadd(gadd(gsub(gadd(gsub(gsub(gadd(gadd(gadd(gsub(gadd(gsub(gadd(gsub(gsub(gadd(gadd(gsub(gsub(gsub(gadd(gadd(gmulsg(1600,
63 gadd(gmul(gmul(gsqr(a0), gsqr(a4)), gsqr(a5)), gmul(gmul(gsqr(a1), gsqr(a2)),
64 gsqr(a6)))), gmulsg(1600, gadd(gmul(gmul(gmul(a0, a1), a2), gpowgs(a5, 3)),
65 gmul(gmul(gmul(gpowgs(a1, 3), a4), a5), a6)))), gmulsg(640,
66 gadd(gmul(gmul(gmul(gmul(a0, a1), a3), a4), gsqr(a5)),
67 gmul(gmul(gmul(gmul(gsqr(a1), a2), a3), a5), a6)))), gmulsg(4000,
68 gadd(gmul(gmul(gsqr(a0), a3), gpowgs(a5, 3)), gmul(gmul(gpowgs(a1, 3), a3),
69 gsqr(a6))))), gmulsg(384, gadd(gmul(gmul(gmul(a0, a1), gpowgs(a4, 3)), a5),
70 gmul(gmul(gmul(a1, gpowgs(a2, 3)), a5), a6)))), gmulsg(640,
71 gadd(gmul(gmul(gmul(a0, gsqr(a2)), a4), gsqr(a5)), gmul(gmul(gmul(gsqr(a1),
72 a2), gsqr(a4)), a6)))), gmulsg(80, gadd(gmul(gmul(gmul(a0, a2), gsqr(a3)),
73 gsqr(a5)), gmul(gmul(gmul(gsqr(a1), gsqr(a3)), a4), a6)))), gmulsg(192,
74 gadd(gmul(gmul(gmul(gmul(a0, a2), a3), gsqr(a4)), a5), gmul(gmul(gmul(gmul(a1,
75 gsqr(a2)), a3), a4), a6)))), gmulsg(48, gadd(gmul(gmul(gmul(a0, gpowgs(a3, 3)),
76 a4), a5), gmul(gmul(gmul(a1, a2), gpowgs(a3, 3)), a6)))), gmulsg(224,
77 gadd(gmul(gmul(gmul(gsqr(a1), a3), gsqr(a4)), a5), gmul(gmul(gmul(a1,
78 gsqr(a2)), a3), gsqr(a5))))), gmulsg(64, gadd(gmul(gsqr(a1), gpowgs(a4, 4)),
79 gmul(gpowgs(a2, 4), gsqr(a5))))), gmulsg(64, gadd(gmul(gmul(gmul(a1, a2), a3),
80 gpowgs(a4, 3)), gmul(gmul(gmul(gpowgs(a2, 3), a3), a4), a5)))), gmulsg(16,
81 gadd(gmul(gmul(a1, gpowgs(a3, 3)), gsqr(a4)), gmul(gmul(gsqr(a2), gpowgs(a3,
82 3)), a5)))), gmulsg(4096, gadd(gmul(gmul(gsqr(a0), gpowgs(a4, 3)), a6),
83 gmul(gmul(a0, gpowgs(a2, 3)), gsqr(a6))))), gmulsg(6400,
84 gadd(gmul(gmul(gmul(gsqr(a0), a2), gsqr(a5)), a6), gmul(gmul(gmul(a0,
85 gsqr(a1)), a4), gsqr(a6))))), gmulsg(10560, gadd(gmul(gmul(gmul(gmul(gsqr(a0),
86 a3), a4), a5), a6), gmul(gmul(gmul(gmul(a0, a1), a2), a3), gsqr(a6))))),
87 gmulsg(2624, gadd(gmul(gmul(gmul(gmul(a0, a1), a3), gsqr(a4)), a6),
88 gmul(gmul(gmul(gmul(a0, gsqr(a2)), a3), a5), a6)))),
89 gmul(gmul(gmul(gmul(gmulsg(4432, a0), a1), gsqr(a3)), a5), a6)),
90 gmul(gmul(gmulsg(8, a2), gpowgs(a3, 4)), a4)), gpowgs(a3, 6)), gmul(gmulsg(320,
91 gpowgs(a1, 3)), gpowgs(a5, 3))), gmul(gmul(gmul(gmulsg(64, gsqr(a1)), a2), a4),
92 gsqr(a5))), gmul(gmul(gmulsg(176, gsqr(a1)), gsqr(a3)), gsqr(a5))),
93 gmul(gmul(gmul(gmulsg(128, a1), gsqr(a2)), gsqr(a4)), a5)),
94 gmul(gmul(gmul(gmul(gmulsg(112, a1), a2), gsqr(a3)), a4), a5)),
95 gmul(gmul(gmulsg(28, a1), gpowgs(a3, 4)), a5)), gmul(gmul(gmulsg(16, gsqr(a2)),
96 gsqr(a3)), gsqr(a4))), gmul(gmulsg(5120, gpowgs(a0, 3)), gpowgs(a6, 3))),
97 gmul(gmul(gmulsg(2544, gsqr(a0)), gsqr(a3)), gsqr(a6))), gmul(gmul(gmulsg(312,
98 a0), gpowgs(a3, 4)), a6)), gmul(gmul(gmul(gmulsg(14336, gsqr(a0)), a2), a4),
99 gsqr(a6))), gmul(gmul(gmul(gmulsg(1024, a0), gsqr(a2)), gsqr(a4)), a6)),
100 gmul(gmul(gmul(gmulsg(2560, gsqr(a0)), a1), a5), gsqr(a6))),
101 gmul(gmul(gmul(gmulsg(2240, a0), gsqr(a1)), gsqr(a5)), a6)),
102 gmul(gmul(gmul(gmul(gmul(gmulsg(6528, a0), a1), a2), a4), a5), a6)),
103 gmul(gmul(gmul(gmul(gmulsg(1568, a0), a2), gsqr(a3)), a4), a6)), -10));
104 }
105 
106 /********************************************************************/
107 /**                                                                **/
108 /**   A REDUCTION ALGORITHM "A LA TATE" FOR CURVES OF GENUS 2      **/
109 /**                                                                **/
110 /********************************************************************/
111 /* Based on genus2reduction-0.3, http://www.math.u-bordeaux.fr/~liu/G2R/
112  * by Qing Liu <liu@math.u-bordeaux.fr>
113  * and Henri Cohen <cohen@math.u-bordeaux.fr>
114 
115  * Qing Liu: Modeles minimaux des courbes de genre deux
116  * J. fuer die Reine und Angew. Math., 453 (1994), 137-164.
117  * http://www.math.u-bordeaux.fr/~liu/articles/modregE.ps */
118 
119 /* some auxiliary polynomials, gp2c-generated */
120 
121 /*
122 apol2(a0,a1,a2) = -5*a1^2+12*a0*a2;
123 */
124 static GEN
apol2(GEN a0,GEN a1,GEN a2)125 apol2(GEN a0, GEN a1, GEN a2)
126 {
127   return gadd(gmulsg(-5, gsqr(a1)), gmul(gmulsg(12, a0), a2));
128 }
129 
130 /*
131 apol3(a0,a1,a2,a3) = 5*a1^3+9*a0*(-2*a1*a2+3*a0*a3);
132 */
133 static GEN
apol3(GEN a0,GEN a1,GEN a2,GEN a3)134 apol3(GEN a0, GEN a1, GEN a2, GEN a3)
135 {
136   return gadd(gmulsg(5, gpowgs(a1, 3)), gmul(gmulsg(9, a0), gadd(gmul(gmulsg(-2, a1), a2), gmul(gmulsg(3, a0), a3))));
137 }
138 
139 /*
140 apol5(a0,a1,a2,a3,a4,a5) = a1^5+3*a0*(-2*a1^3*a2+9*a0*a1^2*a3-36*a0^2*a1*a4+108*a0^3*a5);
141 */
142 static GEN
apol5(GEN a0,GEN a1,GEN a2,GEN a3,GEN a4,GEN a5)143 apol5(GEN a0, GEN a1, GEN a2, GEN a3, GEN a4, GEN a5)
144 {
145   return gadd(gpowgs(a1, 5), gmul(gmulsg(3, a0), gadd(gsub(gadd(gmul(gmulsg(-2, gpowgs(a1, 3)), a2), gmul(gmul(gmulsg(9, a0), gsqr(a1)), a3)), gmul(gmul(gmulsg(36, gsqr(a0)), a1), a4)), gmul(gmulsg(108, gpowgs(a0, 3)), a5))));
146 }
147 
148 /*
149 bpol2(a0,a1,a2,a3,a4) = 2*a2^2-5*a1*a3+10*a0*a4;
150 */
151 static GEN
bpol2(GEN a0,GEN a1,GEN a2,GEN a3,GEN a4)152 bpol2(GEN a0, GEN a1, GEN a2, GEN a3, GEN a4)
153 {
154   return gadd(gsub(gmulsg(2, gsqr(a2)), gmul(gmulsg(5, a1), a3)), gmul(gmulsg(10, a0), a4));
155 }
156 
157 static const long VERYBIG = (1L<<20);
158 static long
myval(GEN x,GEN p)159 myval(GEN x, GEN p) { return signe(x)? Z_pval(x,p): VERYBIG; }
160 static long
my3val(GEN x)161 my3val(GEN x) { return signe(x)? Z_lval(x,3): VERYBIG; }
162 /* b in Z[i], return v_3(b) */
163 static long
myval_zi(GEN b)164 myval_zi(GEN b) { return minss(my3val(real_i(b)), my3val(imag_i(b))); }
165 /* b in Z[i, Y]/(Y^2-3), return v_Y(b) */
166 static long
myval_zi2(GEN b)167 myval_zi2(GEN b)
168 {
169   long v0, v1;
170   b = lift_shallow(b);
171   v0 = myval_zi(RgX_coeff(b,0));
172   v1 = myval_zi(RgX_coeff(b,1));
173   return minss(2*v0, 2*v1+1);
174 }
175 
176 /* min(a,b,c) */
177 static long
min3(long a,long b,long c)178 min3(long a, long b, long c)
179 {
180   long m = a;
181   if (b < m) m = b;
182   if (c < m) m = c;
183   return m;
184 }
185 
186 /* Vector of p-adic factors (over Q_p) to accuracy r of pol. */
187 static GEN
padicfactors(GEN pol,GEN p,long r)188 padicfactors(GEN pol, GEN p, long r) { return gel(factorpadic(pol,p,r),1); }
189 
190 /* x(1/t)*t^6, deg x <= 6 */
191 static GEN
RgX_recip6(GEN x)192 RgX_recip6(GEN x)
193 {
194   long lx = lg(x), i, j;
195   GEN y = cgetg(9, t_POL);
196   y[1] = x[1];
197   for (i=8,j=2; j < lx; i--,j++) gel(y,i) = gel(x,j);
198   for (       ; j <  9; i--,j++) gel(y,i) = gen_0;
199   return normalizepol_lg(y, 9);
200 }
201 /* extract coefficients of a polynomial a0 X^6 + ... + a6, of degree <= 6 */
202 static void
RgX_to_06(GEN q,GEN * a0,GEN * a1,GEN * a2,GEN * a3,GEN * a4,GEN * a5,GEN * a6)203 RgX_to_06(GEN q, GEN *a0, GEN *a1, GEN *a2, GEN *a3, GEN *a4, GEN *a5, GEN *a6)
204 {
205   *a0 = gen_0;
206   *a1 = gen_0;
207   *a2 = gen_0;
208   *a3 = gen_0;
209   *a4 = gen_0;
210   *a5 = gen_0;
211   *a6 = gen_0;
212   switch(degpol(q))
213   {
214     case 6: *a0 = gel(q,8); /*fall through*/
215     case 5: *a1 = gel(q,7); /*fall through*/
216     case 4: *a2 = gel(q,6); /*fall through*/
217     case 3: *a3 = gel(q,5); /*fall through*/
218     case 2: *a4 = gel(q,4); /*fall through*/
219     case 1: *a5 = gel(q,3); /*fall through*/
220     case 0: *a6 = gel(q,2); /*fall through*/
221   }
222 }
223 /* extract coefficients a0,...a3 of a polynomial a0 X^6 + ... + a6 */
224 static void
RgX_to_03(GEN q,GEN * a0,GEN * a1,GEN * a2,GEN * a3)225 RgX_to_03(GEN q, GEN *a0, GEN *a1, GEN *a2, GEN *a3)
226 {
227   *a0 = gen_0;
228   *a1 = gen_0;
229   *a2 = gen_0;
230   *a3 = gen_0;
231   switch(degpol(q))
232   {
233     case 6: *a0 = gel(q,8); /*fall through*/
234     case 5: *a1 = gel(q,7); /*fall through*/
235     case 4: *a2 = gel(q,6); /*fall through*/
236     case 3: *a3 = gel(q,5); /*fall through*/
237   }
238 }
239 
240 /* deg(H mod p) = 3, return v_p( disc(correspondig p-adic factor) ) */
241 static long
discpart(GEN H,GEN p,long prec)242 discpart(GEN H, GEN p, long prec)
243 {
244   GEN list, prod, dis;
245   long i, j;
246 
247   if (degpol(FpX_red(H,p)) != 3)
248     pari_err_BUG("discpart [must not reach]"); /* LCOV_EXCL_LINE */
249   list = padicfactors(H,p,prec);
250   prod = pol_1(varn(H));
251   for(i = 1; i < lg(list); i++)
252   {
253     GEN t = gel(list,i);
254     for(j = 3; j < lg(t); j++) /* include if nonconstant mod p */
255       if (!valp(gel(t,j))) { prod = RgX_mul(prod,t); break; }
256   }
257   if (degpol(prod) != 3) pari_err_BUG("discpart [prod degree]");
258   dis = RgX_disc(prod);
259   return gequal0(dis)? prec+1: valp(dis);
260 }
261 
262 /* B = b0 X^6 + ... + b6 a ZX, 0 <= j <= 3.
263  * Let theta_j(H) := min { v_p(b_i) / (i - j), j < i <= 6 } >= 0.
264  * Return 60 theta \in Z */
265 static long
theta_j(GEN B,GEN p,long j)266 theta_j(GEN B, GEN p, long j)
267 {
268   long i, t = 60*myval(RgX_coeff(B,5-j), p);
269   for(i = 2+j; i <= 6; i++)
270     t = minss(t, myval(RgX_coeff(B,6-i), p) * (60 / (i-j)));
271   return t;
272 }
273 /* compute 6 * theta_3 for B in Z[i][X], p = 3 */
274 static long
theta_3_zi(GEN B)275 theta_3_zi(GEN B)
276 {
277   long v2 = myval_zi(RgX_coeff(B,2));
278   long v1 = myval_zi(RgX_coeff(B,1));
279   long v0 = myval_zi(RgX_coeff(B,0));
280   return min3(6*v2, 3*v1, 2*v0);
281 }
282 /* compute 6 * theta_3 for B in (Z[i,Y]/(Y^2-3))[X], p = 3 */
283 static long
theta_3_zi2(GEN B)284 theta_3_zi2(GEN B)
285 {
286   long v2 = myval_zi2(RgX_coeff(B,2));
287   long v1 = myval_zi2(RgX_coeff(B,1));
288   long v0 = myval_zi2(RgX_coeff(B,0));
289   return min3(6*v2, 3*v1, 2*v0);
290 }
291 
292 /* Set maxord to the maximal multiplicity of a factor. If there is at least
293  * a triple root (=> maxord >= 3) return it, else return NULL */
294 static GEN
factmz(GEN Q,GEN p,long * maxord)295 factmz(GEN Q, GEN p, long *maxord)
296 {
297   GEN z = FpX_factor_squarefree(Q, p);
298   long m = lg(z)-1; /* maximal multiplicity */
299   *maxord = m;
300   return (m >= 3)? FpX_oneroot(gel(z,m), p): NULL;
301 }
302 
303 /* H integral ZX of degree 5 or 6, p > 2. Modify until
304  *   y^2 = p^alpha H is minimal over Z_p, alpha = 0,1
305  * Return [H,lambda,60*theta,alpha,quad,beta], where
306  *   - quad = 1 if H has a root of order 3 in F_p^2 \ F_p, 0 otherwise
307  *   - 0 <= lambda <= 3, index of a coefficient with valuation 0
308  *   - theta = theta_j(H(x + r), p, lambda), 60*theta in Z, where r is a root
309  *   of H mod p
310  *   - beta >= -1 s.t. H = p^n H0(r + p^beta * X) for some n, r in Z, where
311  *   H0 is the initial H or polrecip(H) */
312 static GEN
polymini(GEN H,GEN p)313 polymini(GEN H, GEN p)
314 {
315   GEN a0, a1, a2, a3, Hp, rac;
316   long t60, alpha, lambda, maxord, quad = 0, beta = 0;
317 
318   alpha = ZX_pvalrem(H, p, &H) & 1;
319   RgX_to_03(H, &a0,&a1,&a2,&a3);
320   if (dvdii(a0,p) && dvdii(a1,p) && dvdii(a2,p) && dvdii(a3,p))
321   {
322     H = RgX_recip6(H);
323     RgX_to_03(H, &a0,&a1,&a2,&a3);
324   }
325   if (!dvdii(a3,p)) lambda = 3;
326   else if (!dvdii(a2,p)) lambda = 2;
327   else if (!dvdii(a1,p)) lambda = 1;
328   else lambda = 0;
329 
330   for(;;)
331   { /* lambda <= 3, t60 = 60*theta */
332     long e;
333     t60 = theta_j(H,p,lambda); e = t60 / 60;
334     if (e)
335     {
336       GEN pe = powiu(p,e);
337       /* H <- H(p^e X) / p^(e(6-lambda)) */
338       H = ZX_Z_divexact(ZX_unscale_div(H,pe), powiu(pe,5-lambda));
339       alpha = (alpha + lambda*e)&1;
340       beta += e;
341       t60 -= 60*e;
342     }
343     /* 0 <= t < 60 */
344     Hp = FpX_red(H, p);
345     if (t60) break;
346 
347     rac = factmz(Hp,p, &maxord);
348     if (maxord <= 2)
349     {
350       if (degpol(Hp) <= 3) break;
351       goto end;
352     }
353     else
354     { /* maxord >= 3 */
355       if (!rac) { quad = 1; goto end; }
356       if (signe(rac)) H = ZX_translate(H, rac);
357       lambda = 6-maxord;
358     }
359   }
360 
361   if (lambda <= 2)
362   {
363     if (myval(RgX_coeff(H,2),p) > 1-alpha &&
364         myval(RgX_coeff(H,1),p) > 2-alpha &&
365         myval(RgX_coeff(H,0),p) > 3-alpha)
366     {
367       H = ZX_unscale(H, p);
368       if (alpha) H = ZX_Z_mul(H, p);
369       return polymini(H, p);
370     }
371   }
372   else if (lambda == 3 && alpha == 1)
373   {
374     if (degpol(Hp) == 3)
375     {
376       if (myval(RgX_coeff(H,6),p) >= 3 &&
377           myval(RgX_coeff(H,5),p) >= 2)
378       { /* too close to root [Kodaira symbol for y^2 = p^alpha*H not
379            implemented when alpha = 1]: go back one step */
380         H = ZX_rescale(H, p); /* H(x/p)p^(deg H) */
381         H = ZX_Z_divexact(H, powiu(p, degpol(H)-3)); /* H(x/p)p^3 */
382         t60 += 60; alpha = 0; beta--;
383       }
384     }
385     else if (degpol(Hp) == 6 && t60)
386     {
387       rac = factmz(RgX_mulXn(Hp, -3), p, &maxord);
388       if (maxord == 3)
389       {
390         GEN T = ZX_unscale(ZX_translate(H,rac),p); /* H(rac + px) */
391         if (ZX_pval(T,p)>= 3)
392         {
393           H = RgX_Rg_div(T, powiu(p,3));
394           alpha = 0; beta++;
395           t60 = theta_j(H,p,3);
396           if (!t60)
397           {
398             Hp = FpX_red(H, p);
399             if (!signe(FpX_disc(Hp,p)))
400             {
401               rac = FpX_oneroot(Hp, p);
402               t60 = theta_j(ZX_translate(H,rac),p,3);
403             }
404           }
405         }
406       }
407     }
408   }
409 end:
410   return mkvec2(H, mkvecsmall5(lambda,t60,alpha,quad,beta));
411 }
412 
413 /* a in Q[i], return a^3 mod 3 */
414 static GEN
zi_pow3mod(GEN a)415 zi_pow3mod(GEN a)
416 {
417   GEN x, y;
418   if (typ(a) != t_COMPLEX) return gmodgs(a,3);
419   x = gmodgs(gel(a,1), 3);
420   y = gmodgs(gel(a,2), 3);
421   return mkcomplex(x, negi(y));
422 }
423 static GEN
polymini_zi(GEN pol)424 polymini_zi(GEN pol) /* polynome minimal dans Z[i] */
425 {
426   GEN polh, rac, a0, a1, a2, a3, a4, a5, a6, p = utoipos(3);
427   long alpha, beta = 0, t6;
428 
429   alpha = ZX_pval(pol,p) & 1;
430   polh = alpha? RgX_Rg_div(pol, p): pol;
431   rac = mkcomplex(Fp_div(RgX_coeff(polh,3), RgX_coeff(polh,6), p), gen_1);
432   for(;;)
433   {
434     long e;
435     polh = RgX_translate(polh, rac);
436     t6 = theta_3_zi(polh); e = t6 / 6;
437     if (e)
438     {
439       GEN pe = powiu(p,e);
440       polh = RgX_Rg_div(RgX_unscale(polh,pe), powiu(pe,3));
441       alpha = (alpha+e)&1;
442       t6 -= e * 6; beta += e;
443     }
444     RgX_to_06(polh, &a0,&a1,&a2,&a3,&a4,&a5,&a6);
445     if (t6 || !myval_zi(a4) || !myval_zi(a5)) break;
446     rac = zi_pow3mod(gdiv(a6, gneg(a3)));
447   }
448   if (alpha && myval_zi(a0) >= 3 && myval_zi(a1) >= 2 && myval_zi(a2) >= 1)
449   {
450     t6 += 6; beta--; alpha = 0;
451   }
452   if (alpha && beta >= 1) pari_err_BUG("quadratic");
453   return mkvecsmall3(t6, alpha, beta);
454 }
455 
456 /* pol is a ZX, minimal polynomial over Z_3[i,Y]/(Y^2-3) */
457 static GEN
polymini_zi2(GEN pol)458 polymini_zi2(GEN pol)
459 {
460   long alpha, beta, t6;
461   GEN a0, a1, a2, a3, a4, a5, a6;
462   GEN polh, rac, y = pol_x(fetch_var()), p = utoipos(3);
463 
464   if (ZX_pval(pol,p)) pari_err_BUG("polymini_zi2 [polynomial not minimal]");
465   y = mkpolmod(y, gsubgs(gsqr(y), 3)); /* mod(y,y^2-3) */
466   polh = gdivgs(RgX_unscale(pol, y),27); /* H(y*x) / 27 */
467   if (myval_zi2(RgX_coeff(polh,4)) <= 0 ||
468       myval_zi2(RgX_coeff(polh,2)) <= 0)
469   {
470     (void)delete_var();
471     return mkvecsmall2(0,0);
472   }
473 
474   if (myval_zi2(gsub(RgX_coeff(polh,6), RgX_coeff(polh,0))) > 0)
475     rac = gen_I();
476   else
477     rac = gen_1;
478   alpha = 0;
479   beta  = 0;
480   for(;;)
481   {
482     long e;
483     polh = RgX_translate(polh, rac);
484     t6 = theta_3_zi2(polh); e = t6 / 6;
485     if (e)
486     {
487       GEN pent = gpowgs(y, e);
488       polh = RgX_Rg_div(RgX_unscale(polh, pent), gpowgs(pent,3));
489       alpha = (alpha+e)&1;
490       t6 -= 6*e; beta += e;
491     }
492     RgX_to_06(polh, &a0,&a1,&a2,&a3,&a4,&a5,&a6);
493     if (t6 || !myval_zi2(a4) || !myval_zi2(a5)) break;
494     a3 = liftpol_shallow(a3); if (typ(a3)==t_POL) a3 = RgX_coeff(a3,0);
495     a6 = liftpol_shallow(a6); if (typ(a6)==t_POL) a6 = RgX_coeff(a6,0);
496     rac = zi_pow3mod(gdiv(a6,gneg(a3)));
497   }
498   if (alpha)
499   {
500     if (myval_zi2(a0) < 3 || myval_zi2(a1) < 2 || myval_zi2(a2) < 1)
501       pari_err_BUG("polymini_zi2 [alpha]");
502     t6 += 6; beta--;
503   }
504   (void)delete_var();
505   if (odd(beta)) pari_err_BUG("quartic [type over Z[i] must be [K-K-(2*m)]]");
506   return mkvecsmall2(t6, beta);
507 }
508 
509 struct igusa {
510   GEN j2, i4, j4, j6, j8, j10, i12;
511   GEN a0, A2, A3, A5, B2;
512 };
513 struct igusa_p {
514   long eps, tt, r1, r2, tame;
515   GEN p, stable, val, neron;
516   const char *type;
517 };
518 
519 /* initialize Ip */
520 static void
stable_reduction(struct igusa * I,struct igusa_p * Ip,GEN p)521 stable_reduction(struct igusa *I, struct igusa_p *Ip, GEN p)
522 {
523   static const long d[9] = { 0,60,30,30,20,15,12,10 }; /* 120 / deg(X) */
524   GEN j2 = I->j2, i4 = I->i4, j4 = I->j4, j6 = I->j6, j8 = I->j8;
525   GEN val, J, v, Ieps, j10 = I->j10, i12 = I->i12;
526   long s, r1, r2, r3, r4, i, eps;
527 
528   Ip->tame = 0;
529   Ip->neron = NULL;
530   Ip->type = NULL;
531   Ip->p = p;
532   Ip->val = val = cgetg(9, t_VECSMALL);
533   val[1] = myval(j2,p);
534   val[2] = myval(j4,p);
535   val[3] = myval(i4,p);
536   val[4] = myval(j6,p);
537   val[5] = myval(j8,p);
538   val[6] = myval(j10,p);
539   val[7] = myval(i12,p);
540   switch(itos_or_0(p))
541   {
542     case 2:  eps = 4; val[8] = val[5]; Ieps = j8; break;
543     case 3:  eps = 3; val[8] = val[4]; Ieps = j6; break;
544     default: eps = 1; val[8] = val[1]; Ieps = gdivgs(j2,12); break;
545   }
546 
547   v = cgetg(8,t_VECSMALL);
548   for(i = 1; i <= 7; i++) v[i] = val[i] * d[i];
549   s = vecsmall_min(v);
550   Ip->eps  = eps;
551 
552   r1 = 3*eps*val[3];
553   r3 = eps*val[6] + val[8];
554   r2 = eps*val[7];
555   r4 = min3(r1, r2, r3);
556 
557   /* s = max(v_p(X) / deg(X)) */
558   J = cgetg(1, t_VEC);
559   if (s == v[6])
560     Ip->tt = 1;
561   else if (s == v[7])
562   {
563     J = mkvec( Fp_to_mod(gmod(gdiv(gpowgs(i4,3),i12), p), p) );
564     Ip->tt = 2;
565   }
566   else if (s == v[3])
567     Ip->tt = (val[2] == val[3] || 2*val[4] == 3*val[3])? 3: 4;
568   else if (r3 == r4)
569   {
570     GEN a,b, P, sj, pj, t = gmul(gpowgs(j10,eps),Ieps);
571     sj = gaddsg(1728, gdiv(gpowgs(i12,eps), t));
572     pj = gdiv(gpowgs(i4,3*eps), t);
573     a = gmod(sj, p);
574     b = gmod(pj, p);
575     P = mkpoln(3, gen_1, Fp_neg(a,p), b, 0); /* X^2 - SX + P: roots j1,j2 */
576     J = FpX_roots(P, p);
577     switch(lg(J)-1)
578     {
579       case 0:
580         P = FpX_to_mod(P, p);
581         a = FpX_to_mod(pol_x(0), p);
582         b = FpX_to_mod(deg1pol_shallow(b, gen_m1,0), p);
583         J = mkvec2(mkpolmod(a,P), mkpolmod(b,P)); break;
584       case 1:
585         a = Fp_to_mod(gel(J,1), p);
586         J = mkvec2(a, a); break;
587       case 2:
588         settyp(J, t_VEC);
589         J = FpV_to_mod(J, p); break;
590     }
591     Ip->tt = 5;
592   }
593   else if (r2 == r4)
594   {
595     J = mkvec( Fp_to_mod(gmod(gdiv(gpowgs(i4,3),i12), p), p) );
596     Ip->tt = 6;
597   }
598   else
599     Ip->tt = 7; /* r1 == r4 */
600   Ip->stable = mkvec2(stoi(Ip->tt), J);
601 }
602 
603 struct red {
604   const char *t, *pages;
605   double tnum;
606   GEN g;
607 };
608 
609 /* destroy v */
610 static GEN
zv_snf(GEN v)611 zv_snf(GEN v)
612 {
613   long i, l = lg(v);
614   for (i = 1; i < l; i++)
615   {
616     long j, a = v[i];
617     for (j = i+1; j < l; j++)
618     {
619       long b = v[j], d = ugcd(a,b);
620       v[i] = a = a*(b/d);
621       v[j] = d;
622     }
623   }
624   for (i = l-1; i > 0; i--)
625     if (v[i] != 1) { setlg(v, i+1); break; }
626   return zv_to_ZV(v);
627 }
628 
629 static GEN
cyclic(long n)630 cyclic(long n)
631 { return (n <= 1)? cgetg(1, t_VECSMALL): mkvecsmall(n); }
632 static GEN
dicyclic(long a,long b)633 dicyclic(long a, long b)
634 {
635   long d;
636   if (!a) a = 1;
637   if (!b) b = 1;
638   if (a < b) lswap(a,b);
639   d = ugcd(a,b);
640   if (d == 1) return cyclic(a*b);
641   return mkvecsmall2(a*b/d, d);
642 }
643 /* Z/2xZ/2, resp Z/4 for n even, resp. odd */
644 static GEN
groupH(long n)645 groupH(long n) { return odd(n)? cyclic(4): dicyclic(2,2); }
646 
647 static long
get_red(struct red * S,struct igusa_p * Ip,GEN polh,GEN p,long alpha,long r)648 get_red(struct red *S, struct igusa_p *Ip, GEN polh, GEN p, long alpha, long r)
649 {
650   GEN val = Ip->val;
651   long indice;
652   switch(r)
653   {
654     case 0:
655       indice = FpX_is_squarefree(FpX_red(polh,p), p)
656                ? 0
657                : val[6] - val[7] + val[8]/Ip->eps;
658       S->t = stack_sprintf("I{%ld}", indice);
659       S->tnum = 1;
660       S->pages = "159-177";
661       S->g = cyclic(indice);
662       return indice ? indice: 1;
663     case 6:
664       if (alpha == 0) /* H(px) /p^3 */
665         polh = ZX_Z_divexact(ZX_unscale_div(polh,p), sqri(p));
666       indice = FpX_is_squarefree(FpX_red(polh,p), p)
667                ? 0
668                : val[6] - val[7] + val[8]/Ip->eps;
669       S->t = stack_sprintf("I*{%ld}", indice);
670       S->tnum = 1.5;
671       S->pages = "159-177";
672       S->g = groupH(indice);
673       return indice + 5;
674     case 3:
675       S->t = "III";
676       S->tnum = 3;
677       S->pages = "161-177";
678       S->g = cyclic(2);
679       return 2;
680     case 9:
681       S->t = "III*";
682       S->tnum = 3.5;
683       S->pages = "162-177";
684       S->g = cyclic(2);
685       return 8;
686     case 2:
687       S->t = "II";
688       S->tnum = 2;
689       S->pages = "159-174";
690       S->g = cyclic(1);
691       return 1;
692     case 8:
693       S->t = "IV*";
694       S->tnum = 4.5;
695       S->pages = "160-175";
696       S->g = cyclic(3);
697       return 7;
698     case 4:
699       S->t = "IV";
700       S->tnum = 4;
701       S->pages = "160-174";
702       S->g = cyclic(3);
703       return 3;
704     case 10:
705       S->t = "II*";
706       S->tnum = 2.5;
707       S->pages = "160-174";
708       S->g = cyclic(1);
709       return 9;
710     default: pari_err_BUG("get_red [type]");
711       S->t = "";
712       S->tnum = 0;
713       S->pages = ""; /* gcc -Wall */
714       S->g = NULL;
715       return -1; /*LCOV_EXCL_LINE*/
716   }
717 }
718 
719 /* reduce a/b; assume b > 0 */
720 static void
ssQ_red(long a,long b,long * n,long * d)721 ssQ_red(long a, long b, long *n, long *d)
722 {
723   long g = ugcd(labs(a), b);
724   if (g > 1) { a /= g; b /= g; }
725   *n = a; *d = b;
726 }
727 /* denom(a/b); assume b > 0 */
728 static long
ssQ_denom(long a,long b)729 ssQ_denom(long a, long b)
730 {
731   long g = ugcd(labs(a), b);
732   return g == 1? b: b / g;
733 }
734 /* n = lcm(d, denom(a/b)); r = (a/b * n mod n); assume b > 0 and d > 0 */
735 static void
get_nr(long d,long a,long b,long * n,long * r)736 get_nr(long d, long a, long b, long *n, long *r)
737 {
738   long c, A, B;
739   ssQ_red(a, b, &A,&B);
740   c = d / ugcd(d, B);
741   *n = B * c;
742   *r = umodsu(A * c, *n);
743 }
744 /* n = lcm(denom(a/b), denom(c/d)); r = (a/b * n mod n); q = (c/d * n mod n);
745  * assume b > 0 and d > 0 */
746 static void
get_nrq(long a,long b,long c,long d,long * n,long * r,long * q)747 get_nrq(long a, long b, long c, long d, long *n, long *r, long *q)
748 {
749   long g, A, B, C, D;
750   ssQ_red(a, b, &A,&B);
751   ssQ_red(c, d, &C,&D);
752   g = ugcd(B,D);
753   *n = B * (D/g);
754   *r = umodsu(A * (D/g), *n);
755   *q = umodsu(C * (B/g), *n);
756 }
757 
758 /* Ip->tt = 1 */
759 static long
tame_1(struct igusa * I,struct igusa_p * Ip)760 tame_1(struct igusa *I, struct igusa_p *Ip)
761 {
762   GEN p = Ip->p, val = Ip->val;
763   long condp = -1, va0, va5, r, n;
764   va0 = myval(I->a0,p);
765   va5 = myval(I->A5,p);
766   if (!gequal0(I->A5) && 20*va0+val[6] > 6*va5)
767     get_nr(ssQ_denom(5*val[6]-6*va5, 40), val[6]-2*va5, 20, &n,&r);
768   else
769     get_nr(ssQ_denom(5*va0-val[6], 10), 10*va0-val[6], 30, &n,&r);
770   switch(n)
771   {
772     case 1:
773       condp = 0;
774       Ip->type = "[I{0-0-0}] page 155";
775       Ip->neron = cyclic(1); break;
776     case 2:
777       switch(r)
778       {
779         case 0:
780           condp = 4;
781           Ip->type = "[I*{0-0-0}] page 155";
782           Ip->neron = mkvecsmall4(2,2,2,2); break;
783         case 1:
784           condp = 2;
785           Ip->type = "[II] page 155";
786           Ip->neron = cyclic(1); break;
787         default: pari_err_BUG("tame_1 [bug1]");
788       }
789       break;
790     case 4:
791       condp = 4;
792       Ip->type = "[VI] page 156";
793       Ip->neron = dicyclic(2,2); break;
794     default: pari_err_BUG("tame_1 [bug8]");
795   }
796   return condp;
797 }
798 
799 /* (4.2) */
800 static long
tame_234_init(struct igusa * I,struct igusa_p * Ip,long * n,long * q,long * r)801 tame_234_init(struct igusa *I, struct igusa_p *Ip, long *n, long *q, long *r)
802 {
803   long va0, va5, vb2, v12 = -1, flc = 1;
804   GEN p = Ip->p;
805   switch(Ip->tt)
806   {
807     case 2: v12 = myval(I->i12,  Ip->p); break;
808     case 3: v12 = 3*myval(I->i4, Ip->p); break;
809     case 4: v12 = 6*myval(I->j2, Ip->p); break;
810   }
811   va0 = myval(I->a0,p);
812   va5 = myval(I->A5,p);
813   vb2 = myval(I->B2,p);
814   if (9*vb2 >= 6*va0+v12 && 36*va5 >= 120*va0+5*v12)
815   {
816     get_nrq(12*va0-v12,36, 6*va0-v12,12, n, r, q);
817   }
818   else if (120*va0+5*v12 > 36*va5 && 60*vb2 >= 12*va5+5*v12)
819   {
820     ssQ_red(36*va5-25*v12,240, q,n);
821     *r = umodsu(-2* *q, *n);
822   }
823   else /* 6*va0+v12 > 9*vb2 && 12*va5+5*v12 > 60*vb2 */
824   {
825     get_nrq(v12-6*vb2,12, v12-9*vb2,12, n,r,q);
826     flc = 0;
827   }
828   return flc;
829 }
830 
831 /* Ip->tt = 2 */
832 static long
tame_2(struct igusa * I,struct igusa_p * Ip)833 tame_2(struct igusa *I, struct igusa_p *Ip)
834 {
835   long condp = -1, d, n, q, r;
836   GEN val = Ip->val;
837   (void)tame_234_init(I, Ip, &n, &q, &r);
838   d = n * (6*val[6]-5*val[7]) / 6;
839   switch(n)
840   {
841     case 1: condp = 1;
842       Ip->type = stack_sprintf("[I{%ld-0-0}] page 170", d);
843       Ip->neron = cyclic(d); break;
844     case 2:
845       switch(r)
846       {
847         case 0: condp = 4;
848           Ip->type = stack_sprintf("[I*{%ld-0-0}] page 171",d/2);
849           Ip->neron = shallowconcat(dicyclic(2,2),groupH(d/2)); break;
850         case 1:
851           switch(q)
852           {
853             case 0: condp = 2;
854               Ip->type = stack_sprintf("[II*{%ld-0}] page 172",d/2);
855               Ip->neron = cyclic(1); break;
856             case 1: condp = 3;
857               Ip->type = stack_sprintf("[II{%ld-0}] page 171",d/2);
858               Ip->neron = cyclic(2*d); break;
859             default: pari_err_BUG("tame2 [bug10]");
860           }
861           break;
862         default: pari_err_BUG("tame2 [bug11]");
863       }
864       break;
865     case 3: condp = 3;
866       Ip->neron = cyclic(d);
867       switch(r)
868       {
869         case 1:
870           Ip->type = stack_sprintf("[II{%ld}-IV] page 175", (d-2)/3);
871           break;
872         case 2:
873           Ip->type = stack_sprintf("[II{%ld}-IV*] page 175", (d-1)/3);
874           break;
875         default: pari_err_BUG("tame2 [bug12]");
876       }
877       break;
878     case 4:
879       switch(r)
880       {
881         case 1:
882           switch(q)
883           {
884             case 1: condp = 3;
885               Ip->type = stack_sprintf("[II{%ld}-III] page 177",(d-2)/4);
886               Ip->neron = cyclic(d/2); break;
887             case 3: condp = 4;
888               Ip->type = stack_sprintf("[II*{%ld}-III*] page 178",(d-2)/4);
889               Ip->neron = cyclic(8); break;
890             default: pari_err_BUG("tame2 [bug13]");
891           }
892           break;
893         case 3:
894           switch(q)
895           {
896             case 1: condp = 4;
897               Ip->type = stack_sprintf("[II*{%ld}-III] page 178",(d-2)/4);
898               Ip->neron = cyclic(8); break;
899             case 3: condp = 3;
900               Ip->type = stack_sprintf("[II{%ld}-III*] page 178",(d-2)/4);
901               Ip->neron = cyclic(d/2); break;
902             default: pari_err_BUG("tame2 [bug14]");
903           }
904           break;
905         default: pari_err_BUG("tame2 [bug15]");
906       }
907       break;
908     case 6:
909       switch(r)
910       {
911         case 2: condp = 4;
912           Ip->type = stack_sprintf("[II*-II*{%ld}] page 176", (d-4)/6);
913           Ip->neron = groupH((d+2)/6); break;
914         case 4: condp = 4;
915           Ip->type = stack_sprintf("[II-II*{%ld}] page 176", (d-2)/6);
916           Ip->neron = groupH((d+4)/6); break;
917         default: pari_err_BUG("tame2 [bug16]");
918       }
919       break;
920     default: pari_err_BUG("tame2 [bug17]");
921   }
922   return condp;
923 }
924 
925 /* Ip->tt = 3 */
926 static long
tame_3(struct igusa * I,struct igusa_p * Ip)927 tame_3(struct igusa *I, struct igusa_p *Ip)
928 {
929   long condp = -1, n, q, r, va5, d1, d2;
930   long flc = tame_234_init(I, Ip, &n, &q, &r);
931   GEN val = Ip->val;
932 
933   va5 = 2*val[6]-5*val[3];
934   d1 = minss(n * (val[7]-3*val[3]), n * va5 / 4);
935   d2 = n * va5 / 2 - d1;
936   switch(n)
937   {
938     case 1: condp = 2;
939       Ip->type = stack_sprintf("[I{%ld-%ld-0}] page 179", d1,d2);
940       Ip->neron = dicyclic(d1,d2); break;
941     case 2:
942       switch(r)
943       {
944         case 0: condp = 4;
945           Ip->type = stack_sprintf("[I*{%ld-%ld-0}] page 180", d1/2,d2/2);
946           Ip->neron = shallowconcat(groupH(d1/2),groupH(d2/2)); break;
947         case 1: condp = 3;
948           if (flc)
949           {
950             Ip->type = stack_sprintf("[2I{%ld}-0] page 181", d1);
951             Ip->neron = cyclic(d1);
952           }
953           else
954           { /* FIXME: "or" same with d1<->d2 */
955             Ip->type = stack_sprintf("[II{%ld-%ld}] page 182",d1/2,d2/2);
956             Ip->neron = ((d1*d2-4)&7)? cyclic(2*d1): dicyclic(d1,2);
957           }
958           break;
959         default: pari_err_BUG("tame3 [bug20]");
960       }
961       break;
962     case 4: condp = 4;
963       Ip->type = stack_sprintf("[III{%ld}] page 182", d1/2);
964       Ip->neron = groupH(d1/2); break;
965     default: pari_err_BUG("tame3 [bug21]");
966   }
967   return condp;
968 }
969 
970 /* Ip->tt = 4 */
971 static long
tame_4(struct igusa * I,struct igusa_p * Ip)972 tame_4(struct igusa *I, struct igusa_p *Ip)
973 {
974   long condp = -1, d1,d2,d3, f1,f2, g, h, n, q, r, vl,vn,vm, e1,e2,e3;
975   GEN val = Ip->val;
976   (void)tame_234_init(I, Ip, &n, &q, &r);
977   vl = val[6]-5*val[1];
978   vn = val[7]-6*val[1];
979   vm = val[2]-2*val[1]; /* all >= 0 */
980   e1 = min3(2*vl, 3*vn, 6*vm);
981   e2 = minss(6*vl - e1, 12*vn - 2*e1); /* >= 0 */
982   e3 = 12*vl - (2*e1+e2); /* >= 0 */
983   d1 = e1*n / 6;
984   d2 = e2*n / 12;
985   d3 = e3*n / 12;
986   g = d1*d2 + d1*d3 + d2*d3;
987   h = ugcd(ugcd(d1,d2),d3);
988   switch(n)
989   {
990     case 1: condp = 2;
991       Ip->type = stack_sprintf("[I{%ld-%ld-%ld}] page 182",d1,d2,d3);
992       Ip->neron = dicyclic(h,g/h); break;
993     case 2:
994       switch(r)
995       {
996         case 0: condp = 4;
997           Ip->type = stack_sprintf("[I*{%ld-%ld-%ld}] page 183",d1/2,d2/2,d3/2);
998           Ip->neron = shallowconcat(groupH(g/4), groupH(2-((h&2)>>1))); break;
999         case 1:
1000           if      (d1 == d2 || d1 == d3) f2 = d1;
1001           else if (d2 == d3) f2 = d2;
1002           else {
1003             pari_err_BUG("tame4 [bug23]");
1004             return -1; /*LCOV_EXCL_LINE*/
1005           }
1006           f1 = d1+d2+d3-2*f2;
1007           switch(q)
1008           {
1009             case 0: condp = 3;
1010               Ip->type = stack_sprintf("[II*{%ld-%ld}] page 184", f1/2,f2);
1011               Ip->neron = cyclic(f2); break;
1012             case 1: condp = 3;
1013               Ip->type = stack_sprintf("[II{%ld-%ld}] page 183", f1/2,f2);
1014               Ip->neron = cyclic(2*f1+f2); break;
1015             default: pari_err_BUG("tame4 [bug24]");
1016           }
1017           break;
1018         default: pari_err_BUG("tame4 [bug25]");
1019       }
1020       break;
1021     case 3: condp = 4;
1022       Ip->type = stack_sprintf("[III{%ld}] page 184",d1);
1023       Ip->neron = (d1%3)? cyclic(9): dicyclic(3,3); break;
1024     case 6: condp = 4;
1025       Ip->type = stack_sprintf("[III*{%ld}] page 184",d1/2);
1026       Ip->neron = cyclic(1); break;
1027     default: pari_err_BUG("tame4 [bug26]");
1028   }
1029   return condp;
1030 }
1031 
1032 /* p = 3 */
1033 static void
tame_567_init_3(struct igusa_p * Ip,long dk,long * pd,long * pn,long * pdm,long * pr)1034 tame_567_init_3(struct igusa_p *Ip, long dk,
1035                 long *pd, long *pn, long *pdm, long *pr)
1036 {
1037   long n = 1 + Ip->r1/6;
1038   *pd = n * dk / 36; /* / (12*Ip->eps) */
1039   *pn = n;
1040   *pr = -1; /* unused */
1041   *pdm = 0;
1042 }
1043 
1044 /* (4.3) */
1045 static void
tame_567_init(struct igusa * I,struct igusa_p * Ip,long dk,long * pd,long * pn,long * pdm,long * pr)1046 tame_567_init(struct igusa *I, struct igusa_p *Ip, long dk,
1047               long *pd, long *pn, long *pdm, long *pr)
1048 {
1049   long ndk, ddk;
1050   GEN p = Ip->p, val = Ip->val;
1051 
1052   if (equaliu(p,3)) { tame_567_init_3(Ip, dk, pd, pn, pdm, pr); return; }
1053   /* assume p > 3, Ip->eps = 1 */
1054   ssQ_red(dk, 12, &ndk, &ddk);
1055   if (! odd(val[8]))
1056   {
1057     long va0 = myval(I->a0,p), va2 = myval(I->A2,p), va3 = myval(I->A3,p);
1058     long va5 = myval(I->A5,p), vb2 = myval(I->B2,p);
1059     long v1 = 2*va3-4*va0-val[1],   v2 = 6*va5-20*va0-5*val[1];
1060     long v3 = 3*vb2-2*va0-2*val[1], v4 = 10*vb2-2*va5-5*val[1];
1061     if (v3 >= 0 && v2 >= 0 && v1 >= 0)
1062     {
1063       if (v1==0 || v2==0) get_nr(ddk, va0+val[1], 6,pn,pr); /* Prop 4.3.1 (a) */
1064       else
1065       { /* Prop 4.3.1 (d) */
1066         long v5 = myval(subii(mulii(I->A2,I->A3),mului(3,I->A5)),p);
1067         if (gequal0(I->A2)) pari_err_BUG("tame567 [bug27]");
1068         get_nr(ddk, 12*va0 + min3(dk, 6*va3-9*va2, 4*v5 - 10*va2), 24, pn,pr);
1069       }
1070     }
1071     else if (v2 < 0 && v4 >= 0)
1072       get_nr(ddk, 2*va5+val[1], 8, pn,pr); /* Prop 4.3.1 (b) */
1073     else /* (v3 < 0 && v4 < 0) */
1074       get_nr(ddk, vb2, 4, pn,pr); /* Prop 4.3.1 (c) */
1075     *pd = (*pn/ddk) * ndk;
1076   }
1077   else
1078   {
1079     *pr = ndk;
1080     *pn = 2*ddk;
1081     *pd = 2*ndk;
1082   }
1083   *pdm = umodsu(*pd, *pn);
1084 }
1085 
1086 static long
tame_5(struct igusa * I,struct igusa_p * Ip)1087 tame_5(struct igusa *I, struct igusa_p *Ip)
1088 {
1089   long condp = -1, d, n, dm, r, dk;
1090   GEN val = Ip->val;
1091 
1092   dk = Ip->eps*val[6]-5*val[8];
1093   tame_567_init(I, Ip, dk, &d, &n, &dm, &r);
1094   if (! odd(val[8]))
1095   {
1096     switch(n)
1097     {
1098       case 1: condp = 0;
1099         Ip->type = stack_sprintf("[I{0}-I{0}-%ld] page 158", d);
1100         Ip->neron = cyclic(1); break;
1101       case 2:
1102         switch(dm)
1103         {
1104           case 0: condp = 4;
1105             Ip->type = stack_sprintf("[I*{0}-I*{0}-%ld] page 158",(d-2)/2);
1106             Ip->neron = mkvecsmall4(2,2,2,2); break;
1107           case 1: condp = 2;
1108             Ip->type = stack_sprintf("[I{0}-I*{0}-%ld] page 159",(d-1)/2);
1109             Ip->neron = dicyclic(2,2); break;
1110         }
1111         break;
1112       case 3:
1113         switch(dm)
1114         {
1115           case 0: condp = 4;
1116             Ip->type = stack_sprintf("[IV-IV*-%ld] page 165",(d-3)/3);
1117             Ip->neron = dicyclic(3,3); break;
1118           case 1:
1119             switch(r)
1120             {
1121               case 0: case 1: condp = 2;
1122                 Ip->type = stack_sprintf("[I{0}-IV-%ld] page 160",(d-1)/3);
1123                 Ip->neron = cyclic(3); break;
1124               case 2: condp = 4;
1125                 Ip->type = stack_sprintf("[IV*-IV*-%ld] page 166",(d-4)/3);
1126                 Ip->neron = dicyclic(3,3); break;
1127             }
1128             break;
1129           case 2:
1130             switch(r)
1131             {
1132               case 0: case 2: condp = 2;
1133                 Ip->type = stack_sprintf("[I{0}-IV*-%ld] page 160",(d-2)/3);
1134                 Ip->neron = cyclic(3); break;
1135               case 1: condp = 4;
1136                 Ip->type = stack_sprintf("[IV-IV-%ld] page 165",(d-2)/3);
1137                 Ip->neron = dicyclic(3,3); break;
1138             }
1139             break;
1140         }
1141         break;
1142       case 4:
1143         switch(dm)
1144         {
1145           case 0: condp = 4;
1146             Ip->type = stack_sprintf("[III-III*-%ld] page 169",(d-4)/4);
1147             Ip->neron = dicyclic(2,2); break;
1148           case 1:
1149             switch(r)
1150             {
1151               case 0: case 1: condp = 2;
1152                 Ip->type = stack_sprintf("[I{0}-III-%ld] page 161",(d-1)/4);
1153                 Ip->neron = cyclic(2); break;
1154               case 2: case 3: condp = 4;
1155                 Ip->type = stack_sprintf("[I*{0}-III*-%ld] page 162",(d-5)/4);
1156                 Ip->neron = mkvecsmall3(2,2,2); break;
1157             }
1158             break;
1159           case 2: condp = 4;
1160             Ip->neron = dicyclic(2,2);
1161             switch(r)
1162             {
1163               case 1:
1164                 Ip->type = stack_sprintf("[III-III-%ld] page 169",(d-2)/4);
1165                 break;
1166               case 3:
1167                 Ip->type = stack_sprintf("[III*-III*-%ld] page 169",(d-6)/4);
1168                 break;
1169               default: pari_err_BUG("tame5 [bug29]");
1170             }
1171             break;
1172           case 3:
1173             switch(r)
1174             {
1175               case 0: case 3: condp = 2;
1176                 Ip->type = stack_sprintf("[I{0}-III*-%ld] page 162",(d-3)/4);
1177                 Ip->neron = cyclic(2); break;
1178               case 1: case 2: condp = 4;
1179                 Ip->type = stack_sprintf("[I*{0}-III-%ld] page 162",(d-3)/4);
1180                 Ip->neron = mkvecsmall3(2,2,2); break;
1181             }
1182             break;
1183         }
1184         break;
1185       case 6:
1186         switch(dm)
1187         {
1188           case 0: condp = 4;
1189             Ip->type = stack_sprintf("[II-II*-%ld] page 163",(d-6)/6);
1190             Ip->neron = cyclic(1); break;
1191           case 1:
1192             switch(r)
1193             {
1194               case 0: case 1: condp = 2;
1195                 Ip->type = stack_sprintf("[I{0}-II-%ld] page 159",(d-1)/6);
1196                 Ip->neron = cyclic(1); break;
1197               case 2: case 5: condp = 4;
1198                 Ip->type = stack_sprintf("[II*-IV-%ld] page 164",(d-7)/6);
1199                 Ip->neron = cyclic(3); break;
1200               case 3: case 4: condp = 4;
1201                 Ip->type = stack_sprintf("[I*{0}-IV*-%ld] page 161",(d-7)/6);
1202                 Ip->neron = mkvecsmall2(6,2); break;
1203             }
1204             break;
1205           case 2:
1206             switch(r)
1207             {
1208               case 1: condp = 4;
1209                 Ip->type = stack_sprintf("[II-II-%ld] page 163",(d-2)/6);
1210                 Ip->neron = cyclic(1); break;
1211               case 3: case 5: condp = 4;
1212                 Ip->type = stack_sprintf("[I*{0}-II*-%ld] page 160",(d-8)/6);
1213                 Ip->neron = dicyclic(2,2); break;
1214               default: pari_err_BUG("tame5 [bug30]");
1215             }
1216             break;
1217           case 3:
1218             Ip->neron = cyclic(3);
1219             switch(r)
1220             {
1221               case 1: case 2: condp = 4;
1222                 Ip->type = stack_sprintf("[II-IV-%ld] page 164",(d-3)/6);
1223                 break;
1224               case 4: case 5: condp = 4;
1225                 Ip->type = stack_sprintf("[II*-IV*-%ld] page 164",(d-9)/6);
1226                 break;
1227               default: pari_err_BUG("tame5 [bug31]");
1228             }
1229             break;
1230           case 4:
1231             switch(r)
1232             {
1233               case 1: case 3: condp = 4;
1234                 Ip->type = stack_sprintf("[I*{0}-II-%ld] page 160",(d-4)/6);
1235                 Ip->neron = dicyclic(2,2); break;
1236               case 5: condp = 4;
1237                 Ip->type = stack_sprintf("[II*-II*-%ld] page 163",(d-10)/6);
1238                 Ip->neron = cyclic(1); break;
1239               default: pari_err_BUG("tame5 [bug32]");
1240             }
1241             break;
1242           case 5:
1243             switch(r)
1244             {
1245               case 0: case 5: condp = 2;
1246                 Ip->type = stack_sprintf("[I{0}-II*-%ld] page 160",(d-5)/6);
1247                 Ip->neron = cyclic(1); break;
1248               case 1: case 4: condp = 4;
1249                 Ip->type = stack_sprintf("[II-IV*-%ld] page 164",(d-5)/6);
1250                 Ip->neron = cyclic(3); break;
1251               case 2: case 3: condp = 4;
1252                 Ip->type = stack_sprintf("[I*{0}-IV-%ld] page 161",(d-5)/6);
1253                 Ip->neron = mkvecsmall2(6,2); break;
1254             }
1255             break;
1256           default: pari_err_BUG("tame5 [bug33]");
1257         }
1258         break;
1259       case 12:
1260         condp = 4;
1261         switch(dm)
1262         {
1263           case 1:
1264             switch(r)
1265             {
1266               case 3: case 10:
1267                 Ip->type = stack_sprintf("[II*-III-%ld] page 166",(d-13)/12);
1268                 Ip->neron = cyclic(2); break;
1269               case 4: case 9:
1270                 Ip->type = stack_sprintf("[III*-IV-%ld] page 167",(d-13)/12);
1271                 Ip->neron = cyclic(6); break;
1272               default: pari_err_BUG("tame5 [bug34]");
1273             }
1274             break;
1275           case 5:
1276             switch(r)
1277             {
1278               case 2: case 3:
1279                 Ip->type = stack_sprintf("[II-III-%ld] page 166",(d-5)/12);
1280                 Ip->neron = cyclic(2); break;
1281               case 8: case 9:
1282                 Ip->type = stack_sprintf("[III*-IV*-%ld] page 168",(d-17)/12);
1283                 Ip->neron = cyclic(6); break;
1284               default: pari_err_BUG("tame5 [bug35]");
1285             }
1286             break;
1287           case 7:
1288             switch(r)
1289             {
1290               case 3: case 4:
1291                 Ip->type = stack_sprintf("[III-IV-%ld] page 167",(d-7)/12);
1292                 Ip->neron = cyclic(6); break;
1293               case 9: case 10:
1294                 Ip->type = stack_sprintf("[II*-III*-%ld] page 167",(d-19)/12);
1295                 Ip->neron = cyclic(2); break;
1296               default: pari_err_BUG("tame5 [bug36]");
1297             }
1298             break;
1299           case 11:
1300             switch(r)
1301             {
1302               case 3: case 8:
1303                 Ip->type = stack_sprintf("[III-IV*-%ld] page 168",(d-11)/12);
1304                 Ip->neron = cyclic(6); break;
1305               case 2: case 9:
1306                 Ip->type = stack_sprintf("[II-III*-%ld] page 166",(d-11)/12);
1307                 Ip->neron = cyclic(2); break;
1308               default: pari_err_BUG("tame5 [bug37]");
1309             }
1310             break;
1311           default: pari_err_BUG("tame5 [bug38]");
1312         }
1313         break;
1314       default: pari_err_BUG("tame5 [bug39]");
1315     }
1316   }
1317   else
1318   {
1319     r %= (n >> 1);
1320     switch(n)
1321     {
1322       case 2: condp = 2;
1323         Ip->type = stack_sprintf("[2I{0}-%ld] page 159",(d/2));
1324         Ip->neron = cyclic(1); break;
1325       case 4: condp = 4;
1326         Ip->type = stack_sprintf("[2I*{0}-%ld] page 159",(d/2-1)/2);
1327         Ip->neron = dicyclic(2,2); break;
1328       case 6: condp = 4;
1329         Ip->neron = cyclic(3);
1330         switch(r)
1331           {
1332           case 1:
1333             Ip->type = stack_sprintf("[2IV-%ld] page 165",(d/2-1)/3);
1334             break;
1335           case 2:
1336             Ip->type = stack_sprintf("[2IV*-%ld] page 165",(d/2-2)/3);
1337             break;
1338           default: pari_err_BUG("tame5 [bug40]");
1339           }
1340         break;
1341       case 8: condp = 4;
1342         Ip->neron = cyclic(2);
1343         switch(r)
1344         {
1345           case 1:
1346             Ip->type = stack_sprintf("[2III-%ld] page 168",(d/2-1)/4);
1347             break;
1348           case 3:
1349             Ip->type = stack_sprintf("[2III*-%ld] page 168",(d/2-3)/4);
1350             break;
1351           default: pari_err_BUG("tame5 [bug41]");
1352         }
1353         break;
1354       case 12: condp = 4;
1355         Ip->neron = cyclic(1);
1356         switch(r)
1357         {
1358           case 1:
1359             Ip->type = stack_sprintf("[2II-%ld] page 162",(d/2-1)/6);
1360             break;
1361           case 5:
1362             Ip->type = stack_sprintf("[2II*-%ld] page 163",(d/2-5)/6);
1363             break;
1364           default: pari_err_BUG("tame5 [bug42]");
1365         }
1366         break;
1367       default: pari_err_BUG("tame5 [bug43]");
1368     }
1369   }
1370   return condp;
1371 }
1372 
1373 static long
tame_6(struct igusa * I,struct igusa_p * Ip)1374 tame_6(struct igusa *I, struct igusa_p *Ip)
1375 {
1376   long condp = -1, d, d1, n, dm, r, dk;
1377   GEN val = Ip->val;
1378 
1379   dk = Ip->eps*val[7]-6*val[8];
1380   tame_567_init(I, Ip, dk, &d, &n, &dm, &r);
1381   d1 = n * (Ip->eps*(val[6]-val[7])+val[8]) / Ip->eps;
1382   switch(n)
1383   {
1384     case 1: condp = 1;
1385       Ip->type = stack_sprintf("[I{0}-I{%ld}-%ld] page 170",d1,d);
1386       Ip->neron = cyclic(d1); break;
1387     case 2:
1388       switch(dm)
1389       {
1390         case 0: condp = 4;
1391           Ip->type=stack_sprintf("[I*{0}-I*{%ld}-%ld] page 171", d1/2,(d-2)/2);
1392           Ip->neron = shallowconcat(groupH(d1/2), dicyclic(2,2)); break;
1393         case 1: return -1;
1394         default: pari_err_BUG("tame6 [bug44]");
1395       }
1396       break;
1397     case 3: condp = 3;
1398       Ip->neron = dicyclic(3,d1/3);
1399       switch(dm)
1400       {
1401         case 1:
1402           Ip->type = stack_sprintf("[I{%ld}-IV-%ld] page 173",d1/3,(d-1)/3);
1403           break;
1404         case 2:
1405           Ip->type = stack_sprintf("[I{%ld}-IV*-%ld] page 173",d1/3,(d-2)/3);
1406           break;
1407         default: pari_err_BUG("tame6 [bug45]");
1408       }
1409       break;
1410     case 4:
1411       switch(dm)
1412       {
1413         case 1:
1414           switch(r)
1415           {
1416             case 0: case 1: condp = 3;
1417               Ip->type=stack_sprintf("[I{%ld}-III-%ld] page 176",d1/4,(d-1)/4);
1418               Ip->neron = dicyclic(2,d1/4); break;
1419             case 2: case 3: condp = 4;
1420               Ip->type=stack_sprintf("[I*{%ld}-III*-%ld] page 177",d1/4,(d-5)/4);
1421               Ip->neron = shallowconcat(groupH(d1/4), cyclic(2)); break;
1422             default: pari_err_BUG("tame6 [bug46]");
1423           }
1424           break;
1425         case 3:
1426           switch(r)
1427           {
1428             case 0: case 3: condp = 3;
1429               Ip->type=stack_sprintf("[I{%ld}-III*-%ld] page 176",d1/4,(d-3)/4);
1430               Ip->neron = dicyclic(2,d1/4); break;
1431             case 1: case 2: condp = 4;
1432               Ip->type=stack_sprintf("[I*{%ld}-III-%ld] page 177",d1/4,(d-3)/4);
1433               Ip->neron = shallowconcat(groupH(d1/4), cyclic(2)); break;
1434             default: pari_err_BUG("tame6 [bug47]");
1435           }
1436           break;
1437         default: pari_err_BUG("tame6 [bug48]");
1438       }
1439       break;
1440     case 6:
1441       switch(dm)
1442       {
1443         case 1:
1444           switch(r)
1445           {
1446             case 0: case 1: condp = 3;
1447               Ip->type = stack_sprintf("[I{%ld}-II-%ld] page 172",d1/6,(d-1)/6);
1448               Ip->neron = cyclic(d1/6); break;
1449             case 3: case 4: condp = 4;
1450               Ip->type=stack_sprintf("[I*{%ld}-IV*-%ld] page 174",d1/6,(d-7)/6);
1451               Ip->neron = shallowconcat(groupH(d1/6), cyclic(3)); break;
1452             default: pari_err_BUG("tame6 [bug49]");
1453           }
1454           break;
1455         case 2: condp = 4;
1456           Ip->type = stack_sprintf("[I*{%ld}-II*-%ld] page 174",d1/6,(d-8)/6);
1457           Ip->neron = groupH(d1/6); break;
1458         case 4: condp = 4;
1459           Ip->type = stack_sprintf("[I*{%ld}-II-%ld] page 173",d1/6,(d-4)/6);
1460           Ip->neron = groupH(d1/6); break;
1461         case 5:
1462           switch(r)
1463           {
1464             case 0: case 5: condp = 3;
1465               Ip->type=stack_sprintf("[I{%ld}-II*-%ld] page 172",d1/6,(d-5)/6);
1466               Ip->neron = cyclic(d1/6); break;
1467             case 2: case 3: condp = 4;
1468               Ip->type=stack_sprintf("[I*{%ld}-IV-%ld] page 174",d1/6,(d-5)/6);
1469               Ip->neron = shallowconcat(groupH(d1/6), cyclic(3)); break;
1470             default: pari_err_BUG("tame6 [bug50]");
1471           }
1472           break;
1473         default: pari_err_BUG("tame6 [bug51]");
1474       }
1475       break;
1476     default: pari_err_BUG("tame6 [bug52]");
1477   }
1478   return condp;
1479 }
1480 
1481 static long
tame_7(struct igusa * I,struct igusa_p * Ip)1482 tame_7(struct igusa *I, struct igusa_p *Ip)
1483 {
1484   long condp = -1, d, D, d1, d2, n, dm, r, dk;
1485   GEN val = Ip->val;
1486 
1487   dk = 3*(Ip->eps*val[3]-2*val[8]);
1488   tame_567_init(I, Ip, dk, &d, &n, &dm, &r);
1489   D = n * (Ip->eps*(val[6]-3*val[3])+val[8]) / Ip->eps;
1490   d1 = minss(n * (val[7]-3*val[3]), D/2);
1491   d2 = D - d1;
1492   /* d1 <= d2 */
1493   switch(n)
1494   {
1495     case 1: condp = 2;
1496       Ip->type = stack_sprintf("[I{%ld}-I{%ld}-%ld] page 179",d1,d2,d);
1497       Ip->neron = dicyclic(d1,d2); break;
1498     case 2:
1499       if (odd(val[8]))
1500       {
1501         condp = 3;
1502         Ip->type = stack_sprintf("[2I{%ld}-%ld] page 181",d1,d/2);
1503         Ip->neron = cyclic(d1);
1504       }
1505       else if (dm == 0)
1506       {
1507         condp = 4;
1508         Ip->type = stack_sprintf("[I*{%ld}-I*{%ld}-%ld] page 180", d1/2,d2/2,(d-2)/2);
1509         Ip->neron = shallowconcat(groupH(d1/2),groupH(d2/2));
1510       }
1511       else
1512       {
1513         GEN H;
1514         if (d1 != d2) return -1;
1515         condp = 3; H = groupH(d1/2);
1516         Ip->type = stack_sprintf("[I{%ld}-I*{%ld}-%ld] page 180", d1/2,d1/2,(d-1)/2);
1517         Ip->neron = shallowconcat(H, H);
1518       }
1519       break;
1520     case 4: condp = 4;
1521       Ip->type = stack_sprintf("[2I*{%ld}-%ld] page 181",d1/2,(d-2)/4);
1522       Ip->neron = groupH(d1/2); break;
1523     default: pari_err_BUG("tame7 [bug55]");
1524   }
1525   return condp;
1526 }
1527 
1528 static long labelm3(GEN polh, long t60, long alpha, long Dmin, struct igusa *I, struct igusa_p *Ip);
1529 static long
tame(GEN polh,long t60,long alpha,long Dmin,struct igusa * I,struct igusa_p * Ip)1530 tame(GEN polh, long t60, long alpha, long Dmin, struct igusa *I, struct igusa_p *Ip)
1531 {
1532   long d;
1533   Ip->tame = 1;
1534   switch(Ip->tt)
1535   {
1536     case 1: return tame_1(I,Ip);
1537     case 2: return tame_2(I,Ip);
1538     case 3: return tame_3(I,Ip);
1539     case 4: return tame_4(I,Ip);
1540     case 5: return tame_5(I,Ip);
1541     case 6: d = tame_6(I,Ip); break;
1542     default:d = tame_7(I,Ip); break;
1543   }
1544   if (d < 0) d = labelm3(polh,t60,alpha,Dmin,I,Ip); /* => tt=6 or 7 */
1545   return d;
1546 }
1547 
1548 /* maxc = maximum conductor valuation at p */
1549 static long
get_maxc(GEN p)1550 get_maxc(GEN p)
1551 {
1552   switch (itos_or_0(p))
1553   {
1554     case 2:  return 20; break;
1555     case 3:  return 10; break;
1556     case 5:  return 9; break;
1557     default: return 4; break; /* p > 5 */
1558   }
1559 }
1560 
1561 /* p = 3 */
1562 static long
quartic(GEN polh,long alpha,long Dmin,struct igusa_p * Ip)1563 quartic(GEN polh, long alpha, long Dmin, struct igusa_p *Ip)
1564 {
1565   GEN val = Ip->val, p = Ip->p;
1566   GEN polf = polymini_zi2(ZX_Z_mul(polh, powiu(p, alpha)));
1567   long condp = -1, d, R, r1, beta;
1568   r1 = polf[1];
1569   beta = polf[2];
1570   R = beta/2;
1571   switch(Ip->tt)
1572   {
1573     case 1: case 5: d = 0;break;
1574     case 3: d = val[6] - 5*val[3]/2;break;
1575     case 7: d = val[6] - 3*val[3] + val[8]/Ip->eps;break;
1576     default: pari_err_BUG("quartic [type choices]");
1577              d = 0; /*LCOV_EXCL_LINE*/
1578   }
1579   switch(r1)
1580   {
1581     case 0:
1582       if (d)
1583       {
1584         condp = 3;
1585         Ip->type = stack_sprintf("[2I{%ld}-%ld] page 181",d,R);
1586         Ip->neron = cyclic(d);
1587       }
1588       else
1589       {
1590         condp = 2;
1591         Ip->neron = cyclic(1);
1592         if (R) Ip->type = stack_sprintf("[2I{0}-%ld] page 159",R);
1593         else   Ip->type = "[II] page 155";
1594       }
1595       break;
1596     case 6: condp = 4;
1597       Ip->type = stack_sprintf("[2I*{%ld}-%ld] pages 159, 181",d,R);
1598       Ip->neron = dicyclic(2,2); break;
1599     case 3: condp = 4;
1600       Ip->type = stack_sprintf("[2III-%ld] page 168",R);
1601       Ip->neron = cyclic(2); break;
1602     case 9: condp = 4;
1603       Ip->type = stack_sprintf("[2III*-%ld] page 168",R);
1604       Ip->neron = cyclic(2); break;
1605     case 2: condp = Dmin-12*R-13;
1606       Ip->type = stack_sprintf("[2II-%ld] page 162",R);
1607       Ip->neron = cyclic(1); break;
1608     case 8: condp = Dmin-12*R-19;
1609       Ip->type = stack_sprintf("[2IV*-%ld] page 165",R);
1610       Ip->neron = cyclic(3); break;
1611     case 4: condp = Dmin-12*R-15;
1612       Ip->type = stack_sprintf("[2IV-%ld] page 165",R);
1613       Ip->neron = cyclic(3); break;
1614     case 10: condp = Dmin-12*R-21;
1615       Ip->type = stack_sprintf("[2II*-%ld] page 163",R);
1616       Ip->neron = cyclic(1); break;
1617     default: pari_err_BUG("quartic [type1]");
1618   }
1619   if (condp > get_maxc(p) || condp < 0) pari_err_BUG("quartic [conductor]");
1620   return condp;
1621 }
1622 
1623 static long
litredtp(long alpha,long alpha1,long t60,long t60_1,GEN polh,GEN polh1,long Dmin,long R,struct igusa * I,struct igusa_p * Ip)1624 litredtp(long alpha, long alpha1, long t60, long t60_1, GEN polh, GEN polh1,
1625          long Dmin, long R, struct igusa *I, struct igusa_p *Ip)
1626 {
1627   GEN val = Ip->val, p = Ip->p;
1628   long condp = -1, indice, d;
1629 
1630   if ((Ip->r1 == 0||Ip->r1 == 6) && (Ip->r2 == 0||Ip->r2 == 6))
1631   { /* (r1,r2) = (0,0), (0,6), (6,0) or (6,6) */
1632     if (Ip->tt == 5)
1633     {
1634       switch(Ip->r1 + Ip->r2)
1635       {
1636       case 0: /* (0,0) */
1637         condp = 0;
1638         Ip->type = stack_sprintf("[I{0}-I{0}-%ld] page 158",R);
1639         Ip->neron = cyclic(1); break;
1640       case 6: /* (0,6) or (6,0) */
1641         condp = 2;
1642         Ip->type = stack_sprintf("[I{0}-I*{0}-%ld] page 159",R);
1643         Ip->neron = dicyclic(2,2); break;
1644       case 12: /* (6,6) */
1645         condp = 4;
1646         Ip->type = stack_sprintf("[I*{0}-I*{0}-%ld] page 158",R);
1647         Ip->neron = mkvecsmall4(2,2,2,2); break;
1648       }
1649       return condp;
1650     }
1651     if (Ip->r1 == Ip->r2) return tame(polh, t60, alpha, Dmin, I, Ip);
1652     if (Ip->tt == 6)
1653     {
1654       d = val[6] - val[7] + val[8]/Ip->eps;
1655       if (Ip->r1 && alpha1 == 0) /* H(px) / p^3 */
1656         polh1 = ZX_Z_divexact(ZX_unscale_div(polh1,p), sqri(p));
1657       if (FpX_is_squarefree(FpX_red(polh1,p),p))
1658       { indice = 0; condp = 3-Ip->r2/6; }
1659       else
1660       { indice = d; condp = 3-Ip->r1/6; }
1661     }
1662     else
1663     { /* Ip->tt == 7 */
1664       long d1;
1665       d = val[6] - 3*val[3] + val[8]/Ip->eps;
1666       if (t60_1 == 60) /* H(px) / p^3 */
1667         polh1 = ZX_Z_divexact(ZX_unscale_div(polh1,p), sqri(p));
1668       d1 = minss(val[7]-3*val[3],d/2);
1669       if (d == 2*d1) indice = d1;
1670       else
1671       {
1672         indice = discpart(polh1,p,d1+1);
1673         if (indice>= d1+1) indice = d-d1; else indice = d1;
1674       }
1675       condp = 3;
1676     }
1677     if (Ip->r1) indice = d - indice; /* (r1,r2) = (6,0) */
1678     Ip->neron = shallowconcat(cyclic(indice),groupH(d-indice));
1679     Ip->type = stack_sprintf("[I{%ld}-I*{%ld}-%ld] page %ld",
1680                              indice,d-indice,R, (Ip->tt==6)? 170L: 180L);
1681     return condp;
1682   }
1683   if (Ip->tt == 7) pari_err_BUG("litredtp [switch ri]");
1684   {
1685     struct red __S1, __S2, *S1 = &__S1, *S2 = &__S2;
1686     long f1 = get_red(S1, Ip, polh1, p, alpha1, Ip->r1);
1687     long f2 = get_red(S2, Ip, polh,  p, alpha,  Ip->r2);
1688     /* reorder to normalize representation */
1689     if (S1->tnum > S2->tnum || (S1->tnum == S2->tnum && f1 > f2))
1690     { struct red *S = S1; S1 = S2; S2 = S; }
1691     Ip->type = stack_sprintf("[%s-%s-%ld] pages %s", S1->t,S2->t, R, S1->pages);
1692     Ip->neron = shallowconcat(S1->g, S2->g);
1693     condp = Dmin - (f1 + f2) + ((R >= 0)? 2-12*R: 4);
1694   }
1695   if (condp > get_maxc(p)) pari_err_BUG("litredtp [conductor]");
1696   return condp;
1697 }
1698 
1699 static long
labelm3(GEN h1,long t60_1,long alpha1,long Dmin,struct igusa * I,struct igusa_p * Ip)1700 labelm3(GEN h1, long t60_1, long alpha1, long Dmin, struct igusa *I, struct igusa_p *Ip)
1701 {
1702   GEN h, pm, vs, val = Ip->val, p = Ip->p;
1703   long alpha, t60, lambda, beta, R;
1704 
1705   pm = polymini(ZX_Z_mul(RgX_recip6(h1), powiu(p,alpha1)), p);
1706   h  = gel(pm,1); vs = gel(pm,2);
1707   lambda= vs[1];
1708   t60   = vs[2];
1709   alpha = vs[3];
1710   beta  = vs[5];
1711   if (lambda != 3) pari_err_BUG("labelm3 [lambda != 3]");
1712   R = beta-(alpha1+alpha);
1713   if (odd(R)) pari_err_BUG("labelm3 [R odd]");
1714   R /= 2;
1715   if (R <= -2) pari_err_BUG("labelm3 [R <= -2]");
1716   if (val[8] % (2*Ip->eps)) pari_err_BUG("labelm3 [val(eps2)]");
1717   if (R >= 0 && (alpha+alpha1) >= 1) pari_err_BUG("labelm3 [minimal equation]");
1718   Ip->r1 = t60_1 / 10 + 6*alpha1;
1719   Ip->r2 = t60 / 10 + 6*alpha;
1720   return litredtp(alpha, alpha1, t60, t60_1, h, h1, Dmin, R, I, Ip);
1721 }
1722 
1723 /* p = 3 */
1724 static long
quadratic(GEN polh,long alpha,long Dmin,struct igusa * I,struct igusa_p * Ip)1725 quadratic(GEN polh, long alpha, long Dmin, struct igusa *I, struct igusa_p *Ip)
1726 {
1727   long alpha1 = alpha, beta, t6, R;
1728   GEN vs = polymini_zi(ZX_Z_mul(polh, powiu(Ip->p,alpha)));
1729   t6 = vs[1];
1730   alpha = vs[2];
1731   beta  = vs[3];
1732   R = beta-alpha;
1733   if (R >= 0 && alpha1)
1734   {
1735     Dmin -= 10;
1736     if (DEBUGLEVEL)
1737       err_printf("(Care: minimal discriminant over Z[i] smaller than over Z)\n");
1738   }
1739   Ip->r2 = Ip->r1 = t6 + 6*alpha;
1740   return litredtp(alpha, alpha, t6*10, t6*10, polh, polh, Dmin, R, I, Ip);
1741 }
1742 
1743 static long
genus2localred(struct igusa * I,struct igusa_p * Ip,GEN p,GEN polmini)1744 genus2localred(struct igusa *I, struct igusa_p *Ip, GEN p, GEN polmini)
1745 {
1746   GEN val, vs, polh, list, c1, c2, c3, c4, c5, c6, prod;
1747   long i, vb5, vb6, d, Dmin, alpha, lambda, t60;
1748   long condp = -1, indice, vc6, mm, nb, dism;
1749 
1750   stable_reduction(I, Ip, p);
1751   val = Ip->val; Dmin = val[6];
1752   if (Dmin == 0)
1753   {
1754     Ip->tame = 1;
1755     Ip->type = "[I{0-0-0}] page 155";
1756     Ip->neron = cyclic(1); return 0;
1757   }
1758   if (Dmin == 1)
1759   {
1760     Ip->type = "[I{1-0-0}] page 170";
1761     Ip->neron = cyclic(1); return 1;
1762   }
1763   if (Dmin == 2) switch(Ip->tt)
1764   {
1765     case 2:
1766       Ip->type = "[I{2-0-0}] page 170";
1767       Ip->neron = cyclic(2); return 1;
1768     case 3:
1769       Ip->type = "[I{1-1-0}] page 179";
1770       Ip->neron = cyclic(1); return 2;
1771     case 5:
1772       if (cmpis(p,3) <= 0) pari_err_BUG("genus2localred [tt 1]");
1773       Ip->type = "[I{0}-II-0] page 159";
1774       Ip->neron = cyclic(1); return 2;
1775     default: pari_err_BUG("genus2localred [tt 2]");
1776   }
1777   if (absequaliu(p,2)) return -1;
1778   polh = gel(polmini,1); vs = gel(polmini,2);
1779   lambda = vs[1];
1780   t60    = vs[2];
1781   alpha  = vs[3];
1782   if (vs[4]) return equaliu(p,3)? quadratic(polh, alpha, Dmin, I, Ip):
1783                                   tame(polh, t60, alpha, Dmin, I, Ip);
1784   if (!t60 && lambda<= 2)
1785   {
1786     if (Ip->tt >= 5) pari_err_BUG("genus2localred [tt 3]");
1787     return tame(polh, t60, alpha, Dmin, I, Ip);
1788   }
1789   if (Dmin == 3)
1790   {
1791     switch(Ip->tt)
1792     {
1793       case 2: return tame(polh, t60, alpha, Dmin, I, Ip);
1794       case 3: Ip->type = "[I{2-1-0}] page 179"; Ip->neron = cyclic(2); return 2;
1795       case 4: Ip->type = "[I{1-1-1}] page 182"; Ip->neron = cyclic(3); return 2;
1796       case 5:
1797         if (equaliu(p,3) && t60 != 30)
1798           return labelm3(polh,t60,alpha,Dmin,I,Ip);
1799         Ip->type = "[I{0}-III-0] page 161"; Ip->neron = cyclic(2); return 2;
1800       case 6:
1801         if (equaliu(p,3)) pari_err_BUG("genus2localred [conductor]");
1802         Ip->type = "[I{1}-II-0] page 172"; Ip->neron = cyclic(1); return 3;
1803     }
1804     pari_err_BUG("genus2localred [switch tt 4]");
1805     return -1; /* LCOV_EXCL_LINE */
1806   }
1807   switch(lambda)
1808   {
1809     case 0:
1810       switch(t60+alpha)
1811       {
1812         case 10:
1813           condp = Dmin-1;
1814           Ip->type = "[V] page 156";
1815           Ip->neron = cyclic(3); break;
1816         case 11:
1817           condp = Dmin-11;
1818           Ip->type = "[V*] page 156";
1819           Ip->neron = cyclic(3); break;
1820         case 12:
1821           condp = Dmin-2;
1822           Ip->type = "[IX-2] page 157";
1823           Ip->neron = cyclic(5); break;
1824         case 13:
1825           condp = Dmin-12;
1826           Ip->type = "[VIII-4] page 157";
1827           Ip->neron = cyclic(1); break;
1828         case 24:
1829           condp = Dmin-8;
1830           Ip->type = "[IX-4] page 158";
1831           Ip->neron = cyclic(5);
1832           break;
1833         case 15: case 16:
1834           if (Ip->tt>= 5) pari_err_BUG("genus2localred [tt 6]");
1835           return tame(polh, t60, alpha, Dmin, I, Ip);
1836         case 20: case 21:
1837           {
1838             GEN b0, b1, b2, b3, b4, b5, b6, b02, b03, b04, b05;
1839             RgX_to_06(polh, &b0,&b1,&b2,&b3,&b4,&b5,&b6);
1840             vb5 = myval(b5,p);
1841             vb6 = myval(b6,p);
1842             if (vb6 >= 3)
1843             {
1844               if (vb5 < 2) pari_err_BUG("genus2localred [red1]");
1845               if (vb5 >= 3)
1846               {
1847                 condp = Dmin-8;
1848                 Ip->type = "[II*-IV-(-1)] page 164";
1849                 Ip->neron = cyclic(3);
1850               }
1851               else
1852               {
1853                 condp = Dmin-7;
1854                 Ip->type = "[IV-III*-(-1)] page 167";
1855                 Ip->neron = cyclic(6);
1856               }
1857               break;
1858             }
1859             if (dvdii(b0,p)) pari_err_BUG("genus2localred [b0]");
1860             b02 = gsqr(b0);
1861             b03 = gmul(b02, b0);
1862             b04 = gmul(b03, b0);
1863             b05 = gmul(b04, b0);
1864             c1 = gmul2n(b1,-1);
1865             c2 = gmul2n(gsub(gmul(b0,b2), gsqr(c1)),-1);
1866             c3 = gmul2n(gsub(gmul(b02,b3), gmul2n(gmul(c1,c2),1)),-1);
1867             c4 = gsub(gmul(b03,b4), gadd(gmul2n(gmul(c1,c3),1),gsqr(c2)));
1868             c5 = gsub(gmul(b04,b5), gmul2n(gmul(c2,c3),1));
1869             c6 = gsub(gmul(b05,b6), gsqr(c3));
1870             /* b0^5*H(x/b0) = (x^3+c1*x^2+c2*x+c3)^2+c4*x^2+c5*x+c6 */
1871             vc6 = myval(c6,p);
1872             if (vc6 == 2)
1873             {
1874               if (alpha)
1875               {
1876                 condp = Dmin-16;
1877                 Ip->type = "[IV] page 155";
1878                 Ip->neron = cyclic(1);
1879               }
1880               else
1881               {
1882                 condp = Dmin-6;
1883                 Ip->type = "[III] page 155";
1884                 Ip->neron = dicyclic(3,3);
1885               }
1886             }
1887             else
1888             {
1889               if (myval(c3,p) > 1) pari_err_BUG("genus2localred [c3]");
1890               mm = min3(3*myval(c4,p)-4, 3*myval(c5,p)-5, 3*vc6-6);
1891               if (alpha)
1892               {
1893                 condp = Dmin-mm-16;
1894                 Ip->type = stack_sprintf("[III*{%ld}] page 184", mm);
1895                 Ip->neron = cyclic(1);
1896               }
1897               else
1898               {
1899                 condp = Dmin-mm-6;
1900                 Ip->type = stack_sprintf("[III{%ld}] page 184", mm);
1901                 Ip->neron = (mm%3)? cyclic(9): dicyclic(3,3);
1902               }
1903             }
1904           }
1905           break;
1906         case 30:
1907           return equaliu(p,3)? quartic(polh, alpha, Dmin, Ip)
1908                              : tame(polh, t60, alpha, Dmin, I, Ip);
1909         default: pari_err_BUG("genus2localred [red2]");
1910       }
1911       break;
1912     case 1:
1913       switch(t60+alpha)
1914       {
1915         case 12:
1916           condp = Dmin;
1917           Ip->type = "[VIII-1] page 156";
1918           Ip->neron = cyclic(1); break;
1919         case 13:
1920           condp = Dmin-10;
1921           Ip->type = "[IX-3] page 157";
1922           Ip->neron = cyclic(5); break;
1923         case 24:
1924           condp = Dmin-4;
1925           Ip->type = "[IX-1] page 157";
1926           Ip->neron = cyclic(5); break;
1927         case 25:
1928           condp = Dmin-14;
1929           Ip->type = "[VIII-3] page 157";
1930           Ip->neron = cyclic(1); break;
1931         case 36:
1932           condp = Dmin-8;
1933           Ip->type = "[VIII-2] page 157";
1934           Ip->neron = cyclic(1); break;
1935         case 15:
1936           condp = Dmin-1;
1937           Ip->type = "[VII] page 156";
1938           Ip->neron = cyclic(2); break;
1939         case 16:
1940           condp = Dmin-11;
1941           Ip->type = "[VII*] page 156";
1942           Ip->neron = cyclic(2); break;
1943         case 20:
1944           if (cmpis(p,3))
1945           {
1946             d = 6*val[6]-5*val[7]-2;
1947             if (d%6) pari_err_BUG("genus2localred [index]");
1948             dism = (d/6);
1949           }
1950           else
1951           {
1952             list = padicfactors(polh,p,Dmin-5);
1953             nb = lg(list);
1954             prod = pol_1(varn(polh));
1955             for(i = 1;i<nb;i++)
1956             {
1957               GEN c = gel(list,i);
1958               if (valp(gel(c,2)) && degpol(c)<= 2) prod = RgX_mul(prod,c);
1959             }
1960             if (degpol(prod) > 2) pari_err_BUG("genus2localred [padicfactors]");
1961             dism = valp(RgX_disc(prod)) - 1;
1962           }
1963           condp = Dmin-dism-3;
1964           Ip->type = stack_sprintf("[II-II*{%ld}] page 176", dism);
1965           Ip->neron = groupH(dism+1); break;
1966         case 21:
1967           vb6 = myval(RgX_coeff(polh,0),p);
1968           if (vb6<2) pari_err_BUG("genus2localred [red3]");
1969           condp = Dmin-14;
1970           Ip->type = "[IV*-II{0}] page 175";
1971           Ip->neron = cyclic(1); break;
1972         case 30:
1973           vb5 = myval(RgX_coeff(polh,1),p);
1974           if (vb5 == 2)
1975           {
1976             if (Ip->tt >= 5) pari_err_BUG("genus2localred [tt 6]");
1977             return tame(polh, t60, alpha, Dmin, I, Ip);
1978           }
1979           condp = Dmin-7;
1980           Ip->type = "[II*-III-(-1)] page 167";
1981           Ip->neron = cyclic(2); break;
1982       }
1983       break;
1984     case 2:
1985       if (ugcd(t60, 60) == 15) /* denom(theta) = 4 */
1986       {
1987         if (Ip->tt>4) pari_err_BUG("genus2localred [tt 5]");
1988         return tame(polh, t60, alpha, Dmin, I, Ip);
1989       }
1990       if (!equaliu(p,3) && ugcd(t60, 60) == 20) /* denom(theta) = 3 */
1991         return tame(polh, t60, alpha, Dmin, I, Ip);
1992       list = padicfactors(polh,p,Dmin-10*alpha);
1993       nb = lg(list); prod = pol_1(varn(polh));
1994       for(i = 1;i<nb;i++)
1995       {
1996         GEN c = gel(list,i);
1997         if (!valp(gel(c,2))) prod = RgX_mul(prod,c);
1998       }
1999       switch(degpol(prod))
2000       {
2001         GEN e0, e1, e2;
2002         case 0:
2003           dism = 0; break;
2004         case 1:
2005           e1 = gel(prod,3);
2006           dism = 2*valp(e1); break;
2007         case 2:
2008           e0 = gel(prod,2);
2009           e1 = gel(prod,3);
2010           e2 = gel(prod,4);
2011           dism = valp(gsub(gsqr(e1),gmul2n(gmul(e0,e2),2))); break;
2012         default:
2013           pari_err_BUG("genus2localred [padicfactors 2]");
2014           dism = 0;
2015       }
2016       switch(t60/5+alpha-4)
2017       {
2018         case 0:
2019           condp = Dmin-dism-1;
2020           Ip->type = stack_sprintf("[IV-II{%ld}] page 175", dism);
2021           Ip->neron = cyclic(3*dism+2); break;
2022         case 1:
2023           condp = Dmin-dism-10;
2024           Ip->type = stack_sprintf("[II*-II*{%ld}] page 176",dism);
2025           Ip->neron = groupH(dism+1); break;
2026         case 2: case 3:
2027           if (myval(RgX_coeff(polh,0),p) == 2)
2028           {
2029             if (Ip->tt>4) pari_err_BUG("genus2localred [tt 5]");
2030             return tame(polh, t60, alpha, Dmin, I, Ip);
2031           }
2032           dism++;
2033           indice = val[6]-(5*val[3]/2)-dism;
2034           condp = Dmin-dism-indice-2;
2035           Ip->type = stack_sprintf("[II{%ld-%ld}] page 182", dism,indice);
2036           Ip->neron = both_odd(dism,indice)? dicyclic(2,2*dism): cyclic(4*dism);
2037           break;
2038         case 4:
2039           condp = Dmin-dism-5;
2040           Ip->type = stack_sprintf("[IV*-II{%ld}] page 175",dism+1);
2041           Ip->neron = cyclic(3*dism+4); break;
2042       }
2043       break;
2044     case 3:
2045       if (!equaliu(p,3) || Ip->tt <= 4)
2046         return tame(polh, t60, alpha, Dmin, I, Ip);
2047       return labelm3(polh,t60,alpha,Dmin,I,Ip); /* p = 3 */
2048     default: pari_err_BUG("genus2localred [switch lambda]");
2049   }
2050   if (condp < 2 || condp > get_maxc(p))
2051     pari_err_BUG("genus2localred [conductor]");
2052   return condp;
2053 }
2054 
2055 static long
chk_pol(GEN P)2056 chk_pol(GEN P) {
2057   switch(typ(P))
2058   {
2059     case t_INT: break;
2060     case t_POL: RgX_check_ZX(P,"genus2red"); return varn(P); break;
2061     default: pari_err_TYPE("genus2red", P);
2062   }
2063   return -1;
2064 }
2065 
2066 /* P,Q are ZX, study Y^2 + Q(X) Y = P(X) */
2067 GEN
genus2red(GEN PQ,GEN p)2068 genus2red(GEN PQ, GEN p)
2069 {
2070   pari_sp av = avma;
2071   struct igusa I;
2072   GEN P, Q, D;
2073   GEN j22, j42, j2j6, a0,a1,a2,a3,a4,a5,a6, V,polr,facto,factp, vecmini, cond;
2074   long i, l, dd, vP,vQ;
2075 
2076   PQ = Q_remove_denom(PQ, &D);
2077   if (typ(PQ) == t_VEC && lg(PQ) == 3)
2078   {
2079     P = gel(PQ,1);
2080     Q = gel(PQ,2);
2081   }
2082   else
2083   {
2084     P = PQ;
2085     Q = gen_0;
2086   }
2087 
2088   vP = chk_pol(P);
2089   vQ = chk_pol(Q);
2090   if (vP < 0)
2091   {
2092     if (vQ < 0) pari_err_TYPE("genus2red",mkvec2(P,Q));
2093     P = scalarpol(P,vQ);
2094   }
2095   else if (vQ < 0) Q = scalarpol(Q,vP);
2096   if (p && typ(p) != t_INT) pari_err_TYPE("genus2red", p);
2097   if (D) P = ZX_Z_mul(P,D);
2098 
2099   polr = ZX_add(ZX_sqr(Q), gmul2n(P,2)); /* ZX */
2100   switch(degpol(polr))
2101   {
2102     case 5: case 6: break;
2103     default: pari_err_DOMAIN("genus2red","genus","!=", gen_2,mkvec2(P,Q));
2104   }
2105 
2106   RgX_to_03(polr, &a0,&a1,&a2,&a3);
2107   I.j10 = !signe(a0)? mulii(sqri(a1), ZX_disc(polr)): ZX_disc(polr);
2108   if (!signe(I.j10))
2109     pari_err_DOMAIN("genus2red","genus","<",gen_2,mkvec2(P,Q));
2110   I.j10 = gmul2n(I.j10, -12); /* t_INT */
2111 
2112   if (p == NULL)
2113   {
2114     facto = absZ_factor(I.j10);
2115     factp = gel(facto,1);
2116   }
2117   else
2118   {
2119     factp = mkcol(p);
2120     facto = mkmat2(factp, mkcol(gen_1));
2121   }
2122   l = lg(factp);
2123   vecmini = cgetg(l, t_COL);
2124   for(i = 1; i<l; i++)
2125   {
2126     GEN l = gel(factp,i), pm;
2127     if (i == 1 && absequaliu(l, 2)) { gel(vecmini,1) = gen_0; continue; }
2128     gel(vecmini,i) = pm = polymini(polr, l);
2129     polr = ZX_Q_mul(gel(pm,1), powiu(l, gel(pm,2)[3]));
2130   }
2131   RgX_to_06(polr, &a0,&a1,&a2,&a3,&a4,&a5,&a6);
2132   I.j10 = !signe(a0)? mulii(sqri(a1), ZX_disc(polr)): ZX_disc(polr);
2133   I.j10 = gmul2n(I.j10,-12);
2134 
2135   I.a0 = a0;
2136   I.A2 = apol2(a0,a1,a2);
2137   I.A3 = apol3(a0,a1,a2,a3);
2138   I.A5 = apol5(a0,a1,a2,a3,a4,a5);
2139   I.B2 = bpol2(a0,a1,a2,a3,a4);
2140 
2141   I.j2 = igusaj2(a0,a1,a2,a3,a4,a5,a6);
2142   I.j4 = igusaj4(a0,a1,a2,a3,a4,a5,a6);
2143   I.i4 = gsub(gsqr(I.j2), gmulsg(24,I.j4));
2144   I.j6 = igusaj6(a0,a1,a2,a3,a4,a5,a6);
2145   j42 = gsqr(I.j4);
2146   j22 = gsqr(I.j2);
2147   j2j6 = gmul(I.j2,I.j6);
2148   I.j8 = gmul2n(gsub(j2j6,j42), -2);
2149   I.i12= gmul2n(gsub(gadd(gmul(j22,j42),gmulsg(36,gmul(j2j6,I.j4))),
2150                      gadd(gadd(gmulsg(32,gmul(j42,I.j4)),gmul(j2j6,j22)),gmulsg(108,gsqr(I.j6)))),-2);
2151 
2152   for(i = 1; i < l; i++)
2153     gcoeff(facto,i,2) = stoi(Q_pval(I.j10, gel(factp,i)));
2154   dd = ZX_pval(polr,gen_2) & (~1); /* = 2 floor(val/2) */
2155   polr = gmul2n(polr, -dd);
2156 
2157   V = cgetg(l, t_VEC);
2158   for (i = 1; i < l; i++)
2159   {
2160     GEN q = gel(factp,i), red, N = NULL;
2161     struct igusa_p Ip;
2162     long f = genus2localred(&I, &Ip, q, gel(vecmini,i));
2163     gcoeff(facto,i,2) = stoi(f);
2164     if (Ip.tame) Ip.type = stack_strcat("(tame) ", Ip.type);
2165     if (f >= 0)
2166       N = zv_snf(Ip.neron);
2167     if (DEBUGLEVEL)
2168     {
2169       if (!p) err_printf("p = %Ps\n", q);
2170       err_printf("(potential) stable reduction: %Ps\n", Ip.stable);
2171       if (f >= 0) {
2172         err_printf("reduction at p: %s, %Ps", Ip.type, N);
2173         err_printf(", f = %ld\n", f);
2174       }
2175     }
2176     red = f >= 0? mkvec2(strtoGENstr(Ip.type), N): cgetg(1, t_VEC);
2177     gel(V, i) = mkvec3(q, Ip.stable, red);
2178   }
2179   if (p) V = gel(V,1);
2180   cond = factorback(facto);
2181   /* remove denominator 2 coming from f = -1 in genuslocalred(, p = 2) */
2182   if (typ(cond) != t_INT) cond = gel(cond,1);
2183   return gerepilecopy(av, mkvec4(cond, facto, polr, V));
2184 }
2185