Line data Source code
1 : /* Copyright (C) 2000, 2012 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 : /*******************************************************************/
18 : /* */
19 : /* Conversion --> t_SER */
20 : /* */
21 : /*******************************************************************/
22 : /* x RgX in variable t, return x * (1 + O(t^(l-2))) assuming l >= 2 and
23 : * e = v_t(x). Length and valuation changed by normalization, stripping
24 : * leading zero terms. Shallow if copy = 0, else gc clean */
25 : static GEN
26 2592888 : RgX_to_ser_i(GEN x, long l, long e, int copy)
27 : {
28 2592888 : long i = 2, lx = lg(x), vx = varn(x);
29 : GEN y;
30 2592888 : if (lx == 2) return zeroser(vx, minss(l - 2, e));
31 2592790 : if (l <= 2)
32 : {
33 14 : if (l == 2 && e != LONG_MAX) return zeroser(vx, e);
34 0 : pari_err_BUG("RgX_to_ser (l < 2)");
35 : }
36 2592776 : y = cgetg(l,t_SER);
37 2592776 : if (e == LONG_MAX) { e = 1; lx = 3; } /* e.g. Mod(0,3) * x^0 */
38 2592755 : else if (e)
39 : {
40 17199 : long w = e;
41 486353 : while (isrationalzero(gel(x,2))) { x++; w--; }
42 17199 : lx -= e;
43 17199 : if (w)
44 : { /* keep type information, e.g. Mod(0,3) + x */
45 49 : GEN z = gel(x,2); /* = 0 */
46 49 : if (lx <= 2) gel(y,2) = copy? gcopy(z): z;
47 42 : else { x += w; gel(y,2) = gadd(gel(x,2), z); }
48 49 : i++;
49 : }
50 : }
51 2592776 : y[1] = evalvarn(vx) | evalvalser(e); /* must come here because of LONG_MAX */
52 2592776 : if (lx > l) lx = l;
53 : /* 2 <= lx <= l */
54 2592776 : if (copy)
55 210 : for (; i <lx; i++) gel(y,i) = gcopy(gel(x,i));
56 : else
57 10536934 : for (; i <lx; i++) gel(y,i) = gel(x,i);
58 11838715 : for ( ; i < l; i++) gel(y,i) = gen_0;
59 2592776 : return normalizeser(y);
60 : }
61 : /* enlarge/truncate t_POL x to a t_SER with lg l */
62 : GEN
63 2590312 : RgX_to_ser(GEN x, long l) { return RgX_to_ser_i(x, l, RgX_val(x), 0); }
64 : GEN
65 1673 : RgX_to_ser_inexact(GEN x, long l)
66 : {
67 1673 : long i, lx = lg(x);
68 1673 : int first = 1;
69 2898 : for (i = 2; i < lx && gequal0(gel(x,i)); i++) /* RgX_valrem + normalizeser */
70 1225 : if (first && !isexactzero(gel(x,i)))
71 : {
72 14 : pari_warn(warner,"normalizing a series with 0 leading term");
73 14 : first = 0;
74 : }
75 1673 : return RgX_to_ser_i(x, l, i - 2, 0);
76 : }
77 : /* *pd t_POL normalized wrt exact zeros; normalize fully, keeping type
78 : * information */
79 : static long
80 1722 : RgX_valrem_type(GEN *pd, long *warn)
81 : {
82 1722 : GEN d = *pd, z = gel(d,2);
83 : long v;
84 1722 : if (!gequal0(z)) return 0;
85 56 : *warn = 1;
86 56 : if (!signe(d)) { *pd = scalarpol_shallow(z, varn(d)); return degpol(d); }
87 49 : v = RgX_valrem_inexact(d, &d);
88 49 : if (lg(d) > 2)
89 49 : gel(d,2) = gadd(gel(d,2), z);
90 : else
91 0 : d = scalarpol_shallow(z, varn(d));
92 49 : *pd = d; return v;
93 : }
94 : static GEN
95 875 : _rfrac_to_ser(GEN x, long l, long copy)
96 : {
97 875 : GEN a = gel(x,1), d = gel(x,2);
98 875 : long warn = 0, v = varn(d), e;
99 875 : if (l == 2) return zeroser(v, gvaluation(x, pol_x(v)));
100 861 : e = - RgX_valrem(d, &d);
101 861 : e -= RgX_valrem_type(&d, &warn);
102 861 : if (!signe(d)) pari_err_INV("rfrac_to_ser", gel(x,2));
103 861 : if (typ(a) != t_POL || varn(a) != v)
104 : {
105 623 : a = RgX_Rg_mul(RgXn_inv(d, l - 2), a);
106 623 : e += RgX_valrem_type(&a, &warn);
107 : }
108 : else
109 : {
110 238 : e += RgX_valrem(a, &a);
111 238 : e += RgX_valrem_type(&a, &warn);
112 238 : a = RgXn_div(a, d, l - 2);
113 : }
114 861 : if (warn) pari_warn(warner,"normalizing a series with 0 leading term");
115 861 : a = RgX_to_ser_i(a, l, 0, copy);
116 861 : setvalser(a, valser(a) + e); return a;
117 : }
118 : GEN
119 35 : rfrac_to_ser(GEN x, long l) { return _rfrac_to_ser(x, l, 1); }
120 : GEN
121 840 : rfrac_to_ser_i(GEN x, long l) { return _rfrac_to_ser(x, l, 0); }
122 :
123 : static GEN
124 6090 : RgV_to_ser_i(GEN x, long v, long l, int copy)
125 : {
126 6090 : long j, lx = minss(lg(x), l-1);
127 : GEN y;
128 6090 : if (lx == 1) return zeroser(v, l-2);
129 6083 : y = cgetg(l, t_SER); y[1] = evalsigne(1)|evalvarn(v)|evalvalser(0);
130 6083 : x--;
131 6083 : if (copy)
132 507157 : for (j = 2; j <= lx; j++) gel(y,j) = gcopy(gel(x,j));
133 : else
134 287259 : for (j = 2; j <= lx; j++) gel(y,j) = gel(x,j);
135 6202 : for ( ; j < l; j++) gel(y,j) = gen_0;
136 6083 : return normalizeser(y);
137 : }
138 : GEN
139 5817 : RgV_to_ser(GEN x, long v, long l) { return RgV_to_ser_i(x, v, l, 0); }
140 :
141 : /* x a t_SER, prec >= 0 */
142 : GEN
143 1484 : sertoser(GEN x, long prec)
144 : {
145 1484 : long i, lx = lg(x), l;
146 : GEN y;
147 1484 : if (lx == 2) return zeroser(varn(x), prec);
148 1456 : l = prec+2; lx = minss(lx, l);
149 1456 : y = cgetg(l,t_SER); y[1] = x[1];
150 103824 : for (i = 2; i < lx; i++) gel(y,i) = gel(x,i);
151 2009 : for ( ; i < l; i++) gel(y,i) = gen_0;
152 1456 : return y;
153 : }
154 :
155 : /* R(1/x) = x^v * n/d, val(n) = val(d) = 0 */
156 : long
157 147 : rfracrecip(GEN *pn, GEN *pd)
158 : {
159 147 : long v = degpol(*pd);
160 147 : if (typ(*pn) == t_POL && varn(*pn) == varn(*pd))
161 : {
162 70 : v -= degpol(*pn);
163 70 : (void)RgX_valrem(*pn, pn); *pn = RgX_recip(*pn);
164 : }
165 147 : (void)RgX_valrem(*pd, pd); *pd = RgX_recip(*pd);
166 147 : return v;
167 : }
168 :
169 : /* R(1/x) + O(x^N) */
170 : GEN
171 84 : rfracrecip_to_ser_absolute(GEN R, long N)
172 : {
173 84 : GEN n = gel(R,1), d = gel(R,2);
174 84 : long v = rfracrecip(&n, &d); /* R(1/x) = x^v * n/d, val(n) = val(d) = 0 */
175 84 : if (N <= v) return zeroser(varn(d), N);
176 84 : R = rfrac_to_ser_i(mkrfrac(n, d), N-v+2);
177 84 : setvalser(R, v); return R;
178 : }
179 :
180 : /* assume prec >= 0 */
181 : GEN
182 34083 : scalarser(GEN x, long v, long prec)
183 : {
184 : long i, l, s;
185 : GEN y;
186 :
187 34083 : if (isexactzero(x))
188 : {
189 1554 : if (isrationalzero(x)) return zeroser(v, prec);
190 483 : y = cgetg(3, t_SER);
191 483 : y[1] = evalsigne(0) | _evalvalser(prec) | evalvarn(v);
192 483 : gel(y,2) = gcopy(x); return y;
193 : }
194 32529 : l = prec + 2; y = cgetg(l, t_SER); s = !gequal0(x);
195 32529 : y[1] = evalsigne(s) | _evalvalser(0) | evalvarn(v);
196 81683 : gel(y,2) = gcopy(x); for (i=3; i<l; i++) gel(y,i) = gen_0;
197 32529 : return y;
198 : }
199 :
200 : /* assume x a t_[SER|POL], apply gtoser to all coeffs */
201 : static GEN
202 7 : coefstoser(GEN x, long v, long prec)
203 : {
204 : long i, lx;
205 7 : GEN y = cgetg_copy(x, &lx); y[1] = x[1];
206 21 : for (i=2; i<lx; i++) gel(y,i) = gtoser(gel(x,i), v, prec);
207 7 : return y;
208 : }
209 :
210 : static void
211 14 : err_ser_priority(GEN x, long v) { pari_err_PRIORITY("Ser", x, "<", v); }
212 : /* x a t_POL */
213 : static GEN
214 63 : poltoser(GEN x, long v, long prec)
215 : {
216 63 : long s = varncmp(varn(x), v);
217 63 : if (s < 0) err_ser_priority(x,v);
218 56 : if (s > 0) return scalarser(x, v, prec);
219 42 : return RgX_to_ser_i(x, prec+2, RgX_val(x), 1);
220 : }
221 : /* x a t_RFRAC */
222 : static GEN
223 77 : rfractoser(GEN x, long v, long prec)
224 : {
225 77 : long s = varncmp(varn(gel(x,2)), v);
226 : pari_sp av;
227 77 : if (s < 0) err_ser_priority(x,v);
228 70 : if (s > 0) return scalarser(x, v, prec);
229 35 : av = avma; return gc_upto(av, rfrac_to_ser(x, prec+2));
230 : }
231 : GEN
232 4230971 : toser_i(GEN x)
233 : {
234 4230971 : switch(typ(x))
235 : {
236 156037 : case t_SER: return x;
237 1673 : case t_POL: return RgX_to_ser_inexact(x, precdl+2);
238 140 : case t_RFRAC: return rfrac_to_ser_i(x, precdl+2);
239 : }
240 4073121 : return NULL;
241 : }
242 :
243 : /* conversion: prec ignored if t_VEC or t_SER in variable v */
244 : GEN
245 434 : gtoser(GEN x, long v, long prec)
246 : {
247 434 : long tx = typ(x);
248 :
249 434 : if (v < 0) v = 0;
250 434 : if (prec < 0) pari_err_DOMAIN("Ser", "precision", "<", gen_0, stoi(prec));
251 434 : if (tx == t_SER)
252 : {
253 35 : long s = varncmp(varn(x), v);
254 35 : if (s < 0) return coefstoser(x, v, prec);
255 28 : else if (s > 0) return scalarser(x, v, prec);
256 14 : return gcopy(x);
257 : }
258 399 : if (is_scalar_t(tx)) return scalarser(x,v,prec);
259 357 : switch(tx)
260 : {
261 63 : case t_POL: return poltoser(x, v, prec);
262 77 : case t_RFRAC: return rfractoser(x, v, prec);
263 42 : case t_QFB: return RgV_to_ser_i(x, v, 4+1, 1);
264 14 : case t_VECSMALL: x = zv_to_ZV(x);/*fall through*/
265 168 : case t_VEC: case t_COL:
266 168 : if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("Ser", x, "<=", v);
267 161 : return RgV_to_ser_i(x, v, lg(x)+1, 1);
268 : }
269 7 : pari_err_TYPE("Ser",x);
270 : return NULL; /* LCOV_EXCL_LINE */
271 : }
272 : /* impose prec */
273 : GEN
274 175 : gtoser_prec(GEN x, long v, long prec)
275 : {
276 175 : pari_sp av = avma;
277 175 : if (v < 0) v = 0;
278 175 : if (prec < 0) pari_err_DOMAIN("Ser", "precision", "<", gen_0, stoi(prec));
279 168 : switch(typ(x))
280 : {
281 28 : case t_SER: if (varn(x) != v) break;
282 21 : return gc_GEN(av, sertoser(x, prec));
283 28 : case t_QFB:
284 28 : x = RgV_to_ser_i(mkvec3(gel(x,1),gel(x,2),gel(x,3)), v, prec+2, 1);
285 28 : return gc_upto(av, x);
286 14 : case t_VECSMALL: x = zv_to_ZV(x);/*fall through*/
287 42 : case t_VEC: case t_COL:
288 42 : if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("Ser", x, "<=", v);
289 42 : return RgV_to_ser_i(x, v, prec+2, 1);
290 : }
291 77 : return gtoser(x, v, prec);
292 : }
293 : GEN
294 504 : Ser0(GEN x, long v, GEN d, long prec)
295 : {
296 504 : if (!d) return gtoser(x, v, prec);
297 175 : if (typ(d) != t_INT)
298 : {
299 7 : d = gceil(d);
300 7 : if (typ(d) != t_INT) pari_err_TYPE("Ser [precision]",d);
301 : }
302 175 : return gtoser_prec(x, v, itos(d));
303 : }
|
