1 /*
2 Copyright (C) 2008, 2009 William Hart
3 Copyright (C) 2010, 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 <http://www.gnu.org/licenses/>.
12 */
13
14 #ifdef T
15
16 #include "templates.h"
17
18 static void
__TEMPLATE(T,poly_divrem_divconquer)19 __TEMPLATE(T, poly_divrem_divconquer) (TEMPLATE(T, struct) * Q,
20 TEMPLATE(T, struct) * R,
21 const TEMPLATE(T, struct) * A, slong lenA,
22 const TEMPLATE(T, struct) * B, slong lenB,
23 const TEMPLATE(T, t) invB,
24 const TEMPLATE(T, ctx_t) ctx)
25 {
26 if (lenA < 2 * lenB - 1)
27 {
28 /*
29 Convert unbalanced division into a 2 n1 - 1 by n1 division
30 */
31
32 const slong n1 = lenA - lenB + 1;
33 const slong n2 = lenB - n1;
34
35 const TEMPLATE(T, struct) * p1 = A + n2;
36 const TEMPLATE(T, struct) * d1 = B + n2;
37 const TEMPLATE(T, struct) * d2 = B;
38
39 TEMPLATE(T, struct) * W =
40 _TEMPLATE(T, vec_init) ((2 * n1 - 1) + lenB - 1, ctx);
41
42 TEMPLATE(T, struct) * d1q1 = R + n2;
43 TEMPLATE(T, struct) * d2q1 = W + (2 * n1 - 1);
44
45 _TEMPLATE(T, poly_divrem_divconquer_recursive) (Q, d1q1, W, p1, d1, n1,
46 invB, ctx);
47
48 /*
49 Compute d2q1 = Q d2, of length lenB - 1
50 */
51
52 if (n1 >= n2)
53 _TEMPLATE(T, poly_mul) (d2q1, Q, n1, d2, n2, ctx);
54 else
55 _TEMPLATE(T, poly_mul) (d2q1, d2, n2, Q, n1, ctx);
56
57 /*
58 Compute BQ = d1q1 * x^n1 + d2q1, of length lenB - 1;
59 then compute R = A - BQ
60 */
61
62 _TEMPLATE(T, vec_swap) (R, d2q1, n2, ctx);
63 _TEMPLATE(T, poly_add) (R + n2, R + n2, n1 - 1, d2q1 + n2, n1 - 1,
64 ctx);
65 _TEMPLATE(T, poly_sub) (R, A, lenA, R, lenA, ctx);
66
67 _TEMPLATE(T, vec_clear) (W, (2 * n1 - 1) + lenB - 1, ctx);
68 }
69 else /* lenA = 2 * lenB - 1 */
70 {
71 TEMPLATE(T, struct) * W = _TEMPLATE(T, vec_init) (lenA, ctx);
72
73 _TEMPLATE(T, poly_divrem_divconquer_recursive) (Q, R, W, A, B, lenB,
74 invB, ctx);
75
76 _TEMPLATE(T, poly_sub) (R, A, lenB - 1, R, lenB - 1, ctx);
77
78 _TEMPLATE(T, vec_clear) (W, lenA, ctx);
79 }
80 }
81
82 void
_TEMPLATE(T,poly_divrem_divconquer)83 _TEMPLATE(T, poly_divrem_divconquer) (TEMPLATE(T, struct) * Q,
84 TEMPLATE(T, struct) * R,
85 const TEMPLATE(T, struct) * A,
86 slong lenA, const TEMPLATE(T,
87 struct) * B,
88 slong lenB, const TEMPLATE(T, t) invB,
89 const TEMPLATE(T, ctx_t) ctx)
90 {
91 if (lenA <= 2 * lenB - 1)
92 {
93 __TEMPLATE(T, poly_divrem_divconquer) (Q, R, A, lenA, B, lenB, invB,
94 ctx);
95 }
96 else /* lenA > 2 * lenB - 1 */
97 {
98 slong shift, n = 2 * lenB - 1;
99 TEMPLATE(T, struct) * QB, *W;
100
101 _TEMPLATE(T, vec_set) (R, A, lenA, ctx);
102 W = _TEMPLATE(T, vec_init) (2 * n, ctx);
103 QB = W + n;
104
105 while (lenA >= n)
106 {
107 shift = lenA - n;
108 _TEMPLATE(T, poly_divrem_divconquer_recursive) (Q + shift, QB,
109 W, R + shift, B,
110 lenB, invB, ctx);
111 _TEMPLATE(T, poly_sub) (R + shift, R + shift, n, QB, n, ctx);
112 lenA -= lenB;
113 }
114
115 if (lenA >= lenB)
116 {
117 __TEMPLATE(T, poly_divrem_divconquer) (Q, W, R, lenA, B, lenB,
118 invB, ctx);
119 _TEMPLATE(T, vec_swap) (W, R, lenA, ctx);
120 }
121
122 _TEMPLATE(T, vec_clear) (W, 2 * n, ctx);
123 }
124 }
125
126 void
TEMPLATE(T,poly_divrem_divconquer)127 TEMPLATE(T, poly_divrem_divconquer) (TEMPLATE(T, poly_t) Q,
128 TEMPLATE(T, poly_t) R,
129 const TEMPLATE(T, poly_t) A,
130 const TEMPLATE(T, poly_t) B,
131 const TEMPLATE(T, ctx_t) ctx)
132 {
133 const slong lenA = A->length;
134 const slong lenB = B->length;
135 const slong lenQ = lenA - lenB + 1;
136
137 TEMPLATE(T, struct) * q, *r;
138 TEMPLATE(T, t) invB;
139
140 if (lenA < lenB)
141 {
142 TEMPLATE(T, poly_set) (R, A, ctx);
143 TEMPLATE(T, poly_zero) (Q, ctx);
144 return;
145 }
146
147 TEMPLATE(T, init) (invB, ctx);
148 TEMPLATE(T, inv) (invB, TEMPLATE(T, poly_lead) (B, ctx), ctx);
149
150 if (Q == A || Q == B)
151 {
152 q = _TEMPLATE(T, vec_init) (lenQ, ctx);
153 }
154 else
155 {
156 TEMPLATE(T, poly_fit_length) (Q, lenQ, ctx);
157 q = Q->coeffs;
158 }
159
160 if (R == A || R == B)
161 {
162 r = _TEMPLATE(T, vec_init) (lenA, ctx);
163 }
164 else
165 {
166 TEMPLATE(T, poly_fit_length) (R, lenA, ctx);
167 r = R->coeffs;
168 }
169
170 _TEMPLATE(T, poly_divrem_divconquer) (q, r, A->coeffs, lenA,
171 B->coeffs, lenB, invB, ctx);
172
173 if (Q == A || Q == B)
174 {
175 _TEMPLATE(T, vec_clear) (Q->coeffs, Q->alloc, ctx);
176 Q->coeffs = q;
177 Q->alloc = lenQ;
178 Q->length = lenQ;
179 }
180 else
181 {
182 _TEMPLATE(T, poly_set_length) (Q, lenQ, ctx);
183 }
184
185 if (R == A || R == B)
186 {
187 _TEMPLATE(T, vec_clear) (R->coeffs, R->alloc, ctx);
188 R->coeffs = r;
189 R->alloc = lenA;
190 R->length = lenA;
191 }
192 _TEMPLATE(T, poly_set_length) (R, lenB - 1, ctx);
193 _TEMPLATE(T, poly_normalise) (R, ctx);
194
195 TEMPLATE(T, clear) (invB, ctx);
196 }
197
198
199 #endif
200