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
_fq_zech_mpoly_scalar_addmul_fq_zech(fq_zech_struct * Acoeffs,ulong * Aexps,fq_zech_struct * Bcoeffs,const ulong * Bexps,slong Blen,fq_zech_struct * Ccoeffs,const ulong * Cexps,slong Clen,const fq_zech_t d,slong N,const ulong * cmpmask,const fq_zech_ctx_t fqctx)14 static slong _fq_zech_mpoly_scalar_addmul_fq_zech(
15 fq_zech_struct * Acoeffs, ulong * Aexps,
16 fq_zech_struct * Bcoeffs, const ulong * Bexps, slong Blen,
17 fq_zech_struct * Ccoeffs, const ulong * Cexps, slong Clen,
18 const fq_zech_t d,
19 slong N,
20 const ulong * cmpmask,
21 const fq_zech_ctx_t fqctx)
22 {
23 slong i = 0, j = 0, k = 0;
24 fq_zech_t p;
25
26 FLINT_ASSERT(!fq_zech_is_zero(d, fqctx));
27
28 fq_zech_init(p, fqctx);
29
30 while (i < Blen && j < Clen)
31 {
32 int cmp = mpoly_monomial_cmp(Bexps + i*N, Cexps + j*N, N, cmpmask);
33
34 if (cmp > 0)
35 {
36 mpoly_monomial_set(Aexps + k*N, Bexps + i*N, N);
37 fq_zech_set(Acoeffs + k, Bcoeffs + i, fqctx);
38 i++;
39 k++;
40 }
41 else if (cmp == 0)
42 {
43 mpoly_monomial_set(Aexps + k*N, Bexps + i*N, N);
44 fq_zech_mul(p, Ccoeffs + j, d, fqctx);
45 fq_zech_add(Acoeffs + k, Bcoeffs + i, p, fqctx);
46 k += !fq_zech_is_zero(Acoeffs + k, fqctx);
47 i++;
48 j++;
49 }
50 else
51 {
52 mpoly_monomial_set(Aexps + k*N, Cexps + j*N, N);
53 fq_zech_mul(Acoeffs + k, Ccoeffs + j, d, fqctx);
54 j++;
55 k++;
56 }
57 }
58
59 while (i < Blen)
60 {
61 mpoly_monomial_set(Aexps + k*N, Bexps + i*N, N);
62 fq_zech_set(Acoeffs + k, Bcoeffs + i, fqctx);
63 i++;
64 k++;
65 }
66
67 while (j < Clen)
68 {
69 mpoly_monomial_set(Aexps + k*N, Cexps + j*N, N);
70 fq_zech_mul(Acoeffs + k, Ccoeffs + j, d, fqctx);
71 j++;
72 k++;
73 }
74
75 fq_zech_clear(p, fqctx);
76
77 return k;
78 }
79
fq_zech_mpoly_scalar_addmul_fq_zech(fq_zech_mpoly_t A,const fq_zech_mpoly_t B,const fq_zech_mpoly_t C,const fq_zech_t d,const fq_zech_mpoly_ctx_t ctx)80 void fq_zech_mpoly_scalar_addmul_fq_zech(
81 fq_zech_mpoly_t A,
82 const fq_zech_mpoly_t B,
83 const fq_zech_mpoly_t C,
84 const fq_zech_t d,
85 const fq_zech_mpoly_ctx_t ctx)
86 {
87 flint_bitcnt_t Abits;
88 slong N;
89 ulong * Bexps = B->exps, * Cexps = C->exps;
90 ulong * cmpmask;
91 int freeBexps = 0, freeCexps = 0;
92 TMP_INIT;
93
94 if (fq_zech_mpoly_is_zero(B, ctx))
95 {
96 fq_zech_mpoly_scalar_mul_fq_zech(A, C, d, ctx);
97 return;
98 }
99 else if (fq_zech_mpoly_is_zero(C, ctx) || fq_zech_is_zero(d, ctx->fqctx))
100 {
101 fq_zech_mpoly_set(A, B, ctx);
102 return;
103 }
104
105 TMP_START;
106 Abits = FLINT_MAX(B->bits, C->bits);
107 N = mpoly_words_per_exp(Abits, ctx->minfo);
108 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
109 mpoly_get_cmpmask(cmpmask, N, Abits, ctx->minfo);
110
111 if (Abits != B->bits)
112 {
113 freeBexps = 1;
114 Bexps = (ulong *) flint_malloc(N*B->length*sizeof(ulong));
115 mpoly_repack_monomials(Bexps, Abits, B->exps, B->bits, B->length, ctx->minfo);
116 }
117
118 if (Abits != C->bits)
119 {
120 freeCexps = 1;
121 Cexps = (ulong *) flint_malloc(N*C->length*sizeof(ulong));
122 mpoly_repack_monomials(Cexps, Abits, C->exps, C->bits, C->length, ctx->minfo);
123 }
124
125 if (A == B || A == C)
126 {
127 fq_zech_mpoly_t T;
128 fq_zech_mpoly_init3(T, B->length + C->length, Abits, ctx);
129 T->length = _fq_zech_mpoly_scalar_addmul_fq_zech(T->coeffs, T->exps,
130 B->coeffs, Bexps, B->length,
131 C->coeffs, Cexps, C->length,
132 d, N, cmpmask, ctx->fqctx);
133 fq_zech_mpoly_swap(A, T, ctx);
134 fq_zech_mpoly_clear(T, ctx);
135 }
136 else
137 {
138 fq_zech_mpoly_fit_length_reset_bits(A, B->length + C->length, Abits, ctx);
139 A->length = _fq_zech_mpoly_scalar_addmul_fq_zech(A->coeffs, A->exps,
140 B->coeffs, Bexps, B->length,
141 C->coeffs, Cexps, C->length,
142 d, N, cmpmask, ctx->fqctx);
143 }
144
145 if (freeBexps)
146 flint_free(Bexps);
147
148 if (freeCexps)
149 flint_free(Cexps);
150
151 TMP_END;
152 }
153