1 /*
2 Copyright 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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include "fmpz_mod_mpoly.h"
15 #include "profiler.h"
16
17
18 slong count = 0;
19 slong total_super = 0;
20
profile_gcd(const fmpz_mod_mpoly_t realG,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx,unsigned int algo)21 void profile_gcd(
22 const fmpz_mod_mpoly_t realG,
23 const fmpz_mod_mpoly_t A,
24 const fmpz_mod_mpoly_t B,
25 const fmpz_mod_mpoly_ctx_t ctx,
26 unsigned int algo)
27 {
28 fmpz_mod_mpoly_t G;
29 timeit_t timer;
30 slong hensel = -1, brown = -1, zippel = -1, zippel2 = -1, super = -1;
31
32 fmpz_mod_mpoly_init(G, ctx);
33
34 if (algo & MPOLY_GCD_USE_BROWN)
35 {
36 fmpz_mod_mpoly_add(G, A, B, ctx);
37 timeit_start(timer);
38 fmpz_mod_mpoly_gcd_brown(G, A, B, ctx);
39 timeit_stop(timer);
40 brown = timer->wall;
41 if (!fmpz_mod_mpoly_equal(G, realG, ctx))
42 {
43 flint_printf("brown is wrong\n");
44 flint_abort();
45 }
46 }
47 flint_printf("%10wd ", brown);
48 fflush(stdout);
49
50 if (algo & MPOLY_GCD_USE_HENSEL)
51 {
52 fmpz_mod_mpoly_add(G, A, B, ctx);
53 timeit_start(timer);
54 fmpz_mod_mpoly_gcd_hensel(G, A, B, ctx);
55 timeit_stop(timer);
56 hensel = timer->wall;
57 if (!fmpz_mod_mpoly_equal(G, realG, ctx))
58 {
59 flint_printf("hensel is wrong\n");
60 flint_abort();
61 }
62 }
63 flint_printf("%10wd ", hensel);
64 fflush(stdout);
65
66 if (algo & MPOLY_GCD_USE_ZIPPEL2)
67 {
68 fmpz_mod_mpoly_add(G, A, B, ctx);
69 timeit_start(timer);
70 fmpz_mod_mpoly_gcd_zippel2(G, A, B, ctx);
71 timeit_stop(timer);
72 zippel2 = timer->wall;
73 if (!fmpz_mod_mpoly_equal(G, realG, ctx))
74 {
75 flint_printf("zippel2 is wrong\n");
76 flint_abort();
77 }
78 }
79 flint_printf("%10wd ", zippel2);
80 fflush(stdout);
81
82 if (0 && algo & MPOLY_GCD_USE_ZIPPEL)
83 {
84 fmpz_mod_mpoly_add(G, A, B, ctx);
85 timeit_start(timer);
86 fmpz_mod_mpoly_gcd_zippel(G, A, B, ctx);
87 timeit_stop(timer);
88 zippel = timer->wall;
89 if (!fmpz_mod_mpoly_equal(G, realG, ctx))
90 {
91 flint_printf("zippel is wrong\n");
92 flint_abort();
93 }
94 }
95 flint_printf("%10wd ", zippel);
96 fflush(stdout);
97
98 {
99 fmpz_mod_mpoly_add(G, A, B, ctx);
100 timeit_start(timer);
101 fmpz_mod_mpoly_gcd(G, A, B, ctx);
102 timeit_stop(timer);
103 super = timer->wall;
104 if (!fmpz_mod_mpoly_equal(G, realG, ctx))
105 {
106 flint_printf("super is wrong\n");
107 flint_abort();
108 }
109 }
110
111 count++;
112
113 flint_printf("%10wd #%wd\n", super, count);
114 fflush(stdout);
115
116 total_super += super;
117
118 fmpz_mod_mpoly_clear(G, ctx);
119 }
120
121
print_banner()122 void print_banner()
123 {
124 flint_printf("| brown | hensel | zippel2 | zippel | super |\n");
125 flint_printf("+----------+----------+----------+----------+----------+\n");
126 }
127
main(int argc,char * argv[])128 int main(int argc, char *argv[])
129 {
130 slong i, j;
131 const char * vars[] = {"x", "y", "z", "t" ,"u", "v", "w", "s", "p"};
132 fmpz_t p;
133
134 fmpz_init(p);
135
136 fmpz_set_ui(p, 3);
137 fmpz_pow_ui(p, p, 70);
138 fmpz_nextprime(p, p, 1);
139
140 print_banner();
141 for (i = 3; i <= 10; i++)
142 {
143 fmpz_mod_mpoly_ctx_t ctx;
144 fmpz_mod_mpoly_t g, a, b, t;
145
146 fmpz_mod_mpoly_ctx_init(ctx, i, ORD_DEGREVLEX, p);
147 fmpz_mod_mpoly_init(a, ctx);
148 fmpz_mod_mpoly_init(b, ctx);
149 fmpz_mod_mpoly_init(g, ctx);
150 fmpz_mod_mpoly_init(t, ctx);
151
152 fmpz_mod_mpoly_one(g, ctx);
153 fmpz_mod_mpoly_one(a, ctx);
154 fmpz_mod_mpoly_one(b, ctx);
155 for (j = 0; j < i; j++)
156 {
157 fmpz_mod_mpoly_gen(t, j, ctx);
158 fmpz_mod_mpoly_add_si(t, t, 1, ctx);
159 fmpz_mod_mpoly_mul(g, g, t, ctx);
160 fmpz_mod_mpoly_gen(t, j, ctx);
161 fmpz_mod_mpoly_sub_si(t, t, 2, ctx);
162 fmpz_mod_mpoly_mul(a, a, t, ctx);
163 fmpz_mod_mpoly_gen(t, j, ctx);
164 fmpz_mod_mpoly_add_si(t, t, 2, ctx);
165 fmpz_mod_mpoly_mul(b, b, t, ctx);
166 }
167 fmpz_mod_mpoly_sub_si(g, g, 2, ctx);
168 fmpz_mod_mpoly_add_si(a, a, 2, ctx);
169 fmpz_mod_mpoly_sub_si(b, b, 2, ctx);
170
171 fmpz_mod_mpoly_mul(a, a, g, ctx);
172 fmpz_mod_mpoly_mul(b, b, g, ctx);
173 fmpz_mod_mpoly_make_monic(g, g, ctx);
174
175 profile_gcd(g, a, b, ctx, MPOLY_GCD_USE_BROWN | MPOLY_GCD_USE_HENSEL);
176
177 fmpz_mod_mpoly_clear(a, ctx);
178 fmpz_mod_mpoly_clear(b, ctx);
179 fmpz_mod_mpoly_clear(g, ctx);
180 fmpz_mod_mpoly_clear(t, ctx);
181 fmpz_mod_mpoly_ctx_clear(ctx);
182 }
183
184 print_banner();
185 for (i = 50; i < 100; i += 4)
186 {
187 fmpz_mod_mpoly_ctx_t ctx;
188 fmpz_mod_mpoly_t a, b, t, m;
189
190 fmpz_mod_mpoly_ctx_init(ctx, 2, ORD_LEX, p);
191 fmpz_mod_mpoly_init(a, ctx);
192 fmpz_mod_mpoly_init(b, ctx);
193 fmpz_mod_mpoly_init(t, ctx);
194 fmpz_mod_mpoly_init(m, ctx);
195
196 fmpz_mod_mpoly_set_str_pretty(a, "1+x+y", vars, ctx);
197 fmpz_mod_mpoly_pow_ui(a, a, i, ctx);
198 fmpz_mod_mpoly_set_str_pretty(m, "x", vars, ctx);
199 fmpz_mod_mpoly_add(a, a, m, ctx);
200
201 fmpz_mod_mpoly_set_str_pretty(b, "1-2*x-y", vars, ctx);
202 fmpz_mod_mpoly_pow_ui(b, b, i, ctx);
203 fmpz_mod_mpoly_set_str_pretty(m, "y", vars, ctx);
204 fmpz_mod_mpoly_add(b, b, m, ctx);
205
206 fmpz_mod_mpoly_set_str_pretty(t, "2-x+y", vars, ctx);
207 fmpz_mod_mpoly_pow_ui(t, t, i, ctx);
208 fmpz_mod_mpoly_set_str_pretty(m, "x-y", vars, ctx);
209 fmpz_mod_mpoly_add(t, t, m, ctx);
210
211 fmpz_mod_mpoly_mul(a, a, t, ctx);
212 fmpz_mod_mpoly_mul(b, b, t, ctx);
213 fmpz_mod_mpoly_make_monic(t, t, ctx);
214
215 profile_gcd(t, a, b, ctx, MPOLY_GCD_USE_BROWN | MPOLY_GCD_USE_HENSEL);
216
217 fmpz_mod_mpoly_clear(a, ctx);
218 fmpz_mod_mpoly_clear(b, ctx);
219 fmpz_mod_mpoly_clear(t, ctx);
220 fmpz_mod_mpoly_clear(m, ctx);
221 fmpz_mod_mpoly_ctx_clear(ctx);
222 }
223
224 print_banner();
225 for (i = 5; i < 15; i += 1)
226 {
227 fmpz_mod_mpoly_ctx_t ctx;
228 fmpz_mod_mpoly_t a, b, t, m;
229
230 fmpz_mod_mpoly_ctx_init(ctx, 2, ORD_LEX, p);
231 fmpz_mod_mpoly_init(a, ctx);
232 fmpz_mod_mpoly_init(b, ctx);
233 fmpz_mod_mpoly_init(t, ctx);
234 fmpz_mod_mpoly_init(m, ctx);
235
236 fmpz_mod_mpoly_set_str_pretty(a, "1+x^31+y^51", vars, ctx);
237 fmpz_mod_mpoly_pow_ui(a, a, i, ctx);
238 fmpz_mod_mpoly_set_str_pretty(m, "x^7", vars, ctx);
239 fmpz_mod_mpoly_add(a, a, m, ctx);
240
241 fmpz_mod_mpoly_set_str_pretty(b, "1-2*x^23-y^47", vars, ctx);
242 fmpz_mod_mpoly_pow_ui(b, b, i, ctx);
243 fmpz_mod_mpoly_set_str_pretty(m, "y^9", vars, ctx);
244 fmpz_mod_mpoly_add(b, b, m, ctx);
245
246 fmpz_mod_mpoly_set_str_pretty(t, "2-x^39+y^24", vars, ctx);
247 fmpz_mod_mpoly_pow_ui(t, t, i, ctx);
248 fmpz_mod_mpoly_set_str_pretty(m, "x^6*y^7", vars, ctx);
249 fmpz_mod_mpoly_add(t, t, m, ctx);
250
251 fmpz_mod_mpoly_mul(a, a, t, ctx);
252 fmpz_mod_mpoly_mul(b, b, t, ctx);
253 fmpz_mod_mpoly_make_monic(t, t, ctx);
254
255 profile_gcd(t, a, b, ctx, MPOLY_GCD_USE_BROWN | MPOLY_GCD_USE_HENSEL);
256
257 fmpz_mod_mpoly_clear(a, ctx);
258 fmpz_mod_mpoly_clear(b, ctx);
259 fmpz_mod_mpoly_clear(t, ctx);
260 fmpz_mod_mpoly_clear(m, ctx);
261 fmpz_mod_mpoly_ctx_clear(ctx);
262 }
263
264 print_banner();
265 for (i = 15; i < 30; i++)
266 {
267 fmpz_mod_mpoly_ctx_t ctx;
268 fmpz_mod_mpoly_t a, b, t, m;
269
270 fmpz_mod_mpoly_ctx_init(ctx, 3, ORD_LEX, p);
271 fmpz_mod_mpoly_init(a, ctx);
272 fmpz_mod_mpoly_init(b, ctx);
273 fmpz_mod_mpoly_init(t, ctx);
274 fmpz_mod_mpoly_init(m, ctx);
275
276 fmpz_mod_mpoly_set_str_pretty(a, "1+x+y+z", vars, ctx);
277 fmpz_mod_mpoly_pow_ui(a, a, i, ctx);
278 fmpz_mod_mpoly_set_str_pretty(m, "x", vars, ctx);
279 fmpz_mod_mpoly_add(a, a, m, ctx);
280
281 fmpz_mod_mpoly_set_str_pretty(b, "1-2*x-y+z", vars, ctx);
282 fmpz_mod_mpoly_pow_ui(b, b, i, ctx);
283 fmpz_mod_mpoly_set_str_pretty(m, "y", vars, ctx);
284 fmpz_mod_mpoly_add(b, b, m, ctx);
285
286 fmpz_mod_mpoly_set_str_pretty(t, "3+x+y-2*z", vars, ctx);
287 fmpz_mod_mpoly_pow_ui(t, t, i, ctx);
288 fmpz_mod_mpoly_set_str_pretty(m, "z", vars, ctx);
289 fmpz_mod_mpoly_add(t, t, m, ctx);
290
291 fmpz_mod_mpoly_mul(a, a, t, ctx);
292 fmpz_mod_mpoly_mul(b, b, t, ctx);
293 fmpz_mod_mpoly_make_monic(t, t, ctx);
294
295 profile_gcd(t, a, b, ctx, MPOLY_GCD_USE_ALL);
296
297 fmpz_mod_mpoly_clear(a, ctx);
298 fmpz_mod_mpoly_clear(b, ctx);
299 fmpz_mod_mpoly_clear(t, ctx);
300 fmpz_mod_mpoly_clear(m, ctx);
301 fmpz_mod_mpoly_ctx_clear(ctx);
302 }
303
304 print_banner();
305 for (i = 1; i < 7; i++)
306 {
307 fmpz_mod_mpoly_ctx_t ctx;
308 fmpz_mod_mpoly_t a, b, t, m;
309
310 fmpz_mod_mpoly_ctx_init(ctx, 7, ORD_LEX, p);
311 fmpz_mod_mpoly_init(a, ctx);
312 fmpz_mod_mpoly_init(b, ctx);
313 fmpz_mod_mpoly_init(t, ctx);
314 fmpz_mod_mpoly_init(m, ctx);
315
316 fmpz_mod_mpoly_set_str_pretty(a, "1+x+y+z+t+u+2*v+w", vars, ctx);
317 fmpz_mod_mpoly_pow_ui(a, a, i, ctx);
318 fmpz_mod_mpoly_set_str_pretty(m, "x", vars, ctx);
319 fmpz_mod_mpoly_add(a, a, m, ctx);
320
321 fmpz_mod_mpoly_set_str_pretty(b, "1-2*x-y+z+t+u+v+w", vars, ctx);
322 fmpz_mod_mpoly_pow_ui(b, b, i, ctx);
323 fmpz_mod_mpoly_set_str_pretty(m, "y", vars, ctx);
324 fmpz_mod_mpoly_add(b, b, m, ctx);
325
326 fmpz_mod_mpoly_set_str_pretty(t, "1+x+y+z-t-u+v+3*w", vars, ctx);
327 fmpz_mod_mpoly_pow_ui(t, t, i, ctx);
328 fmpz_mod_mpoly_set_str_pretty(m, "z", vars, ctx);
329 fmpz_mod_mpoly_add(t, t, m, ctx);
330
331 fmpz_mod_mpoly_mul(a, a, t, ctx);
332 fmpz_mod_mpoly_mul(b, b, t, ctx);
333 fmpz_mod_mpoly_make_monic(t, t, ctx);
334
335 profile_gcd(t, a, b, ctx, MPOLY_GCD_USE_HENSEL | MPOLY_GCD_USE_ZIPPEL2);
336
337 fmpz_mod_mpoly_clear(a, ctx);
338 fmpz_mod_mpoly_clear(b, ctx);
339 fmpz_mod_mpoly_clear(t, ctx);
340 fmpz_mod_mpoly_clear(m, ctx);
341 fmpz_mod_mpoly_ctx_clear(ctx);
342 }
343
344 print_banner();
345 for (i = 1; i < 5; i++)
346 {
347 fmpz_mod_mpoly_ctx_t ctx;
348 fmpz_mod_mpoly_t a, b, t, m;
349
350 fmpz_mod_mpoly_ctx_init(ctx, 9, ORD_LEX, p);
351 fmpz_mod_mpoly_init(a, ctx);
352 fmpz_mod_mpoly_init(b, ctx);
353 fmpz_mod_mpoly_init(t, ctx);
354 fmpz_mod_mpoly_init(m, ctx);
355
356 fmpz_mod_mpoly_set_str_pretty(a, "1+x+y+z+t+u+2*v+w+s-p", vars, ctx);
357 fmpz_mod_mpoly_pow_ui(a, a, i, ctx);
358 fmpz_mod_mpoly_set_str_pretty(m, "x", vars, ctx);
359 fmpz_mod_mpoly_add(a, a, m, ctx);
360
361 fmpz_mod_mpoly_set_str_pretty(b, "1-2*x-y+z+t+u+v+w-2*s+p", vars, ctx);
362 fmpz_mod_mpoly_pow_ui(b, b, i, ctx);
363 fmpz_mod_mpoly_set_str_pretty(m, "y", vars, ctx);
364 fmpz_mod_mpoly_add(b, b, m, ctx);
365
366 fmpz_mod_mpoly_set_str_pretty(t, "1+x+y+z-t-u+v+3*w+2*s-3*p", vars, ctx);
367 fmpz_mod_mpoly_pow_ui(t, t, i, ctx);
368 fmpz_mod_mpoly_set_str_pretty(m, "z", vars, ctx);
369 fmpz_mod_mpoly_add(t, t, m, ctx);
370
371 fmpz_mod_mpoly_mul(a, a, t, ctx);
372 fmpz_mod_mpoly_mul(b, b, t, ctx);
373 fmpz_mod_mpoly_make_monic(t, t, ctx);
374
375 profile_gcd(t, a, b, ctx, MPOLY_GCD_USE_HENSEL | MPOLY_GCD_USE_ZIPPEL2);
376
377 fmpz_mod_mpoly_clear(a, ctx);
378 fmpz_mod_mpoly_clear(b, ctx);
379 fmpz_mod_mpoly_clear(t, ctx);
380 fmpz_mod_mpoly_clear(m, ctx);
381 fmpz_mod_mpoly_ctx_clear(ctx);
382 }
383
384 flint_printf("--------------------\n");
385 flint_printf("total time: %wd\n", total_super);
386 flint_printf(" record: 25312\n");
387
388 fmpz_clear(p);
389
390 flint_cleanup_master();
391 return 0;
392 }
393
394