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 <https://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 #define FLINT_REVERSE_NEWTON_CUTOFF 4
21
22 void
_fmpq_poly_revert_series_newton(fmpz * Qinv,fmpz_t den,const fmpz * Q,const fmpz_t Qden,slong Qlen,slong n)23 _fmpq_poly_revert_series_newton(fmpz * Qinv, fmpz_t den,
24 const fmpz * Q, const fmpz_t Qden, slong Qlen, slong n)
25 {
26 Qlen = FLINT_MIN(Qlen, n);
27
28 if (Qlen <= 2)
29 {
30 fmpz_zero(Qinv);
31
32 if (Qlen == 2)
33 {
34 fmpz_set(Qinv + 1, Qden);
35 fmpz_set(den, Q + 1);
36 _fmpq_poly_canonicalise(Qinv, den, 2);
37 }
38
39 _fmpz_vec_zero(Qinv + 2, n - 2);
40 }
41 else
42 {
43 fmpz *T, *U, *V;
44 fmpz_t Tden, Uden, Vden;
45
46 T = _fmpz_vec_init(n);
47 U = _fmpz_vec_init(n);
48 V = _fmpz_vec_init(n);
49 fmpz_init(Tden);
50 fmpz_init(Uden);
51 fmpz_init(Vden);
52
53 FLINT_NEWTON_INIT(FLINT_REVERSE_NEWTON_CUTOFF, n)
54
55 FLINT_NEWTON_BASECASE(k)
56 _fmpq_poly_revert_series_lagrange(Qinv, den, Q, Qden, Qlen, k);
57 _fmpz_vec_zero(Qinv + k, n - k);
58 FLINT_NEWTON_END_BASECASE
59
60 FLINT_NEWTON_LOOP(k0, k)
61 _fmpq_poly_compose_series(T, Tden, Q, Qden, FLINT_MIN(Qlen, k), Qinv, den, k0, k);
62 _fmpq_poly_derivative(U, Uden, T, Tden, k); fmpz_zero(U + k - 1);
63 fmpz_zero(T + 1);
64 _fmpq_poly_div_series(V, Vden, T, Tden, k, U, Uden, k, k);
65 _fmpq_poly_canonicalise(V, Vden, k);
66 _fmpq_poly_derivative(T, Tden, Qinv, den, k);
67 _fmpq_poly_mullow(U, Uden, V, Vden, k, T, Tden, k, k);
68 _fmpq_poly_sub(Qinv, den, Qinv, den, k, U, Uden, k);
69
70 FLINT_NEWTON_END_LOOP
71 FLINT_NEWTON_END
72
73 _fmpq_poly_canonicalise(Qinv, den, n);
74
75 _fmpz_vec_clear(T, n);
76 _fmpz_vec_clear(U, n);
77 _fmpz_vec_clear(V, n);
78 fmpz_clear(Tden);
79 fmpz_clear(Uden);
80 fmpz_clear(Vden);
81 }
82 }
83
84 void
fmpq_poly_revert_series_newton(fmpq_poly_t res,const fmpq_poly_t poly,slong n)85 fmpq_poly_revert_series_newton(fmpq_poly_t res,
86 const fmpq_poly_t poly, slong n)
87 {
88 if (poly->length < 2 || !fmpz_is_zero(poly->coeffs)
89 || fmpz_is_zero(poly->coeffs + 1))
90 {
91 flint_printf("Exception (fmpq_poly_revert_series_newton). Input must have \n"
92 "zero constant term and nonzero coefficient of x^1.\n");
93 flint_abort();
94 }
95
96 if (n < 2)
97 {
98 fmpq_poly_zero(res);
99 return;
100 }
101
102 if (res != poly)
103 {
104 fmpq_poly_fit_length(res, n);
105 _fmpq_poly_revert_series_newton(res->coeffs,
106 res->den, poly->coeffs, poly->den, poly->length, n);
107 }
108 else
109 {
110 fmpq_poly_t t;
111 fmpq_poly_init2(t, n);
112 _fmpq_poly_revert_series_newton(t->coeffs,
113 t->den, poly->coeffs, poly->den, poly->length, n);
114 fmpq_poly_swap(res, t);
115 fmpq_poly_clear(t);
116 }
117
118 _fmpq_poly_set_length(res, n);
119 _fmpq_poly_normalise(res);
120 }
121