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