1 /*
2 Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
10 */
11
12 #include "fq_nmod_mpoly.h"
13
_fq_nmod_mpoly_sub(fq_nmod_struct * coeff1,ulong * exp1,fq_nmod_struct * coeff2,const ulong * exp2,slong len2,fq_nmod_struct * coeff3,const ulong * exp3,slong len3,slong N,const ulong * cmpmask,const fq_nmod_ctx_t fqctx)14 slong _fq_nmod_mpoly_sub(fq_nmod_struct * coeff1, ulong * exp1,
15 fq_nmod_struct * coeff2, const ulong * exp2, slong len2,
16 fq_nmod_struct * coeff3, const ulong * exp3, slong len3,
17 slong N, const ulong * cmpmask, const fq_nmod_ctx_t fqctx)
18 {
19 slong i = 0, j = 0, k = 0;
20
21 while (i < len2 && j < len3)
22 {
23 int cmp = mpoly_monomial_cmp(exp2 + i*N, exp3 + j*N, N, cmpmask);
24
25 if (cmp > 0)
26 {
27 mpoly_monomial_set(exp1 + k*N, exp2 + i*N, N);
28 fq_nmod_set(coeff1 + k, coeff2 + i, fqctx);
29 i++;
30 } else if (cmp == 0)
31 {
32 mpoly_monomial_set(exp1 + k*N, exp2 + i*N, N);
33 fq_nmod_sub(coeff1 + k, coeff2 + i, coeff3 + j, fqctx);
34 k -= fq_nmod_is_zero(coeff1 + k, fqctx);
35 i++;
36 j++;
37 } else
38 {
39 mpoly_monomial_set(exp1 + k*N, exp3 + j*N, N);
40 fq_nmod_neg(coeff1 + k, coeff3 + j, fqctx);
41 j++;
42 }
43 k++;
44 }
45
46 while (i < len2)
47 {
48 mpoly_monomial_set(exp1 + k*N, exp2 + i*N, N);
49 fq_nmod_set(coeff1 + k, coeff2 + i, fqctx);
50 i++;
51 k++;
52 }
53
54 while (j < len3)
55 {
56 mpoly_monomial_set(exp1 + k*N, exp3 + j*N, N);
57 fq_nmod_neg(coeff1 + k, coeff3 + j, fqctx);
58 j++;
59 k++;
60 }
61
62 return k;
63 }
64
fq_nmod_mpoly_sub(fq_nmod_mpoly_t poly1,const fq_nmod_mpoly_t poly2,const fq_nmod_mpoly_t poly3,const fq_nmod_mpoly_ctx_t ctx)65 void fq_nmod_mpoly_sub(fq_nmod_mpoly_t poly1, const fq_nmod_mpoly_t poly2,
66 const fq_nmod_mpoly_t poly3, const fq_nmod_mpoly_ctx_t ctx)
67 {
68 slong len1 = 0, max_bits, N;
69 ulong * exp2 = poly2->exps, * exp3 = poly3->exps;
70 ulong * cmpmask;
71 int free2 = 0, free3 = 0;
72 TMP_INIT;
73
74 max_bits = FLINT_MAX(poly2->bits, poly3->bits);
75 N = mpoly_words_per_exp(max_bits, ctx->minfo);
76
77 if (poly2->length == 0)
78 {
79 fq_nmod_mpoly_neg(poly1, poly3, ctx);
80 return;
81 }
82 else if (poly3->length == 0)
83 {
84 fq_nmod_mpoly_set(poly1, poly2, ctx);
85 return;
86 }
87
88 TMP_START;
89 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
90 mpoly_get_cmpmask(cmpmask, N, max_bits, ctx->minfo);
91
92 if (max_bits > poly2->bits)
93 {
94 free2 = 1;
95 exp2 = (ulong *) flint_malloc(N*poly2->length*sizeof(ulong));
96 mpoly_repack_monomials(exp2, max_bits, poly2->exps, poly2->bits,
97 poly2->length, ctx->minfo);
98 }
99
100 if (max_bits > poly3->bits)
101 {
102 free3 = 1;
103 exp3 = (ulong *) flint_malloc(N*poly3->length*sizeof(ulong));
104 mpoly_repack_monomials(exp3, max_bits, poly3->exps, poly3->bits,
105 poly3->length, ctx->minfo);
106 }
107
108 if (poly1 == poly2 || poly1 == poly3)
109 {
110 fq_nmod_mpoly_t temp;
111
112 fq_nmod_mpoly_init2(temp, poly2->length + poly3->length, ctx);
113 fq_nmod_mpoly_fit_bits(temp, max_bits, ctx);
114 temp->bits = max_bits;
115
116 len1 = _fq_nmod_mpoly_sub(temp->coeffs, temp->exps,
117 poly2->coeffs, exp2, poly2->length,
118 poly3->coeffs, exp3, poly3->length,
119 N, cmpmask, ctx->fqctx);
120
121 fq_nmod_mpoly_swap(temp, poly1, ctx);
122 fq_nmod_mpoly_clear(temp, ctx);
123
124 }
125 else
126 {
127 fq_nmod_mpoly_fit_length(poly1, poly2->length + poly3->length, ctx);
128 fq_nmod_mpoly_fit_bits(poly1, max_bits, ctx);
129 poly1->bits = max_bits;
130
131 len1 = _fq_nmod_mpoly_sub(poly1->coeffs, poly1->exps,
132 poly2->coeffs, exp2, poly2->length,
133 poly3->coeffs, exp3, poly3->length,
134 N, cmpmask, ctx->fqctx);
135 }
136
137 if (free2)
138 flint_free(exp2);
139
140 if (free3)
141 flint_free(exp3);
142
143 _fq_nmod_mpoly_set_length(poly1, len1, ctx);
144
145 TMP_END;
146 }
147