1 /*
2 Copyright (C) 2020 Daniel Schultz
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <https://www.gnu.org/licenses/>.
10 */
11
12 #include "fq_zech_mpoly.h"
13
14 /* exponents of B are not multiprecision */
_fq_zech_mpoly_evaluate_one_fq_zech_sp(fq_zech_mpoly_t A,const fq_zech_mpoly_t B,slong var,const fq_zech_t val,const fq_zech_mpoly_ctx_t ctx)15 static void _fq_zech_mpoly_evaluate_one_fq_zech_sp(
16 fq_zech_mpoly_t A,
17 const fq_zech_mpoly_t B,
18 slong var,
19 const fq_zech_t val,
20 const fq_zech_mpoly_ctx_t ctx)
21 {
22 slong i, N, off, shift;
23 ulong * cmpmask, * one;
24 slong Blen = B->length;
25 const fq_zech_struct * Bcoeffs = B->coeffs;
26 const ulong * Bexps = B->exps;
27 flint_bitcnt_t bits = B->bits;
28 slong Alen;
29 fq_zech_struct * Acoeffs;
30 ulong * Aexps;
31 ulong mask, k;
32 int need_sort = 0, cmp;
33 fq_zech_t pp;
34 TMP_INIT;
35
36 FLINT_ASSERT(B->bits <= FLINT_BITS);
37
38 TMP_START;
39
40 fq_zech_init(pp, ctx->fqctx);
41
42 fq_zech_mpoly_fit_length_reset_bits(A, Blen, bits, ctx);
43 Acoeffs = A->coeffs;
44 Aexps = A->exps;
45
46 mask = (-UWORD(1)) >> (FLINT_BITS - bits);
47 N = mpoly_words_per_exp_sp(bits, ctx->minfo);
48 one = (ulong*) TMP_ALLOC(N*sizeof(ulong));
49 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
50 mpoly_gen_monomial_offset_shift_sp(one, &off, &shift, var, bits, ctx->minfo);
51 mpoly_get_cmpmask(cmpmask, N, bits, ctx->minfo);
52
53 Alen = 0;
54 for (i = 0; i < Blen; i++)
55 {
56 k = (Bexps[N*i + off] >> shift) & mask;
57 fq_zech_pow_ui(pp, val, k, ctx->fqctx);
58 fq_zech_mul(Acoeffs + Alen, Bcoeffs + i, pp, ctx->fqctx);
59 if (fq_zech_is_zero(Acoeffs + Alen, ctx->fqctx))
60 continue;
61
62 mpoly_monomial_msub(Aexps + N*Alen, Bexps + N*i, k, one, N);
63 if (Alen < 1)
64 {
65 Alen = 1;
66 continue;
67 }
68 cmp = mpoly_monomial_cmp(Aexps + N*(Alen - 1), Aexps + N*Alen, N, cmpmask);
69 if (cmp != 0)
70 {
71 need_sort |= (cmp < 0);
72 Alen++;
73 continue;
74 }
75 fq_zech_add(Acoeffs + Alen - 1, Acoeffs + Alen - 1, Acoeffs + Alen, ctx->fqctx);
76 Alen -= fq_zech_is_zero(Acoeffs + Alen - 1, ctx->fqctx);
77 }
78 A->length = Alen;
79
80 fq_zech_clear(pp, ctx->fqctx);
81
82 TMP_END;
83
84 if (need_sort)
85 {
86 fq_zech_mpoly_sort_terms(A, ctx);
87 fq_zech_mpoly_combine_like_terms(A, ctx);
88 }
89
90 FLINT_ASSERT(fq_zech_mpoly_is_canonical(A, ctx));
91 }
92
93
94 /* exponents of B are multiprecision */
_fq_zech_mpoly_evaluate_one_fq_zech_mp(fq_zech_mpoly_t A,const fq_zech_mpoly_t B,slong var,const fq_zech_t val,const fq_zech_mpoly_ctx_t ctx)95 static void _fq_zech_mpoly_evaluate_one_fq_zech_mp(
96 fq_zech_mpoly_t A,
97 const fq_zech_mpoly_t B,
98 slong var,
99 const fq_zech_t val,
100 const fq_zech_mpoly_ctx_t ctx)
101 {
102 slong i, N, off;
103 ulong * cmpmask, * one, * tmp;
104 slong Blen = B->length;
105 const fq_zech_struct * Bcoeffs = B->coeffs;
106 const ulong * Bexps = B->exps;
107 flint_bitcnt_t bits = B->bits;
108 slong Alen;
109 fq_zech_struct * Acoeffs;
110 ulong * Aexps;
111 fmpz_t k;
112 int need_sort = 0, cmp;
113 fq_zech_t pp;
114 TMP_INIT;
115
116 FLINT_ASSERT(B->bits > FLINT_BITS);
117 FLINT_ASSERT((B->bits % FLINT_BITS) == 0);
118
119 TMP_START;
120
121 fmpz_init(k);
122 fq_zech_init(pp, ctx->fqctx);
123
124 fq_zech_mpoly_fit_length_reset_bits(A, Blen, bits, ctx);
125 Acoeffs = A->coeffs;
126 Aexps = A->exps;
127
128 N = mpoly_words_per_exp_mp(bits, ctx->minfo);
129 one = (ulong *) TMP_ALLOC(3*N*sizeof(ulong));
130 cmpmask = one + N;
131 tmp = cmpmask + N;
132 off = mpoly_gen_monomial_offset_mp(one, var, bits, ctx->minfo);
133 mpoly_get_cmpmask(cmpmask, N, bits, ctx->minfo);
134
135 Alen = 0;
136 for (i = 0; i < Blen; i++)
137 {
138 fmpz_set_ui_array(k, Bexps + N*i + off, bits/FLINT_BITS);
139 fq_zech_pow(pp, val, k, ctx->fqctx);
140 fq_zech_mul(Acoeffs + Alen, Bcoeffs + i, pp, ctx->fqctx);
141 if (fq_zech_is_zero(Acoeffs + Alen, ctx->fqctx))
142 continue;
143
144 mpoly_monomial_mul_fmpz(tmp, one, N, k);
145 mpoly_monomial_sub_mp(Aexps + N*Alen, Bexps + N*i, tmp, N);
146 if (Alen < 1)
147 {
148 Alen = 1;
149 continue;
150 }
151 cmp = mpoly_monomial_cmp(Aexps + N*(Alen - 1), Aexps + N*Alen, N, cmpmask);
152 if (cmp != 0)
153 {
154 need_sort |= (cmp < 0);
155 Alen++;
156 continue;
157 }
158 fq_zech_add(Acoeffs + Alen - 1, Acoeffs + Alen - 1, Acoeffs + Alen, ctx->fqctx);
159 Alen -= fq_zech_is_zero(Acoeffs + Alen - 1, ctx->fqctx);
160 }
161 A->length = Alen;
162
163 fq_zech_clear(pp, ctx->fqctx);
164 fmpz_clear(k);
165
166 TMP_END;
167
168 if (need_sort)
169 {
170 fq_zech_mpoly_sort_terms(A, ctx);
171 fq_zech_mpoly_combine_like_terms(A, ctx);
172 }
173
174 FLINT_ASSERT(fq_zech_mpoly_is_canonical(A, ctx));
175 }
176
177
fq_zech_mpoly_evaluate_one_fq_zech(fq_zech_mpoly_t A,const fq_zech_mpoly_t B,slong var,const fq_zech_t val,const fq_zech_mpoly_ctx_t ctx)178 void fq_zech_mpoly_evaluate_one_fq_zech(
179 fq_zech_mpoly_t A,
180 const fq_zech_mpoly_t B,
181 slong var,
182 const fq_zech_t val,
183 const fq_zech_mpoly_ctx_t ctx)
184 {
185 if (fq_zech_mpoly_is_zero(B, ctx))
186 {
187 fq_zech_mpoly_zero(A, ctx);
188 return;
189 }
190
191 if (B->bits <= FLINT_BITS)
192 _fq_zech_mpoly_evaluate_one_fq_zech_sp(A, B, var, val, ctx);
193 else
194 _fq_zech_mpoly_evaluate_one_fq_zech_mp(A, B, var, val, ctx);
195 }
196
197