1 /*
2 Copyright (C) 2013 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 /* (a + bx^c)^g where a = f[0] and b = f[flen-1] */
15 void
_acb_poly_binomial_pow_acb_series(acb_ptr h,acb_srcptr f,slong flen,const acb_t g,slong len,slong prec)16 _acb_poly_binomial_pow_acb_series(acb_ptr h, acb_srcptr f, slong flen, const acb_t g, slong len, slong prec)
17 {
18 slong i, j, d;
19 acb_t t;
20
21 acb_init(t);
22
23 d = flen - 1;
24 acb_pow(h, f, g, prec);
25 acb_div(t, f + d, f, prec);
26
27 for (i = 1, j = d; j < len; i++, j += d)
28 {
29 acb_sub_ui(h + j, g, i - 1, prec);
30 acb_mul(h + j, h + j, h + j - d, prec);
31 acb_mul(h + j, h + j, t, prec);
32 acb_div_ui(h + j, h + j, i, prec);
33 }
34
35 if (d > 1)
36 {
37 for (i = 1; i < len; i++)
38 if (i % d != 0)
39 acb_zero(h + i);
40 }
41
42 acb_clear(t);
43 return;
44 }
45
46 void
_acb_poly_pow_acb_series(acb_ptr h,acb_srcptr f,slong flen,const acb_t g,slong len,slong prec)47 _acb_poly_pow_acb_series(acb_ptr h,
48 acb_srcptr f, slong flen, const acb_t g, slong len, slong prec)
49 {
50 int f_binomial, g_exact, g_int;
51
52 while (flen > 0 && acb_is_zero(f + flen - 1))
53 flen--;
54
55 if (flen <= 1)
56 {
57 acb_pow(h, f, g, prec);
58 _acb_vec_zero(h + 1, len - 1);
59 return;
60 }
61
62 g_exact = acb_is_exact(g);
63 g_int = acb_is_real(g) && arb_is_int(acb_realref(g));
64 f_binomial = _acb_vec_is_zero(f + 1, flen - 2);
65
66 /* g = small integer */
67 if (g_exact && g_int &&
68 arf_cmpabs_2exp_si(arb_midref(acb_realref(g)), FLINT_BITS - 1) < 0)
69 {
70 slong e, hlen;
71
72 e = arf_get_si(arb_midref(acb_realref(g)), ARF_RND_DOWN);
73 hlen = poly_pow_length(flen, FLINT_ABS(e), len);
74
75 if (e >= 0)
76 {
77 _acb_poly_pow_ui_trunc_binexp(h, f, flen, e, hlen, prec);
78 _acb_vec_zero(h + hlen, len - hlen);
79 return;
80 }
81 else if (!f_binomial)
82 {
83 acb_ptr t;
84 t = _acb_vec_init(hlen);
85 _acb_poly_pow_ui_trunc_binexp(t, f, flen, -e, hlen, prec);
86 _acb_poly_inv_series(h, t, hlen, len, prec);
87 _acb_vec_clear(t, hlen);
88 return;
89 }
90 }
91
92 /* (a + bx^c)^g */
93 if (f_binomial)
94 {
95 _acb_poly_binomial_pow_acb_series(h, f, flen, g, len, prec);
96 return;
97 }
98
99 /* g = +/- 1/2 */
100 if (g_exact && acb_is_real(g) && arf_cmpabs_2exp_si(arb_midref(acb_realref(g)), -1) == 0)
101 {
102 if (arf_sgn(arb_midref(acb_realref(g))) > 0)
103 _acb_poly_sqrt_series(h, f, flen, len, prec);
104 else
105 _acb_poly_rsqrt_series(h, f, flen, len, prec);
106 return;
107 }
108
109 /* f^g = exp(g*log(f)) */
110 _acb_poly_log_series(h, f, flen, len, prec);
111 _acb_vec_scalar_mul(h, h, len, g, prec);
112 _acb_poly_exp_series(h, h, len, len, prec);
113
114 }
115
116 void
acb_poly_pow_acb_series(acb_poly_t h,const acb_poly_t f,const acb_t g,slong len,slong prec)117 acb_poly_pow_acb_series(acb_poly_t h,
118 const acb_poly_t f, const acb_t g, slong len, slong prec)
119 {
120 slong flen;
121
122 flen = f->length;
123 flen = FLINT_MIN(flen, len);
124
125 if (len == 0)
126 {
127 acb_poly_zero(h);
128 return;
129 }
130
131 if (acb_is_zero(g))
132 {
133 acb_poly_one(h);
134 return;
135 }
136
137 if (flen == 0)
138 {
139 acb_poly_zero(h);
140 return;
141 }
142
143 if (f == h)
144 {
145 acb_poly_t t;
146 acb_poly_init2(t, len);
147 _acb_poly_pow_acb_series(t->coeffs, f->coeffs, flen, g, len, prec);
148 _acb_poly_set_length(t, len);
149 _acb_poly_normalise(t);
150 acb_poly_swap(t, h);
151 acb_poly_clear(t);
152 }
153 else
154 {
155 acb_poly_fit_length(h, len);
156 _acb_poly_pow_acb_series(h->coeffs, f->coeffs, flen, g, len, prec);
157 _acb_poly_set_length(h, len);
158 _acb_poly_normalise(h);
159 }
160 }
161
162