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