1 /*
2 Copyright (C) 2011 William Hart
3 Copyright (C) 2011 Sebastian Pancratz
4 Copyright (C) 2013 Mike Hansen
5
6 This file is part of FLINT.
7
8 FLINT is free software: you can redistribute it and/or modify it under
9 the terms of the GNU Lesser General Public License (LGPL) as published
10 by the Free Software Foundation; either version 2.1 of the License, or
11 (at your option) any later version. See <https://www.gnu.org/licenses/>.
12 */
13
14 #ifdef T
15
16 #include "templates.h"
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <gmp.h>
21 #include "flint.h"
22 #include "ulong_extras.h"
23 #include "mpn_extras.h"
24
25 #define __mul(C, lenC, A, lenA, B, lenB) \
26 do { \
27 if ((lenA) != 0 && (lenB) != 0) \
28 { \
29 if ((lenA) >= (lenB)) \
30 _TEMPLATE(T, poly_mul)((C), (A), (lenA), (B), (lenB), ctx); \
31 else \
32 _TEMPLATE(T, poly_mul)((C), (B), (lenB), (A), (lenA), ctx); \
33 (lenC) = (lenA) + (lenB) - 1; \
34 } \
35 else \
36 { \
37 (lenC) = 0; \
38 } \
39 } while (0)
40
41
42 #define __swap(u, l, v, m) \
43 do { \
44 { TEMPLATE(T, struct)* _; _ = (u), (u) = (v), (v) = _;} \
45 { slong _; _ = (l), (l) = (m), (m) = _;} \
46 } while (0)
47
48 int
main(void)49 main(void)
50 {
51 int i, result;
52 flint_rand_t state;
53 flint_randinit(state);
54
55 flint_printf("hgcd....");
56 fflush(stdout);
57
58 /*
59 Find coprime polys, multiply by another poly
60 and check the GCD is that poly
61 */
62 for (i = 0; i < 10 * flint_test_multiplier(); i++)
63 {
64 TEMPLATE(T, poly_t) a, b, c, d, c1, d1, s, t;
65 TEMPLATE(T, ctx_t) ctx;
66 TEMPLATE(T, struct) * M[4];
67 slong lenM[4];
68 slong sgnM;
69
70 TEMPLATE(T, ctx_randtest) (ctx, state);
71
72 TEMPLATE(T, poly_init) (a, ctx);
73 TEMPLATE(T, poly_init) (b, ctx);
74 TEMPLATE(T, poly_init) (c, ctx);
75 TEMPLATE(T, poly_init) (d, ctx);
76 TEMPLATE(T, poly_init) (c1, ctx);
77 TEMPLATE(T, poly_init) (d1, ctx);
78 TEMPLATE(T, poly_init) (s, ctx);
79 TEMPLATE(T, poly_init) (t, ctx);
80
81 do
82 {
83 TEMPLATE(T, poly_randtest_not_zero) (a, state,
84 n_randint(state, 800) + 1,
85 ctx);
86 TEMPLATE(T, poly_randtest_not_zero) (b, state,
87 n_randint(state, 800) + 1,
88 ctx);
89 } while (a->length == b->length);
90
91 if (a->length < b->length)
92 TEMPLATE(T, poly_swap) (a, b, ctx);
93
94 M[0] = _TEMPLATE(T, vec_init) (a->length, ctx);
95 M[1] = _TEMPLATE(T, vec_init) (a->length, ctx);
96 M[2] = _TEMPLATE(T, vec_init) (a->length, ctx);
97 M[3] = _TEMPLATE(T, vec_init) (a->length, ctx);
98
99 TEMPLATE(T, poly_fit_length) (c, a->length, ctx);
100 TEMPLATE(T, poly_fit_length) (d, b->length, ctx);
101
102 sgnM = _TEMPLATE(T, poly_hgcd) (M, lenM,
103 c->coeffs, &(c->length), d->coeffs,
104 &(d->length), a->coeffs, a->length,
105 b->coeffs, b->length, ctx);
106
107 TEMPLATE(T, poly_fit_length) (s, 2 * a->length, ctx);
108 TEMPLATE(T, poly_fit_length) (t, 2 * a->length, ctx);
109
110 /* [c1,d1] := sgnM * M^{-1} [a,b] */
111 {
112 __swap(M[0], lenM[0], M[3], lenM[3]);
113 _TEMPLATE(T, vec_neg) (M[1], M[1], lenM[1], ctx);
114 _TEMPLATE(T, vec_neg) (M[2], M[2], lenM[2], ctx);
115
116 __mul(s->coeffs, s->length, M[0], lenM[0], a->coeffs, a->length);
117 __mul(t->coeffs, t->length, M[1], lenM[1], b->coeffs, b->length);
118 TEMPLATE(T, poly_add) (c1, s, t, ctx);
119 __mul(s->coeffs, s->length, M[2], lenM[2], a->coeffs, a->length);
120 __mul(t->coeffs, t->length, M[3], lenM[3], b->coeffs, b->length);
121 TEMPLATE(T, poly_add) (d1, s, t, ctx);
122 }
123
124 if (sgnM < 0)
125 {
126 TEMPLATE(T, poly_neg) (c1, c1, ctx);
127 TEMPLATE(T, poly_neg) (d1, d1, ctx);
128 }
129
130 result = (TEMPLATE(T, poly_equal) (c, c1, ctx)
131 && TEMPLATE(T, poly_equal) (d, d1, ctx));
132 if (!result)
133 {
134 flint_printf("FAIL:\n");
135 flint_printf("a = "), TEMPLATE(T, poly_print_pretty) (a, "x",
136 ctx),
137 flint_printf("\n\n");
138 flint_printf("b = "), TEMPLATE(T, poly_print_pretty) (b, "x",
139 ctx),
140 flint_printf("\n\n");
141 flint_printf("c = "), TEMPLATE(T, poly_print_pretty) (c, "x",
142 ctx),
143 flint_printf("\n\n");
144 flint_printf("d = "), TEMPLATE(T, poly_print_pretty) (d, "x",
145 ctx),
146 flint_printf("\n\n");
147 flint_printf("c1 = "), TEMPLATE(T, poly_print_pretty) (c1, "x",
148 ctx),
149 flint_printf("\n\n");
150 flint_printf("d1 = "), TEMPLATE(T, poly_print_pretty) (d1, "x",
151 ctx),
152 flint_printf("\n\n");
153 abort();
154 }
155
156 TEMPLATE(T, poly_clear) (a, ctx);
157 TEMPLATE(T, poly_clear) (b, ctx);
158 TEMPLATE(T, poly_clear) (c, ctx);
159 TEMPLATE(T, poly_clear) (d, ctx);
160 TEMPLATE(T, poly_clear) (c1, ctx);
161 TEMPLATE(T, poly_clear) (d1, ctx);
162 TEMPLATE(T, poly_clear) (s, ctx);
163 TEMPLATE(T, poly_clear) (t, ctx);
164
165 _TEMPLATE(T, vec_clear) (M[0], a->length, ctx);
166 _TEMPLATE(T, vec_clear) (M[1], a->length, ctx);
167 _TEMPLATE(T, vec_clear) (M[2], a->length, ctx);
168 _TEMPLATE(T, vec_clear) (M[3], a->length, ctx);
169
170 TEMPLATE(T, ctx_clear) (ctx);
171 }
172
173 flint_randclear(state);
174
175 flint_printf("PASS\n");
176 return 0;
177 }
178
179 #undef __mul
180 #undef __swap
181
182 #endif
183