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 : static GEN
23 2162670 : RgX_to_ser_i(GEN x, long l, long v, int copy)
24 : {
25 2162670 : long i = 2, lx = lg(x), vx = varn(x);
26 : GEN y;
27 2162670 : if (lx == 2) return zeroser(vx, minss(l - 2, v));
28 2162572 : if (l <= 2)
29 : {
30 14 : if (l == 2 && v != LONG_MAX) return zeroser(vx, v);
31 0 : pari_err_BUG("RgX_to_ser (l < 2)");
32 : }
33 2162558 : y = cgetg(l,t_SER);
34 2162558 : if (v == LONG_MAX) { v = 1; lx = 3; } /* e.g. Mod(0,3) * x^0 */
35 2162537 : else if (v)
36 : {
37 15505 : long w = v;
38 482958 : while (isrationalzero(gel(x,2))) { x++; w--; }
39 15505 : lx -= v;
40 15505 : if (w)
41 : { /* keep type information, e.g. Mod(0,3) + x */
42 49 : GEN z = gel(x,2); /* = 0 */
43 49 : if (lx <= 2) gel(y,2) = copy? gcopy(z): z;
44 42 : else { x += w; gel(y,2) = gadd(gel(x,2), z); }
45 49 : i++;
46 : }
47 : }
48 2162558 : y[1] = evalvarn(vx) | evalvalser(v); /* must come here because of LONG_MAX */
49 2162558 : if (lx > l) lx = l;
50 : /* 2 <= lx <= l */
51 2162558 : if (copy)
52 210 : for (; i <lx; i++) gel(y,i) = gcopy(gel(x,i));
53 : else
54 8827780 : for (; i <lx; i++) gel(y,i) = gel(x,i);
55 9845099 : for ( ; i < l; i++) gel(y,i) = gen_0;
56 2162558 : return normalizeser(y);
57 : }
58 : /* enlarge/truncate t_POL x to a t_SER with lg l */
59 : GEN
60 2160129 : RgX_to_ser(GEN x, long l) { return RgX_to_ser_i(x, l, RgX_val(x), 0); }
61 : GEN
62 1680 : RgX_to_ser_inexact(GEN x, long l)
63 : {
64 1680 : long i, lx = lg(x);
65 1680 : int first = 1;
66 2912 : for (i = 2; i < lx && gequal0(gel(x,i)); i++) /* ~ RgX_valrem + normalizeser */
67 1232 : if (first && !isexactzero(gel(x,i)))
68 : {
69 14 : pari_warn(warner,"normalizing a series with 0 leading term");
70 14 : first = 0;
71 : }
72 1680 : return RgX_to_ser_i(x, l, i - 2, 0);
73 : }
74 : /* *pd t_POL normalized wrt exact zeros; normalize fully, keeping type
75 : * information */
76 : static long
77 1638 : RgX_valrem_type(GEN *pd, long *warn)
78 : {
79 1638 : GEN d = *pd, z = gel(d,2);
80 : long v;
81 1638 : if (!gequal0(z)) return 0;
82 56 : *warn = 1;
83 56 : if (!signe(d)) { *pd = scalarpol_shallow(z, varn(d)); return degpol(d); }
84 49 : v = RgX_valrem_inexact(d, &d);
85 49 : if (lg(d) > 2)
86 49 : gel(d,2) = gadd(gel(d,2), z);
87 : else
88 0 : d = scalarpol_shallow(z, varn(d));
89 49 : *pd = d; return v;
90 : }
91 : static GEN
92 833 : _rfrac_to_ser(GEN x, long l, long copy)
93 : {
94 833 : GEN a = gel(x,1), d = gel(x,2);
95 833 : long warn = 0, v = varn(d), e;
96 833 : if (l == 2) return zeroser(v, gvaluation(x, pol_x(v)));
97 819 : e = - RgX_valrem(d, &d);
98 819 : e -= RgX_valrem_type(&d, &warn);
99 819 : if (!signe(d)) pari_err_INV("rfrac_to_ser", gel(x,2));
100 819 : if (typ(a) != t_POL || varn(a) != v)
101 : {
102 588 : a = RgX_Rg_mul(RgXn_inv(d, l - 2), a);
103 588 : e += RgX_valrem_type(&a, &warn);
104 : }
105 : else
106 : {
107 231 : e += RgX_valrem(a, &a);
108 231 : e += RgX_valrem_type(&a, &warn);
109 231 : a = RgXn_div(a, d, l - 2);
110 : }
111 819 : if (warn) pari_warn(warner,"normalizing a series with 0 leading term");
112 819 : a = RgX_to_ser_i(a, l, 0, copy);
113 819 : setvalser(a, valser(a) + e); return a;
114 : }
115 : GEN
116 35 : rfrac_to_ser(GEN x, long l) { return _rfrac_to_ser(x, l, 1); }
117 : GEN
118 798 : rfrac_to_ser_i(GEN x, long l) { return _rfrac_to_ser(x, l, 0); }
119 :
120 : static GEN
121 5852 : RgV_to_ser_i(GEN x, long v, long l, int copy)
122 : {
123 5852 : long j, lx = minss(lg(x), l-1);
124 : GEN y;
125 5852 : if (lx == 1) return zeroser(v, l-2);
126 5845 : y = cgetg(l, t_SER); y[1] = evalsigne(1)|evalvarn(v)|evalvalser(0);
127 5845 : x--;
128 5845 : if (copy)
129 507157 : for (j = 2; j <= lx; j++) gel(y,j) = gcopy(gel(x,j));
130 : else
131 277011 : for (j = 2; j <= lx; j++) gel(y,j) = gel(x,j);
132 5964 : for ( ; j < l; j++) gel(y,j) = gen_0;
133 5845 : return normalizeser(y);
134 : }
135 : GEN
136 5579 : RgV_to_ser(GEN x, long v, long l) { return RgV_to_ser_i(x, v, l, 0); }
137 :
138 : /* x a t_SER, prec >= 0 */
139 : GEN
140 1414 : sertoser(GEN x, long prec)
141 : {
142 1414 : long i, lx = lg(x), l;
143 : GEN y;
144 1414 : if (lx == 2) return zeroser(varn(x), prec);
145 1386 : l = prec+2; lx = minss(lx, l);
146 1386 : y = cgetg(l,t_SER); y[1] = x[1];
147 103565 : for (i = 2; i < lx; i++) gel(y,i) = gel(x,i);
148 1876 : for ( ; i < l; i++) gel(y,i) = gen_0;
149 1386 : return y;
150 : }
151 :
152 : /* R(1/x) = x^v * n/d, val(n) = val(d) = 0 */
153 : long
154 112 : rfracrecip(GEN *pn, GEN *pd)
155 : {
156 112 : long v = degpol(*pd);
157 112 : if (typ(*pn) == t_POL && varn(*pn) == varn(*pd))
158 : {
159 63 : v -= degpol(*pn);
160 63 : (void)RgX_valrem(*pn, pn); *pn = RgX_recip(*pn);
161 : }
162 112 : (void)RgX_valrem(*pd, pd); *pd = RgX_recip(*pd);
163 112 : return v;
164 : }
165 :
166 : /* R(1/x) + O(x^N) */
167 : GEN
168 56 : rfracrecip_to_ser_absolute(GEN R, long N)
169 : {
170 56 : GEN n = gel(R,1), d = gel(R,2);
171 56 : long v = rfracrecip(&n, &d); /* R(1/x) = x^v * n/d, val(n) = val(d) = 0 */
172 56 : if (N <= v) return zeroser(varn(d), N);
173 56 : R = rfrac_to_ser_i(mkrfrac(n, d), N-v+2);
174 56 : setvalser(R, v); return R;
175 : }
176 :
177 : /* assume prec >= 0 */
178 : GEN
179 30478 : scalarser(GEN x, long v, long prec)
180 : {
181 : long i, l, s;
182 : GEN y;
183 :
184 30478 : if (isexactzero(x))
185 : {
186 1701 : if (isrationalzero(x)) return zeroser(v, prec);
187 490 : y = cgetg(3, t_SER);
188 490 : y[1] = evalsigne(0) | _evalvalser(prec) | evalvarn(v);
189 490 : gel(y,2) = gcopy(x); return y;
190 : }
191 28777 : l = prec + 2; y = cgetg(l, t_SER); s = !gequal0(x);
192 28777 : y[1] = evalsigne(s) | _evalvalser(0) | evalvarn(v);
193 71386 : gel(y,2) = gcopy(x); for (i=3; i<l; i++) gel(y,i) = gen_0;
194 28777 : return y;
195 : }
196 :
197 : /* assume x a t_[SER|POL], apply gtoser to all coeffs */
198 : static GEN
199 7 : coefstoser(GEN x, long v, long prec)
200 : {
201 : long i, lx;
202 7 : GEN y = cgetg_copy(x, &lx); y[1] = x[1];
203 21 : for (i=2; i<lx; i++) gel(y,i) = gtoser(gel(x,i), v, prec);
204 7 : return y;
205 : }
206 :
207 : static void
208 14 : err_ser_priority(GEN x, long v) { pari_err_PRIORITY("Ser", x, "<", v); }
209 : /* x a t_POL */
210 : static GEN
211 63 : poltoser(GEN x, long v, long prec)
212 : {
213 63 : long s = varncmp(varn(x), v);
214 63 : if (s < 0) err_ser_priority(x,v);
215 56 : if (s > 0) return scalarser(x, v, prec);
216 42 : return RgX_to_ser_i(x, prec+2, RgX_val(x), 1);
217 : }
218 : /* x a t_RFRAC */
219 : static GEN
220 77 : rfractoser(GEN x, long v, long prec)
221 : {
222 77 : long s = varncmp(varn(gel(x,2)), v);
223 : pari_sp av;
224 77 : if (s < 0) err_ser_priority(x,v);
225 70 : if (s > 0) return scalarser(x, v, prec);
226 35 : av = avma; return gerepileupto(av, rfrac_to_ser(x, prec+2));
227 : }
228 : GEN
229 4057648 : toser_i(GEN x)
230 : {
231 4057648 : switch(typ(x))
232 : {
233 140532 : case t_SER: return x;
234 1680 : case t_POL: return RgX_to_ser_inexact(x, precdl+2);
235 140 : case t_RFRAC: return rfrac_to_ser_i(x, precdl+2);
236 : }
237 3915296 : return NULL;
238 : }
239 :
240 : /* conversion: prec ignored if t_VEC or t_SER in variable v */
241 : GEN
242 434 : gtoser(GEN x, long v, long prec)
243 : {
244 434 : long tx = typ(x);
245 :
246 434 : if (v < 0) v = 0;
247 434 : if (prec < 0) pari_err_DOMAIN("Ser", "precision", "<", gen_0, stoi(prec));
248 434 : if (tx == t_SER)
249 : {
250 35 : long s = varncmp(varn(x), v);
251 35 : if (s < 0) return coefstoser(x, v, prec);
252 28 : else if (s > 0) return scalarser(x, v, prec);
253 14 : return gcopy(x);
254 : }
255 399 : if (is_scalar_t(tx)) return scalarser(x,v,prec);
256 357 : switch(tx)
257 : {
258 63 : case t_POL: return poltoser(x, v, prec);
259 77 : case t_RFRAC: return rfractoser(x, v, prec);
260 42 : case t_QFB: return RgV_to_ser_i(x, v, 4+1, 1);
261 14 : case t_VECSMALL: x = zv_to_ZV(x);/*fall through*/
262 168 : case t_VEC: case t_COL:
263 168 : if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("Ser", x, "<=", v);
264 161 : return RgV_to_ser_i(x, v, lg(x)+1, 1);
265 : }
266 7 : pari_err_TYPE("Ser",x);
267 : return NULL; /* LCOV_EXCL_LINE */
268 : }
269 : /* impose prec */
270 : GEN
271 175 : gtoser_prec(GEN x, long v, long prec)
272 : {
273 175 : pari_sp av = avma;
274 175 : if (v < 0) v = 0;
275 175 : if (prec < 0) pari_err_DOMAIN("Ser", "precision", "<", gen_0, stoi(prec));
276 168 : switch(typ(x))
277 : {
278 28 : case t_SER: if (varn(x) != v) break;
279 21 : return gerepilecopy(av, sertoser(x, prec));
280 28 : case t_QFB:
281 28 : x = RgV_to_ser_i(mkvec3(gel(x,1),gel(x,2),gel(x,3)), v, prec+2, 1);
282 28 : return gerepileupto(av, x);
283 14 : case t_VECSMALL: x = zv_to_ZV(x);/*fall through*/
284 42 : case t_VEC: case t_COL:
285 42 : if (varncmp(gvar(x), v) <= 0) pari_err_PRIORITY("Ser", x, "<=", v);
286 42 : return RgV_to_ser_i(x, v, prec+2, 1);
287 : }
288 77 : return gtoser(x, v, prec);
289 : }
290 : GEN
291 504 : Ser0(GEN x, long v, GEN d, long prec)
292 : {
293 504 : if (!d) return gtoser(x, v, prec);
294 175 : if (typ(d) != t_INT)
295 : {
296 7 : d = gceil(d);
297 7 : if (typ(d) != t_INT) pari_err_TYPE("Ser [precision]",d);
298 : }
299 175 : return gtoser_prec(x, v, itos(d));
300 : }
|
