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 637 : compute_Lcyc(GEN Lsprk, GEN moo)
47 : {
48 637 : long i, l = lg(Lsprk), len = l+lg(moo)-1;
49 637 : GEN Lcyc = cgetg(len,t_VEC);
50 1260 : for (i = 1; i < l; i++) gel(Lcyc,i) = sprk_get_cyc(gel(Lsprk,i));
51 756 : for ( ; i < len; i++) gel(Lcyc,i) = mkvec(gen_2);
52 637 : return Lcyc;
53 : }
54 :
55 : static long
56 2051 : 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 637 : localstar(GEN nf, GEN famod, GEN moo)
61 : {
62 637 : GEN Lcyc, cyc, Lsprk, Lgenfil, P = gel(famod,1), E = gel(famod,2);
63 637 : long i, l = lg(P);
64 :
65 637 : Lsprk = cgetg(l, t_VEC);
66 637 : Lgenfil = cgetg(l, t_VEC);
67 1260 : for (i = 1; i < l; i++)
68 : {
69 623 : long e, k = itos(gel(E,i));
70 623 : GEN v, sprk = log_prk_init(nf, gel(P,i), k, NULL);
71 623 : gel(Lsprk,i) = sprk;
72 623 : gel(Lgenfil,i) = v = cgetg(k+1, t_VEC);
73 : /* log on sprk of generators of U_{e-1}/U_e(pr) */
74 623 : gel(v, 1) = col_ei(sprk_get_ncp(sprk), 1);
75 1288 : for (e = 2; e <= k; e++) gel(v, e) = sprk_log_gen_pr2(nf, sprk, e);
76 : }
77 637 : Lcyc = compute_Lcyc(Lsprk, moo);
78 637 : if (lg(Lcyc) > 1)
79 315 : cyc = shallowconcat1(Lcyc);
80 : else
81 322 : cyc = cgetg(1, t_VEC);
82 637 : 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 199199 : nfembedlog(GEN *pnf, GEN x, long prec)
90 : {
91 199199 : pari_sp av = avma;
92 199199 : GEN logs, cxlogs, nf = *pnf;
93 199199 : long k, r1, r2, n, extrabit, extranfbit = 0, nfprec, nfprec0, logprec;
94 :
95 199199 : nfprec0 = nf_get_prec(nf);
96 199199 : nf_get_sign(nf, &r1, &r2);
97 199199 : n = r1 + 2*r2;
98 199199 : logprec = prec;
99 199199 : extrabit = expu(n) + gexpo(nf_get_M(nf)) + 10;
100 199199 : if (typ(x) == t_MAT)
101 : {
102 198562 : long l = lg(gel(x,2));
103 198562 : if (l > 1)
104 : {
105 198522 : extrabit += expu(l) + gexpo(gel(x,2));
106 198522 : extranfbit = gexpo(gel(x,1));
107 : }
108 : }
109 : else
110 : {
111 637 : x = nf_to_scalar_or_basis(nf,x);
112 637 : extranfbit = gexpo(x);
113 : }
114 199199 : 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 199199 : nfprec = prec + nbits2extraprec(extrabit + extranfbit);
119 199199 : logprec = prec + nbits2extraprec(extrabit);
120 199199 : if (nfprec0 < nfprec)
121 : {
122 3316 : if (DEBUGLEVEL)
123 0 : err_printf(" nfembedlog: increasing prec %d -> %d\n", nfprec0, nfprec);
124 3316 : *pnf = nf = nfnewprec_shallow(nf, nfprec);
125 3316 : av = avma;
126 : }
127 199199 : if (!(cxlogs = nf_cxlog(nf, x, logprec))) return gc_NULL(av);
128 199199 : if (!(cxlogs = nf_cxlog_normalize(nf, cxlogs, logprec))) return gc_NULL(av);
129 199199 : logs = cgetg(n+1,t_COL);
130 731225 : for (k = 1; k <= r1+r2; k++) gel(logs,k) = real_i(gel(cxlogs,k));
131 417772 : for ( ; k <= n; k++) gel(logs,k) = gmul2n(imag_i(gel(cxlogs,k-r2)), -1);
132 199199 : extrabit = gexpo(logs);
133 199199 : if (extrabit < 0) extrabit = 0;
134 199199 : prec += nbits2extraprec(extrabit);
135 199199 : return gc_upto(av, gdiv(logs, Pi2n(1,prec)));
136 : }
137 :
138 : static GEN
139 1925 : gchar_Sval(GEN nf, GEN S, GEN x)
140 : {
141 1925 : GEN res = cgetg(lg(S),t_COL);
142 : long i;
143 1925 : if (typ(x)==t_MAT)
144 3507 : for (i=1; i<lg(S); i++)
145 2219 : gel(res, i) = famat_nfvalrem(nf, x, gel(S,i), NULL);
146 : else
147 1015 : for (i=1; i<lg(S); i++)
148 378 : gel(res, i) = stoi(nfval(nf, x, gel(S,i)));
149 1925 : return res;
150 : }
151 :
152 : /* true nf; log_prk(x*pi_pr^{-v_pr(x)}), sign(sigma(x)) */
153 : static GEN
154 194404 : gchar_logm(GEN nf, GEN locs, GEN x)
155 : {
156 194404 : GEN moo, loga, Lsprk = locs_get_Lsprk(locs);
157 194404 : long i, l = lg(Lsprk);
158 :
159 194404 : moo = locs_get_m_infty(locs);
160 194404 : if (typ(x) != t_MAT) x = to_famat_shallow(x, gen_1);
161 194404 : loga = cgetg(l+1, t_VEC);
162 614096 : for (i = 1; i < l; i++)
163 : {
164 419692 : GEN sprk = gel(Lsprk, i), pr = sprk_get_pr(sprk);
165 419692 : GEN g = vec_append(gel(x,1), pr_get_gen(pr));
166 419692 : GEN v = famat_nfvalrem(nf, x, pr, NULL);
167 419692 : GEN e = vec_append(gel(x,2), gneg(v));
168 419692 : gel(loga, i) = famat_zlog_pr(nf, g, e, sprk, NULL);
169 : }
170 194404 : gel(loga, i) = zc_to_ZC( nfsign_arch(nf, x, moo) );
171 194404 : return shallowconcat1(loga);
172 : }
173 :
174 : static GEN
175 1925 : gchar_nflog(GEN *pnf, GEN zm, GEN S, GEN x, long prec)
176 : {
177 1925 : GEN emb = nfembedlog(pnf, x, prec);
178 1925 : if (!emb) return NULL;
179 1925 : 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 6713 : embedcol(GEN v, long size, long shift)
192 : {
193 : long k;
194 6713 : GEN w = zerocol(size);
195 51436 : for (k = 1; k < lg(v); k++) gel(w, shift+k) = gel(v,k);
196 6713 : return w;
197 : }
198 :
199 : /* write exact rationals from real approximation, in place. */
200 : static void
201 5048 : shallow_clean_rat(GEN v, long k0, long k1, GEN den, long prec)
202 : {
203 5048 : long k, e, bit = -prec2nbits(prec)/2;
204 45261 : for (k = k0; k <= k1; k++)
205 : {
206 40213 : GEN t = gel(v,k);
207 40213 : if (den) t = gmul(t, den);
208 40213 : t = grndtoi(t, &e);
209 40213 : 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 40213 : if (e > bit)
212 : pari_err_BUG("gcharinit, non rational entry"); /*LCOV_EXCL_LINE*/
213 40213 : gel(v, k) = den? gdiv(t, den): t;
214 : }
215 5048 : }
216 :
217 : /* FIXME: export ? */
218 : static GEN
219 1239 : rowreverse(GEN m)
220 : {
221 : long i, l;
222 : GEN v;
223 1239 : if (lg(m) == 1) return m;
224 1239 : l = lgcols(m); v = cgetg(l, t_VECSMALL);
225 4557 : for (i = 1; i < l; i++) v[i] = l - i;
226 1239 : 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 1239 : hnf_block(GEN m, long r0, long nr, long c0, long nc)
233 : {
234 : GEN block, u, uu;
235 1239 : long nm = lg(m)-1, k;
236 1239 : pari_sp av = avma;
237 :
238 1239 : block = matslice(m, r0+1, r0+nr, c0+1, c0+nc);
239 1239 : block = Q_remove_denom(block, NULL);
240 1239 : block = rowreverse(block);
241 :
242 1239 : (void)ZM_hnfall(block, &u, 1);
243 1239 : vecreverse_inplace(u); uu = matid(nm); /* embed in matid */
244 6188 : for (k = 1; k <= nc; k++) gel(uu,c0+k) = embedcol(gel(u,k),nm,c0);
245 1239 : return gc_GEN(av, uu);
246 : }
247 :
248 : /* (matrix, starting row, nb rows, starting column, nb columns) */
249 : static GEN
250 740 : lll_block(GEN m, long r0, long nr, long c0, long nc)
251 : {
252 : GEN block, u, uu;
253 740 : long nm = lg(m)-1, k;
254 740 : pari_sp av = avma;
255 :
256 740 : block = matslice(m, r0+1, r0+nr, c0+1, c0+nc);
257 740 : u = lll(block); if (lg(u) <= nc) return NULL;
258 728 : vecreverse_inplace(u); uu = matid(nm); /* embed in matid */
259 2492 : for (k = 1; k <= nc; k++) gel(uu,c0+k) = embedcol(gel(u,k),nm,c0);
260 728 : return gc_GEN(av, uu);
261 : }
262 :
263 : /* insert a column at index i, no copy */
264 : static GEN
265 2586 : shallowmatinsert(GEN m, GEN x, long i)
266 : {
267 2586 : long k, n = lg(m);
268 2586 : GEN mm = cgetg(n+1,t_MAT);
269 15095 : for (k=1; k < i; k++) gel(mm, k) = gel(m, k);
270 2586 : gel(mm, i) = x;
271 3940 : for (k=i; k < n; k++) gel(mm, k+1) = gel(m, k);
272 2586 : return mm;
273 : }
274 :
275 : static GEN
276 2586 : vec_v0(long n, long n0, long r1, long r2)
277 : {
278 : long k;
279 2586 : GEN C = zerocol(n);
280 5995 : for (k = 1; k <= r1; k++) gel(C, n0++) = gen_1;
281 5915 : for (k = 1; k <= r2; k++) gel(C, n0++) = gen_2;
282 2586 : return C;
283 : }
284 :
285 : /* select cm embeddings; return a matrix */
286 : /* TODO detect if precision was insufficient */
287 : static GEN
288 231 : cm_select(GEN nf, GEN cm, long prec)
289 : {
290 : GEN emb, keys, v, m_sel, imag_emb;
291 231 : long nalg, d_cm, r_cm, c, i, j, r2 = nf_get_r2(nf);
292 : pari_sp av;
293 :
294 231 : d_cm = degpol(gel(cm, 1)); /* degree of the cm field; even */
295 231 : nalg = d_cm / 2; /* nb of clusters */
296 231 : r_cm = nf_get_degree(nf) / d_cm; /* nb by cluster; nalg * r_cm = r2 */
297 231 : m_sel = zeromatcopy(nalg, r2); /* selection matrix */
298 231 : av = avma;
299 : /* group complex embeddings */
300 231 : emb = nfeltembed(nf, gel(cm, 2), NULL, prec);
301 : /* sort */
302 231 : imag_emb = imag_i(emb);
303 231 : keys = gadd(gmul(mppi(prec), real_i(emb)), gabs(imag_emb, prec));
304 231 : v = indexsort(keys);
305 :
306 637 : for (j = c = 1; c <= nalg; c++)
307 : {
308 406 : int ref = gsigne(gel(imag_emb, v[j]));
309 406 : gcoeff(m_sel, c, v[j]) = gen_1;
310 406 : j++;
311 518 : 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 231 : 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 637 : bestS(GEN bnf)
329 : {
330 637 : GEN v, S, hw, hv = bnf_get_no(bnf), DL, dl;
331 : long i, lS;
332 : ulong l;
333 : forprime_t P;
334 :
335 637 : 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 637 : gcharDLdata(GEN bnf, GEN S, GEN DL)
365 : {
366 637 : GEN M, h, Minv, Lalpha, t, dl, alpha, gen, cyc = bnf_get_cyc(bnf);
367 : long i;
368 637 : M = shallowmatconcat(DL);
369 637 : h = bnf_get_no(bnf);
370 637 : gen = bnf_get_gen(bnf);
371 :
372 : /* compute right inverse of M modulo cyc */
373 637 : M = shallowtrans(M);
374 637 : M = shallowmatconcat(mkcol2(M,diagonal_shallow(cyc)));
375 637 : Minv = matinvmod(M,h);
376 637 : Minv = vecslice(Minv,1,lg(Minv)-lg(cyc));
377 637 : Minv = shallowtrans(Minv);
378 :
379 637 : Lalpha = cgetg(lg(Minv),t_VEC);
380 966 : 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 637 : return mkvec2(Minv, Lalpha);
389 : }
390 :
391 : /* compute basis of characters; gc[1] generating family, as rows */
392 : GEN
393 658 : gcharinit(GEN bnf, GEN mod, long prec)
394 : {
395 658 : 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 658 : long evalprec = prec, nfprec, mprec, embprec;
400 :
401 658 : 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 658 : 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 567 : GEN fu = bnf_get_sunits(bnf);
441 567 : if (!fu) fu = bnf_get_fu(bnf); /* impose fundamental units */
442 567 : nf = shallowcopy(bnf_get_nf(bnf));
443 567 : nfprec = nf_get_prec(nf);
444 : }
445 :
446 : /* Dirichlet group + make sure mod contains archimedean places */
447 651 : mod = check_mod_factored(nf,mod,NULL,&fa2,&archp,NULL);
448 637 : sort_factor(fa2, (void*)&cmp_prime_ideal, &cmp_nodata);
449 637 : zm = localstar(nf, fa2, archp);
450 637 : zmcyc = locs_get_cyc(zm);
451 :
452 : /* set of primes S and valuations of generators */
453 637 : S = bestS(bnf);
454 637 : DLdata = gel(S,2);
455 637 : S = gel(S,1);
456 637 : DLdata = gcharDLdata(bnf, S, DLdata);
457 :
458 637 : nf_get_sign(nf, &r1, &r2);
459 637 : n = r1+2*r2;
460 637 : ns = lg(S) - 1;
461 637 : nu = r1+r2-1 + ns;
462 637 : nc = lg(zmcyc) - 1;
463 637 : nm = ns+nc+n; /* number of parameters = ns + nc + r1 + r2 + r2 */
464 :
465 : /* units and S-units */
466 637 : sfu = gel(bnfunits(bnf,S), 1);
467 637 : sfu = vec_shorten(sfu, nu); /* remove torsion */
468 :
469 : /* root of unity */
470 637 : order = bnf_get_tuN(bnf);
471 637 : z = bnf_get_tuU(bnf);
472 :
473 : /* maximal cm subfield */
474 637 : 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 637 : if (DEBUGLEVEL>2) err_printf("start matrix m\n");
496 637 : m = cgetg(nm + 1, t_MAT);
497 637 : mprec = nbits2prec(nm+10) + EXTRAPREC64;
498 637 : embprec = mprec;
499 : for(;;)
500 : {
501 1925 : for (k = 1; k <= nu; k++)
502 : { /* Lambda_S (S-units) then Lambda_f, fund. units */
503 1288 : logx = gchar_nflog(&nf, zm, S, gel(sfu,k), embprec);
504 1288 : if (!logx) break;
505 1288 : gel(m, k) = logx;
506 : }
507 637 : 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 1666 : for (k = 1; k <= nc; k++) /* Gamma, structure of (Z/m)* */
513 : {
514 1029 : C = zerocol(nm);
515 1029 : gel(C, ns+k) = gel(zmcyc, k);
516 1029 : gel(m, nu+k) = C;
517 : }
518 : /* zeta, root of unity */
519 637 : gel(m, nu+nc+1) = gchar_nflog(&nf, zm, S, z, mprec);
520 637 : shallow_clean_rat(gel(m, nu+nc+1), 1, nm, stoi(order), mprec);
521 1491 : for (k = 1; k <= r2; k++) /* embed Z^r_2 */
522 : {
523 854 : C = zerocol(nm);
524 854 : gel(C, ns+nc+r1+r2+k) = gen_1;
525 854 : gel(m, nu+nc+1+k) = C;
526 : }
527 637 : if (DEBUGLEVEL>1) err_printf("matrix m = %Ps\n", m);
528 :
529 637 : m0 = m;
530 637 : u0 = gen_0;
531 637 : rel = U = Ui = gen_0;
532 637 : cyc = gen_0;
533 637 : m_inv = gen_0;
534 :
535 : /* only m and m_inv depend on prec and are recomputed under gcharnewprec. */
536 637 : 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 637 : m = gchar_hnfreduce_shallow(gc, cm);
553 :
554 : /* C) compute snf basis of torsion subgroup */
555 637 : rel = shallowtrans(matslice(m, 1, ns+nc, 1, ns+nc));
556 637 : gchar_snfbasis_shallow(gc, rel);
557 :
558 : /* D) transpose inverse m_inv = (m0*u)~^-1 (may increase precision) */
559 637 : gcharmat_tinverse(gc, m, prec);
560 637 : return gc_GEN(av, gc);
561 : }
562 :
563 : /* b) do HNF reductions + LLL, keep base change u0 */
564 : static GEN
565 637 : gchar_hnfreduce_shallow(GEN gc, GEN cm)
566 : {
567 637 : GEN bnf = gchar_get_bnf(gc), nf = gchar_get_nf(gc), u, u0, m;
568 637 : long order, r1, r2, ns, nc, n, nu, nm, nalg = 0, mprec;
569 :
570 637 : nf_get_sign(nf, &r1, &r2);
571 637 : n = r1 + 2*r2;
572 637 : nu = r1 + r2 - 1;
573 637 : ns = gchar_get_ns(gc);
574 637 : nc = gchar_get_nc(gc);
575 637 : nm = ns+nc+n; /* ns + nc + r1 + r2 + r2 */
576 637 : order = 2*bnf_get_tuN(bnf);
577 637 : u0 = matid(nm);
578 637 : m = shallowcopy(gchar_get_m0(gc)); /* keep m0 unchanged */
579 637 : mprec = gprecision(m);
580 637 : if (DEBUGLEVEL>1) err_printf("matrix m = %Ps\n", m);
581 637 : if (nc)
582 : { /* keep steps 1&2 to make sure we have zeta_m */
583 315 : u = hnf_block(m, ns,nc, ns+nu,nc+1);
584 315 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
585 315 : if (DEBUGLEVEL>2) err_printf("step 1 -> %Ps\n", m);
586 315 : u = hnf_block(m, ns,nc, ns,nu+nc);
587 315 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
588 315 : if (DEBUGLEVEL>2) err_printf("step 2 -> %Ps\n", m);
589 : }
590 637 : if (r2)
591 : {
592 378 : u = hnf_block(m, nm-r2,r2, nm-r2-1,r2+1);
593 378 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
594 378 : if (DEBUGLEVEL>2) err_printf("step 3 -> %Ps\n", m);
595 : }
596 : /* remove last column */
597 637 : setlg(u0, nm); setlg(m, nm);
598 637 : if (DEBUGLEVEL>2) err_printf("remove last col -> %Ps\n", m);
599 :
600 637 : if (!gequal0(cm))
601 : {
602 : GEN v, Nargs;
603 231 : long bit = - prec2nbits(mprec) + 16 + expu(order);
604 : /* reduce on Norm arguments */
605 231 : v = cm_select(nf, cm, gchar_get_nfprec(gc));
606 231 : if (DEBUGLEVEL>2) err_printf("cm_select -> %Ps\n", v);
607 231 : nalg = nbrows(v);
608 231 : gchar_set_u0(gc, u0);
609 : for(;;)
610 14 : {
611 : long e, emax, i;
612 245 : Nargs = gmul(v, rowslice(m, nm-r2+1, nm));
613 245 : if (DEBUGLEVEL>2) err_printf("Nargs -> %Ps\n", Nargs);
614 245 : emax = bit-1;
615 1134 : for (i = ns+nc+1; i < lg(Nargs); i++)
616 : {
617 889 : gel(Nargs,i) = grndtoi(gmulgs(gel(Nargs,i), order), &e);
618 889 : emax = maxss(emax,e);
619 : }
620 245 : 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 231 : if (DEBUGLEVEL>2) err_printf("rounded Nargs -> %Ps\n", Nargs);
625 231 : u = hnf_block(Nargs, 0, nalg, ns+nc, n-1);
626 231 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
627 231 : 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 637 : gchar_set_nalg(gc, nalg);
632 637 : gchar_set_u0(gc, u0);
633 :
634 : /* TODO factor these two LLL reduction codes in a function? */
635 637 : if (nalg > 0)
636 : {
637 231 : GEN u = NULL;
638 : while (1)
639 : {
640 231 : if ((u = lll_block(m, ns+nc, n, ns+nc, nalg))) break;
641 0 : mprec = precdbl(mprec); m = gcharmatnewprec_shallow(gc, mprec);
642 : }
643 231 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
644 231 : if (DEBUGLEVEL>1) err_printf("after LLL reduction (CM block) -> %Ps\n", m);
645 : }
646 637 : gchar_set_u0(gc, u0);
647 :
648 637 : if (nu > 0)
649 : {
650 497 : GEN u = NULL;
651 : while (1)
652 : {
653 509 : 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 497 : u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
657 497 : if (DEBUGLEVEL>1) err_printf("after LLL reduction (trans block) -> %Ps\n", m);
658 : }
659 637 : 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 637 : gchar_snfbasis_shallow(GEN gc, GEN rel)
665 : {
666 : long n, r1, r2, lU, lUi;
667 : GEN nf, cyc, U, Ui;
668 :
669 637 : nf = gchar_get_nf(gc);
670 637 : nf_get_sign(nf, &r1, &r2);
671 637 : n = r1+2*r2;
672 :
673 637 : rel = ZM_hnf(rel);
674 637 : if (DEBUGLEVEL>1) err_printf("relations after hnf: %Ps\n", rel);
675 637 : cyc = ZM_snf_group(rel, &U, &Ui);
676 637 : if (lg(cyc)==1)
677 : {
678 294 : cyc = zerovec(n-1);
679 294 : U = shallowconcat(zeromat(n-1,lg(rel)-1),matid(n-1));
680 294 : 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 637 : rel = shallowtrans(rel);
690 637 : lU = lg(U);
691 637 : lUi = lg(Ui);
692 637 : U = shallowtrans(U);
693 637 : if (lg(U)==1 && lUi>1) U = zeromat(0,lUi-1);
694 637 : Ui = shallowtrans(Ui);
695 637 : if (lg(Ui)==1 && lU>1) Ui = zeromat(0,lU-1);
696 :
697 637 : gchar_set_nfree(gc, n-1);
698 637 : gchar_set_ntors(gc, (lg(cyc)-1) - (n-1));
699 637 : gchar_set_cyc(gc, shallowconcat(cyc, real_0(gchar_get_prec(gc))));
700 637 : gchar_set_HUUi(gc, rel, U, Ui);
701 637 : }
702 :
703 : static long
704 3006 : mextraprec(GEN m_inv, GEN m, GEN gc)
705 : {
706 3006 : return nbits2extraprec(2*maxss(gexpo(m_inv),1) + expu(lg(m))
707 3006 : + 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 1827 : gcharmat_tinverse(GEN gc, GEN m, long prec)
715 : {
716 : GEN m_inv;
717 1827 : long k, r1, r2, ns, nc, nalg, nm, mprec, bitprec = prec2nbits(prec);
718 :
719 1827 : nf_get_sign(gchar_get_nf(gc), &r1, &r2);
720 1827 : ns = gchar_get_ns(gc);
721 1827 : nc = gchar_get_nc(gc);
722 1827 : nalg = gchar_get_nalg(gc);
723 1827 : nm = ns + nc + r1 + 2*r2;
724 1827 : if (lg(m)==1) { gchar_set_basis(gc,zeromat(0,nm)); return; }
725 :
726 1806 : mprec = gprecision(m); /* possibly 0, if m is exact */
727 : while (1)
728 780 : {
729 2586 : long targetmprec = 0;
730 : GEN v0, mm;
731 : /* insert at column ns+nc+r1+r2, or last column if cm */
732 2586 : v0 = vec_v0(nm, ns+nc+1, r1, r2);
733 2586 : mm = shallowmatinsert(m, v0, nalg? nm: nm-r2);
734 2586 : if (DEBUGLEVEL>1) err_printf("add v0 -> %Ps\n", mm);
735 2586 : mm = shallowtrans(mm);
736 2586 : m_inv = RgM_inv(mm); /* invert matrix, may need to increase prec */
737 2586 : if (m_inv)
738 : {
739 2586 : if (DEBUGLEVEL>1) err_printf("inverse: %Ps\n",m_inv);
740 2586 : m_inv = vecsplice(m_inv, nalg? nm: nm-r2); /* remove v0 */
741 2586 : if (DEBUGLEVEL>1) err_printf("v0 removed: %Ps\n", m_inv);
742 2586 : m_inv = shallowtrans(m_inv);
743 2586 : if (!mprec) break;
744 : /* enough precision? */
745 : /* |B - A^(-1)| << |B|.|Id-B*A| */
746 2502 : if (gexpo(m_inv) + gexpo(gsub(RgM_mul(m_inv, m), gen_1)) + expu(lg(m))
747 2502 : <= -bitprec)
748 : {
749 : /* |A^(-1) - (A+H)^(-1)| << |H|.|A^(-1)|^2 */
750 1816 : targetmprec = prec + mextraprec(m_inv,m,gc);
751 1816 : if (mprec >= targetmprec) break;
752 : }
753 686 : else targetmprec = 0;
754 : }
755 780 : mprec = maxss(precdbl(mprec), targetmprec);
756 780 : if (mprec < DEFAULTPREC) mprec = DEFAULTPREC;
757 780 : m = gcharmatnewprec_shallow(gc, mprec); /* m0 * u0 to higher prec */
758 : }
759 : /* clean rationals */
760 1806 : if (nc)
761 : { /* reduce mod exponent of the group, TODO reduce on each component */
762 987 : GEN zmcyc = locs_get_cyc(gchar_get_zm(gc));
763 987 : GEN e = ZV_lcm(zmcyc);
764 3693 : for (k = 1; k <= nc; k++)
765 2706 : shallow_clean_rat(gel(m_inv, ns+k), 1, nm - 1, /*zmcyc[k]*/e, prec);
766 : }
767 1806 : if (r2)
768 : {
769 2750 : for (k = 1; k <= r2; k++)
770 1705 : shallow_clean_rat(gel(m_inv,nm-k+1), 1, nm-1, NULL, prec);
771 : }
772 1806 : if (nalg)
773 : {
774 : long i, j, e;
775 1651 : for (i = 1; i<=r1+r2; i++)
776 3190 : for (j = 1; j <= nalg; j++)
777 : {
778 2129 : e = gexpo(gcoeff(m_inv, ns+nc+j, ns+nc+i));
779 2129 : if (e > -20) /* TODO is this bound too permissive? */
780 : pari_err_BUG("gcharinit (nonzero entry)"); /*LCOV_EXCL_LINE*/
781 2129 : gcoeff(m_inv, ns+nc+j, ns+nc+i) = gen_0;
782 : }
783 : }
784 1806 : if (DEBUGLEVEL>1) err_printf("cyc/cm cleaned: %Ps", shallowtrans(m_inv));
785 : /* normalize characters, parameters mod Z */
786 5443 : for (k = 1; k <= ns+nc; k++) gel(m_inv, k) = gfrac(gel(m_inv, k));
787 : /* increase relative prec of real values */
788 1806 : gchar_set_basis(gc, gprec_w(m_inv, prec));
789 : }
790 :
791 : /* make sure same argument determination is chosen */
792 : static void
793 4795 : same_arg(GEN v1, GEN v2, long s1, long s2)
794 : {
795 4795 : long i1, i2, l = lg(v1);
796 23476 : for (i1 = s1, i2 = s2; i1 < l; i1++,i2++)
797 : {
798 18681 : GEN d = grndtoi(gsub(gel(v2,i2),gel(v1,i1)), NULL);
799 18681 : if (signe(d)) gel(v1,i1) = gadd(gel(v1,i1), d);
800 : }
801 4795 : }
802 :
803 : static void
804 4795 : vaffect_shallow(GEN x, long i0, GEN y)
805 : {
806 : long i;
807 40603 : for (i = 1; i < lg(y); i++)
808 35808 : gel(x, i+i0) = gel(y, i);
809 4795 : }
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 1996 : 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 1996 : ns = gchar_get_ns(gc);
821 1996 : nc = gchar_get_nc(gc);
822 1996 : nf = gchar_get_nf(gc);
823 1996 : sfu = gchar_get_sfu(gc);
824 1996 : nf_get_sign(nf, &r1, &r2);
825 1996 : nfprec = nf_get_prec(gchar_get_nf(gc));
826 :
827 1996 : m0 = gchar_get_m0(gc);
828 1996 : u0 = gchar_get_u0(gc);
829 :
830 1996 : if (DEBUGLEVEL>3) err_printf("gcharmatnewprec_shallow mprec=%d nfprec=%d\n", mprec, nfprec);
831 :
832 1996 : embprec = mprec; l = lg(sfu);
833 : while(1)
834 : { /* recompute the nfembedlogs of s-units and fundamental units */
835 6791 : for (k = 1; k < l; k++) /* ns s-units, then r1+r2-1 fundamental units */
836 : {
837 4795 : GEN emb = nfembedlog(&nf, gel(sfu,k), embprec);
838 4795 : if (!emb) break;
839 4795 : same_arg(emb, gel(m0,k), r1+r2, ns+nc+r1+r2);
840 4795 : vaffect_shallow(gel(m0, k), ns+nc, emb);
841 : }
842 1996 : 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 1996 : gchar_set_nf(gc, nf);
849 1996 : gchar_set_nfprec(gc, nfprec);
850 1996 : return RgM_ZM_mul(m0, u0);
851 : }
852 :
853 : static void _check_gchar_group(GEN gc, long flag);
854 : static void
855 11558 : 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 2107 : gcharnewprec_i(GEN gc, long newprec)
860 : {
861 : long prec, prec0, nfprec, nfprec0, mprec;
862 2107 : GEN gc2 = shallowcopy(gc);
863 :
864 2107 : _check_gchar_group(gc2, 1); /* ignore illegal prec */
865 2086 : prec = gchar_get_prec(gc2);
866 2086 : nfprec = gchar_get_nfprec(gc2);
867 :
868 2086 : if (newprec > prec)
869 : { /* increase precision */
870 1470 : if (DEBUGLEVEL) pari_warn(warnprec,"gcharnewprec",newprec);
871 1470 : nfprec += newprec - prec;
872 1470 : prec = newprec;
873 1470 : gchar_copy_precs(gc, gc2);
874 1470 : gchar_set_prec(gc2, prec);
875 1470 : gchar_set_nfprec(gc2, nfprec);
876 : }
877 :
878 2086 : nfprec0 = nf_get_prec(gchar_get_nf(gc2));
879 2086 : if (nfprec0 && nfprec > nfprec0)
880 : {
881 367 : if (DEBUGLEVEL) pari_warn(warnprec,"gcharnewprec (nf)", nfprec);
882 367 : gchar_set_nf(gc2, nfnewprec_shallow(gchar_get_nf(gc2), nfprec));
883 : }
884 :
885 2086 : prec0 = gprecision(gchar_get_basis(gc2));
886 2086 : if (prec0 && prec > prec0)
887 : {
888 1190 : GEN cyc, m, m0 = gchar_get_m0(gc);
889 1190 : if (DEBUGLEVEL) pari_warn(warnprec,"gcharnewprec (minv)", prec);
890 1190 : gchar_set_m0(gc2, shallowcopy(m0));
891 1190 : mprec = prec + mextraprec(gchar_get_basis(gc), m0, gc);
892 1190 : m = gcharmatnewprec_shallow(gc2, mprec);
893 1190 : if (DEBUGLEVEL>2) err_printf("m0*u0 recomputed -> %Ps\n", m);
894 1190 : gcharmat_tinverse(gc2, m, prec);
895 1190 : cyc = shallowcopy(gchar_get_cyc(gc2));
896 1190 : gel(cyc, lg(cyc)-1) = real_0(prec);
897 1190 : gchar_set_cyc(gc2, cyc);
898 : }
899 2086 : return gc2;
900 : }
901 :
902 : /* newprec: requested evalprec */
903 : GEN
904 2107 : gcharnewprec(GEN gc, long newprec)
905 : {
906 2107 : pari_sp av = avma;
907 2107 : GEN gc2 = gcharnewprec_i(gc, newprec + EXTRAPREC64);
908 2086 : gchar_set_evalprec(gc2, newprec);
909 2086 : return gc_GEN(av, gc2);
910 : }
911 :
912 : static void
913 13651 : check_localstar(GEN x)
914 : {
915 13651 : 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 13644 : }
919 :
920 : int
921 9038 : is_gchar_group(GEN gc)
922 : {
923 9038 : return (typ(gc) == t_VEC
924 9038 : && lg(gc) == GC_LENGTH + 1
925 519 : && typ(gel(gc, 8)) == t_VEC
926 519 : && lg(gel(gc, 8)) == 3
927 519 : && typ(gmael(gc, 8, 1)) == t_VECSMALL
928 519 : && typ(gmael(gc, 8, 2)) == t_VECSMALL
929 519 : && (checkbnf_i(gchar_get_bnf(gc)) != NULL)
930 18076 : && (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 13665 : _check_gchar_group(GEN gc, long flag)
937 : {
938 : GEN b, bnf, nf;
939 13665 : if (typ(gc) != t_VEC || lg(gc) != GC_LENGTH + 1)
940 14 : pari_err_TYPE("char group", gc);
941 13651 : check_localstar(gchar_get_zm(gc));
942 13644 : if (typ(gchar_get_loccyc(gc)) != t_VEC)
943 0 : pari_err_TYPE("gchar group (loccyc)", gc);
944 13644 : b = gchar_get_basis(gc);
945 13644 : if (typ(b) != t_MAT) pari_err_TYPE("gchar group (basis)", gc);
946 13637 : bnf = gchar_get_bnf(gc); checkbnf(bnf);
947 13637 : nf = gchar_get_nf(gc); checknf(nf);
948 13637 : 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 13637 : 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 13630 : if (!flag)
953 : {
954 11544 : long prec0 = gprecision(b), nfprec0 = nf_get_prec(nf);
955 11544 : if ((prec0 && gchar_get_prec(gc) > prec0)
956 11537 : || (nfprec0 && gchar_get_nfprec(gc) > nfprec0))
957 7 : pari_err_PREC("gchar group, please call gcharnewprec");
958 : }
959 13623 : }
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) retgc_const(av, cgetg(1, t_VEC));
1110 182 : return gc_GEN(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 9683 : check_gchar_i(GEN chi, long l, GEN *s)
1121 : {
1122 9683 : long i, n = lg(chi)-1;
1123 9683 : if (n == l)
1124 : { /* norm component */
1125 2941 : if (s)
1126 : {
1127 2920 : *s = gel(chi,l);
1128 2920 : switch(typ(*s))
1129 : {
1130 2913 : case t_INT:
1131 : case t_FRAC:
1132 : case t_REAL:
1133 2913 : case t_COMPLEX: break;
1134 7 : default: pari_err_TYPE("check_gchar_i [s]", *s);
1135 : }
1136 : }
1137 2934 : chi = vecslice(chi, 1, l-1);
1138 : }
1139 : else
1140 : { /* missing norm component */
1141 6742 : if (n != l-1) pari_err_DIM("check_gchar_i [chi]");
1142 6735 : if (s) *s = gen_0;
1143 : }
1144 46064 : 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 9662 : return chi;
1147 : }
1148 :
1149 : static GEN
1150 7197 : check_gchar(GEN gc, GEN chi, GEN *s)
1151 : {
1152 7197 : if (typ(chi)!=t_COL) pari_err_TYPE("check_gchar [chi]", chi);
1153 7176 : return check_gchar_i(chi, lg(gchar_get_cyc(gc))-1, s);
1154 : }
1155 :
1156 : static long
1157 2507 : safelgcols(GEN m)
1158 : {
1159 2507 : return lg(m)==1 ? 1 : lg(gel(m,1));
1160 : }
1161 :
1162 : static GEN
1163 2507 : check_gchari(GEN gc, GEN chi, GEN *s)
1164 : {
1165 2507 : if (typ(chi)!=t_VEC) pari_err_TYPE("check_gchari [chi]", chi);
1166 2507 : 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 7197 : gchar_internal(GEN gc, GEN chi, GEN *s)
1173 : {
1174 7197 : chi = check_gchar(gc, chi, s);
1175 7155 : return ZV_ZM_mul(chi, gchar_get_Ui(gc));
1176 : }
1177 :
1178 : static GEN
1179 24410 : RgV_frac_inplace(GEN v, long n)
1180 : {
1181 : long i;
1182 54013 : for (i = 1; i <= n; i++) gel(v,i) = gfrac(gel(v,i));
1183 24410 : return v;
1184 : }
1185 : /* from internal basis form, return the R/Z components */
1186 : static GEN
1187 9136 : gchari_duallog(GEN gc, GEN chi)
1188 : {
1189 9136 : chi = RgV_RgM_mul(chi, gchar_get_basis(gc)); /* take exponents mod Z */
1190 9136 : 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 1358 : conductor_expo_pr(GEN gens_fil, GEN chip)
1196 : {
1197 : long i;
1198 2345 : for (i = lg(gens_fil) - 1; i > 0; i--)
1199 : {
1200 1960 : GEN gens = gel(gens_fil, i);
1201 : long j;
1202 3227 : for (j = 1; j < lg(gens); j++)
1203 : {
1204 2240 : GEN v = gmul(chip, gel(gens, j));
1205 2240 : 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 1856 : gcharlog_conductor_f(GEN gc, GEN chi, GEN *faN)
1214 : {
1215 : GEN zm, P, E, Lsprk, ufil;
1216 : long i, l, ic;
1217 :
1218 1856 : if (gchar_get_nc(gc) == 0)
1219 : {
1220 568 : if (faN) *faN = trivial_fact();
1221 568 : return gen_1;
1222 : }
1223 1288 : zm = gchar_get_zm(gc);
1224 1288 : Lsprk = locs_get_Lsprk(zm);
1225 1288 : ufil = locs_get_Lgenfil(zm);
1226 1288 : P = gel(locs_get_famod(zm), 1);
1227 1288 : l = lg(Lsprk); E = cgetg(l, t_COL);
1228 2646 : for (i = 1, ic = gchar_get_ns(gc); i < l ; i++)
1229 : {
1230 1358 : GEN gens = gel(ufil, i), chip;
1231 1358 : long ncp = sprk_get_ncp(gel(Lsprk,i));
1232 :
1233 1358 : chip = vecslice(chi, ic + 1, ic + ncp);
1234 1358 : gel(E, i) = conductor_expo_pr(gens, chip);
1235 1358 : ic += ncp;
1236 : }
1237 1288 : if (faN) *faN = famat_remove_trivial(mkmat2(P,E));
1238 1288 : return idealfactorback(gchar_get_nf(gc), P, E, 0); /* red = 0 */
1239 : }
1240 :
1241 : /* ={sigma} s.t. k_sigma = 1 */
1242 : static GEN
1243 17130 : gcharlog_conductor_oo(GEN gc, GEN chi)
1244 : {
1245 17130 : long noo, i, n0 = gchar_get_ns(gc) + gchar_get_nc(gc);
1246 : GEN k_real, chi_oo, moo, den;
1247 :
1248 17130 : moo = locs_get_m_infty(gchar_get_zm(gc));
1249 17130 : noo = lg(moo)-1;
1250 17130 : k_real = vecslice(chi, n0-noo+1, n0);
1251 17130 : chi_oo = zerovec(gchar_get_r1(gc));
1252 18740 : 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 17130 : 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 gc_GEN(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 = gc_GEN(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 gc_GEN(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 gc_GEN(av, shallowconcat1(mkcol2(logchi,s)));
1415 : }
1416 :
1417 : static GEN
1418 192479 : gcharisprincipal(GEN gc, GEN x, GEN *palpha)
1419 : {
1420 192479 : GEN t, v, bnf = gchar_get_bnf(gc), DLdata = gchar_get_DLdata(gc);
1421 192479 : t = bnfisprincipal0(bnf, x, nf_GENMAT | nf_FORCE); v = gel(t, 1);
1422 192479 : *palpha = famat_reduce(famat_mul(nffactorback(bnf, gel(DLdata,2), v),
1423 192479 : gel(t, 2)));
1424 192479 : 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 192479 : gchar_log(GEN gc, GEN x, GEN logchi, long prec)
1431 : {
1432 192479 : GEN zm, v, alpha, arch_log = NULL, zm_log, nf;
1433 :
1434 192479 : nf = gchar_get_nf(gc);
1435 192479 : zm = gchar_get_zm(gc);
1436 192479 : v = gcharisprincipal(gc, x, &alpha); /* alpha a GENMAT */
1437 192479 : if (DEBUGLEVEL>2) err_printf("v %Ps\n", v);
1438 192479 : zm_log = gchar_logm(nf, zm, alpha);
1439 192479 : if (DEBUGLEVEL>2) err_printf("zm_log(alpha) %Ps\n", zm_log);
1440 :
1441 192479 : if (logchi)
1442 : { /* check if precision is sufficient, take care of gexpo = -infty */
1443 188692 : long e, bit = expu(nf_get_degree(nf) + lg(zm_log)-1) + 3;
1444 188692 : e = gexpo(logchi); if (e > 0) bit += e;
1445 188692 : e = gexpo(gel(alpha,1)); if (e > 0) bit += e;
1446 188692 : prec += nbits2extraprec(bit);
1447 : }
1448 : for(;;)
1449 : {
1450 192479 : arch_log = nfembedlog(&nf, alpha, prec);
1451 192479 : if (arch_log) break;
1452 0 : prec = precdbl(prec);
1453 : }
1454 192479 : if (DEBUGLEVEL>2) err_printf("arch log %Ps\n", arch_log);
1455 192479 : 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 gc_GEN(av, vec_append(gchar_log(gc, x, NULL, prec), norm));
1469 : }
1470 :
1471 : static GEN
1472 203350 : gcharlog_eval_raw(GEN logchi, GEN logx)
1473 203350 : { 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 188692 : gchari_eval(GEN gc, GEN chi, GEN x, long flag, GEN logchi, GEN s, long prec)
1481 : {
1482 : GEN z, norm, logx;
1483 188692 : 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 188692 : logx = gchar_log(gc, x, logchi, prec);
1495 188692 : norm = gequal0(s)? NULL: idealnorm(gchar_get_nf(gc), x);
1496 188692 : z = gcharlog_eval_raw(logchi, logx);
1497 188692 : if (flag)
1498 : {
1499 184898 : z = expIPiC(gmul2n(z, 1), prec);
1500 184898 : 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 188692 : if (DEBUGLEVEL>1) err_printf("char value %Ps\n", z);
1505 188692 : 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 gc_upto(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 gc_GEN(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 2498657 : gaddmul(GEN x, GEN y, GEN z)
1745 : {
1746 : pari_sp av;
1747 2498657 : if (typ(z) == t_INT)
1748 : {
1749 2210828 : if (!signe(z)) return x;
1750 13948 : if (equali1(z)) return gadd(x,y);
1751 : }
1752 296926 : if (isintzero(x)) return gmul(y,z);
1753 118860 : av = avma;
1754 118860 : return gc_upto(av, gadd(x, gmul(y,z)));
1755 : }
1756 :
1757 : GEN
1758 840 : vecan_gchar(GEN an, long n, long prec)
1759 : {
1760 : forprime_t iter;
1761 840 : GEN gc = gel(an,1), chi = gel(an,2), P = gel(an,3), PZ = gel(an,4);
1762 840 : GEN BOUND = stoi(n), v = vec_ei(n, 1);
1763 840 : GEN gp = cgetipos(3), nf, chilog, s;
1764 : ulong p;
1765 :
1766 : /* prec increase: 1/n*log(N(pmax)) < log(pmax) */
1767 840 : if (DEBUGLEVEL > 1)
1768 0 : err_printf("vecan_gchar: need extra prec %ld\n", nbits2extraprec(expu(n)));
1769 840 : gc = gcharnewprec(gc, prec + nbits2extraprec(expu(n)));
1770 840 : chilog = gchari_duallog(gc, check_gchari(gc, chi, &s));
1771 :
1772 840 : 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 840 : u_forprime_init(&iter, 2, n);
1782 186291 : while ((p = u_forprime_next(&iter)))
1783 : {
1784 : GEN L;
1785 : long j;
1786 185451 : int check = !umodiu(PZ,p);
1787 185451 : gp[2] = p;
1788 185451 : L = idealprimedec_limit_norm(nf, gp, BOUND);
1789 370496 : for (j = 1; j < lg(L); j++)
1790 : {
1791 185045 : GEN pr = gel(L, j), ch;
1792 : pari_sp av;
1793 : ulong k, q;
1794 185045 : if (check && gen_search(P, pr, (void*)cmp_prime_ideal, cmp_nodata) > 0)
1795 448 : continue;
1796 : /* TODO: extract code and use precom sprk? */
1797 184597 : av = avma;
1798 184597 : ch = gchari_eval(gc, chi, pr, 1, chilog, gen_0, prec);
1799 184597 : ch = gc_upto(av, ch);
1800 184597 : q = upr_norm(pr);
1801 184597 : gel(v, q) = gadd(gel(v, q), ch);
1802 2683254 : for (k = 2*q; k <= (ulong)n; k += q)
1803 2498657 : gel(v, k) = gaddmul(gel(v, k), ch, gel(v, k/q));
1804 : }
1805 : }
1806 : /* weight, could consider shifting s at eval instead */
1807 840 : if (!gequal0(s))
1808 : {
1809 462 : GEN ns = dirpowers(n, gneg(s), prec);
1810 : long j;
1811 292901 : for (j = 1; j <= n; j++)
1812 292439 : if (gel(v,j) != gen_0) gel(v, j) = gmul(gel(v,j),gel(ns,j));
1813 : }
1814 840 : return v;
1815 : }
1816 :
1817 : GEN
1818 28 : eulerf_gchar(GEN an, GEN p, long prec)
1819 : {
1820 28 : 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 28 : if (DEBUGLEVEL > 1)
1827 0 : err_printf("vecan_gchar: need extra prec %ld\n", nbits2extraprec(expi(p)));
1828 28 : gc = gcharnewprec(gc, prec + nbits2extraprec(expi(p)));
1829 28 : chilog = gchari_duallog(gc, check_gchari(gc, chi, &s));
1830 :
1831 28 : nf = gchar_get_nf(gc);
1832 28 : f = pol_1(0);
1833 28 : check = dvdii(PZ, p);
1834 28 : L = idealprimedec(nf, p); l = lg(L);
1835 63 : for (j = 1; j < l; j++)
1836 : {
1837 35 : GEN pr = gel(L, j), ch;
1838 35 : if (check && gen_search(P, pr, (void*)cmp_prime_ideal, cmp_nodata) > 0)
1839 0 : continue;
1840 35 : ch = gchari_eval(gc, chi, pr, 1, chilog, s, prec);
1841 35 : f = gmul(f, gsub(gen_1, monomial(ch, pr_get_f(pr), 0)));
1842 : }
1843 28 : return mkrfrac(gen_1,f);
1844 : }
1845 :
1846 : static GEN
1847 1590 : cleanup_vga(GEN vga, long prec)
1848 : {
1849 : GEN ind;
1850 : long bitprec, i, l;
1851 1590 : if (!prec) return vga; /* already exact */
1852 1590 : bitprec = prec2nbits(prec);
1853 1590 : vga = shallowcopy(vga); l = lg(vga);
1854 5981 : for (i = 1; i < l; i++)
1855 : {
1856 4391 : GEN z = gel(vga,i);
1857 4391 : if (typ(z) != t_COMPLEX) continue;
1858 2564 : if (gexpo(gel(z,2)) < -bitprec+20) gel(vga,i) = gel(z,1);
1859 : }
1860 1590 : ind = indexsort(imag_i(vga));
1861 4391 : for (i = 2; i < l; i++)
1862 : {
1863 2801 : GEN z = gel(vga,ind[i]), t;
1864 2801 : if (typ(z) != t_COMPLEX) continue;
1865 1695 : t = imag_i(gel(vga, ind[i-1]));
1866 1695 : if (gexpo(gsub(gel(z,2), t)) < -bitprec+20)
1867 364 : gel(vga, ind[i]) = mkcomplex(gel(z,1), t);
1868 : }
1869 5981 : for (i = 1; i < l; i++)
1870 : {
1871 4391 : GEN z = gel(vga,i);
1872 4391 : if (typ(z) != t_COMPLEX) continue;
1873 2564 : gel(vga, i) = mkcomplex(gel(z,1), bestappr(gel(z,2), int2n(bitprec/2)));
1874 : }
1875 1590 : return vga;
1876 : }
1877 :
1878 : /* TODO: move to lfunutils, use lfunzeta and lfunzetak */
1879 : GEN
1880 1639 : 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 1639 : nf = gchar_get_nf(gc);
1887 1639 : ns = gchar_get_ns(gc);
1888 1639 : nc = gchar_get_nc(gc);
1889 1639 : nm = gchar_get_nm(gc);
1890 1639 : nf_get_sign(nf, &r1, &r2);
1891 1639 : chi = check_gchari(gc, chi, &s);
1892 1639 : chilog = gchari_duallog(gc, chi);
1893 1639 : s = gadd(s0,s); chiw = shallowconcat(chi, s);
1894 1639 : if (!gequal0(imag_i(s)))
1895 7 : pari_err_IMPL("lfun for gchar with imaginary norm component");
1896 1632 : cond_f = gcharlog_conductor_f(gc, chilog, &faN);
1897 1632 : P = gel(faN, 1); /* prime ideals dividing cond(chi) */
1898 1632 : cond_oo = gcharlog_conductor_oo(gc, chilog);
1899 :
1900 1632 : NN = mulii(idealnorm(nf, cond_f), absi_shallow(nf_get_disc(nf)));
1901 1632 : if (equali1(NN)) return lfunshift(lfuncreate(gen_1), gneg(s), 0,
1902 : prec2nbits(gchar_get_evalprec(gc)));
1903 1611 : 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 1590 : v_phi = gmul(vecslice(chilog, ns+nc+1, ns+nc+r1+r2), gen_I());
1909 1590 : v_arg = gdivgs(gabs(vecslice(chilog,ns+nc+r1+r2+1,nm),BITS_IN_LONG), 2);
1910 1590 : vga_r = gadd(vecslice(v_phi, 1, r1), cond_oo);
1911 1590 : vga_c = gadd(vecslice(v_phi, r1+1, r1+r2), v_arg);
1912 1590 : sig = shallowconcat1(mkvec3(vga_r,vga_c,gadd(vga_c,const_vec(r2,gen_1))));
1913 : /* TODO: remove cleanup when gammamellinv takes ldata*/
1914 1590 : sig = cleanup_vga(sig, gchar_get_prec(gc));
1915 1590 : k = gen_1;
1916 1590 : if (!gequal0(s))
1917 : {
1918 : long j;
1919 2803 : for (j = 1; j < lg(sig); j++) gel(sig, j) = gadd(gel(sig, j), s);
1920 778 : k = gsub(k, gmulgs(s,2));
1921 : }
1922 :
1923 : #define tag(x,t) mkvec2(mkvecsmall(t),x)
1924 1590 : 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 519 : lfungchar(GEN gc, GEN chi)
1930 : {
1931 519 : pari_sp av = avma;
1932 : GEN s;
1933 519 : check_gchar_group(gc);
1934 519 : chi = gchar_internal(gc, chi, &s);
1935 505 : return gc_GEN(av, gchari_lfun(gc, chi, s));
1936 : }
|
