1 /*
2 Copyright (C) 2010 Sebastian Pancratz
3 Copyright (C) 2014 William Hart
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 <gmp.h>
14 #include "flint.h"
15 #include "fmpz.h"
16 #include "fmpz_vec.h"
17 #include "fmpq_poly.h"
18
_fmpq_poly_add_series_can(fmpz * rpoly,fmpz_t rden,const fmpz * poly1,const fmpz_t den1,slong len1,const fmpz * poly2,const fmpz_t den2,slong len2,slong n,int can)19 void _fmpq_poly_add_series_can(fmpz * rpoly, fmpz_t rden,
20 const fmpz * poly1, const fmpz_t den1, slong len1,
21 const fmpz * poly2, const fmpz_t den2, slong len2, slong n, int can)
22 {
23 int trunc = 0;
24 slong max, min;
25 fmpz_t d;
26
27 if (n < len1)
28 {
29 len1 = n;
30 trunc = 1;
31 }
32
33 if (n < len2)
34 {
35 len2 = n;
36 trunc = 1;
37 }
38
39 max = FLINT_MAX(len1, len2);
40 min = FLINT_MIN(len1, len2);
41
42 if (fmpz_equal(den1, den2))
43 {
44 _fmpz_poly_add(rpoly, poly1, len1, poly2, len2);
45
46 if (fmpz_is_one(den1) || !can)
47 fmpz_set(rden, den1);
48 else
49 {
50 fmpz_init(d);
51 _fmpz_vec_content(d, rpoly, max);
52
53 if (!fmpz_is_one(d))
54 fmpz_gcd(d, d, den1);
55
56 if (fmpz_is_one(d))
57 fmpz_set(rden, den1);
58 else
59 {
60 _fmpz_vec_scalar_divexact_fmpz(rpoly, rpoly, max, d);
61 fmpz_divexact(rden, den1, d);
62 }
63
64 fmpz_clear(d);
65 }
66
67 return;
68 }
69
70 fmpz_init(d);
71 fmpz_one(d);
72 if (!fmpz_is_one(den1) && !fmpz_is_one(den2))
73 fmpz_gcd(d, den1, den2);
74
75 if (fmpz_is_one(d))
76 {
77 _fmpz_vec_scalar_mul_fmpz(rpoly, poly1, len1, den2);
78 _fmpz_vec_scalar_addmul_fmpz(rpoly, poly2, min, den1);
79 if (len1 < len2)
80 _fmpz_vec_scalar_mul_fmpz(rpoly + min, poly2 + min, max - min, den1);
81 fmpz_mul(rden, den1, den2);
82
83 /* may not be canonical if actual truncation happened */
84 if (can && trunc)
85 {
86 if (_fmpz_vec_is_zero(rpoly, max))
87 fmpz_one(rden);
88 else
89 {
90 _fmpz_vec_content(d, rpoly, max);
91 if (!fmpz_is_one(d))
92 fmpz_gcd(d, d, rden);
93
94 if (!fmpz_is_one(d))
95 {
96 _fmpz_vec_scalar_divexact_fmpz(rpoly, rpoly, max, d);
97 fmpz_divexact(rden, rden, d);
98 }
99 }
100 }
101 }
102 else
103 {
104 fmpz_t den11;
105 fmpz_t den22;
106 fmpz_init(den11);
107 fmpz_init(den22);
108 fmpz_divexact(den11, den1, d);
109 fmpz_divexact(den22, den2, d);
110
111 _fmpz_vec_scalar_mul_fmpz(rpoly, poly1, len1, den22);
112 _fmpz_vec_scalar_addmul_fmpz(rpoly, poly2, len2, den11);
113 if (len1 < len2)
114 _fmpz_vec_scalar_mul_fmpz(rpoly + min, poly2 + min, max - min, den11);
115
116 if (_fmpz_vec_is_zero(rpoly, max))
117 fmpz_one(rden);
118 else
119 {
120 if (can)
121 {
122 fmpz_t e;
123 fmpz_init(e);
124 _fmpz_vec_content(e, rpoly, max);
125
126 if (fmpz_is_one(e))
127 fmpz_mul(rden, den1, den22);
128 else
129 {
130 if (trunc) /* there may be extra common factors if truncation occurred */
131 {
132 fmpz_mul(rden, den1, den22);
133 fmpz_gcd(e, e, rden);
134
135 if (!fmpz_is_one(e))
136 {
137 _fmpz_vec_scalar_divexact_fmpz(rpoly, rpoly, max, e);
138 fmpz_divexact(rden, rden, e);
139 }
140 } else
141 {
142 fmpz_gcd(e, e, d);
143
144 if (fmpz_is_one(e))
145 fmpz_mul(rden, den1, den22);
146 else
147 {
148 _fmpz_vec_scalar_divexact_fmpz(rpoly, rpoly, max, e);
149 fmpz_divexact(den11, den1, e);
150 fmpz_mul(rden, den11, den22);
151 }
152 }
153 }
154 fmpz_clear(e);
155 } else
156 fmpz_mul(rden, den1, den22);
157 }
158 fmpz_clear(den11);
159 fmpz_clear(den22);
160 }
161 fmpz_clear(d);
162 }
163
_fmpq_poly_add_series(fmpz * rpoly,fmpz_t rden,const fmpz * poly1,const fmpz_t den1,slong len1,const fmpz * poly2,const fmpz_t den2,slong len2,slong n)164 void _fmpq_poly_add_series(fmpz * rpoly, fmpz_t rden,
165 const fmpz * poly1, const fmpz_t den1, slong len1,
166 const fmpz * poly2, const fmpz_t den2, slong len2, slong n)
167 {
168 _fmpq_poly_add_series_can(rpoly, rden, poly1, den1, len1, poly2, den2, len2, n, 1);
169 }
170
fmpq_poly_add_series_can(fmpq_poly_t res,const fmpq_poly_t poly1,const fmpq_poly_t poly2,slong n,int can)171 void fmpq_poly_add_series_can(fmpq_poly_t res, const fmpq_poly_t poly1,
172 const fmpq_poly_t poly2, slong n, int can)
173 {
174 slong len1 = poly1->length, len2, max;
175
176 if (n < 0)
177 n = 0;
178
179 if (poly1 == poly2) /* Set res = 2 * poly1 */
180 {
181 len1 = FLINT_MIN(len1, n);
182
183 fmpq_poly_fit_length(res, len1);
184 _fmpq_poly_set_length(res, len1);
185
186 if (fmpz_is_even(poly1->den))
187 {
188 _fmpz_vec_set(res->coeffs, poly1->coeffs, len1);
189 fmpz_fdiv_q_2exp(res->den, poly1->den, 1);
190 }
191 else
192 {
193 _fmpz_vec_scalar_mul_2exp(res->coeffs, poly1->coeffs, len1, 1);
194 fmpz_set(res->den, poly1->den);
195 }
196
197 /* may not be canonical if actual truncation happened */
198 if (len1 < poly1->length)
199 {
200 if (can)
201 {
202 fmpz_t e;
203 fmpz_init(e);
204 _fmpz_vec_content(e, res->coeffs, len1);
205 if (!fmpz_is_one(e))
206 fmpz_gcd(e, e, res->den);
207
208 if (!fmpz_is_one(e))
209 {
210 _fmpz_vec_scalar_divexact_fmpz(res->coeffs, res->coeffs, len1, e);
211 fmpz_divexact(res->den, res->den, e);
212 }
213 fmpz_clear(e);
214 }
215
216 _fmpq_poly_normalise(res);
217 }
218
219 return;
220 }
221
222 len2 = poly2->length;
223 max = FLINT_MAX(len1, len2);
224 max = FLINT_MIN(max, n);
225
226 fmpq_poly_fit_length(res, max);
227
228 if (res != poly2)
229 _fmpq_poly_add_series_can(res->coeffs, res->den,
230 poly1->coeffs, poly1->den, len1,
231 poly2->coeffs, poly2->den, len2, n, can);
232 else
233 _fmpq_poly_add_series_can(res->coeffs, res->den,
234 poly2->coeffs, poly2->den, len2,
235 poly1->coeffs, poly1->den, len1, n, can);
236
237 _fmpq_poly_set_length(res, max);
238 _fmpq_poly_normalise(res);
239 }
240
fmpq_poly_add_series(fmpq_poly_t res,const fmpq_poly_t poly1,const fmpq_poly_t poly2,slong n)241 void fmpq_poly_add_series(fmpq_poly_t res, const fmpq_poly_t poly1,
242 const fmpq_poly_t poly2, slong n)
243 {
244 fmpq_poly_add_series_can(res, poly1, poly2, n, 1);
245 }
246
247