1 /*
2 Copyright (C) 2010 Sebastian Pancratz
3 Copyright (C) 2011 Fredrik Johansson
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_poly.h"
17 #include "fmpq_poly.h"
18
19
20 static void
_set_vec(fmpz * rnum,fmpz_t den,const fmpz * xnum,const fmpz * xden,slong len)21 _set_vec(fmpz * rnum, fmpz_t den,
22 const fmpz * xnum, const fmpz * xden, slong len)
23 {
24 slong j;
25 fmpz_t t;
26 fmpz_init(t);
27 fmpz_one(den);
28
29 for (j = 0; j < len; j++)
30 fmpz_lcm(den, den, xden + j);
31
32 for (j = 0; j < len; j++)
33 {
34 fmpz_divexact(t, den, xden + j);
35 fmpz_mul(rnum + j, xnum + j, t);
36 }
37
38 fmpz_clear(t);
39 }
40
41 void
_fmpq_poly_revert_series_lagrange(fmpz * Qinv,fmpz_t den,const fmpz * Q,const fmpz_t Qden,slong Qlen,slong n)42 _fmpq_poly_revert_series_lagrange(fmpz * Qinv, fmpz_t den,
43 const fmpz * Q, const fmpz_t Qden, slong Qlen, slong n)
44 {
45 slong i;
46 fmpz *R, *S, *T, *dens, *tmp;
47 fmpz_t Rden, Sden, Tden;
48
49 Qlen = FLINT_MIN(Qlen, n);
50
51 if (Qlen <= 2)
52 {
53 fmpz_zero(Qinv);
54
55 if (Qlen == 2)
56 {
57 fmpz_set(Qinv + 1, Qden);
58 fmpz_set(den, Q + 1);
59 _fmpq_poly_canonicalise(Qinv, den, 2);
60 }
61
62 _fmpz_vec_zero(Qinv + 2, n - 2);
63 }
64 else
65 {
66 dens = _fmpz_vec_init(n);
67 R = _fmpz_vec_init(n - 1);
68 S = _fmpz_vec_init(n - 1);
69 T = _fmpz_vec_init(n - 1);
70 fmpz_init(Rden);
71 fmpz_init(Sden);
72 fmpz_init(Tden);
73
74 fmpz_zero(Qinv);
75 fmpz_one(dens);
76 fmpz_set(Qinv + 1, Qden);
77 fmpz_set(dens + 1, Q + 1);
78
79 _fmpq_poly_inv_series(R, Rden, Q + 1, Qden, Qlen - 1, n - 1);
80 _fmpq_poly_canonicalise(R, Rden, n - 1);
81
82 _fmpz_vec_set(S, R, n - 1);
83 fmpz_set(Sden, Rden);
84
85 for (i = 2; i < n; i++)
86 {
87 _fmpq_poly_mullow(T, Tden, S, Sden, n - 1, R, Rden, n - 1, n - 1);
88 _fmpq_poly_canonicalise(T, Tden, n - 1);
89 fmpz_set(Qinv + i, T + i - 1);
90 fmpz_mul_ui(dens + i, Tden, i);
91 tmp = S; S = T; T = tmp;
92 fmpz_swap(Sden, Tden);
93 }
94
95 _set_vec(Qinv, den, Qinv, dens, n);
96 _fmpq_poly_canonicalise(Qinv, den, n);
97
98 _fmpz_vec_clear(R, n - 1);
99 _fmpz_vec_clear(S, n - 1);
100 _fmpz_vec_clear(T, n - 1);
101 _fmpz_vec_clear(dens, n);
102 fmpz_clear(Rden);
103 fmpz_clear(Sden);
104 fmpz_clear(Tden);
105 }
106 }
107
108 void
fmpq_poly_revert_series_lagrange(fmpq_poly_t res,const fmpq_poly_t poly,slong n)109 fmpq_poly_revert_series_lagrange(fmpq_poly_t res,
110 const fmpq_poly_t poly, slong n)
111 {
112 if (poly->length < 2 || !fmpz_is_zero(poly->coeffs)
113 || fmpz_is_zero(poly->coeffs + 1))
114 {
115 flint_printf("Exception (fmpq_poly_revert_series_lagrange). Input must have \n"
116 "zero constant term and nonzero coefficient of x^1.\n");
117 flint_abort();
118 }
119
120 if (n < 2)
121 {
122 fmpq_poly_zero(res);
123 return;
124 }
125
126 if (res != poly)
127 {
128 fmpq_poly_fit_length(res, n);
129 _fmpq_poly_revert_series_lagrange(res->coeffs,
130 res->den, poly->coeffs, poly->den, poly->length, n);
131 }
132 else
133 {
134 fmpq_poly_t t;
135 fmpq_poly_init2(t, n);
136 _fmpq_poly_revert_series_lagrange(t->coeffs,
137 t->den, poly->coeffs, poly->den, poly->length, n);
138 fmpq_poly_swap(res, t);
139 fmpq_poly_clear(t);
140 }
141
142 _fmpq_poly_set_length(res, n);
143 _fmpq_poly_normalise(res);
144 }
145