1 /*
2 Copyright (C) 2010, 2011 Sebastian Pancratz
3 Copyright (C) 2013 Mike Hansen
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 #ifdef T
14
15 #include "templates.h"
16
17 #define FQ_POLY_INV_NEWTON_CUTOFF 64
18
19 void
_TEMPLATE(T,poly_inv_series_newton)20 _TEMPLATE(T, poly_inv_series_newton) (TEMPLATE(T, struct) * Qinv,
21 const TEMPLATE(T, struct) * Q, slong n,
22 const TEMPLATE(T, t) cinv,
23 const TEMPLATE(T, ctx_t) ctx)
24 {
25 if (n == 1) /* {Q,1} x* cinv == 1 mod (x) */
26 {
27 TEMPLATE(T, set) (Qinv, cinv, ctx);
28 }
29 else
30 {
31 const slong alloc = FLINT_MAX(n, 3 * FQ_POLY_INV_NEWTON_CUTOFF);
32 slong *a, i, m;
33 TEMPLATE(T, struct) * W;
34
35 W = _TEMPLATE(T, vec_init) (alloc, ctx);
36
37 for (i = 1; (WORD(1) << i) < n; i++) ;
38
39 a = (slong *) flint_malloc(i * sizeof(slong));
40 a[i = 0] = n;
41 while (n >= FQ_POLY_INV_NEWTON_CUTOFF)
42 a[++i] = (n = (n + 1) / 2);
43
44 /* Base case */
45 {
46 TEMPLATE(T, struct) * Qrev = W + 2 * FQ_POLY_INV_NEWTON_CUTOFF;
47
48 _TEMPLATE(T, poly_reverse) (Qrev, Q, n, n, ctx);
49 _TEMPLATE(T, vec_zero) (W, 2 * n - 2, ctx);
50 TEMPLATE(T, one) (W + (2 * n - 2), ctx);
51 _TEMPLATE(T, poly_div_basecase) (Qinv, W, W, 2 * n - 1, Qrev, n,
52 cinv, ctx);
53 _TEMPLATE(T, poly_reverse) (Qinv, Qinv, n, n, ctx);
54 }
55
56 for (i--; i >= 0; i--)
57 {
58 m = n;
59 n = a[i];
60
61 _TEMPLATE(T, poly_mullow) (W, Q, n, Qinv, m, n, ctx);
62 _TEMPLATE(T, poly_mullow) (Qinv + m, Qinv, m, W + m, n - m, n - m,
63 ctx);
64 _TEMPLATE(T, poly_neg) (Qinv + m, Qinv + m, n - m, ctx);
65 }
66
67 _TEMPLATE(T, vec_clear) (W, alloc, ctx);
68 flint_free(a);
69 }
70 }
71
72 void
TEMPLATE(T,poly_inv_series_newton)73 TEMPLATE(T, poly_inv_series_newton) (TEMPLATE(T, poly_t) Qinv,
74 const TEMPLATE(T, poly_t) Q, slong n,
75 const TEMPLATE(T, ctx_t) ctx)
76 {
77 TEMPLATE(T, t) cinv;
78 TEMPLATE(T, struct) * Qcopy;
79 int Qalloc;
80
81 if (Q->length >= n)
82 {
83 Qcopy = Q->coeffs;
84 Qalloc = 0;
85 }
86 else
87 {
88 Qcopy = _TEMPLATE(T, vec_init) (n, ctx);
89 _TEMPLATE(T, vec_set) (Qcopy, Q->coeffs, Q->length, ctx);
90 Qalloc = 1;
91 }
92
93 TEMPLATE(T, init) (cinv, ctx);
94 TEMPLATE(T, inv) (cinv, Q->coeffs, ctx);
95
96 if (Qinv != Q)
97 {
98 TEMPLATE(T, poly_fit_length) (Qinv, n, ctx);
99 _TEMPLATE(T, poly_inv_series_newton) (Qinv->coeffs, Qcopy, n, cinv,
100 ctx);
101 }
102 else
103 {
104 TEMPLATE(T, struct) * t = _TEMPLATE(T, vec_init) (n, ctx);
105
106 _TEMPLATE(T, poly_inv_series_newton) (t, Qcopy, n, cinv, ctx);
107
108 _TEMPLATE(T, vec_clear) (Qinv->coeffs, Qinv->alloc, ctx);
109 Qinv->coeffs = t;
110 Qinv->alloc = n;
111 Qinv->length = n;
112 }
113 _TEMPLATE(T, poly_set_length) (Qinv, n, ctx);
114 _TEMPLATE(T, poly_normalise) (Qinv, ctx);
115
116 if (Qalloc)
117 _TEMPLATE(T, vec_clear) (Qcopy, n, ctx);
118 TEMPLATE(T, clear) (cinv, ctx);
119 }
120
121
122 #endif
123