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