1 /*
2 Copyright (C) 2012 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "acb_poly.h"
13
14 /* allow changing this from the test code */
15 ARB_DLL slong acb_poly_newton_exp_cutoff = 0;
16
17 /* with inverse=1 simultaneously computes g = exp(-x) to length n
18 with inverse=0 uses g as scratch space, computing
19 g = exp(-x) only to length (n+1)/2 */
20 static void
_acb_poly_exp_series_newton(acb_ptr f,acb_ptr g,acb_srcptr h,slong len,slong prec,int inverse,slong cutoff)21 _acb_poly_exp_series_newton(acb_ptr f, acb_ptr g,
22 acb_srcptr h, slong len, slong prec, int inverse, slong cutoff)
23 {
24 slong alloc;
25 acb_ptr T, U, hprime;
26
27 alloc = 3 * len;
28 T = _acb_vec_init(alloc);
29 U = T + len;
30 hprime = U + len;
31
32 _acb_poly_derivative(hprime, h, len, prec);
33 acb_zero(hprime + len - 1);
34
35 NEWTON_INIT(cutoff, len)
36
37 /* f := exp(h) + O(x^m), g := exp(-h) + O(x^m2) */
38 NEWTON_BASECASE(n)
39 _acb_poly_exp_series_basecase(f, h, n, n, prec);
40 _acb_poly_inv_series(g, f, (n + 1) / 2, (n + 1) / 2, prec);
41 NEWTON_END_BASECASE
42
43 /* extend from length m to length n */
44 NEWTON_LOOP(m, n)
45
46 slong m2 = (m + 1) / 2;
47 slong l = m - 1; /* shifted for derivative */
48
49 /* g := exp(-h) + O(x^m) */
50 _acb_poly_mullow(T, f, m, g, m2, m, prec);
51 _acb_poly_mullow(g + m2, g, m2, T + m2, m - m2, m - m2, prec);
52 _acb_vec_neg(g + m2, g + m2, m - m2);
53
54 /* U := h' + g (f' - f h') + O(x^(n-1))
55 Note: should replace h' by h' mod x^(m-1) */
56 _acb_vec_zero(f + m, n - m);
57 _acb_poly_mullow(T, f, n, hprime, n, n, prec); /* should be mulmid */
58 _acb_poly_derivative(U, f, n, prec); acb_zero(U + n - 1); /* should skip low terms */
59 _acb_vec_sub(U + l, U + l, T + l, n - l, prec);
60 _acb_poly_mullow(T + l, g, n - m, U + l, n - m, n - m, prec);
61 _acb_vec_add(U + l, hprime + l, T + l, n - m, prec);
62
63 /* f := f + f * (h - int U) + O(x^n) = exp(h) + O(x^n) */
64 _acb_poly_integral(U, U, n, prec); /* should skip low terms */
65 _acb_vec_sub(U + m, h + m, U + m, n - m, prec);
66 _acb_poly_mullow(f + m, f, n - m, U + m, n - m, n - m, prec);
67
68 /* g := exp(-h) + O(x^n) */
69 /* not needed if we only want exp(x) */
70 if (n == len && inverse)
71 {
72 _acb_poly_mullow(T, f, n, g, m, n, prec);
73 _acb_poly_mullow(g + m, g, m, T + m, n - m, n - m, prec);
74 _acb_vec_neg(g + m, g + m, n - m);
75 }
76
77 NEWTON_END_LOOP
78
79 NEWTON_END
80
81 _acb_vec_clear(T, alloc);
82 }
83
84 void
_acb_poly_exp_series(acb_ptr f,acb_srcptr h,slong hlen,slong n,slong prec)85 _acb_poly_exp_series(acb_ptr f, acb_srcptr h, slong hlen, slong n, slong prec)
86 {
87 hlen = FLINT_MIN(hlen, n);
88
89 if (hlen == 1)
90 {
91 acb_exp(f, h, prec);
92 _acb_vec_zero(f + 1, n - 1);
93 }
94 else if (n == 2)
95 {
96 acb_exp(f, h, prec);
97 acb_mul(f + 1, f, h + 1, prec); /* safe since hlen >= 2 */
98 }
99 else if (_acb_vec_is_zero(h + 1, hlen - 2)) /* h = a + bx^d */
100 {
101 slong i, j, d = hlen - 1;
102 acb_t t;
103 acb_init(t);
104 acb_set(t, h + d);
105 acb_exp(f, h, prec);
106 for (i = 1, j = d; j < n; j += d, i++)
107 {
108 acb_mul(f + j, f + j - d, t, prec);
109 acb_div_ui(f + j, f + j, i, prec);
110 _acb_vec_zero(f + j - d + 1, hlen - 2);
111 }
112 _acb_vec_zero(f + j - d + 1, n - (j - d + 1));
113 acb_clear(t);
114 }
115 else
116 {
117 slong cutoff;
118
119 if (acb_poly_newton_exp_cutoff != 0)
120 cutoff = acb_poly_newton_exp_cutoff;
121 else if (prec <= 256)
122 cutoff = 750;
123 else
124 cutoff = 1e5 / pow(log(prec), 3);
125
126 if (hlen <= cutoff)
127 {
128 _acb_poly_exp_series_basecase(f, h, hlen, n, prec);
129 }
130 else
131 {
132 acb_ptr g, t;
133 acb_t u;
134 int fix;
135
136 g = _acb_vec_init((n + 1) / 2);
137 fix = (hlen < n || h == f || !acb_is_zero(h));
138
139 if (fix)
140 {
141 t = _acb_vec_init(n);
142 _acb_vec_set(t + 1, h + 1, hlen - 1);
143 }
144 else
145 t = (acb_ptr) h;
146
147 acb_init(u);
148 acb_exp(u, h, prec);
149
150 _acb_poly_exp_series_newton(f, g, t, n, prec, 0, cutoff);
151
152 if (!acb_is_one(u))
153 _acb_vec_scalar_mul(f, f, n, u, prec);
154
155 _acb_vec_clear(g, (n + 1) / 2);
156 if (fix)
157 _acb_vec_clear(t, n);
158 acb_clear(u);
159 }
160 }
161 }
162
163 void
acb_poly_exp_series(acb_poly_t f,const acb_poly_t h,slong n,slong prec)164 acb_poly_exp_series(acb_poly_t f, const acb_poly_t h, slong n, slong prec)
165 {
166 slong hlen = h->length;
167
168 if (n == 0)
169 {
170 acb_poly_zero(f);
171 return;
172 }
173
174 if (hlen == 0)
175 {
176 acb_poly_one(f);
177 return;
178 }
179
180 if (hlen == 1)
181 n = 1;
182
183 acb_poly_fit_length(f, n);
184 _acb_poly_exp_series(f->coeffs, h->coeffs, hlen, n, prec);
185 _acb_poly_set_length(f, n);
186 _acb_poly_normalise(f);
187 }
188