Line data Source code
1 : /* Copyright (C) 2020 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 : #include "pari.h"
15 : #include "paripriv.h"
16 :
17 : #define DEBUGLEVEL DEBUGLEVEL_bnf
18 :
19 : /* if x a famat, assume it is a true unit (very costly to check even that
20 : * it's an algebraic integer) */
21 : GEN
22 1645 : bnfisunit(GEN bnf, GEN x)
23 : {
24 1645 : long tx = typ(x), i, r1, RU, e, n, prec;
25 1645 : pari_sp av = avma;
26 : GEN t, logunit, ex, nf, pi2_sur_w, rx, emb, arg;
27 :
28 1645 : bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
29 1645 : RU = lg(bnf_get_logfu(bnf));
30 1645 : n = bnf_get_tuN(bnf); /* # { roots of 1 } */
31 1645 : if (tx == t_MAT)
32 : { /* famat, assumed OK */
33 1309 : if (lg(x) != 3) pari_err_TYPE("bnfisunit [not a factorization]", x);
34 : } else {
35 336 : x = nf_to_scalar_or_basis(nf,x);
36 336 : if (typ(x) != t_COL)
37 : { /* rational unit ? */
38 : GEN v;
39 : long s;
40 140 : if (typ(x) != t_INT || !is_pm1(x)) return cgetg(1,t_COL);
41 133 : s = signe(x); set_avma(av); v = zerocol(RU);
42 133 : gel(v,RU) = utoi((s > 0)? 0: n>>1);
43 133 : return v;
44 : }
45 196 : if (!isint1(Q_denom(x))) { set_avma(av); return cgetg(1,t_COL); }
46 : }
47 :
48 1505 : r1 = nf_get_r1(nf);
49 1505 : prec = nf_get_prec(nf);
50 1505 : for (i = 1;; i++)
51 0 : {
52 : GEN rlog;
53 1505 : nf = bnf_get_nf(bnf);
54 1505 : logunit = bnf_get_logfu(bnf);
55 1505 : rlog = real_i(logunit);
56 1505 : rx = nflogembed(nf,x,&emb, prec);
57 1505 : if (rx)
58 : {
59 1505 : GEN logN = RgV_sum(rx); /* log(Nx), should be ~ 0 */
60 1505 : if (gexpo(logN) > -20)
61 : { /* precision problem ? */
62 14 : if (typ(logN) != t_REAL) { set_avma(av); return cgetg(1,t_COL); } /*no*/
63 28 : if (i == 1 && typ(x) != t_MAT &&
64 28 : !is_pm1(nfnorm(nf, x))) { set_avma(av); return cgetg(1, t_COL); }
65 : }
66 : else
67 : {
68 504 : ex = RU == 1? cgetg(1,t_COL)
69 1491 : : RgM_solve(rlog, rx); /* ~ fundamental units exponents */
70 1491 : if (ex) { ex = grndtoi(ex, &e); if (e < -4) break; }
71 : }
72 : }
73 0 : if (i == 1)
74 0 : prec = nbits2prec(gexpo(x) + 128);
75 : else
76 : {
77 0 : if (i > 4) pari_err_PREC("bnfisunit");
78 0 : prec = precdbl(prec);
79 : }
80 0 : if (DEBUGLEVEL) pari_warn(warnprec,"bnfisunit",prec);
81 0 : bnf = bnfnewprec_shallow(bnf, prec);
82 : }
83 : /* choose a large embedding => small relative error */
84 2069 : for (i = 1; i < RU; i++)
85 1418 : if (signe(gel(rx,i)) > -1) break;
86 1491 : if (RU == 1) t = gen_0;
87 : else
88 : {
89 987 : t = imag_i( row_i(logunit,i, 1,RU-1) );
90 987 : t = RgV_dotproduct(t, ex);
91 987 : if (i > r1) t = gmul2n(t, -1);
92 : }
93 1491 : if (typ(emb) != t_MAT) arg = garg(gel(emb,i), prec);
94 : else
95 : {
96 1309 : GEN p = gel(emb,1), e = gel(emb,2);
97 1309 : long j, l = lg(p);
98 1309 : arg = NULL;
99 62132 : for (j = 1; j < l; j++)
100 : {
101 60823 : GEN a = gmul(gel(e,j), garg(gel(gel(p,j),i), prec));
102 60823 : arg = arg? gadd(arg, a): a;
103 : }
104 : }
105 1491 : t = gsub(arg, t); /* = arg(the missing root of 1) */
106 1491 : pi2_sur_w = divru(mppi(prec), n>>1); /* 2pi / n */
107 1491 : e = umodiu(roundr(divrr(t, pi2_sur_w)), n);
108 1491 : if (n > 2)
109 : {
110 7 : GEN z = algtobasis(nf, bnf_get_tuU(bnf)); /* primitive root of 1 */
111 7 : GEN ro = RgV_dotproduct(row(nf_get_M(nf), i), z);
112 7 : GEN p2 = roundr(divrr(garg(ro, prec), pi2_sur_w));
113 7 : e *= Fl_inv(umodiu(p2,n), n);
114 7 : e %= n;
115 : }
116 1491 : gel(ex,RU) = utoi(e); setlg(ex, RU+1); return gerepilecopy(av, ex);
117 : }
118 :
119 : /* split M a square ZM in HNF as [H, B; 0, Id], H in HNF without 1-eigenvalue */
120 : static GEN
121 1008 : hnfsplit(GEN M, GEN *pB)
122 : {
123 1008 : long i, l = lg(M);
124 3510 : for (i = l-1; i; i--)
125 3209 : if (!equali1(gcoeff(M,i,i))) break;
126 1008 : if (!i) { *pB = zeromat(0, l-1); return cgetg(1, t_MAT); }
127 707 : *pB = matslice(M, 1, i, i+1, l-1); return matslice(M, 1, i, 1, i);
128 : }
129 :
130 : /* S a list of prime ideal in idealprimedec format. If pH != NULL, set it to
131 : * the HNF of the S-class group and return bnfsunit, else return bnfunits */
132 : static GEN
133 1442 : bnfsunit_i(GEN bnf, GEN S, GEN *pH, GEN *pA, GEN *pden)
134 : {
135 1442 : long FLAG, i, lS = lg(S);
136 : GEN M, U1, U2, U, V, H, Sunit, B, g;
137 :
138 1442 : if (!is_vec_t(typ(S))) pari_err_TYPE("bnfsunit",S);
139 1442 : bnf = checkbnf(bnf);
140 1442 : if (lS == 1)
141 : {
142 434 : *pA = cgetg(1,t_MAT);
143 434 : *pden = gen_1; return cgetg(1,t_VEC);
144 : }
145 1008 : M = cgetg(lS,t_MAT); /* relation matrix for the S class group */
146 1008 : g = cgetg(lS,t_MAT); /* principal part */
147 1008 : FLAG = pH ? 0: nf_GENMAT;
148 6083 : for (i = 1; i < lS; i++)
149 : {
150 5075 : GEN pr = gel(S,i);
151 5075 : checkprid(pr);
152 5075 : if (pH)
153 1834 : gel(M,i) = isprincipal(bnf, pr);
154 : else
155 : {
156 3241 : GEN v = bnfisprincipal0(bnf, pr, FLAG);
157 3241 : gel(M,i) = gel(v,1);
158 3241 : gel(g,i) = gel(v,2);
159 : }
160 : }
161 : /* S class group and S units, use ZM_hnflll to get small 'U' */
162 1008 : M = shallowconcat(M, diagonal_shallow(bnf_get_cyc(bnf)));
163 1008 : H = ZM_hnflll(M, &U, 1); setlg(U, lS); if (pH) *pH = H;
164 1008 : U1 = rowslice(U,1, lS-1);
165 1008 : U2 = rowslice(U,lS, lg(M)-1); /* (M | cyc) [U1; U2] = 0 */
166 1008 : H = ZM_hnflll(U1, pH? NULL: &V, 0);
167 : /* U1 = upper left corner of U, invertible. S * U1 = principal ideals
168 : * whose generators generate the S-units */
169 1008 : H = hnfsplit(H, &B);
170 : /* [ H B ] [ H^-1 - H^-1 B ]
171 : * U1 * V = HNF(U1) = [ 0 Id ], inverse = [ 0 Id ]
172 : * S * HNF(U1) = integral generators for S-units = Sunit */
173 1008 : Sunit = cgetg(lS, t_VEC);
174 1008 : if (pH)
175 : {
176 161 : long nH = lg(H) - 1;
177 161 : FLAG = nf_FORCE | nf_GEN;
178 1995 : for (i = 1; i < lS; i++)
179 : {
180 1834 : GEN C = NULL, E;
181 1834 : if (i <= nH) E = gel(H,i); else { C = gel(S,i), E = gel(B,i-nH); }
182 1834 : gel(Sunit,i) = gel(isprincipalfact(bnf, C, S, E, FLAG), 2);
183 : }
184 : }
185 : else
186 : {
187 847 : GEN cycgen = bnf_build_cycgen(bnf);
188 847 : U1 = ZM_mul(U1, V);
189 847 : U2 = ZM_mul(U2, V);
190 847 : FLAG = nf_FORCE | nf_GENMAT;
191 4088 : for (i = 1; i < lS; i++)
192 : {
193 3241 : GEN a = famatV_factorback(g, gel(U1,i));
194 3241 : GEN b = famatV_factorback(cycgen, ZC_neg(gel(U2,i)));
195 3241 : gel(Sunit,i) = famat_reduce(famat_mul(a, b));
196 : }
197 : }
198 1008 : H = ZM_inv(H, pden);
199 1008 : *pA = shallowconcat(H, ZM_neg(ZM_mul(H,B))); /* top inverse * den */
200 1008 : return Sunit;
201 : }
202 : GEN
203 168 : bnfsunit(GEN bnf,GEN S,long prec)
204 : {
205 168 : pari_sp av = avma;
206 168 : long i, l = lg(S);
207 168 : GEN v, R, nf, A, den, U, cl, H = NULL;
208 168 : bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
209 168 : v = cgetg(7, t_VEC);
210 168 : gel(v,1) = U = bnfsunit_i(bnf, S, &H, &A, &den);
211 168 : gel(v,2) = mkvec2(A, den);
212 168 : gel(v,3) = cgetg(1,t_VEC); /* dummy */
213 168 : R = bnf_get_reg(bnf);
214 168 : cl = bnf_get_clgp(bnf);
215 168 : if (l != 1)
216 : {
217 161 : GEN u,A, G = bnf_get_gen(bnf), D = ZM_snf_group(H,NULL,&u), h = ZV_prod(D);
218 161 : long lD = lg(D);
219 161 : A = cgetg(lD, t_VEC);
220 182 : for(i = 1; i < lD; i++) gel(A,i) = idealfactorback(nf, G, gel(u,i), 1);
221 161 : cl = mkvec3(h, D, A);
222 161 : R = mpmul(R, h);
223 1995 : for (i = 1; i < l; i++)
224 : {
225 1834 : GEN pr = gel(S,i), p = pr_get_p(pr);
226 1834 : long f = pr_get_f(pr);
227 1834 : R = mpmul(R, logr_abs(itor(p,prec)));
228 1834 : if (f != 1) R = mulru(R, f);
229 1834 : gel(U,i) = nf_to_scalar_or_alg(nf, gel(U,i));
230 : }
231 : }
232 168 : gel(v,4) = R;
233 168 : gel(v,5) = cl;
234 168 : gel(v,6) = S; return gerepilecopy(av, v);
235 : }
236 : GEN
237 1092 : bnfunits(GEN bnf, GEN S)
238 : {
239 1092 : pari_sp av = avma;
240 : GEN A, den, U, fu, tu;
241 1092 : bnf = checkbnf(bnf);
242 1092 : U = bnfsunit_i(bnf, S? S: cgetg(1,t_VEC), NULL, &A, &den);
243 1092 : if (!S) S = cgetg(1,t_VEC);
244 1092 : fu = bnf_compactfu(bnf);
245 1092 : if (!fu)
246 : {
247 : long i, l;
248 231 : fu = bnf_has_fu(bnf); if (!fu) bnf_build_units(bnf);
249 231 : fu = shallowcopy(fu); l = lg(fu);
250 686 : for (i = 1; i < l; i++) gel(fu,i) = to_famat_shallow(gel(fu,i),gen_1);
251 : }
252 1092 : tu = nf_to_scalar_or_basis(bnf_get_nf(bnf), bnf_get_tuU(bnf));
253 1092 : U = shallowconcat(U, vec_append(fu, to_famat_shallow(tu,gen_1)));
254 1092 : return gerepilecopy(av, mkvec4(U, S, A, den));
255 : }
256 : GEN
257 182 : sunits_mod_units(GEN bnf, GEN S)
258 : {
259 182 : pari_sp av = avma;
260 : GEN A, den;
261 182 : bnf = checkbnf(bnf);
262 182 : return gerepilecopy(av, bnfsunit_i(bnf, S, NULL, &A, &den));
263 : }
264 :
265 : /* v_S(x), x in famat form */
266 : static GEN
267 105 : sunit_famat_val(GEN nf, GEN S, GEN x)
268 : {
269 105 : long i, l = lg(S);
270 105 : GEN v = cgetg(l, t_VEC);
271 1085 : for (i = 1; i < l; i++) gel(v,i) = famat_nfvalrem(nf, x, gel(S,i), NULL);
272 105 : return v;
273 : }
274 : /* v_S(x), x algebraic number */
275 : static GEN
276 1218 : sunit_val(GEN nf, GEN S, GEN x, GEN N)
277 : {
278 1218 : long i, l = lg(S);
279 1218 : GEN v = zero_zv(l-1), N0 = N;
280 37093 : for (i = 1; i < l; i++)
281 : {
282 35875 : GEN P = gel(S,i), p = pr_get_p(P);
283 35875 : if (dvdii(N, p)) { v[i] = nfval(nf,x,P); (void)Z_pvalrem(N0, p, &N0); }
284 : }
285 1218 : return is_pm1(N0)? v: NULL;
286 : }
287 :
288 : /* if *px a famat, assume it's an S-unit */
289 : static GEN
290 1624 : make_unit(GEN nf, GEN U, GEN *px)
291 : {
292 1624 : GEN den, v, w, A, gen = gel(U,1), S = gel(U,2), x = *px;
293 1624 : long cH, i, l = lg(S);
294 :
295 1624 : if (l == 1) return cgetg(1, t_COL);
296 1596 : A = gel(U,3); den = gel(U,4);
297 1596 : cH = nbrows(A);
298 1596 : if (typ(x) == t_MAT && lg(x) == 3)
299 : {
300 105 : w = sunit_famat_val(nf, S, x); /* x = S v */
301 105 : v = ZM_ZC_mul(A, w);
302 105 : w += cH; w[0] = evaltyp(t_COL) | _evallg(lg(A) - cH);
303 : }
304 : else
305 : {
306 : GEN N;
307 1491 : x = nf_to_scalar_or_basis(nf,x);
308 1491 : switch(typ(x))
309 : {
310 497 : case t_INT: N = x; if (!signe(N)) return NULL; break;
311 0 : case t_FRAC: N = mulii(gel(x,1),gel(x,2)); break;
312 994 : default: { GEN d = Q_denom(x); N = mulii(idealnorm(nf,gmul(x,d)), d); }
313 : }
314 : /* relevant primes divide N */
315 1491 : if (is_pm1(N)) return zerocol(l-1);
316 1218 : w = sunit_val(nf, S, x, N);
317 1218 : if (!w) return NULL;
318 1204 : v = ZM_zc_mul(A, w);
319 1204 : w += cH; w[0] = evaltyp(t_VECSMALL) | _evallg(lg(A) - cH);
320 1204 : w = zc_to_ZC(w);
321 : }
322 2519 : if (!is_pm1(den)) for (i = 1; i <= cH; i++)
323 : {
324 : GEN r;
325 1210 : gel(v,i) = dvmdii(gel(v,i), den, &r);
326 1210 : if (r != gen_0) return NULL;
327 : }
328 1309 : v = shallowconcat(v, w); /* append bottom of w (= [0 Id] part) */
329 38094 : for (i = 1; i < l; i++)
330 : {
331 36785 : GEN e = gel(v,i);
332 36785 : if (signe(e)) x = famat_mulpow_shallow(x, gel(gen,i), negi(e));
333 : }
334 1309 : *px = x; return v;
335 : }
336 :
337 : static GEN
338 1624 : bnfissunit_i(GEN bnf, GEN x, GEN U)
339 : {
340 1624 : GEN w, nf, v = NULL;
341 1624 : bnf = checkbnf(bnf);
342 1624 : nf = bnf_get_nf(bnf);
343 1624 : if ( (w = make_unit(nf, U, &x)) ) v = bnfisunit(bnf, x);
344 1624 : if (!v || lg(v) == 1) return NULL;
345 1596 : return mkvec2(v, w);
346 : }
347 : static int
348 126 : checkU(GEN U)
349 : {
350 126 : if (typ(U) != t_VEC || lg(U) != 5) return 0;
351 126 : return typ(gel(U,1)) == t_VEC && is_vec_t(typ(gel(U,2)))
352 252 : && typ(gel(U,4))==t_INT;
353 : }
354 : GEN
355 161 : bnfisunit0(GEN bnf, GEN x, GEN U)
356 : {
357 161 : pari_sp av = avma;
358 : GEN z;
359 161 : if (!U) return bnfisunit(bnf, x);
360 126 : if (!checkU(U)) pari_err_TYPE("bnfisunit",U);
361 126 : z = bnfissunit_i(bnf, x, U);
362 126 : if (!z) { set_avma(av); return cgetg(1,t_COL); }
363 112 : return gerepilecopy(av, shallowconcat(gel(z,2), gel(z,1)));
364 : }
365 :
366 : /* OBSOLETE */
367 : static int
368 1498 : checkbnfS_i(GEN v)
369 : {
370 : GEN S, g, w;
371 1498 : if (typ(v) != t_VEC || lg(v) != 7) return 0;
372 1498 : g = gel(v,1); w = gel(v,2); S = gel(v,6);
373 1498 : if (typ(g) != t_VEC || !is_vec_t(typ(S)) || lg(S) != lg(g)) return 0;
374 1498 : return typ(w) == t_VEC && lg(w) == 3;
375 : }
376 : /* OBSOLETE */
377 : GEN
378 1498 : bnfissunit(GEN bnf, GEN bnfS, GEN x)
379 : {
380 1498 : pari_sp av = avma;
381 : GEN z, U;
382 1498 : if (!checkbnfS_i(bnfS)) pari_err_TYPE("bnfissunit",bnfS);
383 1498 : U = mkvec4(gel(bnfS,1), gel(bnfS,6), gmael(bnfS,2,1), gmael(bnfS,2,2));
384 1498 : z = bnfissunit_i(bnf, x, U);
385 1498 : if (!z) { set_avma(av); return cgetg(1,t_COL); }
386 1484 : return gerepilecopy(av, shallowconcat(gel(z,1), gel(z,2)));
387 : }
|
