1 /*
2 Copyright (C) 2011 Sebastian Pancratz
3 Copyright (C) 2015 Claus Fieker
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_resultant_div(fmpz_t rnum,fmpz_t rden,const fmpz * poly1,const fmpz_t den1,slong len1,const fmpz * poly2,const fmpz_t den2,slong len2,const fmpz_t divisor,slong nbits)19 void _fmpq_poly_resultant_div(fmpz_t rnum, fmpz_t rden,
20 const fmpz *poly1, const fmpz_t den1, slong len1,
21 const fmpz *poly2, const fmpz_t den2, slong len2,
22 const fmpz_t divisor, slong nbits)
23 {
24 fmpz_t div, l;
25
26 if (fmpz_is_zero(divisor))
27 {
28 fmpz_zero(rnum);
29 fmpz_one(rden);
30 return;
31 }
32
33 if (len2 == 1)
34 {
35 if (len1 == 1)
36 {
37 fmpz_one(rnum);
38 fmpz_one(rden);
39 }
40 else if (len1 == 2)
41 {
42 fmpz_set(rnum, poly2);
43 fmpz_set(rden, den2);
44 }
45 else
46 {
47 fmpz_pow_ui(rnum, poly2, len1 - 1);
48 if (fmpz_is_one(den2))
49 {
50 fmpz_one(rden);
51 }
52 else
53 {
54 fmpz_pow_ui(rden, den2, len1 - 1);
55 }
56 }
57 fmpz_divexact(rnum, rnum, divisor);
58 }
59 else /* len1 >= len2 >= 2 */
60 {
61 fmpz_t c1, c2, t, la, lb;
62 fmpz *prim1, *prim2;
63
64 fmpz_init(c1);
65 fmpz_init(c2);
66
67 _fmpz_vec_content(c1, poly1, len1);
68 _fmpz_vec_content(c2, poly2, len2);
69
70 prim1 = _fmpz_vec_init(len1);
71 prim2 = _fmpz_vec_init(len2);
72
73 _fmpz_vec_scalar_divexact_fmpz(prim1, poly1, len1, c1);
74 _fmpz_vec_scalar_divexact_fmpz(prim2, poly2, len2, c2);
75
76 fmpz_init(l);
77 if (!fmpz_is_one(c1))
78 {
79 fmpz_pow_ui(l, c1, len2 - 1);
80 fmpz_init(div);
81 fmpz_init(la);
82 fmpz_gcd(div, l, divisor); /* div = gcd(c1^(len2-1), divisor) */
83 fmpz_divexact(la, l, div); /* la = c1^(len2 -1)/gcd */
84 fmpz_divexact(div, divisor, div); /*div /= gcd */
85 nbits = nbits - fmpz_bits(la) + 1;
86 } else {
87 fmpz_init_set(div, divisor);
88 }
89
90 if (!fmpz_is_one(c2))
91 {
92 fmpz_init(lb);
93 fmpz_pow_ui(lb, c2, len1 - 1);
94 fmpz_gcd(l, lb, div);
95 fmpz_divexact(lb, lb, l);
96 fmpz_divexact(div, div, l);
97 nbits = nbits - fmpz_bits(lb) + 1;
98 }
99
100 _fmpz_poly_resultant_modular_div(rnum, prim1, len1, prim2, len2, div, nbits);
101
102 fmpz_init(t);
103 if (!fmpz_is_one(c1))
104 {
105 fmpz_mul(rnum, rnum, la);
106 fmpz_clear(la);
107 }
108 if (!fmpz_is_one(c2))
109 {
110 fmpz_mul(rnum, rnum, lb);
111 fmpz_clear(lb);
112 }
113
114 if (fmpz_is_one(den1))
115 {
116 if (fmpz_is_one(den2))
117 fmpz_one(rden);
118 else
119 fmpz_pow_ui(rden, den2, len1 - 1);
120 }
121 else
122 {
123 if (fmpz_is_one(den2))
124 fmpz_pow_ui(rden, den1, len2 - 1);
125 else
126 {
127 fmpz_pow_ui(rden, den1, len2 - 1);
128 fmpz_pow_ui(t, den2, len1 - 1);
129 fmpz_mul(rden, rden, t);
130 }
131 }
132 _fmpq_canonicalise(rnum, rden);
133 fmpz_clear(t);
134
135 fmpz_clear(c1);
136 fmpz_clear(c2);
137 fmpz_clear(div);
138 fmpz_clear(l);
139 _fmpz_vec_clear(prim1, len1);
140 _fmpz_vec_clear(prim2, len2);
141 }
142 }
143
fmpq_poly_resultant_div(fmpq_t r,const fmpq_poly_t f,const fmpq_poly_t g,const fmpz_t divisor,slong nbits)144 void fmpq_poly_resultant_div(fmpq_t r, const fmpq_poly_t f, const fmpq_poly_t g, const fmpz_t divisor, slong nbits)
145 {
146 const slong len1 = f->length;
147 const slong len2 = g->length;
148
149 if (len1 == 0 || len2 == 0 || fmpz_is_zero(divisor))
150 {
151 fmpq_zero(r);
152 }
153 else
154 {
155 if (len1 >= len2)
156 {
157 _fmpq_poly_resultant_div(fmpq_numref(r), fmpq_denref(r),
158 f->coeffs, f->den, len1,
159 g->coeffs, g->den, len2,
160 divisor, nbits);
161 }
162 else
163 {
164 _fmpq_poly_resultant_div(fmpq_numref(r), fmpq_denref(r),
165 g->coeffs, g->den, len2,
166 f->coeffs, f->den, len1,
167 divisor, nbits);
168
169 if (((len1 | len2) & WORD(1)) == WORD(0))
170 fmpq_neg(r, r);
171 }
172 }
173 }
174
175