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