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