1 /*
2 Copyright (C) 2018 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 "nmod_mpoly.h"
13 #include "thread_pool.h"
14
_nmod_mpoly_divides_try_dense(slong * Adegs,slong * Bdegs,slong nvars,slong Alen,slong Blen)15 static int _nmod_mpoly_divides_try_dense(
16 slong * Adegs,
17 slong * Bdegs,
18 slong nvars,
19 slong Alen,
20 slong Blen)
21 {
22 slong i, total_dense_size;
23 ulong hi;
24
25 FLINT_ASSERT(Alen > WORD(0));
26 FLINT_ASSERT(Blen > WORD(0));
27
28 total_dense_size = WORD(1);
29 for (i = 0; i < nvars; i++)
30 {
31 umul_ppmm(hi, total_dense_size, total_dense_size, Adegs[i] + 1);
32 if (hi != WORD(0) || total_dense_size <= WORD(0))
33 return 0;
34 }
35
36 return total_dense_size < WORD(5000000)
37 && total_dense_size/Alen < WORD(10);
38 }
39
40
41 typedef struct
42 {
43 slong * degs;
44 ulong * exps;
45 slong length;
46 flint_bitcnt_t bits;
47 const mpoly_ctx_struct * mctx;
48 const thread_pool_handle * handles;
49 slong num_handles;
50 }
51 _degrees_arg_struct;
52
53 typedef _degrees_arg_struct _degrees_arg_t[1];
54
_worker_degrees(void * varg)55 static void _worker_degrees(void * varg)
56 {
57 _degrees_arg_struct * arg = (_degrees_arg_struct *) varg;
58
59 mpoly_degrees_si_threaded(arg->degs, arg->exps, arg->length, arg->bits,
60 arg->mctx, arg->handles, arg->num_handles);
61 }
62
_nmod_mpoly_divides_threaded_pool(nmod_mpoly_t Q,const nmod_mpoly_t A,const nmod_mpoly_t B,const nmod_mpoly_ctx_t ctx,const thread_pool_handle * handles,slong num_handles)63 int _nmod_mpoly_divides_threaded_pool(
64 nmod_mpoly_t Q,
65 const nmod_mpoly_t A,
66 const nmod_mpoly_t B,
67 const nmod_mpoly_ctx_t ctx,
68 const thread_pool_handle * handles,
69 slong num_handles)
70 {
71 slong i, * Adegs, * Bdegs;
72 int divides;
73 TMP_INIT;
74
75 TMP_START;
76 divides = -1;
77
78 if (A->bits <= FLINT_BITS &&
79 B->bits <= FLINT_BITS &&
80 A->length > 50)
81 {
82 Adegs = (slong *) TMP_ALLOC(ctx->minfo->nvars*sizeof(slong));
83 Bdegs = (slong *) TMP_ALLOC(ctx->minfo->nvars*sizeof(slong));
84
85 if (num_handles > 0)
86 {
87 slong m = mpoly_divide_threads(num_handles, A->length, B->length);
88 _degrees_arg_t arg;
89
90 FLINT_ASSERT(m >= 0);
91 FLINT_ASSERT(m < num_handles);
92
93 arg->degs = Bdegs;
94 arg->exps = B->exps;
95 arg->length = B->length;
96 arg->bits = B->bits;
97 arg->mctx = ctx->minfo;
98 arg->handles = handles + (m + 1);
99 arg->num_handles = num_handles - (m + 1);
100 thread_pool_wake(global_thread_pool, handles[m], 0,
101 _worker_degrees, arg);
102 mpoly_degrees_si_threaded(Adegs, A->exps, A->length, A->bits,
103 ctx->minfo, handles + 0, m);
104 thread_pool_wait(global_thread_pool, handles[m]);
105 }
106 else
107 {
108 mpoly_degrees_si(Adegs, A->exps, A->length, A->bits, ctx->minfo);
109 mpoly_degrees_si(Bdegs, B->exps, B->length, B->bits, ctx->minfo);
110 }
111
112 /* quick degree check */
113 for (i = 0; i < ctx->minfo->nvars; i++)
114 {
115 if (Adegs[i] < Bdegs[i])
116 {
117 nmod_mpoly_zero(Q, ctx);
118 divides = 0;
119 goto cleanup;
120 }
121 }
122
123 if (_nmod_mpoly_divides_try_dense(Adegs, Bdegs, ctx->minfo->nvars,
124 A->length, B->length))
125 {
126 divides = nmod_mpoly_divides_dense(Q, A, B, ctx);
127 }
128 }
129
130 if (divides == 0 || divides == 1)
131 {
132 /* have answer */
133 goto cleanup;
134 }
135
136 if (num_handles > 0)
137 {
138 divides = _nmod_mpoly_divides_heap_threaded_pool(Q, A, B, ctx,
139 handles, num_handles);
140 }
141 else
142 {
143 divides = nmod_mpoly_divides_monagan_pearce(Q, A, B, ctx);
144 }
145
146 cleanup:
147
148 TMP_END;
149 return divides;
150 }
151
nmod_mpoly_divides(nmod_mpoly_t Q,const nmod_mpoly_t A,const nmod_mpoly_t B,const nmod_mpoly_ctx_t ctx)152 int nmod_mpoly_divides(
153 nmod_mpoly_t Q,
154 const nmod_mpoly_t A,
155 const nmod_mpoly_t B,
156 const nmod_mpoly_ctx_t ctx)
157 {
158 thread_pool_handle * handles;
159 slong num_handles;
160 slong thread_limit;
161 int divides;
162
163 if (B->length == 0)
164 {
165 if (A->length == 0 || nmod_mpoly_ctx_modulus(ctx) == 1)
166 {
167 nmod_mpoly_set(Q, A, ctx);
168 return 1;
169 } else
170 flint_throw(FLINT_DIVZERO, "Exception in nmod_mpoly_divides_threaded: "
171 "Cannot divide by zero.\n");
172 }
173
174 if (1 != n_gcd(B->coeffs[0], ctx->ffinfo->mod.n))
175 {
176 flint_throw(FLINT_IMPINV, "Exception in nmod_mpoly_divides_threaded: "
177 "Cannot invert leading coefficient.\n");
178 }
179
180 thread_limit = A->length/1024;
181
182 if (A->length <= 50)
183 {
184 return nmod_mpoly_divides_monagan_pearce(Q, A, B, ctx);
185 }
186
187 num_handles = flint_request_threads(&handles, thread_limit);
188 divides = _nmod_mpoly_divides_threaded_pool(Q, A, B, ctx, handles, num_handles);
189 flint_give_back_threads(handles, num_handles);
190
191 return divides;
192 }
193
194