Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
10 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 : #include "pari.h"
14 : #include "paripriv.h"
15 :
16 : #define DEBUGLEVEL DEBUGLEVEL_gchar
17 :
18 : static GEN gchari_eval(GEN gc, GEN chi, GEN x, long flag, GEN logchi, GEN s0, long prec);
19 :
20 : /*********************************************************************/
21 : /** **/
22 : /** HECKE GROSSENCHARACTERS **/
23 : /** contributed by Pascal Molin and Aurel Page (2022) **/
24 : /** **/
25 : /*********************************************************************/
26 : /* characters can be represented by:
27 : * - t_COL of coordinates on the snf basis (mostly for gp use): prefix gchar_
28 : * - t_VEC of coordinates on the internal basis: prefix gchari_
29 : * - t_VEC of R/Z components (logs): prefix gcharlog_
30 : *
31 : * see gchar_internal for SNF -> internal
32 : * and gchari_duallog for internal -> R/Z components */
33 :
34 : /* localstar: represent (Z_F/m)^* . {+-1}^r + generators of U_{i-1}(p)/U_i
35 : * structure:
36 : * - [ sprk for p^k | m ] , size np
37 : * - [ Ufil_p for p | m ], size np
38 : * - m_oo, t_VECSMALL of size nci <= r1 (indices of real places)
39 : *
40 : * where Ufil_p = [ Mat([gen[j], t_COL of size ncp]_j) ]_{1<=i<=k} */
41 :
42 : #define GC_LENGTH 12
43 : #define LOCS_LENGTH 4
44 :
45 : static GEN
46 623 : compute_Lcyc(GEN Lsprk, GEN moo)
47 : {
48 623 : long i, l = lg(Lsprk), len = l+lg(moo)-1;
49 623 : GEN Lcyc = cgetg(len,t_VEC);
50 1239 : for (i = 1; i < l; i++) gel(Lcyc,i) = sprk_get_cyc(gel(Lsprk,i));
51 742 : for ( ; i < len; i++) gel(Lcyc,i) = mkvec(gen_2);
52 623 : return Lcyc;
53 : }
54 :
55 : static long
56 2016 : sprk_get_ncp(GEN sprk) { return lg(sprk_get_cyc(sprk)) - 1; }
57 :
58 : /* true nf; modulus = [ factor(m_f), m_oo ] */
59 : static GEN
60 623 : localstar(GEN nf, GEN famod, GEN moo)
61 : {
62 623 : GEN Lcyc, cyc, Lsprk, Lgenfil, P = gel(famod,1), E = gel(famod,2);
63 623 : long i, l = lg(P);
64 :
65 623 : Lsprk = cgetg(l, t_VEC);
66 623 : Lgenfil = cgetg(l, t_VEC);
67 1239 : for (i = 1; i < l; i++)
68 : {
69 616 : long e, k = itos(gel(E,i));
70 616 : GEN v, sprk = log_prk_init(nf, gel(P,i), k, NULL);
71 616 : gel(Lsprk,i) = sprk;
72 616 : gel(Lgenfil,i) = v = cgetg(k+1, t_VEC);
73 : /* log on sprk of generators of U_{e-1}/U_e(pr) */
74 616 : gel(v, 1) = col_ei(sprk_get_ncp(sprk), 1);
75 1267 : for (e = 2; e <= k; e++) gel(v, e) = sprk_log_gen_pr2(nf, sprk, e);
76 : }
77 623 : Lcyc = compute_Lcyc(Lsprk, moo);
78 623 : if (lg(Lcyc) > 1)
79 308 : cyc = shallowconcat1(Lcyc);
80 : else
81 315 : cyc = cgetg(1, t_VEC);
82 623 : return mkvec4(cyc, Lsprk, Lgenfil, mkvec2(famod,moo));
83 : }
84 :
85 : /* (nv * log|x^sigma|/norm, arg(x^sigma))/2*Pi
86 : * substract norm so that we project to the hyperplane
87 : * H : sum n_s x_s = 0 */
88 : static GEN
89 194923 : nfembedlog(GEN *pnf, GEN x, long prec)
90 : {
91 194923 : pari_sp av = avma;
92 194923 : GEN logs, cxlogs, nf = *pnf;
93 194923 : long k, r1, r2, n, extrabit, extranfbit = 0, nfprec, nfprec0, logprec;
94 :
95 194923 : nfprec0 = nf_get_prec(nf);
96 194923 : nf_get_sign(nf, &r1, &r2);
97 194923 : n = r1 + 2*r2;
98 194923 : logprec = prec;
99 194923 : extrabit = expu(n) + gexpo(nf_get_M(nf)) + 10;
100 194923 : if (typ(x) == t_MAT)
101 : {
102 194300 : long l = lg(gel(x,2));
103 194300 : if (l > 1)
104 : {
105 194259 : extrabit += expu(l) + gexpo(gel(x,2));
106 194259 : extranfbit = gexpo(gel(x,1));
107 : }
108 : }
109 : else
110 : {
111 623 : x = nf_to_scalar_or_basis(nf,x);
112 623 : extranfbit = gexpo(x);
113 : }
114 194923 : if (DEBUGLEVEL>3)
115 0 : err_printf(" nfembedlog: prec=%d extrabit=%d nfprec=%d extralogprec=%d\n",
116 : prec, nbits2extraprec(extrabit + extranfbit), nfprec0,
117 : nbits2extraprec(extrabit));
118 194923 : nfprec = prec + nbits2extraprec(extrabit + extranfbit);
119 194923 : logprec = prec + nbits2extraprec(extrabit);
120 194923 : if (nfprec0 < nfprec)
121 : {
122 3253 : if (DEBUGLEVEL)
123 0 : err_printf(" nfembedlog: increasing prec %d -> %d\n", nfprec0, nfprec);
124 3253 : *pnf = nf = nfnewprec_shallow(nf, nfprec);
125 3253 : av = avma;
126 : }
127 194923 : if (!(cxlogs = nf_cxlog(nf, x, logprec))) return gc_NULL(av);
128 194923 : if (!(cxlogs = nf_cxlog_normalize(nf, cxlogs, logprec))) return gc_NULL(av);
129 194923 : logs = cgetg(n+1,t_COL);
130 720620 : for (k = 1; k <= r1+r2; k++) gel(logs,k) = real_i(gel(cxlogs,k));
131 411283 : for ( ; k <= n; k++) gel(logs,k) = gmul2n(imag_i(gel(cxlogs,k-r2)), -1);
132 194923 : extrabit = gexpo(logs);
133 194923 : if (extrabit < 0) extrabit = 0;
134 194923 : prec += nbits2extraprec(extrabit);
135 194923 : return gerepileupto(av, gdiv(logs, Pi2n(1,prec)));
136 : }
137 :
138 : static GEN
139 1904 : gchar_Sval(GEN nf, GEN S, GEN x)
140 : {
141 1904 : GEN res = cgetg(lg(S),t_COL);
142 : long i;
143 1904 : if (typ(x)==t_MAT)
144 3500 : for (i=1; i<lg(S); i++)
145 2219 : gel(res, i) = famat_nfvalrem(nf, x, gel(S,i), NULL);
146 : else
147 1001 : for (i=1; i<lg(S); i++)
148 378 : gel(res, i) = stoi(nfval(nf, x, gel(S,i)));
149 1904 : return res;
150 : }
151 :
152 : /* true nf; log_prk(x*pi_pr^{-v_pr(x)}), sign(sigma(x)) */
153 : static GEN
154 190372 : gchar_logm(GEN nf, GEN locs, GEN x)
155 : {
156 190372 : GEN moo, loga, Lsprk = locs_get_Lsprk(locs);
157 190372 : long i, l = lg(Lsprk);
158 :
159 190372 : moo = locs_get_m_infty(locs);
160 190372 : if (typ(x) != t_MAT) x = to_famat_shallow(x, gen_1);
161 190372 : loga = cgetg(l+1, t_VEC);
162 607852 : for (i = 1; i < l; i++)
163 : {
164 417480 : GEN sprk = gel(Lsprk, i), pr = sprk_get_pr(sprk);
165 417480 : GEN g = vec_append(gel(x,1), pr_get_gen(pr));
166 417480 : GEN v = famat_nfvalrem(nf, x, pr, NULL);
167 417480 : GEN e = vec_append(gel(x,2), gneg(v));
168 417480 : gel(loga, i) = famat_zlog_pr(nf, g, e, sprk, NULL);
169 : }
170 190372 : gel(loga, i) = zc_to_ZC( nfsign_arch(nf, x, moo) );
171 190372 : return shallowconcat1(loga);
172 : }
173 :
174 : static GEN
175 1904 : gchar_nflog(GEN *pnf, GEN zm, GEN S, GEN x, long prec)
176 : {
177 1904 : GEN emb = nfembedlog(pnf, x, prec);
178 1904 : if (!emb) return NULL;
179 1904 : return shallowconcat1(mkvec3(gchar_Sval(*pnf,S,x),
180 : gchar_logm(*pnf,zm,x), emb));
181 : }
182 :
183 : /*******************************************************************/
184 : /* */
185 : /* CHARACTER GROUP */
186 : /* */
187 : /*******************************************************************/
188 :
189 : /* embed v in vector of length size, shifted to the right */
190 : static GEN
191 6657 : embedcol(GEN v, long size, long shift)
192 : {
193 : long k;
194 6657 : GEN w = zerocol(size);
195 51296 : for (k = 1; k < lg(v); k++) gel(w, shift+k) = gel(v,k);
196 6657 : return w;
197 : }
198 :
199 : /* write exact rationals from real approximation, in place. */
200 : static void
201 5050 : shallow_clean_rat(GEN v, long k0, long k1, GEN den, long prec)
202 : {
203 5050 : long k, e, bit = -prec2nbits(prec)/2;
204 45440 : for (k = k0; k <= k1; k++)
205 : {
206 40390 : GEN t = gel(v,k);
207 40390 : if (den) t = gmul(t, den);
208 40390 : t = grndtoi(t, &e);
209 40390 : if (DEBUGLEVEL>1) err_printf("[%Ps*%Ps=%Ps..e=%ld|prec=%ld]\n",
210 0 : gel(v,k), den? den: gen_1, t, e, prec);
211 40390 : if (e > bit)
212 : pari_err_BUG("gcharinit, non rational entry"); /*LCOV_EXCL_LINE*/
213 40390 : gel(v, k) = den? gdiv(t, den): t;
214 : }
215 5050 : }
216 :
217 : /* FIXME: export ? */
218 : static GEN
219 1211 : rowreverse(GEN m)
220 : {
221 : long i, l;
222 : GEN v;
223 1211 : if (lg(m) == 1) return m;
224 1211 : l = lgcols(m); v = cgetg(l, t_VECSMALL);
225 4501 : for (i = 1; i < l; i++) v[i] = l - i;
226 1211 : return rowpermute(m, v);
227 : }
228 :
229 : /* lower-left hnf on subblock m[r0+1..r0+nr, c0+1..c0+nc]
230 : * return base change matrix U */
231 : static GEN
232 1211 : hnf_block(GEN m, long r0, long nr, long c0, long nc)
233 : {
234 : GEN block, u, uu;
235 1211 : long nm = lg(m)-1, k;
236 1211 : pari_sp av = avma;
237 :
238 1211 : block = matslice(m, r0+1, r0+nr, c0+1, c0+nc);
239 1211 : block = Q_remove_denom(block, NULL);
240 1211 : block = rowreverse(block);
241 :
242 1211 : (void)ZM_hnfall(block, &u, 1);
243 1211 : vecreverse_inplace(u); uu = matid(nm); /* embed in matid */
244 6118 : for (k = 1; k <= nc; k++) gel(uu,c0+k) = embedcol(gel(u,k),nm,c0);
245 1211 : return gerepilecopy(av, uu);
246 : }
247 :
248 : /* (matrix, starting row, nb rows, starting column, nb columns) */
249 : static GEN
250 726 : lll_block(GEN m, long r0, long nr, long c0, long nc)
251 : {
252 : GEN block, u, uu;
253 726 : long nm = lg(m)-1, k;
254 726 : pari_sp av = avma;
255 :
256 726 : block = matslice(m, r0+1, r0+nr, c0+1, c0+nc);
257 726 : u = lll(block); if (lg(u) <= nc) return NULL;
258 714 : vecreverse_inplace(u); uu = matid(nm); /* embed in matid */
259 2464 : for (k = 1; k <= nc; k++) gel(uu,c0+k) = embedcol(gel(u,k),nm,c0);
260 714 : return gerepilecopy(av, uu);
261 : }
262 :
263 : /* insert a column at index i, no copy */
264 : static GEN
265 2323 : shallowmatinsert(GEN m, GEN x, long i)
266 : {
267 2323 : long k, n = lg(m);
268 2323 : GEN mm = cgetg(n+1,t_MAT);
269 14596 : for (k=1; k < i; k++) gel(mm, k) = gel(m, k);
270 2323 : gel(mm, i) = x;
271 3677 : for (k=i; k < n; k++) gel(mm, k+1) = gel(m, k);
272 2323 : return mm;
273 : }
274 :
275 : static GEN
276 2323 : vec_v0(long n, long n0, long r1, long r2)
277 : {
278 : long k;
279 2323 : GEN C = zerocol(n);
280 5227 : for (k = 1; k <= r1; k++) gel(C, n0++) = gen_1;
281 5644 : for (k = 1; k <= r2; k++) gel(C, n0++) = gen_2;
282 2323 : return C;
283 : }
284 :
285 : /* select cm embeddings; return a matrix */
286 : /* TODO detect if precision was insufficient */
287 : static GEN
288 224 : cm_select(GEN nf, GEN cm, long prec)
289 : {
290 : GEN emb, keys, v, m_sel, imag_emb;
291 224 : long nalg, d_cm, r_cm, c, i, j, r2 = nf_get_r2(nf);
292 : pari_sp av;
293 :
294 224 : d_cm = degpol(gel(cm, 1)); /* degree of the cm field; even */
295 224 : nalg = d_cm / 2; /* nb of clusters */
296 224 : r_cm = nf_get_degree(nf) / d_cm; /* nb by cluster; nalg * r_cm = r2 */
297 224 : m_sel = zeromatcopy(nalg, r2); /* selection matrix */
298 224 : av = avma;
299 : /* group complex embeddings */
300 224 : emb = nfeltembed(nf, gel(cm, 2), NULL, prec);
301 : /* sort */
302 224 : imag_emb = imag_i(emb);
303 224 : keys = gadd(gmul(mppi(prec), real_i(emb)), gabs(imag_emb, prec));
304 224 : v = indexsort(keys);
305 :
306 623 : for (j = c = 1; c <= nalg; c++)
307 : {
308 399 : int ref = gsigne(gel(imag_emb, v[j]));
309 399 : gcoeff(m_sel, c, v[j]) = gen_1;
310 399 : j++;
311 511 : for (i = 2; i <= r_cm; i++)
312 : {
313 112 : int s = gsigne(gel(imag_emb, v[j]));
314 112 : gcoeff(m_sel, c, v[j]) = (s == ref) ? gen_1 : gen_m1;
315 112 : j++;
316 : }
317 : }
318 224 : return gc_const(av, m_sel);
319 : }
320 :
321 : static GEN gchar_hnfreduce_shallow(GEN gc, GEN cm);
322 : static void gchar_snfbasis_shallow(GEN gc, GEN rel);
323 : static void gcharmat_tinverse(GEN gc, GEN m, long prec);
324 : static GEN gcharmatnewprec_shallow(GEN gc, long mprec);
325 :
326 : /* return a set S of prime ideals such that Cl_S(K) = 1 */
327 : static GEN
328 623 : bestS(GEN bnf)
329 : {
330 623 : GEN v, S, hw, hv = bnf_get_no(bnf), DL, dl;
331 : long i, lS;
332 : ulong l;
333 : forprime_t P;
334 :
335 623 : if (equali1(hv)) return mkvec2(cgetg(1,t_VEC), cgetg(1,t_VEC));
336 210 : v = diagonal_shallow(bnf_get_cyc(bnf));
337 210 : S = cgetg(expi(hv)+2, t_VEC); lS = 1;
338 210 : DL = cgetg(expi(hv)+2, t_VEC);
339 210 : u_forprime_init(&P,2,ULONG_MAX);
340 702 : while ((l = u_forprime_next(&P)))
341 : {
342 702 : pari_sp av = avma;
343 702 : GEN w, Sl = idealprimedec(bnf, utoi(l));
344 702 : long nSl = lg(Sl)-1;
345 1189 : for (i = 1; i < nSl; i++) /* remove one prime ideal */
346 : {
347 697 : dl = isprincipal(bnf, gel(Sl,i));
348 697 : w = ZM_hnf(shallowconcat(v, dl));
349 697 : hw = ZM_det(w);
350 697 : if (cmpii(hw, hv) < 0)
351 : {
352 378 : gel(DL,lS) = dl;
353 378 : gel(S,lS++) = gel(Sl,i);
354 378 : hv = hw; v = w; av = avma;
355 378 : if (equali1(hv)) { setlg(S, lS); setlg(DL, lS); return mkvec2(S,DL); }
356 : }
357 : }
358 492 : set_avma(av);
359 : }
360 : return NULL;/*LCOV_EXCL_LINE*/
361 : }
362 :
363 : static GEN
364 623 : gcharDLdata(GEN bnf, GEN S, GEN DL)
365 : {
366 623 : GEN M, h, Minv, Lalpha, t, dl, alpha, gen, cyc = bnf_get_cyc(bnf);
367 : long i;
368 623 : M = shallowmatconcat(DL);
369 623 : h = bnf_get_no(bnf);
370 623 : gen = bnf_get_gen(bnf);
371 :
372 : /* compute right inverse of M modulo cyc */
373 623 : M = shallowtrans(M);
374 623 : M = shallowmatconcat(mkcol2(M,diagonal_shallow(cyc)));
375 623 : Minv = matinvmod(M,h);
376 623 : Minv = vecslice(Minv,1,lg(Minv)-lg(cyc));
377 623 : Minv = shallowtrans(Minv);
378 :
379 623 : Lalpha = cgetg(lg(Minv),t_VEC);
380 952 : for (i=1; i<lg(Minv); i++)
381 : {
382 : /* gen[i] = (alpha) * prod_j S[j]^Minv[j,i] */
383 329 : t = isprincipalfact(bnf, gel(gen,i), S, gneg(gel(Minv,i)), nf_GENMAT);
384 329 : dl = gel(t, 1); alpha = gel(t, 2);
385 329 : if (!gequal0(dl)) pari_err_BUG("gcharDLdata (non-principal ideal)");
386 329 : gel(Lalpha,i) = alpha;
387 : }
388 623 : return mkvec2(Minv, Lalpha);
389 : }
390 :
391 : /* compute basis of characters; gc[1] generating family, as rows */
392 : GEN
393 644 : gcharinit(GEN bnf, GEN mod, long prec)
394 : {
395 644 : pari_sp av = avma;
396 : GEN nf, zm, zmcyc, S, DLdata, sfu, logx;
397 : GEN fa2, archp, z, C, gc, cm, cyc, rel, U, Ui, m, m_inv, m0, u0;
398 : long n, k, r1, r2, ns, nc, nu, nm, order;
399 644 : long evalprec = prec, nfprec, mprec, embprec;
400 :
401 644 : prec = evalprec + EXTRAPREC64; /* default 64 extra bits */
402 :
403 : /* note on precision:
404 :
405 : - evalprec = precision requested for evaluation
406 :
407 : - prec = precision available
408 : = (relative) precision of parameters = m_inv
409 : = (relative) precision of evaluation of small chars
410 : if no cancellation
411 :
412 : - nfprec = internal nf precision used for
413 : the embedding matrix m
414 :
415 : In the structure we store [evalprec,prec,nfprec]
416 :
417 : When evaluating chi(x) at evalprec,
418 : we need prec >= max(evalprec + exponent(chi), nfprec+exponent(x))
419 : where exponent(x) is the exponent of the number field element alpha
420 : obtained after principalisation of x.
421 :
422 : If prec is not sufficient, we call gcharnewprec.
423 :
424 : To mitigate precision increase, we always initialize the structure
425 : with 64 extra bits.
426 :
427 : Final remark: a gchar struct may have inconsistent values
428 : for prec and nfprec, for example if it has been saved and loaded at
429 : default prec : one should call gcharnewprec then.
430 : */
431 :
432 644 : if (!checkbnf_i(bnf))
433 : {
434 91 : nfprec = prec;
435 91 : bnf = bnfinit0(bnf, 1, NULL, nfprec);
436 84 : nf = shallowcopy(bnf_get_nf(bnf));
437 : }
438 : else
439 : {
440 553 : GEN fu = bnf_get_sunits(bnf);
441 553 : if (!fu) fu = bnf_get_fu(bnf); /* impose fundamental units */
442 553 : nf = shallowcopy(bnf_get_nf(bnf));
443 553 : nfprec = nf_get_prec(nf);
444 : }
445 :
446 : /* Dirichlet group + make sure mod contains archimedean places */
447 637 : mod = check_mod_factored(nf,mod,NULL,&fa2,&archp,NULL);
448 623 : sort_factor(fa2, (void*)&cmp_prime_ideal, &cmp_nodata);
449 623 : zm = localstar(nf, fa2, archp);
450 623 : zmcyc = locs_get_cyc(zm);
451 :
452 : /* set of primes S and valuations of generators */
453 623 : S = bestS(bnf);
454 623 : DLdata = gel(S,2);
455 623 : S = gel(S,1);
456 623 : DLdata = gcharDLdata(bnf, S, DLdata);
457 :
458 623 : nf_get_sign(nf, &r1, &r2);
459 623 : n = r1+2*r2;
460 623 : ns = lg(S) - 1;
461 623 : nu = r1+r2-1 + ns;
462 623 : nc = lg(zmcyc) - 1;
463 623 : nm = ns+nc+n; /* number of parameters = ns + nc + r1 + r2 + r2 */
464 :
465 : /* units and S-units */
466 623 : sfu = gel(bnfunits(bnf,S), 1);
467 623 : sfu = vec_shorten(sfu, nu); /* remove torsion */
468 :
469 : /* root of unity */
470 623 : order = bnf_get_tuN(bnf);
471 623 : z = bnf_get_tuU(bnf);
472 :
473 : /* maximal cm subfield */
474 623 : cm = nfsubfieldscm(nf, 0);
475 :
476 : /*
477 : Now compute matrix of parameters,
478 : until we obtain the right precision
479 : FIXME: make sure generators, units, and discrete logs
480 : do not depend on precision.
481 :
482 : m0 is the matrix of units embeddings
483 : u is the HNF base change, m = m0*u
484 :
485 : subsequent steps may lead to precision increase, we put everything in gc
486 : struct and modify it in place.
487 :
488 : A) sets m0
489 : B) sets U, cyc, rel, U and Ui
490 : C) sets m_inv
491 : */
492 :
493 : /* A) make big matrix m0 of embeddings of units */
494 :
495 623 : if (DEBUGLEVEL>2) err_printf("start matrix m\n");
496 623 : m = cgetg(nm + 1, t_MAT);
497 623 : mprec = nbits2prec(nm+10) + EXTRAPREC64;
498 623 : embprec = mprec;
499 : for(;;)
500 : {
501 1904 : for (k = 1; k <= nu; k++)
502 : { /* Lambda_S (S-units) then Lambda_f, fund. units */
503 1281 : logx = gchar_nflog(&nf, zm, S, gel(sfu,k), embprec);
504 1281 : if (!logx) break;
505 1281 : gel(m, k) = logx;
506 : }
507 623 : if (k > nu) break;
508 0 : if (DEBUGLEVEL) err_printf("gcharinit: increasing embprec %d -> %d\n",
509 : embprec, precdbl(embprec));
510 0 : embprec = precdbl(embprec);
511 : }
512 1645 : for (k = 1; k <= nc; k++) /* Gamma, structure of (Z/m)* */
513 : {
514 1022 : C = zerocol(nm);
515 1022 : gel(C, ns+k) = gel(zmcyc, k);
516 1022 : gel(m, nu+k) = C;
517 : }
518 : /* zeta, root of unity */
519 623 : gel(m, nu+nc+1) = gchar_nflog(&nf, zm, S, z, mprec);
520 623 : shallow_clean_rat(gel(m, nu+nc+1), 1, nm, stoi(order), mprec);
521 1470 : for (k = 1; k <= r2; k++) /* embed Z^r_2 */
522 : {
523 847 : C = zerocol(nm);
524 847 : gel(C, ns+nc+r1+r2+k) = gen_1;
525 847 : gel(m, nu+nc+1+k) = C;
526 : }
527 623 : if (DEBUGLEVEL>1) err_printf("matrix m = %Ps\n", m);
528 :
529 623 : m0 = m;
530 623 : u0 = gen_0;
531 623 : rel = U = Ui = gen_0;
532 623 : cyc = gen_0;
533 623 : m_inv = gen_0;
534 :
535 : /* only m and m_inv depend on prec and are recomputed under gcharnewprec. */
536 623 : gc = mkvecn(GC_LENGTH,
537 : m_inv, /* internal basis, characters as rows */
538 : bnf,
539 : nf,
540 : zm, /* Zk/mod, nc components */
541 : S, /* generators of clgp, ns components */
542 : DLdata,
543 : sfu,
544 : mkvec2(mkvecsmall3(evalprec,prec,nfprec),
545 : mkvecsmall3(0,0,0)), /* ntors, nfree, nalg */
546 : cyc, /* reduced components */
547 : mkvec3(rel, U, Ui), /* internal / SNF base change */
548 : m0, /* embeddings of units */
549 : u0); /* m_inv = (m0 u0)~^-1 */
550 :
551 : /* B) do HNF reductions + LLL (may increase precision) */
552 623 : m = gchar_hnfreduce_shallow(gc, cm);
553 :
554 : /* C) compute snf basis of torsion subgroup */
555 623 : rel = shallowtrans(matslice(m, 1, ns+nc, 1, ns+nc));
556 623 : gchar_snfbasis_shallow(gc, rel);
557 :
558 : /* D) transpose inverse m_inv = (m0*u)~^-1 (may increase precision) */
559 623 : gcharmat_tinverse(gc, m, prec);
560 623 : return gerepilecopy(av, gc);
561 : }
562 :
563 : /* b) do HNF reductions + LLL, keep base change u0 */
564 : static GEN
565 623 : gchar_hnfreduce_shallow(GEN gc, GEN cm)
566 : {
567 623 : GEN bnf = gchar_get_bnf(gc), nf = gchar_get_nf(gc), u, u0, m;
568 623 : long order, r1, r2, ns, nc, n, nu, nm, nalg = 0, mprec;
569 :
570 623 : nf_get_sign(nf, &r1, &r2);
571 623 : n = r1 + 2*r2;
572 623 : nu = r1 + r2 - 1;
573 623 : ns = gchar_get_ns(gc);
574 623 : nc = gchar_get_nc(gc);
575 623 : nm = ns+nc+n; /* ns + nc + r1 + r2 + r2 */
576 623 : order = 2*bnf_get_tuN(bnf);
577 623 : u0 = matid(nm);
578 623 : m = shallowcopy(gchar_get_m0(gc)); /* keep m0 unchanged */
579 623 : mprec = gprecision(m);
580 623 : if (DEBUGLEVEL>1) err_printf("matrix m = %Ps\n", m);
581 623 : if (nc)
582 : { /* keep steps 1&2 to make sure we have zeta_m */
583 308 : u = hnf_block(m, ns,nc, ns+nu,nc+1);
584 308 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
585 308 : if (DEBUGLEVEL>2) err_printf("step 1 -> %Ps\n", m);
586 308 : u = hnf_block(m, ns,nc, ns,nu+nc);
587 308 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
588 308 : if (DEBUGLEVEL>2) err_printf("step 2 -> %Ps\n", m);
589 : }
590 623 : if (r2)
591 : {
592 371 : u = hnf_block(m, nm-r2,r2, nm-r2-1,r2+1);
593 371 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
594 371 : if (DEBUGLEVEL>2) err_printf("step 3 -> %Ps\n", m);
595 : }
596 : /* remove last column */
597 623 : setlg(u0, nm); setlg(m, nm);
598 623 : if (DEBUGLEVEL>2) err_printf("remove last col -> %Ps\n", m);
599 :
600 623 : if (!gequal0(cm))
601 : {
602 : GEN v, Nargs;
603 224 : long bit = - prec2nbits(mprec) + 16 + expu(order);
604 : /* reduce on Norm arguments */
605 224 : v = cm_select(nf, cm, gchar_get_nfprec(gc));
606 224 : if (DEBUGLEVEL>2) err_printf("cm_select -> %Ps\n", v);
607 224 : nalg = nbrows(v);
608 224 : gchar_set_u0(gc, u0);
609 : for(;;)
610 14 : {
611 : long e, emax, i;
612 238 : Nargs = gmul(v, rowslice(m, nm-r2+1, nm));
613 238 : if (DEBUGLEVEL>2) err_printf("Nargs -> %Ps\n", Nargs);
614 238 : emax = bit-1;
615 1120 : for (i = ns+nc+1; i < lg(Nargs); i++)
616 : {
617 882 : gel(Nargs,i) = grndtoi(gmulgs(gel(Nargs,i), order), &e);
618 882 : emax = maxss(emax,e);
619 : }
620 238 : if (emax < bit) break;
621 14 : if (DEBUGLEVEL>1) err_printf("cm select: doubling prec\n");
622 14 : mprec = precdbl(mprec); m = gcharmatnewprec_shallow(gc, mprec);
623 : }
624 224 : if (DEBUGLEVEL>2) err_printf("rounded Nargs -> %Ps\n", Nargs);
625 224 : u = hnf_block(Nargs, 0, nalg, ns+nc, n-1);
626 224 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
627 224 : if (DEBUGLEVEL>2) err_printf("after cm reduction -> %Ps\n", m);
628 : }
629 :
630 : /* apply LLL on Lambda_m, may need to increase prec */
631 623 : gchar_set_nalg(gc, nalg);
632 623 : gchar_set_u0(gc, u0);
633 :
634 : /* TODO factor these two LLL reduction codes in a function? */
635 623 : if (nalg > 0)
636 : {
637 224 : GEN u = NULL;
638 : while (1)
639 : {
640 224 : if ((u = lll_block(m, ns+nc, n, ns+nc, nalg))) break;
641 0 : mprec = precdbl(mprec); m = gcharmatnewprec_shallow(gc, mprec);
642 : }
643 224 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
644 224 : if (DEBUGLEVEL>1) err_printf("after LLL reduction (CM block) -> %Ps\n", m);
645 : }
646 623 : gchar_set_u0(gc, u0);
647 :
648 623 : if (nu > 0)
649 : {
650 490 : GEN u = NULL;
651 : while (1)
652 : {
653 502 : if ((u = lll_block(m, ns+nc, n, ns+nc+nalg, n-1-nalg))) break;
654 12 : mprec = precdbl(mprec); m = gcharmatnewprec_shallow(gc, mprec);
655 : }
656 490 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
657 490 : if (DEBUGLEVEL>1) err_printf("after LLL reduction (trans block) -> %Ps\n", m);
658 : }
659 623 : gchar_set_u0(gc, u0); return m;
660 : }
661 :
662 : /* convert to snf basis of torsion + Z^(r1+2*r2-1) */
663 : static void
664 623 : gchar_snfbasis_shallow(GEN gc, GEN rel)
665 : {
666 : long n, r1, r2, lU, lUi;
667 : GEN nf, cyc, U, Ui;
668 :
669 623 : nf = gchar_get_nf(gc);
670 623 : nf_get_sign(nf, &r1, &r2);
671 623 : n = r1+2*r2;
672 :
673 623 : rel = ZM_hnf(rel);
674 623 : if (DEBUGLEVEL>1) err_printf("relations after hnf: %Ps\n", rel);
675 623 : cyc = ZM_snf_group(rel, &U, &Ui);
676 623 : if (lg(cyc)==1)
677 : {
678 280 : cyc = zerovec(n-1);
679 280 : U = shallowconcat(zeromat(n-1,lg(rel)-1),matid(n-1));
680 280 : Ui = shallowtrans(U);
681 : }
682 343 : else if (n!=1)
683 : {
684 322 : cyc = shallowconcat(cyc, zerovec(n-1));
685 322 : U = shallowmatconcat(mkmat22(U,gen_0,gen_0,matid(n-1)));
686 322 : Ui = shallowmatconcat(mkmat22(Ui,gen_0,gen_0,matid(n-1)));
687 : }
688 :
689 623 : rel = shallowtrans(rel);
690 623 : lU = lg(U);
691 623 : lUi = lg(Ui);
692 623 : U = shallowtrans(U);
693 623 : if (lg(U)==1 && lUi>1) U = zeromat(0,lUi-1);
694 623 : Ui = shallowtrans(Ui);
695 623 : if (lg(Ui)==1 && lU>1) Ui = zeromat(0,lU-1);
696 :
697 623 : gchar_set_nfree(gc, n-1);
698 623 : gchar_set_ntors(gc, (lg(cyc)-1) - (n-1));
699 623 : gchar_set_cyc(gc, shallowconcat(cyc, real_0(gchar_get_prec(gc))));
700 623 : gchar_set_HUUi(gc, rel, U, Ui);
701 623 : }
702 :
703 : static long
704 2519 : mextraprec(GEN m_inv, GEN m, GEN gc)
705 : {
706 2519 : return nbits2extraprec(2*maxss(gexpo(m_inv),1) + expu(lg(m))
707 2519 : + gexpo(gchar_get_u0(gc)) + 10);
708 : }
709 :
710 : /* c) transpose inverse + clean rationals.
711 : * prec = target prec for m^-1,
712 : * mprec = prec of m */
713 : static void
714 1573 : gcharmat_tinverse(GEN gc, GEN m, long prec)
715 : {
716 : GEN m_inv;
717 1573 : long k, r1, r2, ns, nc, nalg, nm, mprec, bitprec = prec2nbits(prec);
718 :
719 1573 : nf_get_sign(gchar_get_nf(gc), &r1, &r2);
720 1573 : ns = gchar_get_ns(gc);
721 1573 : nc = gchar_get_nc(gc);
722 1573 : nalg = gchar_get_nalg(gc);
723 1573 : nm = ns + nc + r1 + 2*r2;
724 1573 : if (lg(m)==1) { gchar_set_basis(gc,zeromat(0,nm)); return; }
725 :
726 1552 : mprec = gprecision(m); /* possibly 0, if m is exact */
727 : while (1)
728 771 : {
729 2323 : long targetmprec = 0;
730 : GEN v0, mm;
731 : /* insert at column ns+nc+r1+r2, or last column if cm */
732 2323 : v0 = vec_v0(nm, ns+nc+1, r1, r2);
733 2323 : mm = shallowmatinsert(m, v0, nalg? nm: nm-r2);
734 2323 : if (DEBUGLEVEL>1) err_printf("add v0 -> %Ps\n", mm);
735 2323 : mm = shallowtrans(mm);
736 2323 : m_inv = RgM_inv(mm); /* invert matrix, may need to increase prec */
737 2323 : if (m_inv)
738 : {
739 2323 : if (DEBUGLEVEL>1) err_printf("inverse: %Ps\n",m_inv);
740 2323 : m_inv = vecsplice(m_inv, nalg? nm: nm-r2); /* remove v0 */
741 2323 : if (DEBUGLEVEL>1) err_printf("v0 removed: %Ps\n", m_inv);
742 2323 : m_inv = shallowtrans(m_inv);
743 2323 : if (!mprec) break;
744 : /* enough precision? */
745 : /* |B - A^(-1)| << |B|.|Id-B*A| */
746 2246 : if (gexpo(m_inv) + gexpo(gsub(RgM_mul(m_inv, m), gen_1)) + expu(lg(m))
747 2246 : <= -bitprec)
748 : {
749 : /* |A^(-1) - (A+H)^(-1)| << |H|.|A^(-1)|^2 */
750 1569 : targetmprec = prec + mextraprec(m_inv,m,gc);
751 1569 : if (mprec >= targetmprec) break;
752 : }
753 677 : else targetmprec = 0;
754 : }
755 771 : mprec = maxss(precdbl(mprec), targetmprec);
756 771 : if (mprec < DEFAULTPREC) mprec = DEFAULTPREC;
757 771 : m = gcharmatnewprec_shallow(gc, mprec); /* m0 * u0 to higher prec */
758 : }
759 : /* clean rationals */
760 1552 : if (nc)
761 : { /* reduce mod exponent of the group, TODO reduce on each component */
762 985 : GEN zmcyc = locs_get_cyc(gchar_get_zm(gc));
763 985 : GEN e = ZV_lcm(zmcyc);
764 3714 : for (k = 1; k <= nc; k++)
765 2729 : shallow_clean_rat(gel(m_inv, ns+k), 1, nm - 1, /*zmcyc[k]*/e, prec);
766 : }
767 1552 : if (r2)
768 : {
769 2736 : for (k = 1; k <= r2; k++)
770 1698 : shallow_clean_rat(gel(m_inv,nm-k+1), 1, nm-1, NULL, prec);
771 : }
772 1552 : if (nalg)
773 : {
774 : long i, j, e;
775 1637 : for (i = 1; i<=r1+r2; i++)
776 3176 : for (j = 1; j <= nalg; j++)
777 : {
778 2122 : e = gexpo(gcoeff(m_inv, ns+nc+j, ns+nc+i));
779 2122 : if (e > -20) /* TODO is this bound too permissive? */
780 : pari_err_BUG("gcharinit (nonzero entry)"); /*LCOV_EXCL_LINE*/
781 2122 : gcoeff(m_inv, ns+nc+j, ns+nc+i) = gen_0;
782 : }
783 : }
784 1552 : if (DEBUGLEVEL>1) err_printf("cyc/cm cleaned: %Ps", shallowtrans(m_inv));
785 : /* normalize characters, parameters mod Z */
786 5212 : for (k = 1; k <= ns+nc; k++) gel(m_inv, k) = gfrac(gel(m_inv, k));
787 : /* increase relative prec of real values */
788 1552 : gchar_set_basis(gc, gprec_w(m_inv, prec));
789 : }
790 :
791 : /* make sure same argument determination is chosen */
792 : static void
793 4551 : same_arg(GEN v1, GEN v2, long s1, long s2)
794 : {
795 4551 : long i1, i2, l = lg(v1);
796 22987 : for (i1 = s1, i2 = s2; i1 < l; i1++,i2++)
797 : {
798 18436 : GEN d = grndtoi(gsub(gel(v2,i2),gel(v1,i1)), NULL);
799 18436 : if (signe(d)) gel(v1,i1) = gadd(gel(v1,i1), d);
800 : }
801 4551 : }
802 :
803 : static void
804 4551 : vaffect_shallow(GEN x, long i0, GEN y)
805 : {
806 : long i;
807 39881 : for (i = 1; i < lg(y); i++)
808 35330 : gel(x, i+i0) = gel(y, i);
809 4551 : }
810 :
811 : /* recompute matrix with precision increased */
812 : /* u0 the base change, returns m0 * u0 */
813 : /* mprec: requested precision for m0 */
814 : static GEN
815 1747 : gcharmatnewprec_shallow(GEN gc, long mprec)
816 : {
817 : GEN nf, m0, u0, sfu;
818 : long k, l, ns, nc, r1, r2, nfprec, embprec;
819 :
820 1747 : ns = gchar_get_ns(gc);
821 1747 : nc = gchar_get_nc(gc);
822 1747 : nf = gchar_get_nf(gc);
823 1747 : sfu = gchar_get_sfu(gc);
824 1747 : nf_get_sign(nf, &r1, &r2);
825 1747 : nfprec = nf_get_prec(gchar_get_nf(gc));
826 :
827 1747 : m0 = gchar_get_m0(gc);
828 1747 : u0 = gchar_get_u0(gc);
829 :
830 1747 : if (DEBUGLEVEL>3) err_printf("gcharmatnewprec_shallow mprec=%d nfprec=%d\n", mprec, nfprec);
831 :
832 1747 : embprec = mprec; l = lg(sfu);
833 : while(1)
834 : { /* recompute the nfembedlogs of s-units and fundamental units */
835 6298 : for (k = 1; k < l; k++) /* ns s-units, then r1+r2-1 fundamental units */
836 : {
837 4551 : GEN emb = nfembedlog(&nf, gel(sfu,k), embprec);
838 4551 : if (!emb) break;
839 4551 : same_arg(emb, gel(m0,k), r1+r2, ns+nc+r1+r2);
840 4551 : vaffect_shallow(gel(m0, k), ns+nc, emb);
841 : }
842 1747 : if (k == l) break;
843 0 : if (DEBUGLEVEL)
844 0 : err_printf("gcharmatnewprec_shallow: increasing embprec %d -> %d\n",
845 : embprec, precdbl(embprec));
846 0 : embprec = precdbl(embprec);
847 : }
848 1747 : gchar_set_nf(gc, nf);
849 1747 : gchar_set_nfprec(gc, nfprec);
850 1747 : return RgM_ZM_mul(m0, u0);
851 : }
852 :
853 : static void _check_gchar_group(GEN gc, long flag);
854 : static void
855 11431 : check_gchar_group(GEN gc) { _check_gchar_group(gc, 0); }
856 :
857 : /* increase prec if needed. newprec: requested precision for m_inv */
858 : static GEN
859 1806 : gcharnewprec_i(GEN gc, long newprec)
860 : {
861 : long prec, prec0, nfprec, nfprec0, mprec;
862 1806 : GEN gc2 = shallowcopy(gc);
863 :
864 1806 : _check_gchar_group(gc2, 1); /* ignore illegal prec */
865 1785 : prec = gchar_get_prec(gc2);
866 1785 : nfprec = gchar_get_nfprec(gc2);
867 :
868 1785 : if (newprec > prec)
869 : { /* increase precision */
870 1202 : if (DEBUGLEVEL) pari_warn(warnprec,"gcharnewprec",newprec);
871 1202 : nfprec += newprec - prec;
872 1202 : prec = newprec;
873 1202 : gchar_copy_precs(gc, gc2);
874 1202 : gchar_set_prec(gc2, prec);
875 1202 : gchar_set_nfprec(gc2, nfprec);
876 : }
877 :
878 1785 : nfprec0 = nf_get_prec(gchar_get_nf(gc2));
879 1785 : if (nfprec0 && nfprec > nfprec0)
880 : {
881 329 : if (DEBUGLEVEL) pari_warn(warnprec,"gcharnewprec (nf)", nfprec);
882 329 : gchar_set_nf(gc2, nfnewprec_shallow(gchar_get_nf(gc2), nfprec));
883 : }
884 :
885 1785 : prec0 = gprecision(gchar_get_basis(gc2));
886 1785 : if (prec0 && prec > prec0)
887 : {
888 950 : GEN cyc, m, m0 = gchar_get_m0(gc);
889 950 : if (DEBUGLEVEL) pari_warn(warnprec,"gcharnewprec (minv)", prec);
890 950 : gchar_set_m0(gc2, shallowcopy(m0));
891 950 : mprec = prec + mextraprec(gchar_get_basis(gc), m0, gc);
892 950 : m = gcharmatnewprec_shallow(gc2, mprec);
893 950 : if (DEBUGLEVEL>2) err_printf("m0*u0 recomputed -> %Ps\n", m);
894 950 : gcharmat_tinverse(gc2, m, prec);
895 950 : cyc = shallowcopy(gchar_get_cyc(gc2));
896 950 : gel(cyc, lg(cyc)-1) = real_0(prec);
897 950 : gchar_set_cyc(gc2, cyc);
898 : }
899 1785 : return gc2;
900 : }
901 :
902 : /* newprec: requested evalprec */
903 : GEN
904 1806 : gcharnewprec(GEN gc, long newprec)
905 : {
906 1806 : pari_sp av = avma;
907 1806 : GEN gc2 = gcharnewprec_i(gc, newprec + EXTRAPREC64);
908 1785 : gchar_set_evalprec(gc2, newprec);
909 1785 : return gerepilecopy(av, gc2);
910 : }
911 :
912 : static void
913 13223 : check_localstar(GEN x)
914 : {
915 13223 : if (typ(x) != t_VEC || lg(x) != LOCS_LENGTH + 1)
916 7 : pari_err_TYPE("char group", x);
917 : /* FIXME: check components once settled */
918 13216 : }
919 :
920 : int
921 8876 : is_gchar_group(GEN gc)
922 : {
923 8876 : return (typ(gc) == t_VEC
924 8876 : && lg(gc) == GC_LENGTH + 1
925 392 : && typ(gel(gc, 8)) == t_VEC
926 392 : && lg(gel(gc, 8)) == 3
927 392 : && typ(gmael(gc, 8, 1)) == t_VECSMALL
928 392 : && typ(gmael(gc, 8, 2)) == t_VECSMALL
929 392 : && (checkbnf_i(gchar_get_bnf(gc)) != NULL)
930 17752 : && (checknf_i(gchar_get_nf(gc)) != NULL));
931 : }
932 :
933 : /* validates the structure format.
934 : * Raise error if inconsistent precision, unless flag=1. */
935 : static void
936 13237 : _check_gchar_group(GEN gc, long flag)
937 : {
938 : GEN b, bnf, nf;
939 13237 : if (typ(gc) != t_VEC || lg(gc) != GC_LENGTH + 1)
940 14 : pari_err_TYPE("char group", gc);
941 13223 : check_localstar(gchar_get_zm(gc));
942 13216 : if (typ(gchar_get_loccyc(gc)) != t_VEC)
943 0 : pari_err_TYPE("gchar group (loccyc)", gc);
944 13216 : b = gchar_get_basis(gc);
945 13216 : if (typ(b) != t_MAT) pari_err_TYPE("gchar group (basis)", gc);
946 13209 : bnf = gchar_get_bnf(gc); checkbnf(bnf);
947 13209 : nf = gchar_get_nf(gc); checknf(nf);
948 13209 : if (!gequal(nf_get_pol(nf),nf_get_pol(bnf_get_nf(bnf))))
949 0 : pari_err_TYPE("gchar group (nf != bnf.nf)", gc);
950 13209 : if (typ(gel(gc,8)) != t_VEC ||typ(gmael(gc,8,1)) != t_VECSMALL)
951 7 : pari_err_TYPE("gchar group (gc[8])", gc);
952 13202 : if (!flag)
953 : {
954 11417 : long prec0 = gprecision(b), nfprec0 = nf_get_prec(nf);
955 11417 : if ((prec0 && gchar_get_prec(gc) > prec0)
956 11410 : || (nfprec0 && gchar_get_nfprec(gc) > nfprec0))
957 7 : pari_err_PREC("gchar group, please call gcharnewprec");
958 : }
959 13195 : }
960 :
961 : /* basis of algebraic characters + cyc vector */
962 : static GEN
963 70 : gchar_algebraic_basis(GEN gc)
964 : {
965 : long ntors, nfree, nc, nm, r2, nalg, n0, k;
966 : GEN basis, args, m, w, normchar, alg_basis, tors_basis;
967 :
968 : /* in snf basis */
969 70 : ntors = gchar_get_ntors(gc); /* number of torsion generators */
970 70 : nfree = gchar_get_nfree(gc);
971 70 : nc = ntors + nfree;
972 : /* in internal basis */
973 70 : n0 = gchar_get_ns(gc) + gchar_get_nc(gc); /* last index of torsion char */
974 70 : r2 = gchar_get_r2(gc);
975 70 : nm = gchar_get_nm(gc);
976 : /* in both */
977 70 : nalg = gchar_get_nalg(gc); /* number of generators of free algebraic chars */
978 :
979 : /* finite order characters have weight 0 */
980 70 : tors_basis = matid(ntors);
981 :
982 : /* the norm is an algebraic character */
983 70 : normchar = gneg(col_ei(nc+1,nc+1));
984 :
985 : /* add sublattice of algebraic */
986 :
987 70 : if (!nalg)
988 : {
989 28 : if (DEBUGLEVEL>2) err_printf("nalg=0\n");
990 28 : return shallowmatconcat(mkvec2(tors_basis,normchar));
991 : }
992 :
993 : /* block of k_s parameters of free algebraic */
994 42 : args = matslice(gchar_get_basis(gc),n0+1,n0+nalg,nm-r2+1,nm);
995 :
996 42 : if (r2 == 1)
997 : {
998 : /* no parity condition */
999 14 : if (DEBUGLEVEL>2) err_printf("r2 = 1 -> args = %Ps\n", args);
1000 14 : alg_basis = matid(nalg);
1001 14 : w = gel(args,1);
1002 : }
1003 : else
1004 : {
1005 : /* parity condition: w + k_s = 0 mod 2 for all s,
1006 : ie solve x.K constant mod 2, ie solve
1007 : x.K' = 0 mod 2 where K' = [ C-C0 ] (substract first column)
1008 : */
1009 : /* select block k_s in char parameters and */
1010 28 : if (DEBUGLEVEL>2) err_printf("block ks -> %Ps\n", args);
1011 28 : m = cgetg(r2, t_MAT);
1012 77 : for (k = 1; k < r2; k++)
1013 49 : gel(m,k) = gsub(gel(args,k+1),gel(args,1));
1014 28 : if (DEBUGLEVEL>2) err_printf("block ks' -> %Ps", m);
1015 28 : alg_basis = shallowtrans(gel(matsolvemod(shallowtrans(m),gen_2,gen_0,1),2));
1016 28 : if (DEBUGLEVEL>2) err_printf("alg_basis -> %Ps\n", alg_basis);
1017 28 : w = gmul(alg_basis, gel(args,1));
1018 28 : if (DEBUGLEVEL>2) err_printf("w -> %Ps\n", w);
1019 : }
1020 : /* add weight to infinite order characters, at position nc+1 */
1021 42 : w = gneg(gdivgs(gmodgs(w, 2), 2));
1022 42 : alg_basis = shallowmatconcat(mkvec3(alg_basis, zerovec(nfree-nalg), w));
1023 :
1024 42 : if (lg(tors_basis)>1)
1025 28 : basis = shallowmatconcat(mkmat22(tors_basis, gen_0, gen_0, alg_basis));
1026 : else
1027 14 : basis = alg_basis;
1028 42 : return shallowmatconcat(mkvec2(shallowtrans(basis),normchar));
1029 : }
1030 : static GEN
1031 154 : gchar_algebraicnormtype(GEN gc, GEN type)
1032 : {
1033 154 : GEN w = NULL, w1, v;
1034 : long i, r1, r2, n;
1035 154 : r1 = gchar_get_r1(gc);
1036 154 : r2 = gchar_get_r2(gc);
1037 210 : for (i=1; i<=r1; i++)
1038 : {
1039 56 : w1 = addii(gmael(type,i,1),gmael(type,i,2));
1040 56 : if (!w) w = w1;
1041 14 : else if (!equalii(w,w1)) return NULL;
1042 : }
1043 203 : for (i=r1+1; i<=r1+r2; i++)
1044 : {
1045 154 : w1 = gmael(type,i,1);
1046 154 : if (!w) w = w1;
1047 42 : else if (!equalii(w,w1)) return NULL;
1048 140 : if (!equalii(w,gmael(type,i,2))) return NULL;
1049 : }
1050 49 : n = lg(gchar_get_cyc(gc))-1;
1051 49 : v = zerocol(n);
1052 49 : gel(v,n) = negi(w);
1053 49 : return mkvec(v);
1054 : }
1055 : static GEN
1056 210 : gchar_algebraicoftype(GEN gc, GEN type)
1057 : {
1058 : long i, r1, r2, nalg, n0, nm;
1059 : GEN p, q, w, k, matk, chi;
1060 : /* in internal basis */
1061 210 : n0 = gchar_get_ns(gc) + gchar_get_nc(gc); /* last index of torsion chars */
1062 210 : r1 = gchar_get_r1(gc);
1063 210 : r2 = gchar_get_r2(gc);
1064 210 : nm = gchar_get_nm(gc);
1065 : /* in both */
1066 210 : nalg = gchar_get_nalg(gc); /* number of generators of free algebraic chars */
1067 :
1068 210 : if (typ(type)!=t_VEC || lg(type) != r1+r2+1)
1069 28 : pari_err_TYPE("gcharalgebraic (1)", type);
1070 448 : for (i = 1; i<=r1+r2; i++)
1071 294 : if (typ(gel(type,i)) != t_VEC
1072 287 : ||lg(gel(type,i)) != 3
1073 280 : ||typ(gmael(type,i,1)) != t_INT
1074 273 : ||typ(gmael(type,i,2)) != t_INT)
1075 28 : pari_err_TYPE("gcharalgebraic (2)", type);
1076 :
1077 154 : chi = gchar_algebraicnormtype(gc, type);
1078 154 : if (chi) return chi;
1079 :
1080 105 : if (!nalg) return NULL;
1081 91 : k = cgetg(r2+1,t_VEC);
1082 91 : p = gmael(type, 1, 1);
1083 91 : q = gmael(type, 1, 2); w = addii(p, q);
1084 91 : gel(k, 1) = subii(q, p);
1085 140 : for (i = 2; i <= r2; i++)
1086 : {
1087 63 : p = gmael(type, i, 1);
1088 63 : q = gmael(type, i, 2);
1089 63 : if (!equalii(w, addii(p, q))) return NULL;
1090 49 : gel(k, i) = subii(q, p);
1091 : }
1092 : /* block of k_s parameters of free algebraic */
1093 77 : matk = matslice(gchar_get_basis(gc),n0+1,n0+nalg,nm-r2+1,nm);
1094 77 : chi = inverseimage(shallowtrans(matk),shallowtrans(k));
1095 77 : if (lg(chi) == 1) return NULL;
1096 182 : for (i=1; i<lg(chi); i++) if (typ(gel(chi,i)) != t_INT) return NULL;
1097 126 : chi = mkvec4(zerocol(gchar_get_ntors(gc)), chi,
1098 63 : zerocol(gchar_get_nfree(gc) - nalg), gmul2n(negi(w),-1));
1099 63 : return mkvec(shallowconcat1(chi));
1100 : }
1101 :
1102 : GEN
1103 280 : gcharalgebraic(GEN gc, GEN type)
1104 : {
1105 280 : pari_sp av = avma;
1106 : GEN b;
1107 280 : check_gchar_group(gc);
1108 280 : b = type? gchar_algebraicoftype(gc, type): gchar_algebraic_basis(gc);
1109 224 : if (!b) { set_avma(av); return cgetg(1, t_VEC); }
1110 182 : return gerepilecopy(av, b);
1111 : }
1112 :
1113 : /*********************************************************************/
1114 : /* */
1115 : /* CHARACTERS */
1116 : /* */
1117 : /*********************************************************************/
1118 : /* return chi without (optional) norm component; set *s = to the latter */
1119 : static GEN
1120 9128 : check_gchar_i(GEN chi, long l, GEN *s)
1121 : {
1122 9128 : long i, n = lg(chi)-1;
1123 9128 : if (n == l)
1124 : { /* norm component */
1125 2527 : if (s)
1126 : {
1127 2506 : *s = gel(chi,l);
1128 2506 : switch(typ(*s))
1129 : {
1130 2499 : case t_INT:
1131 : case t_FRAC:
1132 : case t_REAL:
1133 2499 : case t_COMPLEX: break;
1134 7 : default: pari_err_TYPE("check_gchar_i [s]", *s);
1135 : }
1136 : }
1137 2520 : chi = vecslice(chi, 1, l-1);
1138 : }
1139 : else
1140 : { /* missing norm component */
1141 6601 : if (n != l-1) pari_err_DIM("check_gchar_i [chi]");
1142 6594 : if (s) *s = gen_0;
1143 : }
1144 44912 : for (i = 1; i < l; i++) if (typ(gel(chi,i))!=t_INT)
1145 7 : pari_err_TYPE("check_gchar_i [coefficient]", gel(chi,i));
1146 9107 : return chi;
1147 : }
1148 :
1149 : static GEN
1150 7070 : check_gchar(GEN gc, GEN chi, GEN *s)
1151 : {
1152 7070 : if (typ(chi)!=t_COL) pari_err_TYPE("check_gchar [chi]", chi);
1153 7049 : return check_gchar_i(chi, lg(gchar_get_cyc(gc))-1, s);
1154 : }
1155 :
1156 : static long
1157 2079 : safelgcols(GEN m)
1158 : {
1159 2079 : return lg(m)==1 ? 1 : lg(gel(m,1));
1160 : }
1161 :
1162 : static GEN
1163 2079 : check_gchari(GEN gc, GEN chi, GEN *s)
1164 : {
1165 2079 : if (typ(chi)!=t_VEC) pari_err_TYPE("check_gchari [chi]", chi);
1166 2079 : return check_gchar_i(chi, safelgcols(gchar_get_basis(gc)), s);
1167 : }
1168 :
1169 : /* from coordinates on snf basis, return coordinates on internal basis, set
1170 : * s to the norm component */
1171 : static GEN
1172 7070 : gchar_internal(GEN gc, GEN chi, GEN *s)
1173 : {
1174 7070 : chi = check_gchar(gc, chi, s);
1175 7028 : return ZV_ZM_mul(chi, gchar_get_Ui(gc));
1176 : }
1177 :
1178 : static GEN
1179 23982 : RgV_frac_inplace(GEN v, long n)
1180 : {
1181 : long i;
1182 53543 : for (i = 1; i <= n; i++) gel(v,i) = gfrac(gel(v,i));
1183 23982 : return v;
1184 : }
1185 : /* from internal basis form, return the R/Z components */
1186 : static GEN
1187 8708 : gchari_duallog(GEN gc, GEN chi)
1188 : {
1189 8708 : chi = RgV_RgM_mul(chi, gchar_get_basis(gc)); /* take exponents mod Z */
1190 8708 : return RgV_frac_inplace(chi, gchar_get_ns(gc) + gchar_get_nc(gc));
1191 : }
1192 :
1193 : /* chip has ncp = #zm[1][i].cyc components */
1194 : static GEN
1195 1330 : conductor_expo_pr(GEN gens_fil, GEN chip)
1196 : {
1197 : long i;
1198 2317 : for (i = lg(gens_fil) - 1; i > 0; i--)
1199 : {
1200 1932 : GEN gens = gel(gens_fil, i);
1201 : long j;
1202 3199 : for (j = 1; j < lg(gens); j++)
1203 : {
1204 2212 : GEN v = gmul(chip, gel(gens, j));
1205 2212 : if (denom_i(v) != gen_1) return stoi(i);
1206 : }
1207 : }
1208 385 : return gen_0;
1209 : }
1210 :
1211 : /* assume chi in log form */
1212 : static GEN
1213 1589 : gcharlog_conductor_f(GEN gc, GEN chi, GEN *faN)
1214 : {
1215 : GEN zm, P, E, Lsprk, ufil;
1216 : long i, l, ic;
1217 :
1218 1589 : if (gchar_get_nc(gc) == 0)
1219 : {
1220 329 : if (faN) *faN = trivial_fact();
1221 329 : return gen_1;
1222 : }
1223 1260 : zm = gchar_get_zm(gc);
1224 1260 : Lsprk = locs_get_Lsprk(zm);
1225 1260 : ufil = locs_get_Lgenfil(zm);
1226 1260 : P = gel(locs_get_famod(zm), 1);
1227 1260 : l = lg(Lsprk); E = cgetg(l, t_COL);
1228 2590 : for (i = 1, ic = gchar_get_ns(gc); i < l ; i++)
1229 : {
1230 1330 : GEN gens = gel(ufil, i), chip;
1231 1330 : long ncp = sprk_get_ncp(gel(Lsprk,i));
1232 :
1233 1330 : chip = vecslice(chi, ic + 1, ic + ncp);
1234 1330 : gel(E, i) = conductor_expo_pr(gens, chip);
1235 1330 : ic += ncp;
1236 : }
1237 1260 : if (faN) *faN = famat_remove_trivial(mkmat2(P,E));
1238 1260 : return idealfactorback(gchar_get_nf(gc), P, E, 0); /* red = 0 */
1239 : }
1240 :
1241 : /* ={sigma} s.t. k_sigma = 1 */
1242 : static GEN
1243 16863 : gcharlog_conductor_oo(GEN gc, GEN chi)
1244 : {
1245 16863 : long noo, i, n0 = gchar_get_ns(gc) + gchar_get_nc(gc);
1246 : GEN k_real, chi_oo, moo, den;
1247 :
1248 16863 : moo = locs_get_m_infty(gchar_get_zm(gc));
1249 16863 : noo = lg(moo)-1;
1250 16863 : k_real = vecslice(chi, n0-noo+1, n0);
1251 16863 : chi_oo = zerovec(gchar_get_r1(gc));
1252 18473 : for (i=1; i<=noo; i++)
1253 : {
1254 1610 : den = Q_denom(gel(k_real,i));
1255 1610 : if (den && !equali1(den)) gel(chi_oo, moo[i]) = gen_1;
1256 : }
1257 16863 : return chi_oo;
1258 : }
1259 :
1260 : static GEN
1261 224 : gchari_conductor(GEN gc, GEN chi)
1262 : {
1263 224 : chi = gchari_duallog(gc, chi);
1264 224 : return mkvec2(gcharlog_conductor_f(gc, chi, NULL),
1265 : gcharlog_conductor_oo(gc, chi));
1266 : }
1267 : GEN
1268 224 : gchar_conductor(GEN gc, GEN chi)
1269 : {
1270 224 : pari_sp av = avma;
1271 224 : check_gchar_group(gc);
1272 224 : return gerepilecopy(av, gchari_conductor(gc, gchar_internal(gc, chi, NULL)));
1273 : }
1274 :
1275 : int
1276 189 : gcharisalgebraic(GEN gc, GEN chi, GEN *pq)
1277 : {
1278 : long i, ntors, nfree, n0, nalg, r1, r2;
1279 : GEN w, chii, v;
1280 189 : pari_sp av = avma;
1281 :
1282 189 : check_gchar_group(gc);
1283 : /* in snf basis */
1284 189 : ntors = gchar_get_ntors(gc);
1285 189 : nfree = gchar_get_nfree(gc);
1286 : /* in internal basis */
1287 189 : r1 = gchar_get_r1(gc);
1288 189 : r2 = gchar_get_r2(gc);
1289 189 : n0 = gchar_get_nm(gc) - r2; /* last index before k_s */
1290 : /* in both */
1291 189 : nalg = gchar_get_nalg(gc); /* number of generators of free algebraic chars */
1292 :
1293 189 : chii = gchar_internal(gc, chi, &w);
1294 : /* check component are on algebraic generators */
1295 399 : for (i = ntors+nalg+1; i <= ntors+nfree; i++)
1296 252 : if (!gequal0(gel(chi,i))) return gc_bool(av, 0);
1297 147 : chii = gchari_duallog(gc, chii);
1298 :
1299 147 : if (r1) /* not totally complex: finite order * integral power of norm */
1300 : {
1301 77 : if (typ(w) != t_INT) return gc_bool(av, 0);
1302 63 : w = negi(w);
1303 63 : if (pq)
1304 : { /* set the infinity type */
1305 63 : v = cgetg(r1+r2+1, t_VEC);
1306 147 : for (i = 1; i <= r1; i++) gel(v, i) = mkvec2(w, gen_0);
1307 91 : for ( ; i <= r1+r2; i++) gel(v, i) = mkvec2(w, w);
1308 : }
1309 : }
1310 : else /* totally complex */
1311 : {
1312 : int w2;
1313 70 : w = gneg(gmul2n(w, 1));
1314 70 : if (typ(w) != t_INT) return gc_bool(av, 0);
1315 63 : w2 = mpodd(w);
1316 189 : for (i = 1; i <= r2; i++) /* need k_s + w = 0 mod 2 for all s */
1317 133 : if (mpodd(gel(chii, n0 + i)) != w2) return gc_bool(av, 0);
1318 56 : if (pq)
1319 : { /* set the infinity type */
1320 42 : v = cgetg(r2+1, t_VEC);
1321 140 : for (i = 1; i <= r2; i++)
1322 : {
1323 98 : GEN p = gmul2n(subii(w, gel(chii, n0+i)), -1);
1324 98 : gel(v, i) = mkvec2(p, subii(w, p));
1325 : }
1326 : }
1327 : }
1328 119 : if (pq) { *pq = gerepilecopy(av, v); av = avma; }
1329 119 : return gc_bool(av, 1);
1330 : }
1331 :
1332 : GEN
1333 497 : gcharlocal(GEN gc, GEN chi, GEN v, long prec, GEN* pbid)
1334 : {
1335 497 : pari_sp av = avma;
1336 497 : GEN nf = gchar_get_nf(gc), s, chiv, logchi;
1337 :
1338 497 : check_gchar_group(gc);
1339 497 : chi = gchar_internal(gc, chi, &s);
1340 490 : logchi = gchari_duallog(gc, chi);
1341 490 : if (typ(v) == t_INT) /* v infinite */
1342 : {
1343 385 : long i, r1, r2, tau = itos(v), n0 = gchar_get_ns(gc) + gchar_get_nc(gc);
1344 : GEN phi, k;
1345 385 : nf_get_sign(nf, &r1, &r2);
1346 385 : if (tau <= 0)
1347 7 : pari_err_DOMAIN("gcharlocal [index of an infinite place]", "v", "<=", gen_0, v);
1348 378 : if (tau > r1+r2)
1349 7 : pari_err_DOMAIN("gcharlocal [index of an infinite place]", "v", ">", stoi(r1+r2), v);
1350 371 : if (r1+r2 == 1) phi = gen_0;
1351 357 : else phi = gel(logchi, n0 + tau);
1352 371 : if (tau <= r1) /* v real */
1353 : {
1354 133 : GEN moo = gel(gchar_get_mod(gc),2);
1355 133 : i = zv_search(moo, tau);
1356 133 : if (i==0) k = gen_0;
1357 : else
1358 : {
1359 21 : k = gel(logchi, n0 - lg(moo) + 1 + i); /* 0 or 1/2 */
1360 21 : if (!gequal0(k)) k = gen_1;
1361 : }
1362 : }
1363 : else /* v complex */
1364 238 : k = gel(logchi, n0 + r2 + tau);
1365 371 : if (s) phi = gsub(phi, mulcxI(s));
1366 371 : chiv = mkvec2(k, phi);
1367 : }
1368 : else /* v finite */
1369 : {
1370 105 : GEN P = gchar_get_modP(gc);
1371 : long iv;
1372 105 : checkprid(v);
1373 105 : iv = gen_search(P, v, (void*)cmp_prime_ideal, cmp_nodata);
1374 105 : chiv = gchari_eval(gc, NULL, v, 0, logchi, s, prec);
1375 105 : if (iv <= 0) chiv = mkvec(chiv);
1376 : else
1377 : {
1378 63 : GEN cyc, w, Lsprk, bid, chip = NULL;
1379 63 : long i, ic, l = lg(P);
1380 63 : Lsprk = locs_get_Lsprk(gchar_get_zm(gc));
1381 70 : for (i = 1, ic = gchar_get_ns(gc); i < l; i++)
1382 : {
1383 70 : long ncp = sprk_get_ncp(gel(Lsprk,i));
1384 70 : if (i == iv) { chip = vecslice(logchi, ic + 1, ic + ncp); break; }
1385 7 : ic += ncp;
1386 : }
1387 63 : if (!chip) pari_err_BUG("gcharlocal (chip not found)");
1388 : /* TODO store bid instead of recomputing? */
1389 63 : bid = sprk_to_bid(nf, gel(Lsprk,i), nf_INIT);
1390 63 : cyc = bid_get_cyc(bid);
1391 63 : w = RgV_RgM_mul(chip, gel(bid_get_U(bid),1));
1392 168 : for (i = 1; i < lg(w); i++)
1393 105 : gel(w,i) = modii(gmul(gel(w,i), gel(cyc,i)), gel(cyc,i));
1394 63 : chiv = vec_append(w, chiv);
1395 63 : if (pbid) { *pbid = bid; return gc_all(av, 2, &chiv, pbid); }
1396 : }
1397 : }
1398 448 : return gerepilecopy(av, chiv);
1399 : }
1400 :
1401 :
1402 : /*******************************************************************/
1403 : /* */
1404 : /* EVALUATION OF CHARACTERS */
1405 : /* */
1406 : /*******************************************************************/
1407 : GEN
1408 1799 : gcharduallog(GEN gc, GEN chi)
1409 : {
1410 1799 : pari_sp av = avma;
1411 : GEN logchi, s;
1412 1799 : check_gchar_group(gc);
1413 1799 : logchi = gchari_duallog(gc, gchar_internal(gc, chi, &s));
1414 1799 : return gerepilecopy(av, shallowconcat1(mkcol2(logchi,s)));
1415 : }
1416 :
1417 : static GEN
1418 188468 : gcharisprincipal(GEN gc, GEN x, GEN *palpha)
1419 : {
1420 188468 : GEN t, v, bnf = gchar_get_bnf(gc), DLdata = gchar_get_DLdata(gc);
1421 188468 : t = bnfisprincipal0(bnf, x, nf_GENMAT | nf_FORCE); v = gel(t, 1);
1422 188468 : *palpha = famat_reduce(famat_mul(nffactorback(bnf, gel(DLdata,2), v),
1423 188468 : gel(t, 2)));
1424 188468 : return ZM_ZC_mul(gel(DLdata,1), v);
1425 : }
1426 :
1427 : /* complete log of ideal; if logchi != NULL make sure precision is
1428 : * sufficient to evaluate gcharlog_eval_raw to attain 'prec' */
1429 : static GEN
1430 188468 : gchar_log(GEN gc, GEN x, GEN logchi, long prec)
1431 : {
1432 188468 : GEN zm, v, alpha, arch_log = NULL, zm_log, nf;
1433 :
1434 188468 : nf = gchar_get_nf(gc);
1435 188468 : zm = gchar_get_zm(gc);
1436 188468 : v = gcharisprincipal(gc, x, &alpha); /* alpha a GENMAT */
1437 188468 : if (DEBUGLEVEL>2) err_printf("v %Ps\n", v);
1438 188468 : zm_log = gchar_logm(nf, zm, alpha);
1439 188468 : if (DEBUGLEVEL>2) err_printf("zm_log(alpha) %Ps\n", zm_log);
1440 :
1441 188468 : if (logchi)
1442 : { /* check if precision is sufficient, take care of gexpo = -infty */
1443 184681 : long e, bit = expu(nf_get_degree(nf) + lg(zm_log)-1) + 3;
1444 184681 : e = gexpo(logchi); if (e > 0) bit += e;
1445 184681 : e = gexpo(gel(alpha,1)); if (e > 0) bit += e;
1446 184681 : prec += nbits2extraprec(bit);
1447 : }
1448 : for(;;)
1449 : {
1450 188468 : arch_log = nfembedlog(&nf, alpha, prec);
1451 188468 : if (arch_log) break;
1452 0 : prec = precdbl(prec);
1453 : }
1454 188468 : if (DEBUGLEVEL>2) err_printf("arch log %Ps\n", arch_log);
1455 188468 : return shallowconcat1(mkvec3(v, gneg(zm_log), gneg(arch_log)));
1456 : }
1457 :
1458 : /* gp version, with norm component */
1459 : GEN
1460 70 : gcharlog(GEN gc, GEN x, long prec)
1461 : {
1462 70 : pari_sp av = avma;
1463 : GEN norm;
1464 :
1465 70 : check_gchar_group(gc);
1466 70 : norm = idealnorm(gchar_get_nf(gc), x);
1467 70 : norm = mkcomplex(gen_0, gdiv(glog(norm,prec), Pi2n(1,prec)));
1468 70 : return gerepilecopy(av, vec_append(gchar_log(gc, x, NULL, prec), norm));
1469 : }
1470 :
1471 : static GEN
1472 199339 : gcharlog_eval_raw(GEN logchi, GEN logx)
1473 199339 : { GEN v = RgV_dotproduct(logchi, logx); return gsub(v, ground(v)); }
1474 :
1475 : /* if flag = 1 -> value in C, flag = 0 -> value in R/Z
1476 : * chi in chari format without norm component (given in s)
1477 : * if logchi is provided, assume it has enough precision
1478 : * TODO check that logchi argument is indeed used correctly by callers */
1479 : static GEN
1480 184681 : gchari_eval(GEN gc, GEN chi, GEN x, long flag, GEN logchi, GEN s, long prec)
1481 : {
1482 : GEN z, norm, logx;
1483 184681 : if (!logchi)
1484 : {
1485 3955 : long e, precgc = prec, prec0 = gchar_get_prec(gc);
1486 3955 : logchi = gchari_duallog(gc, chi);
1487 3955 : e = gexpo(logchi); if (e > 0) precgc += nbits2extraprec(e);
1488 3955 : if (precgc > prec0)
1489 : {
1490 14 : gc = gcharnewprec(gc, precgc);
1491 14 : logchi = gchari_duallog(gc, chi);
1492 : }
1493 : }
1494 184681 : logx = gchar_log(gc, x, logchi, prec);
1495 184681 : norm = gequal0(s)? NULL: idealnorm(gchar_get_nf(gc), x);
1496 184681 : z = gcharlog_eval_raw(logchi, logx);
1497 184681 : if (flag)
1498 : {
1499 180887 : z = expIPiC(gmul2n(z, 1), prec);
1500 180887 : if (norm) z = gmul(z, gpow(norm, gneg(s), prec));
1501 : }
1502 3794 : else if (norm)
1503 70 : z = gadd(z, mulcxI(gdiv(gmul(s, glog(norm,prec)), Pi2n(1,prec))));
1504 184681 : if (DEBUGLEVEL>1) err_printf("char value %Ps\n", z);
1505 184681 : return z;
1506 : }
1507 :
1508 : GEN
1509 3969 : gchareval(GEN gc, GEN chi, GEN x, long flag)
1510 : {
1511 : GEN s;
1512 : long prec;
1513 3969 : pari_sp av = avma;
1514 3969 : check_gchar_group(gc);
1515 3969 : prec = gchar_get_evalprec(gc);
1516 3969 : chi = gchar_internal(gc, chi, &s);
1517 3955 : return gerepileupto(av, gchari_eval(gc, chi, x, flag, NULL, s, prec));
1518 : }
1519 :
1520 : /*******************************************************************/
1521 : /* */
1522 : /* IDENTIFICATION */
1523 : /* */
1524 : /*******************************************************************/
1525 : static GEN
1526 112 : col_2ei(long n, long i) { GEN e = zerocol(n); gel(e,i) = gen_2; return e; }
1527 : static GEN
1528 4011 : gchar_identify_init(GEN gc, GEN Lv, long prec)
1529 : {
1530 : GEN M, cyc, mult, Lpr, Lk1, Lphi1, Lk2, Llog, eps, U, P, nf, moo, vlogchi;
1531 : long beps, r1, r2, d, nmiss, n, npr, nk1, nchi, s, i, j, l, dim, n0, ncol;
1532 :
1533 4011 : check_gchar_group(gc);
1534 3990 : n0 = gchar_get_ns(gc) + gchar_get_nc(gc); /* last index of torsion char */
1535 3990 : cyc = gchar_get_cyc(gc);
1536 3990 : P = gchar_get_modP(gc);
1537 3990 : moo = gel(gchar_get_mod(gc), 2);
1538 3990 : nf = gchar_get_nf(gc); nf_get_sign(nf, &r1, &r2);
1539 3990 : nchi = lg(cyc)-2; /* ignore norm */
1540 3990 : d = r1 + 2*r2;
1541 3990 : mult = (nchi >= d)? gel(cyc,1): gen_1;
1542 3990 : s = (8*prec2nbits(prec))/10; mult = shifti(mult, s);
1543 3990 : beps = -(7*s) / 16; eps = real2n(beps, prec);
1544 3990 : l = lg(Lv);
1545 3990 : if (lg(gen_sort_uniq(Lv, (void*)cmp_universal, cmp_nodata)) != l)
1546 7 : pari_err(e_MISC, "components of Lv must be distinct");
1547 : /* log of prime ideals */
1548 3983 : Llog = cgetg(l, t_VEC);
1549 : /* index in Lchiv corresponding to the places */
1550 3983 : Lpr = cgetg(l, t_VECSMALL);
1551 3983 : Lk1 = cgetg(l, t_VECSMALL);
1552 3983 : Lphi1 = zero_zv(r1);
1553 3983 : Lk2 = zero_zv(r2); nmiss = d;
1554 11865 : for (i = 1, npr = nk1 = 0; i < l; i++)
1555 7924 : if (typ(gel(Lv,i)) == t_INT)
1556 : {
1557 4179 : long v = itos(gel(Lv,i));
1558 4179 : if (v <= 0) pari_err_DOMAIN("gcharidentify", "v", "<=", gen_0, stoi(v));
1559 4172 : if (v > r1+r2)
1560 7 : pari_err_DOMAIN("gcharidentify", "v", ">", stoi(r1+r2), stoi(v));
1561 4165 : if (v <= r1)
1562 : { /* don't put in k1 if not in conductor (but keep as phi) */
1563 245 : if (zv_search(moo, v)) { nk1++; Lk1[nk1] = i; }
1564 245 : Lphi1[v] = i; nmiss--;
1565 : }
1566 : else
1567 : {
1568 3920 : Lk2[v-r1] = i; nmiss -= 2;
1569 : }
1570 : }
1571 : else
1572 : {
1573 3745 : GEN pr = gel(Lv,i);
1574 3745 : long lP = lg(P);
1575 3745 : if (idealtyp(&pr, NULL) != id_PRIME)
1576 7 : pari_err_TYPE("gcharidentify [ideal]", pr);
1577 4046 : for (j = 1; j < lP; j++)
1578 329 : if (pr_equal(pr, gel(P,j)))
1579 7 : pari_err_COPRIME("gcharidentify", pr, gel(gchar_get_mod(gc), 1));
1580 3717 : npr++;
1581 3717 : Lpr[npr] = i;
1582 3717 : gel(Llog,npr) = gchar_log(gc, pr, NULL, prec);
1583 : }
1584 3941 : setlg(Llog, npr+1); setlg(Lpr, npr+1); setlg(Lk1, nk1+1);
1585 : /* build matrix M */
1586 3941 : n = npr + nk1; dim = n + d; ncol = n + 1 + nchi;
1587 3941 : M = cgetg(ncol+1, t_MAT);
1588 7658 : for (j = 1; j <= npr; j++) gel(M,j) = col_ei(dim, j);
1589 4053 : for (; j <= n; j++) gel(M,j) = col_2ei(dim, j);
1590 3941 : gel(M,j) = zerocol(dim);
1591 11858 : for (i = n+1; i <= n+r1+r2; i++) gcoeff(M,i,j) = eps;
1592 3941 : vlogchi = RgM_mul(gchar_get_Ui(gc), gchar_get_basis(gc));
1593 19215 : for (j=1; j<=nchi; j++)
1594 : {
1595 15274 : GEN logchi = RgV_frac_inplace(row(vlogchi, j), n0); /* duallog(e_j) */
1596 15274 : GEN chi_oo = gcharlog_conductor_oo(gc, logchi), C = cgetg(dim+1, t_COL);
1597 29932 : for (i=1; i<=npr; i++) gel(C,i) = gcharlog_eval_raw(logchi, gel(Llog,i));
1598 15666 : for (i=1; i<=nk1; i++) gel(C,npr+i) = gel(chi_oo, itos(gel(Lv,Lk1[i])));
1599 76069 : for (i=1; i<=d; i++) gel(C,n+i) = gel(logchi, n0 + i);
1600 15274 : gel(M,n+1+j) = C;
1601 : }
1602 4382 : for (i = 1; i <= r1; i++)
1603 441 : if (!Lphi1[i])
1604 : {
1605 196 : long a = n + i;
1606 959 : for (j = 1; j <= ncol; j++) gcoeff(M,a,j) = gmul2n(gcoeff(M,a,j), beps);
1607 : }
1608 11417 : for (i = 1; i <= r2; i++)
1609 7476 : if (!Lk2[i])
1610 : {
1611 3584 : long a = n + r1 + i, b = a + r2;
1612 25032 : for (j = 1; j <= ncol; j++)
1613 : {
1614 21448 : gcoeff(M,a,j) = gmul2n(gcoeff(M,a,j), beps);
1615 21448 : gcoeff(M,b,j) = gmul2n(gcoeff(M,b,j), beps);
1616 : }
1617 : }
1618 : /* rescale and truncate M, then LLL-reduce M */
1619 3941 : M = gtrunc(RgM_Rg_mul(M, mult));
1620 3941 : U = ZM_lll(M, 0.99, LLL_IM);
1621 3941 : M = ZM_mul(M, U);
1622 3941 : U = rowslice(U, n + 2, n + 1 + nchi);
1623 3941 : return mkvecn(9, M, U, Lpr, Lk1, Lphi1, Lk2, mult, Lv,
1624 : mkvecsmall3(prec, nmiss, beps));
1625 : }
1626 :
1627 : /* TODO return the distance between the character found and the conditions? */
1628 : static GEN
1629 3941 : gchar_identify_i(GEN gc, GEN idinit, GEN Lchiv)
1630 : {
1631 : GEN M, U, Lpr, Lk1, Lphi1, Lk2, mult, cyc, y, x, sumphi, Lv, Norm, nf;
1632 : long i, l, r1, r2, beps, npr, nk1, n, nmiss, nnorm, prec;
1633 3941 : M = gel(idinit,1);
1634 3941 : U = gel(idinit,2);
1635 3941 : Lpr = gel(idinit,3);
1636 3941 : Lk1 = gel(idinit,4);
1637 3941 : Lphi1 = gel(idinit,5);
1638 3941 : Lk2 = gel(idinit,6);
1639 3941 : mult = gel(idinit,7);
1640 3941 : Lv = gel(idinit,8);
1641 3941 : prec = gel(idinit,9)[1];
1642 3941 : nmiss = gel(idinit,9)[2];
1643 3941 : beps = gel(idinit,9)[3];
1644 3941 : npr = lg(Lpr)-1;
1645 3941 : nk1 = lg(Lk1)-1; n = npr + nk1;
1646 3941 : cyc = gchar_get_cyc(gc);
1647 3941 : nf = gchar_get_nf(gc);
1648 3941 : nf_get_sign(nf, &r1, &r2);
1649 3941 : nnorm = 0;
1650 3941 : Norm = gen_0;
1651 :
1652 3941 : l = lg(Lchiv); Lchiv = shallowcopy(Lchiv);
1653 3941 : if (lg(Lv) != l) pari_err_DIM("gcharidentify [#Lv != #Lchiv]");
1654 11739 : for (i = 1; i < l; i++)
1655 : {
1656 7847 : GEN t = gel(Lchiv,i), u;
1657 7847 : if (typ(gel(Lv,i)) != t_INT)
1658 : {
1659 3717 : if (typ(t) == t_VEC) /* value at last component */
1660 14 : gel(Lchiv,i) = t = gel(t, lg(t)-1);
1661 3717 : if (typ(t) == t_COMPLEX)
1662 : {
1663 21 : nnorm++; /* 2 Pi Im(theta) / log N(pr) */
1664 21 : Norm = gadd(Norm, gdiv(gmul(Pi2n(1,prec), gel(t,2)),
1665 21 : glog(idealnorm(nf,gel(Lv,i)),prec)));
1666 21 : gel(Lchiv,i) = t = gel(t,1);
1667 : }
1668 3717 : if (!is_real_t(typ(t)))
1669 14 : pari_err_TYPE("gcharidentify [character value: should be real or complex]", t);
1670 : }
1671 : else
1672 : {
1673 4130 : if (typ(t) != t_VEC)
1674 7 : pari_err_TYPE("gcharidentify [character at infinite place: should be a t_VEC]", t);
1675 4123 : if (lg(t) != 3)
1676 7 : pari_err_DIM("gcharidentify [character at infinite place: should have two components]");
1677 4116 : if (typ(gel(t,1)) != t_INT)
1678 7 : pari_err_TYPE("gcharidentify [k parameter at infinite place: should be a t_INT]", gel(t,1));
1679 4109 : u = gel(t,2);
1680 4109 : if (typ(u) == t_COMPLEX)
1681 : {
1682 133 : nnorm++;
1683 133 : Norm = gsub(Norm, gel(u,2)); u = gel(u,1);
1684 133 : gel(Lchiv, i) = mkvec2(gel(t,1), u);
1685 : }
1686 4109 : if (!is_real_t(typ(u)))
1687 7 : pari_err_TYPE("gcharidentify [phi parameter at infinite place: should be real or complex]", u);
1688 : }
1689 : }
1690 :
1691 : /* construct vector */
1692 3892 : y = zerocol(n + r1 + 2*r2); sumphi = gen_0;
1693 7595 : for (i=1; i<=npr; i++) gel(y,i) = gel(Lchiv, Lpr[i]);
1694 4004 : for (i=1; i<=nk1; i++) gel(y,npr+i) = gmael(Lchiv,Lk1[i],1);
1695 4333 : for (i=1; i<=r1; i++)
1696 441 : if (Lphi1[i])
1697 : {
1698 245 : gel(y, n+i) = x = gmael(Lchiv,Lphi1[i],2);
1699 245 : sumphi = gadd(sumphi, x);
1700 : }
1701 11270 : for (i=1; i<=r2; i++)
1702 7378 : if (Lk2[i])
1703 : {
1704 3857 : long a = n + r1 + i;
1705 3857 : gel(y, a + r2) = gmael(Lchiv,Lk2[i],1);
1706 3857 : gel(y, a) = x = gmael(Lchiv,Lk2[i],2);
1707 3857 : sumphi = gadd(sumphi, gshift(x,1));
1708 : }
1709 3892 : if (nmiss)
1710 : {
1711 1925 : sumphi = gmul2n(gdivgs(sumphi, -nmiss), beps);
1712 2226 : for (i = 1; i <= r1; i++) if (!Lphi1[i]) gel(y, n + i) = sumphi;
1713 5502 : for (i = 1; i <= r2; i++) if (!Lk2[i]) gel(y, n + r1+i) = sumphi;
1714 : }
1715 3892 : y = gtrunc(RgC_Rg_mul(y, mult));
1716 :
1717 : /* find approximation */
1718 3892 : x = ZM_ZC_mul(U, RgM_Babai(M, y));
1719 18970 : for (i = 1; i < lg(cyc)-1; i++) /* ignore norm */
1720 15078 : if (signe(gel(cyc,i))) gel(x,i) = modii(gel(x,i), gel(cyc,i));
1721 3892 : if (nnorm > 0) x = vec_append(x, gdivgu(Norm, lg(Lv)-1));
1722 3892 : return x;
1723 : }
1724 :
1725 : /* TODO export the init interface */
1726 : GEN
1727 4011 : gchar_identify(GEN gc, GEN Lv, GEN Lchiv, long prec)
1728 : {
1729 4011 : pari_sp av = avma;
1730 4011 : GEN idinit = gchar_identify_init(gc, Lv, prec);
1731 3941 : return gerepilecopy(av, gchar_identify_i(gc,idinit,Lchiv));
1732 : }
1733 :
1734 : /*******************************************************************/
1735 : /* */
1736 : /* L FUNCTION */
1737 : /* */
1738 : /*******************************************************************/
1739 :
1740 : /* TODO: could merge with vecan_chigen */
1741 :
1742 : /* return x + yz; y != 0; z = 0,1 "often"; x = 0 "often" */
1743 : static GEN
1744 2476278 : gaddmul(GEN x, GEN y, GEN z)
1745 : {
1746 : pari_sp av;
1747 2476278 : if (typ(z) == t_INT)
1748 : {
1749 2191739 : if (!signe(z)) return x;
1750 13486 : if (equali1(z)) return gadd(x,y);
1751 : }
1752 293244 : if (isintzero(x)) return gmul(y,z);
1753 117467 : av = avma;
1754 117467 : return gerepileupto(av, gadd(x, gmul(y,z)));
1755 : }
1756 :
1757 : GEN
1758 693 : vecan_gchar(GEN an, long n, long prec)
1759 : {
1760 : forprime_t iter;
1761 693 : GEN gc = gel(an,1), chi = gel(an,2), P = gel(an,3), PZ = gel(an,4);
1762 693 : GEN BOUND = stoi(n), v = vec_ei(n, 1);
1763 693 : GEN gp = cgetipos(3), nf, chilog, s;
1764 : ulong p;
1765 :
1766 : /* prec increase: 1/n*log(N(pmax)) < log(pmax) */
1767 693 : if (DEBUGLEVEL > 1)
1768 0 : err_printf("vecan_gchar: need extra prec %ld\n", nbits2extraprec(expu(n)));
1769 693 : gc = gcharnewprec(gc, prec + nbits2extraprec(expu(n)));
1770 693 : chilog = gchari_duallog(gc, check_gchari(gc, chi, &s));
1771 :
1772 693 : nf = gchar_get_nf(gc);
1773 : /* FIXME: when many log of many primes are computed:
1774 : - bnfisprincipal may not be improved
1775 : - however we can precompute the logs of generators
1776 : for principal part
1777 : - when galois, should compute one ideal by orbit.
1778 : - when real, clear imaginary part
1779 : */
1780 :
1781 693 : u_forprime_init(&iter, 2, n);
1782 181902 : while ((p = u_forprime_next(&iter)))
1783 : {
1784 : GEN L;
1785 : long j;
1786 181209 : int check = !umodiu(PZ,p);
1787 181209 : gp[2] = p;
1788 181209 : L = idealprimedec_limit_norm(nf, gp, BOUND);
1789 362250 : for (j = 1; j < lg(L); j++)
1790 : {
1791 181041 : GEN pr = gel(L, j), ch;
1792 : pari_sp av;
1793 : ulong k, q;
1794 181041 : if (check && gen_search(P, pr, (void*)cmp_prime_ideal, cmp_nodata) > 0)
1795 434 : continue;
1796 : /* TODO: extract code and use precom sprk? */
1797 180607 : av = avma;
1798 180607 : ch = gchari_eval(gc, chi, pr, 1, chilog, gen_0, prec);
1799 180607 : ch = gerepileupto(av, ch);
1800 180607 : q = upr_norm(pr);
1801 180607 : gel(v, q) = gadd(gel(v, q), ch);
1802 2656885 : for (k = 2*q; k <= (ulong)n; k += q)
1803 2476278 : gel(v, k) = gaddmul(gel(v, k), ch, gel(v, k/q));
1804 : }
1805 : }
1806 : /* weight, could consider shifting s at eval instead */
1807 693 : if (!gequal0(s))
1808 : {
1809 329 : GEN ns = dirpowers(n, gneg(s), prec);
1810 : long j;
1811 285964 : for (j = 1; j <= n; j++)
1812 285635 : if (gel(v,j) != gen_0) gel(v, j) = gmul(gel(v,j),gel(ns,j));
1813 : }
1814 693 : return v;
1815 : }
1816 :
1817 : GEN
1818 14 : eulerf_gchar(GEN an, GEN p, long prec)
1819 : {
1820 14 : GEN gc = gel(an,1), chi = gel(an,2), P = gel(an,3), PZ = gel(an,4);
1821 : GEN f, L, nf, chilog, s;
1822 : int check;
1823 : long j, l;
1824 :
1825 : /* prec increase: 1/n*log(N(pmax)) < log(pmax) */
1826 14 : if (DEBUGLEVEL > 1)
1827 0 : err_printf("vecan_gchar: need extra prec %ld\n", nbits2extraprec(expi(p)));
1828 14 : gc = gcharnewprec(gc, prec + nbits2extraprec(expi(p)));
1829 14 : chilog = gchari_duallog(gc, check_gchari(gc, chi, &s));
1830 :
1831 14 : nf = gchar_get_nf(gc);
1832 14 : f = pol_1(0);
1833 14 : check = dvdii(PZ, p);
1834 14 : L = idealprimedec(nf, p); l = lg(L);
1835 28 : for (j = 1; j < l; j++)
1836 : {
1837 14 : GEN pr = gel(L, j), ch;
1838 14 : if (check && gen_search(P, pr, (void*)cmp_prime_ideal, cmp_nodata) > 0)
1839 0 : continue;
1840 14 : ch = gchari_eval(gc, chi, pr, 1, chilog, s, prec);
1841 14 : f = gmul(f, gsub(gen_1, monomial(ch, pr_get_f(pr), 0)));
1842 : }
1843 14 : return mkrfrac(gen_1,f);
1844 : }
1845 :
1846 : static GEN
1847 1323 : cleanup_vga(GEN vga, long prec)
1848 : {
1849 : GEN ind;
1850 : long bitprec, i, l;
1851 1323 : if (!prec) return vga; /* already exact */
1852 1323 : bitprec = prec2nbits(prec);
1853 1323 : vga = shallowcopy(vga); l = lg(vga);
1854 5180 : for (i = 1; i < l; i++)
1855 : {
1856 3857 : GEN z = gel(vga,i);
1857 3857 : if (typ(z) != t_COMPLEX) continue;
1858 2086 : if (gexpo(gel(z,2)) < -bitprec+20) gel(vga,i) = gel(z,1);
1859 : }
1860 1323 : ind = indexsort(imag_i(vga));
1861 3857 : for (i = 2; i < l; i++)
1862 : {
1863 2534 : GEN z = gel(vga,ind[i]), t;
1864 2534 : if (typ(z) != t_COMPLEX) continue;
1865 1456 : t = imag_i(gel(vga, ind[i-1]));
1866 1456 : if (gexpo(gsub(gel(z,2), t)) < -bitprec+20)
1867 364 : gel(vga, ind[i]) = mkcomplex(gel(z,1), t);
1868 : }
1869 5180 : for (i = 1; i < l; i++)
1870 : {
1871 3857 : GEN z = gel(vga,i);
1872 3857 : if (typ(z) != t_COMPLEX) continue;
1873 2086 : gel(vga, i) = mkcomplex(gel(z,1), bestappr(gel(z,2), int2n(bitprec/2)));
1874 : }
1875 1323 : return vga;
1876 : }
1877 :
1878 : /* TODO: move to lfunutils, use lfunzeta and lfunzetak */
1879 : GEN
1880 1372 : gchari_lfun(GEN gc, GEN chi, GEN s0)
1881 : {
1882 : GEN nf, chilog, s, cond_f, cond_oo, vga_r, vga_c, chiw;
1883 : GEN v_phi, v_arg, sig, k, NN, faN, P;
1884 : long ns, nc, nm, r1, r2;
1885 :
1886 1372 : nf = gchar_get_nf(gc);
1887 1372 : ns = gchar_get_ns(gc);
1888 1372 : nc = gchar_get_nc(gc);
1889 1372 : nm = gchar_get_nm(gc);
1890 1372 : nf_get_sign(nf, &r1, &r2);
1891 1372 : chi = check_gchari(gc, chi, &s);
1892 1372 : chilog = gchari_duallog(gc, chi);
1893 1372 : s = gadd(s0,s); chiw = shallowconcat(chi, s);
1894 1372 : if (!gequal0(imag_i(s)))
1895 7 : pari_err_IMPL("lfun for gchar with imaginary norm component");
1896 1365 : cond_f = gcharlog_conductor_f(gc, chilog, &faN);
1897 1365 : P = gel(faN, 1); /* prime ideals dividing cond(chi) */
1898 1365 : cond_oo = gcharlog_conductor_oo(gc, chilog);
1899 :
1900 1365 : NN = mulii(idealnorm(nf, cond_f), absi_shallow(nf_get_disc(nf)));
1901 1365 : if (equali1(NN)) return lfunshift(lfuncreate(gen_1), gneg(s), 0,
1902 : prec2nbits(gchar_get_evalprec(gc)));
1903 1344 : if (ZV_equal0(chi)) return lfunshift(lfuncreate(nf), gneg(s), 0,
1904 : prec2nbits(gchar_get_evalprec(gc)));
1905 :
1906 : /* vga_r = vector(r1,k,I*c[ns+nc+k]-s + cond_oo[k]);
1907 : * vga_c = vector(r2,k,I*c[ns+nc+r1+k]+abs(c[ns+nc+r1+r2+k])/2-s) */
1908 1323 : v_phi = gmul(vecslice(chilog, ns+nc+1, ns+nc+r1+r2), gen_I());
1909 1323 : v_arg = gdivgs(gabs(vecslice(chilog,ns+nc+r1+r2+1,nm),BITS_IN_LONG), 2);
1910 1323 : vga_r = gadd(vecslice(v_phi, 1, r1), cond_oo);
1911 1323 : vga_c = gadd(vecslice(v_phi, r1+1, r1+r2), v_arg);
1912 1323 : sig = shallowconcat1(mkvec3(vga_r,vga_c,gadd(vga_c,const_vec(r2,gen_1))));
1913 : /* TODO: remove cleanup when gammamellinv takes ldata*/
1914 1323 : sig = cleanup_vga(sig, gchar_get_prec(gc));
1915 1323 : k = gen_1;
1916 1323 : if (!gequal0(s))
1917 : {
1918 : long j;
1919 2086 : for (j = 1; j < lg(sig); j++) gel(sig, j) = gadd(gel(sig, j), s);
1920 539 : k = gsub(k, gmulgs(s,2));
1921 : }
1922 :
1923 : #define tag(x,t) mkvec2(mkvecsmall(t),x)
1924 1323 : return mkvecn(6, tag(mkvec4(gc, chiw, P, prV_lcm_capZ(P)), t_LFUN_HECKE),
1925 : gen_1, sig, k, NN, gen_0);
1926 : }
1927 :
1928 : GEN
1929 392 : lfungchar(GEN gc, GEN chi)
1930 : {
1931 392 : pari_sp av = avma;
1932 : GEN s;
1933 392 : check_gchar_group(gc);
1934 392 : chi = gchar_internal(gc, chi, &s);
1935 378 : return gerepilecopy(av, gchari_lfun(gc, chi, s));
1936 : }
|
