1 /*
2 Copyright (C) 2008, 2009, 2019 William Hart
3 Copyright (C) 2010 Sebastian Pancratz
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #include <stdlib.h>
14 #include <gmp.h>
15 #include "flint.h"
16 #include "fmpz.h"
17 #include "fmpz_vec.h"
18 #include "fmpz_poly.h"
19
20 #define FLINT_DIVREM_DIVCONQUER_CUTOFF 16
21
22 int
_fmpz_poly_divrem_divconquer_recursive(fmpz * Q,fmpz * BQ,fmpz * W,const fmpz * A,const fmpz * B,slong lenB,int exact)23 _fmpz_poly_divrem_divconquer_recursive(fmpz * Q, fmpz * BQ, fmpz * W,
24 const fmpz * A, const fmpz * B,
25 slong lenB, int exact)
26 {
27 if (lenB <= FLINT_DIVREM_DIVCONQUER_CUTOFF)
28 {
29 _fmpz_vec_zero(BQ, lenB - 1);
30 _fmpz_vec_set(BQ + (lenB - 1), A + (lenB - 1), lenB);
31
32 if (!_fmpz_poly_divrem_basecase(Q, BQ, BQ, 2 * lenB - 1, B, lenB, exact))
33 return 0;
34
35 _fmpz_vec_neg(BQ, BQ, lenB - 1);
36 _fmpz_vec_sub(BQ + (lenB - 1), A + (lenB - 1), BQ + (lenB - 1), lenB);
37 }
38 else
39 {
40 const slong n2 = lenB / 2;
41 const slong n1 = lenB - n2;
42
43 fmpz * W1 = W;
44 fmpz * W2 = W + lenB;
45
46 const fmpz * p1 = A + 2 * n2;
47 const fmpz * p2;
48 const fmpz * d1 = B + n2;
49 const fmpz * d2 = B;
50 const fmpz * d3 = B + n1;
51 const fmpz * d4 = B;
52
53 fmpz * q1 = Q + n2;
54 fmpz * q2 = Q;
55 fmpz * dq1 = BQ + n2;
56 fmpz * d1q1 = BQ + 2 * n2;
57
58 fmpz *d2q1, *d3q2, *d4q2, *t;
59
60 /*
61 Set q1 to p1 div d1, a 2 n1 - 1 by n1 division so q1 ends up
62 being of length n1; d1q1 = d1 q1 is of length 2 n1 - 1
63 */
64
65 if (!_fmpz_poly_divrem_divconquer_recursive(q1, d1q1, W1, p1, d1, n1, exact))
66 return 0;
67
68 /*
69 Compute d2q1 = d2 q1, of length lenB - 1
70 */
71
72 d2q1 = W1;
73 _fmpz_poly_mul(d2q1, q1, n1, d2, n2);
74
75 /*
76 Compute dq1 = d1 q1 x^n2 + d2 q1, of length 2 n1 + n2 - 1
77 */
78
79 _fmpz_vec_swap(dq1, d2q1, n2);
80 _fmpz_vec_add(dq1 + n2, dq1 + n2, d2q1 + n2, n1 - 1);
81
82 /*
83 Compute t = A/x^n2 - dq1, which has length 2 n1 + n2 - 1, but we
84 are not interested in the top n1 coeffs as they will be zero, so
85 this has effective length n1 + n2 - 1
86
87 For the following division, we want to set {p2, 2 n2 - 1} to the
88 top 2 n2 - 1 coeffs of this
89
90 Since the bottom n2 - 1 coeffs of p2 are irrelevant for the
91 division, we in fact set {t, n2} to the relevant coeffs
92 */
93
94 t = BQ;
95 _fmpz_vec_sub(t, A + n2 + (n1 - 1), dq1 + (n1 - 1), n2);
96 p2 = t - (n2 - 1);
97
98 /*
99 Compute q2 = t div d3, a 2 n2 - 1 by n2 division, so q2 will have
100 length n2; let d3q2 = d3 q2, of length 2 n2 - 1
101 */
102
103 d3q2 = W1;
104
105 if (!_fmpz_poly_divrem_divconquer_recursive(q2, d3q2, W2, p2, d3, n2, exact))
106 return 0;
107
108 /*
109 Compute d4q2 = d4 q2, of length n1 + n2 - 1 = lenB - 1
110 */
111
112 d4q2 = W2;
113 _fmpz_poly_mul(d4q2, d4, n1, q2, n2);
114
115 /*
116 Compute dq2 = d3q2 x^n1 + d4q2, of length n1 + 2 n2 - 1
117 */
118
119 _fmpz_vec_swap(BQ, d4q2, n2);
120 _fmpz_vec_add(BQ + n2, BQ + n2, d4q2 + n2, n1 - 1);
121 _fmpz_vec_add(BQ + n1, BQ + n1, d3q2, 2 * n2 - 1);
122
123 /*
124 Note Q = q1 x^n2 + q2, and BQ = dq1 x^n2 + dq2
125 */
126 }
127
128 return 1;
129 }
130