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