1 /*
2     Copyright (C) 2021 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 
16 int
main(void)17 main(void)
18 {
19     slong i, j;
20     slong tmul = 5;
21     FLINT_TEST_INIT(state);
22 
23     flint_printf("derivative....");
24     fflush(stdout);
25 
26     /* Check d(f*g) = df*g + f*dg */
27     for (i = 0; i < tmul * flint_test_multiplier(); i++)
28     {
29         fmpz_mod_mpoly_ctx_t ctx;
30         fmpz_mod_mpoly_t f, g, h, fp, gp, hp, t1, t2;
31         slong len, len1, len2;
32         flint_bitcnt_t exp_bits, exp_bits1, exp_bits2;
33         slong idx;
34 
35         fmpz_mod_mpoly_ctx_init_rand_bits(ctx, state, 20, 200);
36 
37         fmpz_mod_mpoly_init(f, ctx);
38         fmpz_mod_mpoly_init(g, ctx);
39         fmpz_mod_mpoly_init(h, ctx);
40         fmpz_mod_mpoly_init(fp, ctx);
41         fmpz_mod_mpoly_init(gp, ctx);
42         fmpz_mod_mpoly_init(hp, ctx);
43         fmpz_mod_mpoly_init(t1, ctx);
44         fmpz_mod_mpoly_init(t2, ctx);
45 
46         len = n_randint(state, 100);
47         len1 = n_randint(state, 100);
48         len2 = n_randint(state, 100);
49 
50         exp_bits = n_randint(state, 200) + 1;
51         exp_bits1 = n_randint(state, 200) + 1;
52         exp_bits2 = n_randint(state, 200) + 1;
53 
54         fmpz_mod_mpoly_randtest_bits(hp, state, len, exp_bits, ctx);
55         fmpz_mod_mpoly_randtest_bits(fp, state, len, exp_bits1, ctx);
56         fmpz_mod_mpoly_randtest_bits(gp, state, len, exp_bits2, ctx);
57 
58         for (j = 0; j < 4; j++)
59         {
60             fmpz_mod_mpoly_randtest_bound(f, state, len1, exp_bits1, ctx);
61             fmpz_mod_mpoly_randtest_bound(g, state, len2, exp_bits2, ctx);
62 
63             idx = n_randint(state, ctx->minfo->nvars);
64 
65             fmpz_mod_mpoly_mul_johnson(h, f, g, ctx);
66             fmpz_mod_mpoly_assert_canonical(h, ctx);
67 
68             fmpz_mod_mpoly_derivative(hp, h, idx, ctx);
69             fmpz_mod_mpoly_assert_canonical(hp, ctx);
70 
71             fmpz_mod_mpoly_derivative(fp, f, idx, ctx);
72             fmpz_mod_mpoly_assert_canonical(fp, ctx);
73             fmpz_mod_mpoly_derivative(gp, g, idx, ctx);
74             fmpz_mod_mpoly_assert_canonical(gp, ctx);
75 
76             fmpz_mod_mpoly_mul_johnson(t1, f, gp, ctx);
77             fmpz_mod_mpoly_assert_canonical(t1, ctx);
78             fmpz_mod_mpoly_mul_johnson(t2, g, fp, ctx);
79             fmpz_mod_mpoly_assert_canonical(t2, ctx);
80             fmpz_mod_mpoly_add(t1, t1, t2, ctx);
81             fmpz_mod_mpoly_assert_canonical(t1, ctx);
82 
83             if (!fmpz_mod_mpoly_equal(hp, t1, ctx))
84             {
85                 flint_printf("FAIL: Check d(f*g) = df*g + f*dg\n");
86                 flint_printf("i = %wd, j = %wd\n", i, j);
87                 flint_abort();
88                 flint_abort();
89             }
90         }
91 
92         fmpz_mod_mpoly_clear(f, ctx);
93         fmpz_mod_mpoly_clear(g, ctx);
94         fmpz_mod_mpoly_clear(h, ctx);
95         fmpz_mod_mpoly_clear(fp, ctx);
96         fmpz_mod_mpoly_clear(gp, ctx);
97         fmpz_mod_mpoly_clear(hp, ctx);
98         fmpz_mod_mpoly_clear(t1, ctx);
99         fmpz_mod_mpoly_clear(t2, ctx);
100 
101         fmpz_mod_mpoly_ctx_clear(ctx);
102     }
103 
104     /* Check d(f*g) = df*g + f*dg with aliasing */
105     for (i = 0; i < tmul * flint_test_multiplier(); i++)
106     {
107         fmpz_mod_mpoly_ctx_t ctx;
108         fmpz_mod_mpoly_t f, g, h, fp, gp, t1, t2;
109         slong len, len1, len2;
110         flint_bitcnt_t exp_bits, exp_bits1, exp_bits2;
111         slong idx;
112 
113         fmpz_mod_mpoly_ctx_init_rand_bits(ctx, state, 20, 200);
114 
115         fmpz_mod_mpoly_init(f, ctx);
116         fmpz_mod_mpoly_init(g, ctx);
117         fmpz_mod_mpoly_init(h, ctx);
118         fmpz_mod_mpoly_init(fp, ctx);
119         fmpz_mod_mpoly_init(gp, ctx);
120         fmpz_mod_mpoly_init(t1, ctx);
121         fmpz_mod_mpoly_init(t2, ctx);
122 
123         len = n_randint(state, 100);
124         len1 = n_randint(state, 100);
125         len2 = n_randint(state, 100);
126 
127         exp_bits = n_randint(state, 200) + 1;
128         exp_bits1 = n_randint(state, 200) + 1;
129         exp_bits2 = n_randint(state, 200) + 1;
130 
131         fmpz_mod_mpoly_randtest_bits(h, state, len, exp_bits, ctx);
132         fmpz_mod_mpoly_randtest_bits(fp, state, len, exp_bits, ctx);
133         fmpz_mod_mpoly_randtest_bits(gp, state, len, exp_bits, ctx);
134 
135         for (j = 0; j < 4; j++)
136         {
137             fmpz_mod_mpoly_randtest_bits(f, state, len1, exp_bits1, ctx);
138             fmpz_mod_mpoly_randtest_bits(g, state, len2, exp_bits2, ctx);
139 
140             idx = n_randint(state, ctx->minfo->nvars);
141 
142             fmpz_mod_mpoly_mul_johnson(h, f, g, ctx);
143             fmpz_mod_mpoly_assert_canonical(h, ctx);
144 
145             fmpz_mod_mpoly_derivative(h, h, idx, ctx);
146             fmpz_mod_mpoly_assert_canonical(h, ctx);
147             fmpz_mod_mpoly_set(fp, f, ctx);
148             fmpz_mod_mpoly_derivative(fp, fp, idx, ctx);
149             fmpz_mod_mpoly_assert_canonical(fp, ctx);
150             fmpz_mod_mpoly_set(gp, g, ctx);
151             fmpz_mod_mpoly_derivative(gp, gp, idx, ctx);
152             fmpz_mod_mpoly_assert_canonical(gp, ctx);
153 
154             fmpz_mod_mpoly_mul_johnson(t1, f, gp, ctx);
155             fmpz_mod_mpoly_assert_canonical(t1, ctx);
156             fmpz_mod_mpoly_mul_johnson(t2, g, fp, ctx);
157             fmpz_mod_mpoly_assert_canonical(t2, ctx);
158             fmpz_mod_mpoly_add(t1, t1, t2, ctx);
159             fmpz_mod_mpoly_assert_canonical(t1, ctx);
160 
161             if (!fmpz_mod_mpoly_equal(h, t1, ctx))
162             {
163                 flint_printf("FAIL: Check d(f*g) = df*g + f*dg with aliasing\n");
164                 flint_printf("i = %wd, j = %wd\n", i, j);
165                 flint_abort();
166             }
167         }
168 
169         fmpz_mod_mpoly_clear(f, ctx);
170         fmpz_mod_mpoly_clear(g, ctx);
171         fmpz_mod_mpoly_clear(h, ctx);
172         fmpz_mod_mpoly_clear(fp, ctx);
173         fmpz_mod_mpoly_clear(gp, ctx);
174         fmpz_mod_mpoly_clear(t1, ctx);
175         fmpz_mod_mpoly_clear(t2, ctx);
176 
177         fmpz_mod_mpoly_ctx_clear(ctx);
178     }
179 
180     FLINT_TEST_CLEANUP(state);
181     flint_printf("PASS\n");
182     return 0;
183 }
184 
185